#! /usr/bin/gawk -f
# Last edited on 2015-03-05 17:51:02 by stolfilocal

BEGIN \
  {
    abort = -1;
    
    # Caller must define (with "-v") variables {inFile,parmsFile},
    # set the environment variable "TZ" to "UTC".
    
    # Reads a smoothed price data file from {inFile}.
    # Reads a list of bubble parameters from {parmsFile}
    # (see {buf_read_bubble_parms_file} in "bubble_parms_file_functions.gawk").
    # Writes to {stdout} a decomposition of the price into bubbles
    # according to the additive bubble model in the format
    #
    # "{DATE[i]} {TIME} | {PISM[i]} | {PMOD[i]} | {PQUO[i]} | {PBUB[i,k]}"
    # 
    # where {DATE[i]} is a date like "2014-01-20", {TIME} is "00:00:00",
    # {PISM[i]} is the input smoothed price, {PMOD[i]} is the total 
    # of all bubble components (the bubble model), {PQUO[i]}
    # is the ratio {PMOD[i]/PISM[i]},
    # and each {PBUB[i,k]} is the component of the model price
    # due to bubble {k}.
  
    if (ENVIRON["TZ"] != "UTC") { arg_error(("must set TZ to 'UTC'")); }
    if (parmsFile == "") { arg_error(("must define {parmsFile}")); }
    if (inFile == "") { arg_error(("must define {inFile}")); }
    if (hrad == "") { arg_error(("must define {hrad}")); }
    if (hrad !~ /^[0-9]+$/) { arg_error(("invalid {hrad} = \"" hrad "\"")); }
    hrad += 0; # Make {hrad} numeric.
    
    pi = 3.1415926;

    ndays = 0;     # Number of days in input file.
    dy_first = ""; # First date seen.
    dy_last = "";  # Last date seen.
    
    ulp_price = 0.00001; # Unit in last place for output prices.
    
    # Input price data, indexed with {id} in {0..ndays-1}
    split("", date_id); # Dates in file.
    split("", pism_id);  # Smoothed price.
    ndays = read_smoothed_prices(inFile,date_id,pism_id);
    
    nbubs = -1; # Number of bubbles (from params file)
    
    # Bubble parameters (indexed 0 to {nbubs}):
    buf_initialize_bubble_parms_tables();
    
    # Computed bubble data, indexed with {id} in {0..ndays-1}, {kb} in ,0..nbubs-1}: 
    split("", pbub_id_kb); # Bubble component of price.

    nbubs = buf_read_bubble_parms_file(parmsFile,diup_kb,rtup_kb,dfup_kb,didn_kb,rtdn_kb,dfdn_kb,btag_kb,color_kb); 
  
    adjust_bubbles(ndays,date_id,pism_id,hrad, nbubs,diup_kb,rtup_kb,dfup_kb,didn_kb,rtdn_kb,dfdn_kb,btag_kb, pmid_kb, pbub_id_kb);
    
    output_bubbles(ndays,date_id,pism_id,hrad, nbubs,btag_kb,  pbub_id_kb);

    fflush("/dev/stdout");
    
    exit(0);
  }
  
(abort >= 0) { exit(abort); }

function read_smoothed_prices( \
    inFile,date_id,pism_id, \
    kf,date_dy,vbt_dy,vcr_dy,pmd_dy,rlodate,rhidate, \
    ndays,dy,id \
  )
  {
    # Reads from {inFile} the smoothed price,
    # returns in {
    
    # Input date, smoothed price, and volume tables:
    spf_initialize_smoothed_price_tables();
    
    # Clipping date interval (all data):
    rlodate="2009-01-01";
    rhidate="2999-12-31"; 
    kf = 0; # Index of price data file.
    spf_read_smoothed_price_file(\
      inFile,kf,rlodate,rhidate, \
      date_dy,vbt_dy,vcr_dy,pmd_dy \
    );
    
    # Get a list {date_id} of the dates, in order:
    ndays = asorti(date_dy,date_id); # Note, indexed {1..ndays}
    
    # Reindex by day number {id},  {0..ndays-1}:
    for (id = 0; id < ndays; id++)
      { dy = date_id[id+1];
        if (date_dy[dy] != 1) { prog_error(("inconsistent {date_dy,date_id}")); }
        date_id[id] = dy;
        pism_id[id] = pmd_dy[dy,kf];
      }
    delete date_id[ndays];
    
    # Clear auxiliary tables:
    spf_initialize_smoothed_price_tables();
    
    return ndays;
  }
  
