#! /usr/bin/python3 
# Last edited on 2025-08-10 09:08:41 by stolfi

???

import numpy
from math import sin, cos, log, exp, floor, pi, inf
import sys, re
from error_funcs import prog_error, arg_error
from process_funcs import bash, run_command
import argparser

def main():

  page, clip, mktag = parse_options()

  # Read the raw pixel data:
  M, bands = load_pixels(page, clip, mktag)
  np, nb = numpy.shape(M)
  debug("M[0,:]", M[0,:])
  debug("M[:,0]", M[:,0])
  
  # Compute {N} whose rows are the displacements from barycenter:
  bar = barycenter(M)
  debug("max(abs(bar))", numpy.max(numpy.abs(bar)))
  N = numpy.array([ M[ip] - bar for ip in range(np) ])
  debug("max(abs(N))", numpy.max(numpy.abs(N)))
  
  # Compute the principal directions of {N'*N}:
  C = (N.transpose() @ N) / float(np)
  debug("max(C)", numpy.max(C))
  evv = numpy.linalg.eigh(C)
  sys.stderr.write("eigenvalues:\n")
  print(evv.eigenvalues)

  sys.stderr.write("eigenvectors:\n")
  print(evv.eigenvectors)
  
  # Generate images with last three eigenvectors:
  write_eigens(page, clip, mktag, bands, evv, 10)

  return 0
  # ----------------------------------------------------------------------
   
def write_eigens(page, clip, mktag, bands, evv, ne):
  # Writes a file with the eigenvectors with the {ne} largest eigenvalues.
  # Each line has the name of a band, e.g. "MB570AM_027_F",
  # and the corresponding coordinate of each of the {ne}
  # eigenvectors.
  #
  nb = numpy.size(evv.eigenvalues,0)
  efile = f"deriv/davis/{page}/{clip}_{mktag}.eig"
  with open(efile, "w") as wr:
    wr.write(f"nb = {nb}\n")
    wr.write(f"ne = {ne}\n")
    for ib in range(nb):
      wr.write(f"{bands[ib]}")
      for ie in range(ne):
        wr.write(" %+7.4f" % evv.eigenvectors[ib,nb-1-ie])
      wr.write("\n")
    wr.close()
  return
  # ----------------------------------------------------------------------
      
def load_pixels(page, clip, mktag):
  ifile = f"deriv/davis/{page}/{clip}_{mktag}.npy"
  M = numpy.load(ifile)
  debug("shape(M)", numpy.shape(M))
  np, 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

def barycenter(M):
  # Expects {M} to be a {numpy} matrix. Computes the 
  # average of all rows. 
  bar = numpy.sum(M, 0)
  debug("shape(bar)", numpy.shape(bar))
  bar = bar/numpy.size(M, 0)
  return bar
  # ----------------------------------------------------------------------

def parse_options():
  page = sys.argv[1]
  clip = sys.argv[2]   # Clip from which pixels were extracted.       
  mktag = sys.argv[3]  # Tag of mask used to exatract them.
  return page, clip, mktag
  # ----------------------------------------------------------------------
  
def debug(name, arg):
  sys.stderr.write(f"{name}");
  if isinstance(arg, list) or isinstance(arg, tuple) or isinstance(arg, dict):
    sys.stderr.write(f" (length = {len(arg)}):\n")
    sys.stderr.write(r"  ")
  else:
    sys.stderr.write(r" = ")
  sys.stderr.write(f"{arg}\n")
  return
  # ----------------------------------------------------------------------

main()