#! /usr/bin/gawk -f
# Last edited on 2010-10-07 12:03:17 by stolfilocal

# Applies an orthonormal transformation to a discrete signal.

BEGIN{
  if (n == "") { arg_error("must define the dimension {n}"); }
  if (map == "") { arg_error("must define the matrix filename {map}"); }
  if (transp == "") { arg_error("must define the transpose flag {transp}"); }
  pi = 3.1415926;
  
  # Read the basis matrix {phi[0..n-1,0..n-1]}:
  # Each row is an element of the basis, assumed to be 
  # orthonormal under the inner product 
  #   { <f,g> = (1/n)*SUM{ f[k]*g[k] : k IN {0..n-1} } }
  # Thus the columns are the dual basis.
  
  split("", phi);
  read_matrix(map, phi, n, transp);
  check_basis(phi,n,n);
  
  # Read the vector {u[0..n-1]} to be transformed: 
  split("", u); 
  read_vec("/dev/stdin", u, n);
  
  # Compute the mapped vector {v = u*(phi')}:
  split("", v); 
  for (i = 0; i < n; i++) { v[i] = dot_u(phi, n, u, i); }

  # Write out:
  for (i = 0; i < n; i++) { printf "%4d %+9.6f\n", i, v[i]; }
  fflush();
}

function read_matrix(fname,phi,n,transp,   nrows,nline,raw,lin,fld,nfld,tmp,i,j)
  {
    nrows=0;
    nline=0;
    while((getline raw < fname) > 0) { 
      nline++;
      lin = raw;
      gsub(/[\011\015]/, " ", lin);
      gsub(/[ ]*[\#].*$/, "", lin);
      if (! match(lin, /^[ ]*$/))
        { nfld = split(lin, fld, " ");
          if (nfld != (n+1)) { file_error(fname, nline, raw, ("bad line format")); }
          i = nrows;
          if (nrows >= n) { file_error(fname, nline, raw, ("too many rows")); }
          # If {inv} is true, swap the two columns:
          if ((fld[1] + 0) != i) { file_error(fname, nline, raw, ("unexpected row")); }
          for (j = 0; j < n; j++) 
            { if (transp)
                { phi[j,i] = fld[j+2]; }
              else
                { phi[i,j] = fld[j+2]; }
            }
          nrows++;
        }
    }
    if (nrows != n) { file_error(fname, nline, "", ("wrong number of rows")); }
    if (ERRNO != "0") { file_error(fname, nline, "", ERRNO); }
    close (fname);
    if (nline == 0) { arg_error(("file \"" fname "\" empty or missing")); }
    # printf "loaded %6d matrix rows\n", nrows > "/dev/stderr"
  }

function read_vec(fname,u,n,   nrows,nline,raw,lin,fld,nfld,tmp,i)
  {
    nrows=0;
    nline=0;
    while((getline raw < fname) > 0) { 
      nline++;
      lin = raw;
      gsub(/[\011\015]/, " ", lin);
      gsub(/[ ]*[\#].*$/, "", lin);
      if (! match(lin, /^[ ]*$/))
        { nfld = split(lin, fld, " ");
          if (nfld != 2) { file_error(fname, nline, raw, ("bad line format")); }
          i = nrows;
          if (nrows >= n) { file_error(fname, nline, raw, ("too many elements")); }
          # If {inv} is true, swap the two columns:
          if ((fld[1] + 0) != i) { file_error(fname, nline, raw, ("unexpected element")); }
          u[i] = fld[2];
          nrows++;
        }
    }
    if (nrows != n) { file_error(fname, nline, "", ("wrong number of elements")); }
    if (ERRNO != "0") { file_error(fname, nline, "", ERRNO); }
    close (fname);
    if (nline == 0) { arg_error(("file \"" fname "\" empty or missing")); }
    # printf "loaded %6d vector elements\n", nrows > "/dev/stderr"
  }

function check_basis(phi,m,n,  M,i,j,bug,tmp)
  {
    # Checks whether the basis {phi[0..m-1]} with {n} columns is orthonormal.
    bug = 0;
    split("",M);  # Scalar product matrix:
    printf "Scalar product matrix:\n" > "/dev/stderr";
    for (i = 0; i < n; i++)
      { for (j = 0; j < n; j++) 
          { M[i,j] = dot(phi,n,i,j);
            tmp = (i == j);
            if (abs(M[i,j] - tmp) > 1.0e-6) 
              { printf "** <phi[%d],phi[%d]> = %+9.6f\n", M[i,j] > "/dev/stderr";
                bug=1;
              }
          }
      }

    if (bug)
      { 
        printf "Scalar product matrix:\n" > "/dev/stderr";
        for (i = 0; i < n; i++)
          { for (j = 0; j < n; j++) 
              { M[i,j] = dot(phi,n,i,j);
                printf " %+9.6f", M[i,j] > "/dev/stderr";
              }
            printf "\n" > "/dev/stderr";
          }
      }
  }

function dot(phi,n,i,j,   sum,k)
  { 
    # Dot product of rows {i,j} of the matrix {phi} with {n} columns.
    sum = 0;
    for (k = 0; k < n; k++) { sum += phi[i,k]*phi[j,k]; }
    return sum/n;
  }

function dot_u(phi,n,u,i,   sum,k)
  { 
    # Dot product of {u} and row {i} of the matrix {phi} with {n} columns.
    sum = 0;
    for (k = 0; k < n; k++) { sum += phi[i,k]*u[k]; }
    return sum/n;
  }

function abs(x)
  { 
    if (x < 0) { return -x; } else { return x; }
  }

function num_error(i,msg)
  { 
    printf "** element %d: %s\n", i, msg > "/dev/stderr";
    exit 1;
  }

function arg_error(msg)
  { 
    printf "** argument error: %s\n", msg > "/dev/stderr";
    exit 1;
  }

function file_error(fname,nline,raw,msg)
  { 
    printf "%s:%d: ** %s\n", fname, nline, msg > "/dev/stderr";
    if (raw != "") { printf "  line = «%s»\n", raw > "/dev/stderr"; }
    abort = 1;
    exit 1
  }