/* rtfactors.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          5/20/98
 *                *
 * COMPILE AS: gccemuo rtfactors
 * RUN AS: rtfactors <Name of output NC/binary file>  
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: rtfactors \n\
Program to calculate temporal correlation maps (Rt) between a Dependent Variable\n\
and a number of factors (Independent Variables).\n\
\n\
"};

/* local functions */
void Setup(void);
void Process(void);
void Output(void);

/* local variables */
float **depvarTsXYblock, **TsXYblock, ***out;
int errnc;

#define  NUMPARAM 3
#define  DEP      0
#define  FAC      1
#define  R        0


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Temporal Correlation Maps");
  AssignString(&pszNChist, "correlation maps created by rtfactors");

  Initialize(&argc, argv, DOMASK);

  Setup();
  Process();
  Output();

  CleanUp();

  return(0);
}


/* ================================================================= */
void Setup(void)
{
  int v, i;

  /* allocate memory and initialize out[][][] TO 0.0  */
  out = alloc3d_f(0,1-1,0,nVars-1,0,nRows*nCols-1);
  for (v=0; v < nVars;v++)
    for (i=0; i < nRows*nCols;i++)
      out[R][v][i] = 0.0;           R="#009900">/* correlation, R   */

  TsXYblock       = alloc2d_f(0,nTimeSteps-1,0,nBlockRows*nCols-1);
  depvarTsXYblock = alloc2d_f(0,nTimeSteps-1,0,nBlockRows*nCols-1);
}

/* ================================================================= */
/* ======== READ DATA AND CALCULATE CORRELATION COEFFICIENT ======== */
void Process(void)
{
  int v, i, iout, col, row, n, i0, nyblock;   int bug=0;
  int l1, l2, tvar = 0;
  float fill_val;
  double stats_i[2][NUMPARAM], cov_i;
  float ts, txyval, depvartxyval;


  for (v = 1; v < nVars; v++)
  {
    printf("\n v=%d *** \n",v);
    for (nyblock = 0; nyblock < nYBlocks; nyblock++)
    {
      ReadTXYbl(TsXYblock, pstVars[v], nyblock);
      ReadTXYbl(depvarTsXYblock, pstVars[DEP], nyblock);
      i0 = nyblock*nBlockRows*nCols;
      if (bDebug)  printf("\n      nyblock=%d, i0=%d \n", nyblock,i0);

      for (i = i0; i < i0+nBlockRows*nCols; i++)
      {
        iout = (bOutBinary) ? i : getRowRevCol(i, &col, &row);
        if (masktest(pMask[i])) {
          if (i%50000 == 0 && bDebug)  printf("\n pstVars[%d]=%s: i=%d  *  ", v,pstVars[v].param,i);

          for (l1=0; l1<2; l1++)
            for (l2=0; l2<NUMPARAM; l2++)
              stats_i[l1][l2]  = 0.0;

          for (tvar = 0; tvar < nTimeSteps;tvar++)
          {
            ts = TsXYblock[tvar][i-i0];       /* INDEPENDENT VAR */
            /* catch floating-point underflow errors w/ tolerance "TOL" */
            txyval   = (ts > TOL || ts < -TOL) ? ts : 0.0;
            ts = depvarTsXYblock[tvar][i-i0];     /* DEPENDENT VAR */
            depvartxyval = (ts > TOL || ts < -TOL) ? ts : 0.0;

            stats_i[FAC][SUM]    += txyval;
            stats_i[FAC][SUM2]   += txyval * txyval;
            stats_i[DEP][SUM]    += depvartxyval;
            stats_i[DEP][SUM2]   += depvartxyval * depvartxyval;
            stats_i[FAC][SUM2XY] += txyval * depvartxyval;
          }  /* END tvar FOR-LOOP (loop thru each time step) */

          stats_i[FAC][MEAN] = stats_i[FAC][SUM] / (double)nTimeSteps;
          stats_i[FAC][SD]   = sqrt(stats_i[FAC][SUM2]/(double)nTimeSteps - DSQR(stats_i[FAC][MEAN]));
          stats_i[DEP][MEAN] = stats_i[DEP][SUM] / (double)nTimeSteps;
          stats_i[DEP][SD]   = sqrt(stats_i[DEP][SUM2]/(double)nTimeSteps - DSQR(stats_i[DEP][MEAN]));
          cov_i  = stats_i[FAC][SUM2XY]/(double)nTimeSteps - stats_i[FAC][MEAN]*stats_i[DEP][MEAN];

          out[R][v][iout] = (float) cov_i/(stats_i[FAC][SD] * stats_i[DEP][SD]);
        } /* END masktest() IF-BLOCK */
        else
        {         /* outside the pMask: FillValue */
          if (pstVars[v].tEmuType == EMU_SHORT)
            fill_val = (float) pstVars[v].nFillValue;
          else if (pstVars[v].tEmuType == EMU_FLOAT)
            fill_val = pstVars[v].fFillValue;

          out[R][v][iout] = fill_val;
        }

        R="#FF0000">if (IsBad((R="#FF0000">double)out[R][v][iout]))
          out[R][v][iout] = fill_val;   R="#009900">/* BUT FILL_VAL MAY NOT HAVE BEEN READ, IF THE ELSE-BLOCK
                                        ABOVE WASN'T EXECUTED !!! */

      } /* END i-loop (read each pixel within nyblock) */
    } /* END nyblock-loop (loop thru each nyblock) */
  } /* END v FOR-LOOP (loop thru each variable) */
}

/* ================================================================= */
/* ==============  WRITE STATS TO FILE(S)  ========================= */
void Output(void)
{
  int v;
  int ncID;               /* netCDF file ID */
  int tID;                /* time axis ID   */
  int varID;
  char varName[NAMELEN];      /* Name of NC variable(s) to be written      */
  char varLongName[NAMELEN];  /* Long Name of NC variable(s) to be written */
  FILE *fout;


  if (bOutBinary)    /* write output to binary files */
    for (v = 1; v < nVars; v++)
    {
      fout = OUTPUTfile(pszOutputFile, ".correl");
      fwrite(out[R][v], sizeof(float), nRows*nCols, fout);
      fclose(fout);
    }
  else
  {         /* write output to a single netCDF file */
    /* setup new NC file from NC headers */
    SetupDefaultAxes(TRUE, FALSE);
    ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
    for (v = 1; v < nVars; v++)
    {
      sprintf(varName,"%s_Rt", pstVars[v].param);
      sprintf(varLongName,"Rt: %s vs. %s", pstVars[DEP].param,pstVars[v].param);
      varID = DefNCvar(ncID, varName, EMU_FLOAT, "xy", varLongName, "");
      errnc = nc_put_var_float(ncID, varID, out[R][v]);
    }

    errnc = nc_close(ncID);
  }
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}

/* ================================================================= */
void LocalCleanUp(void)
{
  free3d_f(out,0,1-1,0,nVars-1,0,nRows*nCols-1);
  free2d_f(TsXYblock,0,nTimeSteps-1,0,nBlockRows*nCols-1);
  free2d_f(depvarTsXYblock,0,nTimeSteps-1,0,nBlockRows*nCols-1);
}