#! /usr/bin/python3
# Last edited on 2023-04-18 06:01:05 by stolfi

# Writes to stdout a product-to-term (prit) tabe that lists all products
# of pairs of coefficients fronm the canonical basis (CANC) for a 3x3
# window, mapping to the same term those pairs that should have same
# weight in a sharpness estimator because of rotation and flip
# symmetries.

def main():
  ns = 9; # Number of samples in window.
  wrer("generating terms...\n");
  prit = make_all_CANC_terms(ns); 
  nt = wrer("writint prit table...\n");
  for kp in range(prit):
    prk = prit[kp];
    ib1k, ib2k, itk, pnamek = prk;
    wrot("%4d %4d  %4d %s\s" % (ib1k, ib2k, itk, pnamek));
  # ......................................................................

def make_all_CANC_terms(ns):
  # Returns the product-to-term table {prit} for
  # the CANC basis for a window with {ns} samples.
  # Currently only works for {ns==3}.
  #
  # Each entry of the table is a tuple {(ib1k, ib2k, itk, pnamek)}.
  # See {multifok_term_prod_t} in {multifok_term.h}.

  assert(ns == 9); # For the time being.
  
  nb = ns; # Number of basis elements.
  #           0      1      2      3      4      5      6      7      8 
  bname = ["Soo", "Som", "Sop", "Smo", "Spo", "Smm", "Spp", "Smp", "Spm"]; # Basis element names, indexed {0..nb-1}.
  assert(len(bname) == nb);
  
  # Dyadic permutation group generators, indexed {0..nb-1}: */
  # An up-down flip of the window sends sample {i} to sample {flip[i]}:
  flip = [0, 2, 1, 3, 4, 7, 8, 5, 6]; 
  # A CCW  rotation of the window sends sample {i} to sample {turn[i]}.
  turn = [0, 4, 3, 1, 2, 8, 7, 5, 6]; 
  
  nt = 0; # Number of terms found so far.
  prit = [].copy();  # The product-to-term table, indexed {0..np-1}.

  for ib1 in range(nb):
    for ib2 in range(ib1+1):
      # Check if this product dyadic-equivalent to some other product:
      it = nt;
      for jp in range(prit):
        prj = prit[jp];
        jb1, jb2, jt, pnamej = prj;
        if equivalent_prod(ib1,ib2,jb1,jb2,flip,turn):
          # Copy the term index:
          it = jt;
          break
      if it == nt:
        # New term:
        nt++;
    pnamei = (bname[ib1] + "*" + bname[ib2]);
    prit.append((ib1, ib2, it, pnamei));
  # Append the constant term {1}:
  prit.append((-1, -1, nt, "1"));
  nt = nt+1;
  return prit;
  # ......................................................................

def equivalent_prod(ib1,ib2,jb1,jb2,flip,turn):
  # Checks if the product of basis elements {ib1,ib2} 
  # is equivalent to {jb1,jb2} by dyadic operations and 
  # commutativity.
  assert ib1 <= ib2, "bug 1";
  nb = len(flip);
  assert nb == len(turn);
  for kt in range(4):
    for kf in range(2):
      assert jb1 >= 0 and jb1 < nb, "bug 2"
      assert jb2 >= 0 and jb2 < nb, "bug 3"
      # Apply communtativity:
      if jb1 > jb2: 
        jb2, jb1 = jb1, jb2
      if ib1 == ib2 and jb1 == jb2:
        return True;
      # Apply vertical flip:
      jb1 = flip[jb1];
      jb2 = flip[jb2];
    # Apply CCW turn:
    jb1 = turn[jb1];
    jb2 = turn[jb2];
  return False;
  # ......................................................................

def wrot(s):
  sys.stdout.write(s);
  # ......................................................................
  
def wrer(s)
  sys.stderr.write(s);
  # ......................................................................

main()