#! /usr/bin/gawk -f
# Last edited on 2014-05-09 19:51:00 by stolfilocal

BEGIN {
  abort = -1;
  Qmin = get_num_arg(Qmin, "Qmin", 0.0001, 9.9999);
  Qmax = get_num_arg(Qmax, "Qmax", Qmin, 9.9999);
  Qstep = get_num_arg(Qstep, "Qstep", 0.0001, 9.9999);

  split("", D);  # Days since arbitrary date
  split("", P);  # Prices
  split("", W);  # Weights
  np = 0;  # Number of data points.
  debug = 1;
}

(abort >= 0) { exit abort; }

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

/^[ ]*[0-9]+[ ]+[-+]?[0-9]*([.][0-9]*|[0-9])[ ]+[-+]?[0-9]*([.][0-9]*|[0-9])*[ ]*$/ {
  D[np] = $1;
  P[np] = $2;
  W[np] = $3;
  np++
  next;
}

// { data_error(("bad line format")); }

END {
  if (abort >= 0) { exit abort; }
  tmp = ( "/tmp/" PROCINFO["pid"] );
  dfile = ( tmp ".dat" ); # Input data file for {linear_fit}.
  efile = ( tmp ".exp" ); # Output coefficient file from {linear_fit}.
  Qbest = 0;
  Abest = 0;
  Bbest = 0;
  Sbest = 1.0e100;
  for (Q = Qmin; Q <= Qmax + 1e-6*Qstep; Q += Qstep) { 
    # Write data with exponential term using current {Q}:
    if (debug) { printf "=== trying Q = %6.4f ===\n", Q > "/dev/stderr"; }
    system(("rm -f " dfile));
    if (debug) { printf "deleted %s\n", dfile > "/dev/stderr"; }
    for (ip = 0; ip < np; ip++)
      { X = exp(log(Q)*(D[ip] - D[0]));
        printf "%5d %8.2f %8.6f  %3.1f %8.6f\n", D[ip], P[ip], W[ip], 1.0, X > dfile;
      }
    close(dfile)
    # Call {linear_fit} on that data:
    status = system(("linear_fit -weighted 1 -terms 2 -writeFormula " efile " < " dfile));
    if (status != 0) { end_error(("linear_fit failed, status = " status)); }
    # Recover the coefficients:
    split("", C)
    linear_fit_read_coeffs(efile, 2, C);
    M = C["Avg"];
    S = C["Dev"];
    if (S < Sbest) { Qbest = Q; Sbest = S; Abest = C[0]; Bbest = C[1]; }
    if (debug) { printf "Q = %6.4f S = %12.5f Qbest = %6.4f\n", Q, S, Qbest > "/dev/stderr"; }
  }

  # Output result
  printf "best fit: A = %8.2f B = %+8.2f Q = %6.4f\n", Abest, Bbest, Qbest;
  system(("rm -f " efile));
  system(("rm -f " dfile));
}

function linear_fit_read_coeffs(fname,nc,C,   ncr,nlin,lin,fld,nfld) {
  # Reads the "-writeFormula" output {fname} of {linear_fit}.  
  # It must have {nc} coefficients. 
  # Returns the coefficients in {C[0..nc-1]}, the average residual in {C["Avg"]} and 
  # the deviation of the residual in {C["Dev"]}.
  ncr=0;   # Num coeffs read.
  nlin=0;
  while((getline lin < fname) > 0) { 
    nlin++;
    # Remove spurious spaces and comments:
    gsub(/[\#].*$/, "", lin);
    lin = clean_spaces(lin);
    if (lin != "") {
      if (! match(lin, /^[-+ .0-9e]+$/)) { tbl_error(fname, nlin, ("bad line format = \"" lin "\"")); }
      nfld = split(lin, fld, " ");
      if (nlin == 1) { 
        # Number of coeffs:
        if ((nfld != 1) || (fld[1]+0 != nc)) { tbl_error(fname, nlin, ("bad num coeffs = \"" lin "\"")); }
      } else if (ncr < nc) {
        # Coefficient:
        if (nfld != 2) { tbl_error(fname, nlin, ("bad coeff line = \"" lin "\"")); }
        if (fld[1]+0 != ncr) { tbl_error(fname, nlin, ("wrong coeff index = \"" lin "\" - expected " ncr)); }
        C[ncr] = fld[2]+0;
        ncr++;
      } else { 
        # Average or deviation of residual:
        if (nfld != 1) { tbl_error(fname, nlin, ("bad Avg/Dev line format = \"" lin "\"")); }
        if (ncr == nc) { C["Avg"] = fld[1]+0; }
        else if (ncr == nc+1) { C["Dev"] = fld[1]+0; }
        else { tbl_error(fname, nlin, ("spurious line = \"" lin "\"")); }
        ncr++;
      }
    }
  }
  if ((ERRNO != "0") && (ERRNO != "")) { tbl_error(fname, nlin, ERRNO); }
  close (fname);
  if (nlin == 0) { arg_error(("file \"" fname "\" empty or missing")); }
  # printf "loaded %6d map pairs\n", ntbl > "/dev/stderr"
}

function get_num_arg(x,name,xmin,xmax) {
  # Checks whether {x} is defined, numeric (but not "%e"), and in {[xmin _ xmax]}:
  if (x == "") { arg_error(("must define {" name "}")); }
  x = clean_spaces(x);
  if (! match(x, /^[-+]?[0-9]*([.][0-9]*|[0-9])$/)) { arg_error(("argument {" name "} must be numeric")); }
  x = x + 0;
  if ((x < xmin) || (x > xmax)) { arg_error(("argument {" name "} must be in [" xmin " _ " xmax "]")); }
  return x;
}

function clean_spaces(x) {
  gsub(/[\011]/, " ", x);
  gsub(/^[  ]+/, "", x);
  gsub(/[  ]+$/, "", x);
  gsub(/[  ]+/, " ", x);
  return x;
}

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

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

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

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