#! /usr/bin/gawk -f
# Last edited on 2002-02-25 00:11:53 by stolfi

BEGIN{
  abort = -1;
  usage = ( ARGV[0] " [-v avg=1 | -v dif=1] [ -v width=WIDTH ] < INFILE > OUTFILE" );
  
  # Reads a file with one line per recipe, in the format
  # RXNUM RXID RXVAL  where RXVAL is any additive attribute of the recipe.
  # Writes a file in the same format where the RXVAL is 
  # computed from the RXVAL's of WIDTH consecutive
  # recipes, centered on the nominal one. The output RXVAL may be fractional.
  # Treats lines with RXNUM = 000 as breaks in the sequence.
  # 
  # If "avg" is true, the output RXVAL is simply the weighted average of the 
  # RXVAL's in the window.  If "dif" is true, the output RXVAL is 
  # a half-width filtered sequence divided by the full-width one.
  
  if (width == "") { width = 1; }
  if (avg == "") { avg = 0; }
  if (dif == "") { dif = 0; }
  if ((avg && dif) || (! (avg || dif))) 
    { arg_error("must set exactly one of \"avg\" and \"dif\""); }
  # Buffer:
  # Last nbuf lines are in buf[ibuf-nbuf..ibuf-1];
  # Indices are computed modulo `width'.
  split("", bufrn); # Recipe numbers
  split("", bufid); # Recipe ids
  split("", bufrv); # Recipe sizes (tokens)
  nbuf = 0;
  ibuf = 0;
  # Average weight tables:
  split("", avg_wt);
  avg_nw = width;
  compute_weights(avg_wt, avg_nw);
  if (dif)
    { if ((width % 4) != 1) 
        { arg_error("\"dif\" option requires width = 4*k + 1"); }
      split("", dif_wt); 
      dif_nw = (width-1)/2 + 1;
      compute_weights(dif_wt, dif_nw);
    }
}

(abort >= 0) { exit abort; }

/^ *([#]|$)/ { next; }

/^000+[ ]/ { nbuf = 0; print; next; }

/./{ 
  # Append new entry to buffer
  bufrn[ibuf] = $1;
  bufid[ibuf] = $2;
  bufrv[ibuf] = $3;
  ibuf = (ibuf + 1) % width;
  if (nbuf < width) { nbuf++; }
  # Compute running average, if enough entries
  if (nbuf == width)
    { avgrv = weighted_sum(avg_wt, avg_nw);
      if (avg) 
        { outrv = avgrv; }
      else if (dif)
        { 
          midrv = weighted_sum(dif_wt, dif_nw);
          outrv = (avgrv == 0 ? 0 : midrv/avgrv);
        }
      else
        { prog_error("no output option"); }
      jmid = (ibuf + 1 + int((width-1)/2)) % width;
      printf "%s %s %6.2f\n", bufrn[jmid], bufid[jmid], outrv;
    }
  next;
}

END { }

function compute_weights(wt,nw,  j,k)
{
  wt[0] = 1;
  for (k = 1; k < nw; k++)
    { wt[k] = 0;
      for (j = k; j >= 0; j--)
        { wt[j] = (wt[j] + (j > 0 ? wt[j-1] : 0))/2; }
    }
}

function weighted_sum(wt,nw,   jdel,j,k,sum)
{
  jdel = int((width - nw)/2);
  sum = 0;
  for (k = 0; k < nw; k++)
    { j = (ibuf + 1 + jdel + k) % width;
      sum += wt[k] * bufrv[j];
    }
  return sum;
}

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