# Last edited on 2000-01-30 18:42:44 by stolfi

# To be included in gawk scripts

function fctb_factor(a,u,v,eps,   i,j,n,m,ea,la,sla)
  {
    # Factors the nonnegative matrix a[i,j] into a product
    # os nonnegative matrices u[i]*d[i,j]*v[j] 
    # such that the product of the elements along any row
    # or column of "d" is 1.

    m = 0; for (i in u) { u[i] = 0; m++; }
    n = 0; for (j in v) { v[j] = 0; n++; }
    sla = 0;
    for (i in u) 
      { for (j in v) { 
          ea = a[i,j];
          la = log(sqrt(ea*ea + eps*eps)); 
          sla += la; 
          u[i] += la;
          v[j] += la;
        }
      }
    for (i in u) { u[i] = exp((u[i] - sla/2/m)/n); }
    for (j in v) { v[j] = exp((v[j] - sla/2/n)/m); }

    # fctb_print(a,u,v,eps);
  }
        
function fctb_print(a,u,v,eps,   i,j,k,eij,aij,dij,bij,si,sj,ti,tj)
  {
    for (k=1; k<=2; k++)
      { 
        printf "--- %s -------------------------------\n", 
          ( k == 1 ? "anomaly matrix" : "tensor approximation" );

        printf "    ";
        for (j in v) { printf " %7d", j; }
        printf "   %7s", "tot";
        if (k == 1) { printf " %7s", "u_i"; }
        printf "\n";

        printf "    ";
        for (j in v) { printf " %7.7s", "-------------"; }
        printf "   %7.7s", "-------------";
        if (k == 1) { printf " %7.7s", "-------------"; }
        printf "\n";

        split("", si); for (i in u) { si[i] = 0; }
        split("", sj); for (j in v) { sj[j] = 0; }

        for (i in u)
          { printf "%2d |", i;
            for (j in v) 
              { 
                aij = a[i,j];
                bij = u[i]*v[j];
                dij = sqrt(aij*aij + eps*eps)/bij;
                eij = ( k == 1 ? dij : bij );
                si[i] += eij; sj[j] += eij;
                printf " %7s", fctb_format(eij,7,3);
              }
            printf " | %7s", fctb_format(si[i],7,3);
            if (k == 1) { printf " %7s", fctb_format(u[i],7,3); }
            printf "\n";
          }

        printf "    ";
        for (j in v) { printf " %7.7s", "-------------"; }
        printf "   %7.7s", "-------------";
        if (k == 1) { printf " %7.7s", "-------------"; }
        printf "\n";

        printf "tot ";
        for (j in v) { 
          printf " %7s", fctb_format(sj[j],7,3);
        }
        printf "\n";
        
        if (k==1)
          { printf "v_j ";
            for (j in v) { 
              printf " %7s", fctb_format(v[j],7,3);
            }
            printf "\n";
          }
        
        fflush(stdout);
      }
    printf "\n";
  }

function fctb_format(e,w,p,  f)
  {
     f = sprintf("%*.*f", w,p,e);
     if (length(f) > w)
       { f = substr(f, length(f)-w+1, w);
         gsub(/[0-9]/, "*", f);
       }
     else if (e == 0)
       { gsub(/[0]/, " ", f); }
     return(f);
  }