#! /usr/bin/gawk -f
# Last edited on 2014-02-06 00:38:59 by stolfilocal

# Reads a file with raw increment data as written by
# {compute_future_uncertainty.gawk}.
# 
# The client must load (with "-f") the library {write_fni_image.gawk}.
# It must also define (with "-v") the variables
#
#   {rmax}  max time step expected in input.
#   {nh}    number of {d} histogram bins on each side of zero.
#   {dzmax} max absolute increment to include in histogram.
#
# From each input line reads "{i} {i+r} {r} {duvar[i,r]} {dtime[i,r]} {dvbtc[i,r]} {dvnat[i,r]} {dzvar[i,r]}"
# where 
#
#   {i}          the index of a sampling time.
#   {r}          an time index step.
#   {duvar[i,r]} the increment in the independent variable {u} between mids of intervals {i} and {i+r}. 
#   {dtime[i,r]} the difference between times at mids of intervals {i} and {i+r}. 
#   {dvbtc[i,r]} the accumulated BTC volume between mids of intervals {i} and {i+r}. 
#   {dvnat[i,r]} the accumulated national currency volume between mids of intervals {i} and {i+r}. 
#   {dzvar[i,r]} the log-price increment from intervals {i} to {i+r}, that is, {zvar[i+r] - zvar[i]}.
# 
# Writes to standard output a ".fni" float image where the horizontal axis is the time span {r}
# and the vertical axis is the increment {dzvar} in {log10(price)} over that time span, divided by 
# {sqrt(r)} and quantized into {2*nh+1} intervals. 

BEGIN \
  { 
    if (rmax == "") { arg_error(("must define {rmax}")); }
    rmax += 0; 
    if (rmax < 1) { arg_error(("invalid {rmax}")); } 
    
    if (nh == "") { arg_error(("must define {nh}")); }
    nh += 0;
    if (nh < 1) { arg_error(("invalid {nh}")); } 

    if (dzmax == "") { arg_error(("must define {dzmax}")); }
    dzmax += 0;
    if (dzmax <= 0) { arg_error(("invalid {dzmax}")); } 

    nbins = 2*nh + 1;              # Total number of bins in histogram.
    ppbin = int(rmax/nbins + 1.0); # Histogram bin width in pixels.
    
    NC = 1;           # Image channels.
    NX = rmax;        # Image columns {0..NX-1} correspond to steps {1..rmax}. 
    NY = nbins*ppbin; # Image rows {0..NY-1} corresponds to bins {-nh..+nh} of the histogram.
    
    n = 0;          # Number of entries read.
    tav_out = 0;    # Total volume for data points that fell out of the histogram.
    split("", hav); # Element {hav[r,-nh..+nh]} is the total volume for step {r} and each {dz} bin.
    split("", tav); # Element {tav[r]} is the total volume for step {r}.
    for (r = 1; r <= rmax; r++)
      { tav[r] = 0;
        for (k = -nh; k <= +nh; k++) { hav[r,k] = 0; }
      }
  } 

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

/^[ ]*[0-9]+[ ]+[0-9]+[ ]+[0-9]+[ ]/ \
  { 
    if (NF != 8) { data_error(("wrong {NF}")); }
    i = 0+$1; # Start time index.
    j = 0+$2; # Final time index.
    r = 0+$3; # Index difference.
    duvar = 0+$4; # Increment in independent variable {u}.
    dtime = 0+$5; # Increment in time.
    dvbtc = 0+$6; # Increment in BTC trade volume.
    dvnat = 0+$7; # Increment in national currency volume.
    dzvar = 0+$8; # Log price increment.
    
    if ((r < 1) || (r > rmax)) { data_error(("bad time step {r} " r)); }
    if ((i < 0) || (j != i + r)) { data_error(("bad time indices {i,j} " i " " j)); }
    if (duvar <= 0) { data_error(("bad indep val increment {duvar} " duvar)); } 
    if (dtime <= 0) { data_error(("bad time increment {dtime} " dtime)); } 
    if (dvbtc <= 0) { data_error(("bad BTC trade volume {dvbtc} " dvbtc)); } 
    if (dvnat <= 0) { data_error(("bad NAT trade volume {dvnat} " dvnat)); } 
   
    n++;
    
    # Compute the bin index {k}:
    ak = int(0.5 + (nh + 0.49999999)*fabs(dzvar)/sqrt(r)/dzmax);
    if (ak > nh) 
      { dvnat_out += dvnat; }
    else
      { if (dzvar < 0) { k = -ak; } else {k = +ak; }
        # Accumulate in histogram:
        # !!! Should antialias? !!!
        hav[r,k] += dvnat;
        tav[r] += dvnat;
      }
    next;
  }
  
// \
  { data_error("invalid line format"); }
  
END \
  { printf "read %d valid data points\n", n > "/dev/stderr";
    if (dvnat_out > 0)
      { printf "!! total volume outside histogram = %.7f\n", dvnat_out > "/dev/stderr"; }
    # Build image:
    split("", img);
    for (x = 0; x < NX; x++)
      { r = x + 1;
        if (r > rmax) { data_error(("boh?")); }
        for (k = -nh; k <= +nh; k++)
          { for (dy = 0; dy < ppbin; dy++)
              { y = (k + nh)*ppbin + dy;
                if ((y < 0) || (y >= NY)) { data_error(("beuh? k = " k " y = " y)); }
                img[0,x,y] = (tav[r] == 0 ? 0.0 : hav[r,k]/tav[r]);
              }
          }
      }
    write_fni_image(NC, NX, NY, img, "/dev/stdout");
  }
 
function fabs(x)
  { if (x < 0) { x = -x; }
    return x;
  }
         
function data_error(msg)
  { printf "%s:%s: ** %s\n", FILENAME, FNR, msg > "/dev/stderr"; 
    printf "  «%s»\n", $0 > "/dev/stderr"; 
    abort = 1;
    exit(abort);
  } 
          
function arg_error(msg)
  { printf "** %s\n", msg > "/dev/stderr"; 
    abort = 1;
    exit(abort);
  }