#! /usr/bin/gawk -f 
# Last edited on 2002-02-24 12:49:46 by stolfi

BEGIN {
  abort = -1;
  usage = ( ARGV[0] " \\\n" \
    "  -v xFile=XFILE [-v xField=XFIELDNUM] \\\n" \
    "  -v yFile=YFILE [-v yField=YFIELDNUM] \\\n" \
    "  [-v absErr=ABSERR] [-v relErr=RELERR] \\\n" \
    "  > OUT.pgm" \
  );
  
  # Reads two files XFILE and YFILE containing two numeric 
  # sequences xv[0..nx-1] and yv[0..ny-1]. 
  # Writes a PGM image where the pixel (i,j) is the probability 
  # of xv[i] = yv[j] (counting (0,0) from the bottom left).
  
  if (xFile == "") { arg_error("must define \"xFile\""); }
  if (xField == "") { xField = 1; }
  split("", xv);
  nx = read_seq(xFile,xv,xField);
  
  if (yFile == "") { arg_error("must define \"yFile\""); }
  if (yField == "") { yField = 1; }
  split("", yv);
  ny = read_seq(yFile,yv,yField);
  
  # Defaults should be good for integers with 10% acuracy:
  if (absErr == "") { absErr = 0.5 }
  if (relErr == "") { relErr = 0.1 }

  write_image(xv,nx, yv,ny,yField)
}

function read_seq(fname,z,field,   nz,i,nlin,lin,fld,nfld)
{
  nz=0;
  nlin=0;
  while((getline lin < fname) > 0) { 
    nlin++;
    if (! match(lin, /^[ \011]*([#]|$)/))
      { gsub(/[#].*$/, "", lin);
        nfld = split(lin, fld, " ");
        if (nfld < field) 
          { seq_error(fname, nlin, ("bad sequence file entry = \"" lin "\"")); }
        if (fld[field] !~ /^[-+]?[0-9]*([.][0-9]*)?/) 
          { seq_error(fname, nlin, ("bad sequence value = \"" lin "\"")); }
        z[nz] = fld[field];
        nz++;
      }
  }
  if (ERRNO != "0") { seq_error(fname, nlin, ERRNO); }
  close (fname);
  if (nlin == 0) { arg_error(("file \"" fname "\" empty or missing")); }
  # printf "loaded %6d map pairs\n", nz > "/dev/stderr"
  return nz;
}
  
function write_image(xv,nx,yv,ny,  x,y,r,c,v)
{
  write_pgm_header(nx,ny,255);
  for (r = ny-1; r >= 0; r--)
    { for (c = 0; c < nx; c++)
        { x = xv[c]; y = yv[r];
          v = ((x == 0) || (y == 0) ? 1 : prob_eq(x,y,absErr,relErr));
          printf "%3d%s", int(255*v + 0.5), (c == nx-1 ? "\n" : " "); 
        }
    }
}

function write_pgm_header(nx,ny,maxval)
{
  printf "P2\n%d %d\n%d\n", nx, ny, maxval;
}

function prob_eq(x,y,  d2,s2,e2,chi2,prob)
{
  # Temporary definition
  s2 = (x*x + y*y);
  d2 = (x - y)*(x - y);
  e2 = s2*relErr*relErr + absErr*absErr;
  chi2 = d2/e2/2;
  prob = (chi2 > 50 ? 0 : exp(-chi2));
  # printf "%s %s %5.3f\n", x, y, prob > "/dev/stderr";
  return prob;
}

function arg_error(msg)
{ 
  printf "%s\n", msg >> "/dev/stderr";
  abort = 1;
  exit 1
}

function seq_error(f,n,msg)
{ 
  printf "file %s, line %d: %s\n", f, n, msg > "/dev/stderr";
  abort = 1;
  exit 1
}

function data_error(msg)
{ 
  printf "line %d: %s\n", FNR, msg > "/dev/stderr";
  abort = 1;
  exit 1
}