#! /usr/bin/python3
# Last edited on 2025-08-23 16:17:46 by stolfi

  # Collect the relevant spectra and average them:
  smin = numpy.zeros((nb,)); nmin = 0
  smax = numpy.zeros((nb,)); nmax = 0
  for ks in range(ns):
    spix = S[ks,:]
    tpix = numpy.sum(spix)
    if tpix >= talb_min_clip and tpix <= talb_min_take:
      smin += spix; nmin += 1
    if tpix >= talb_max_take and tpix <= talb_max_clip:
      smax += spix; nmax += 1
  debug("nmin", nmin)
  debug("nmax", nmax)
  assert nmin >= 1 and nmax >= 1, "no pixels in range"
  smin /= nmin;
  smax /= nmax;



def compute_deviations_of_analysis_data(masks, bands, savg, smin, smax, rdev, P):
  # For each {mask} in {masks[0..nm]}, writes a file with the analysis
  # results for the subset of samples of clip {clip} in page {page} selected by the 
  # mask {mask}, namely {savg,smin,smax,rdev,P} and spectral bands with names
  # {bands[0..nb-1]}. See the top of this module for the format of this
  # file.
  #
  nb, ne, nm = numpy.shape(P)
  assert len(masks) == nm, "inconsitent {masks}"
  assert len(bands) == nb, "inconsitent {bands}"
  assert numpy.shape(savg) == (nb,nm,)
  assert numpy.shape(smin) == (nb,nm,)
  assert numpy.shape(smax) == (nb,nm,)
  assert numpy.shape(rdev) == (nb,nm,)

  for im in range(nm):
    mask = masks[im]
    sum2_savg = numpy.zeros((nm,))
    sum2_smin = numpy.zeros((nm,))
    sum2_smax = numpy.zeros((nm,))
    sum2_P = numpy.zeros((ne,nm,))
    sum2_rdev = numpy.zeros((nm,))
    for ib in range(nb+1):
      band_tag = bands[ib] if ib < nb else "DEVIATIONS"
      efmt = " %+7.4f" if ib < nb else " %7.4f"

      for svar, sum2_svar in \
        ( (savg, sum2_savg),
          (smin, sum2_smin),
          (smax, sum2_smax),
          (rdev, sum2_rdev)
        ):
        svar_val = svar[ib,im] if ib < nb else sqrt(sum2_svar[im])
        wr.write(efmt % svar_val)
        if ib < nb: sum2_svar[im] += svar_val**2



      for ie in range(ne):
        P_val = P[ib,ie,im] if ib < nb else sqrt(sum2_P[ie,im])
        if ib < nb: sum2_P[ie,im] += P_val**2

  return
  # ----------------------------------------------------------------------


########################################################################

def read_normalization_mask_image(page, clip, mask):
  mfile = f"pages/{page}/clips/{clip}/masks/{mask}.png"
  tfile = ".temp.pgm"
  bash(f"convert {mfile} -depth 16 -colorspace LinearGray {tfile}")
  M = read_pgm_image(tfile, 65535)
  return M
  # ----------------------------------------------------------------------

########################################################################
  
  #  Full color range minus zero:
  hblack = "#000100010001"
  hwhite = "#ffffffffffff"
  
  # Multiply raw clip image by mask, eliminate other zeros:
  pgdir_in = f"MS/davis/{page}"
  cfile = f"{pgdir_in}/{band}/clips/{clip}.pgm"
  mfile = f"pages/{page}/clips/{clip}/masks/{mask}.png"
  xfile = f".{clip}.pgm"
  vfile = f".{clip}.vals"
  bash \
    ( \
      f"convert " + \
      f" \\( {cfile} -set colorspace LinearGray +level-colors '{hblack},{hwhite}' \\)" + \
      f" \\( {mfile} -set colorspace LinearGray -depth 16 -colorspace LinearGray \\)" + \
      f" -compose multiply" + \
      f" -composite" + \
      f" -depth 16 -colorspace LinearGray" + \
      f" {xfile}"
    )
    
  # Extract nonzero sample values in scanline order:
  bash \
    ( \
      f"cat {xfile}" + \
      r" | pnmtoplainpnm" + \
      r" | tr ' ' '\012'" + \
      r" | egrep -e '.'" + \
      r" | gawk '// { if ((FNR > 4) && ($1 > 0)) { print $1; } }'" + \
      f" > {vfile}"
    )
    
  # Read samples, convert to float in {[0 _ 1]}, make a list:
  samples = []
  with open(vfile) as rd:
    nbig = 0; ismp_max = 0
    for line in rd:
      line = re.sub(r"[#].*$", "", line)
      if line.strip() != "":
        ismp = int(line.strip())
        if ismp != 0:
          if ismp > maxval:
            nbig += 1
            if ismp > ismp_max: ismp_max = ismp
            ismp = maxval
          fval = float(ismp)/float(maxval)
          samples.append(fval)
    rd.close()
    if nbig > 0:
      sys.stderr.write(f"{nbig} samples above {maxval}, max = {ismp_max}\n")


# ######################################################################

def save_pixels(M, bands, page, vset):
  odir = f"pages/{page}/vsets/{vset}"
  bash(f"mkdir -p {odir}")
  ofile = f"{odir}/spectra.npy"
  numpy.save(ofile, M)
  bfile = f"odir/bands.txt"
  nb = len(bands)
  assert nb == numpy.size(M,1)
  with open(bfile, "w") as wr:
    for ib in range(nb):
      wr.write(f"{bands[ib]}\n")
    wr.close()
  return 
  # ----------------------------------------------------------------------
      
def load_pixels(page, clip,mask):
  ifile = f"deriv/davis/{page}/{clip}_{mktag}.npy"
  M = numpy.load(ifile)
  debug("shape(M)", numpy.shape(M))
  ns, nb = numpy.shape(M)
  bfile = f"deriv/davis/{page}/{clip}_{mktag}.bands"
  bands = []
  with open(bfile, "r") as rd:
    for line in rd:
      bands.append(line.strip())
  if len(bands) != nb: prog_error("M/bands mismatch")
  return M, bands