function adjust_bubbles( \
    ndays,date_id,pism_id,hrad, nbubs,diup_kb,rtup_kb,dfup_kb,didn_kb,rtdn_kb,dfdn_kb,btag_kb, pmid_kb,\
    pbub_id_kb, \
    id,kb,dy,pini,pfin,ppro,ppsm,diup,rtup,dfup,didn,rtdn,dfdn,pres_id,ppro_id,ppsm_id,idiup,idfup,ididn,idfdn, \
    log_rtup,log_rtdn,wini,wfin,arg,wdot,s_ppro_ppro,s_ppro_pres,pmid \
  )
  {
    # Adjusts bubble component peak values {pmid_kb} to smoothed prices {pism_id}.
    # Assumes given the bubble parameters {nbubs,diup_kb,rtup_kb,rtdn_kb,dfdn_kb,dfup|didn_kb}.
    # The {date_id} is used for error reporting.
    # Saves bubble values to {pbub_id_kb[0..ndays-1,0..nbubs-1]}.
    
    split("", pres_id);    # Residual component of price ({pism_id} minus sum of prev bubbles).
    split("", ppro_id);    # Prototype of next bubble.
    split("", ppsm_id);    # Smoothed prototype.
    
    # Set {pres_id} to the input prices:
    for (id = 0; id < ndays; id++) { pres_id[id] = pism_id[id]; }
    
    # Adjust each bubble:
    for (kb = 0; kb < nbubs; kb++)
      { # Adjust the bubble {kb} to the residual {pres_id}:
        printf "adjusting bubble %d (%s) ...\n", kb, btag_kb[kb] > "/dev/stderr";
        
        # Bubble parameters:
        diup = diup_kb[kb];
        rtup = rtup_kb[kb];
        dfup = dfup_kb[kb];
        didn = didn_kb[kb];
        rtdn = rtdn_kb[kb];
        dfdn = dfdn_kb[kb];
        printf "  parms = %s %8.5f %s %s %8.5f %s\n", diup, rtup, dfup, didn, rtdn, dfdn > "/dev/stderr";
        
        log_rtup = log(rtup);
        log_rtdn = log(rtdn);
        
        # Find the day indices {idiup,idfup,ididn,idfdn} of the critical dates:
        for (id = 0; id < ndays; id++)
          { if (date_id[id] == diup) { idiup = id; }
            if (date_id[id] == dfup) { idfup = id; }
            if (date_id[id] == didn) { ididn = id; }
            if (date_id[id] == dfdn) { idfdn = id; }
          }
        
        # Build the bubble prototype {ppro_id}:
        for (id = 0; id < ndays; id++)
          { if (id <  idfup)
              { ppro = exp(log_rtup*(id - idfup)); }
            else if (id > ididn)
              { ppro = exp(log_rtdn*(id - ididn)); }
            else
              { ppro = 1.0; }
            ppro_id[id] = ppro;
          }
        
        # Smooth the prototype like the input price, yelding {ppsm_id}:
        spf_smooth_series(0, ndays-1, ppro_id, hrad, ppsm_id);
        
        # Compute scalar products:
        s_ppsm_pres = 0;
        s_ppsm_ppsm = 0;
        for (id = 0; id < ndays; id++)
          { dy = date_id[id];
            ppsm = ppsm_id[id];
            if ((id >= idiup) && (id <= idfdn))
              { # Accumulate the scalar products:
                arg = pi*(id - idiup + 0.5)/(idfdn - idiup + 1);
                wdot = 0.5*(1 - cos(arg));
                s_ppsm_pres += wdot*ppsm*pres_id[id];
                s_ppsm_ppsm += wdot*ppsm*ppsm;
              }
          }
          
        # Compute the scale factor:
        pmid = s_ppsm_pres/s_ppsm_ppsm;
        printf "  s_ppsm_pres = %.5f s_ppsm_ppsm = %.5f\n", s_ppsm_pres, s_ppsm_ppsm > "/dev/stderr";
        printf "  pmid = %.5f\n", pmid > "/dev/stderr";
        pmid_kb[kb] = pmid;
        
        # Scale the prototype, save in {pbub_id_kb}, adjust {pres_id}:
        for (id = 0; id < ndays; id++)
          { dy = date_id[id];
            ppsm = ppsm_id[id];
            pbub = pmid*ppsm;
            # if ((id >= idfup|didn - 5) && (id >= idfup|didn + 5))
            #   { printf "  %s %18.5f %10.7f %18.5f\n", dy, pres_id[id], ppsm, pbub > "/dev/stderr"; }
            # Subtract from residual before rounding:
            pres_id[id] -= pbub;
            # Round it and save:
            pbub = (pbub == 0 ? 0.0 : usf_round_value(pbub,ulp_price));
            pbub_id_kb[id,kb] = pbub;
          }
      }
  }

function output_bubbles( \
    ndays,date_id,pism_id, hrad, \
    nbubs,btag_kb, pbub_id_kb, \
    id,kb,dy,tm,ptot_id,s_pbub,pquo \
  )
  {
    
    printf "writing to stdout...\n" > "/dev/stderr";
    # Sets {ptot_id[0..ndays-1]} to the sum of all bubbles:
    split("", ptot_id);
    for (id = 0; id < ndays; id++)
      { s_pbub = 0;
        for (kb = 0; kb < nbubs; kb++) { s_pbub += pbub_id_kb[id,kb]; }
        ptot_id[id] = (s_pbub == 0 ? 0.0 : usf_round_value(s_pbub,ulp_price));
      }
    
    for (id = 0; id < ndays; id++)
      { # Output data for day {id}
        dy = date_id[id];
        tm = "00:00:00";
        printf "%s %s", dy, tm;
        printf " | %.5f | %.5f", pism_id[id], ptot_id[id];
        pquo = (ptot_id[id] == 0 ? 0.0 : pism_id[id]/ptot_id[id]);
        printf " | %.6f", pquo;
        for (kb = 0; kb < nbubs; kb++)
          { printf " | %.5f", pbub_id_kb[id,kb]; }
        printf "\n";
      }
    fflush("/dev/stdout");
  }