/* covarmatrix.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          4/22/98
 *                *
 * COMPILE AS: gccemu covarmatrix
 * RUN AS: covmatrix <parameter>
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: covarmatrix \n\
.\n\
\n\
Command-line switches specific to this program are:\n\
\n\
"};

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

/* local variables */
float **x;
double **sum2_xy, *sum_x, **covtt, *meant;
int tk, tj;
int errnc;

#define  LOADLIMIT  40


int main(int argc, char *argv[])
{
  Initialize(&argc, argv, DOMASK);

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

  CleanUp();

  return(0);
 }

/* ================================================================= */
void Setup(void)
{
  /* allocate memory space and initialize arrays */
  meant   = alloc1d_d(1,nTimeSteps);
  sum_x   = alloc1d_d(1,nTimeSteps);
  x       = alloc2d_f(1,2*LOADLIMIT,0,nRows*nCols-1);
  sum2_xy = alloc2d_d(1,nTimeSteps,1,nTimeSteps);
  covtt   = alloc2d_d(1,nTimeSteps,1,nTimeSteps);

  for (tk=1; tk <= nTimeSteps; tk++)
  {
    sum2_xy[tk][tk] = sum_x[tk] = 0.0;
    covtt[tk][tk]   = meant[tk] = 0.0;
    for (tj=tk+1; tj <= nTimeSteps; tj++)
    {
      sum2_xy[tk][tj] = sum2_xy[tj][tk] = 0.0;
      covtt[tk][tj]   = covtt[tj][tk]   = 0.0;
    }
  }
}

/* ================================================================= */
/* ======== READ DATA AND CALCULATE STATISTICAL PARAMETERS ========= */
void Process(void)
{
  int yr, mo, i;     int bug=0;
  int trow, tcol;
  int l, lcol, lrow, lrow_o, tx, loadtk, loadtj, lyoffset;
  float *y[LOADLIMIT+1];


  /* size of load block, in number of time steps (months) */
  tx = (int) floor((double)nTimeSteps/(double)LOADLIMIT);

  /* ====== Load Data, Calculate Sum2_XY And Sum_X Arrays ====== */
      /* the tk and tj indices will be base *1*, not *0* */
  for(tk=1; tk <= nTimeSteps; tk+=LOADLIMIT)
  {
    /* LOAD data for loadtk? block */
    loadtk = (tx*LOADLIMIT/tk >= 1) ? LOADLIMIT:(nTimeSteps % tk)+1;
    getyrmonth(tk-1, &yr, &mo);
    if (bDebug)  printf("\n %d/%d, tk=%d, loadtk=%d \n", mo+1,yr,tk,loadtk);
    for (l=1; l <= loadtk; l++)
    {
      ReadXY(x[l], pstVars[0], yr, mo);
      if (bDebug)  printf("  x[%d][10200]=%f  ", l,x[l][10200]);
      NEXTyrmo(&yr, &mo);
    }

    lyoffset = 0;
    for(tj=tk; tj <= nTimeSteps; tj+=LOADLIMIT)
    {
      loadtj = (tx*LOADLIMIT/tj >= 1) ? LOADLIMIT:(nTimeSteps % tj)+1;
      if (bDebug)  printf("\n IN tj LOOP, %d/%d, tj=%d, loadtj=%d \n", mo+1,yr,tj,loadtj);

      if (tj > tk)
      {
        for (l=LOADLIMIT+1; l <= LOADLIMIT+loadtj; l++)
        {
          ReadXY(x[l], pstVars[0], yr, mo);
          NEXTyrmo(&yr, &mo);
          if (bDebug)  printf("\n   IN if-block, %d > %d, l=%d   ", tj,tk,l);
        }
        /* lyoffset = loadtj; */
        lyoffset = LOADLIMIT;
      }   /* END (tj > tk) if-block */

      if (bDebug)  printf("\n IN tj LOOP %d/%d, > 2nd readbin block; lyoffset=%d \n", mo+1,yr,lyoffset);

      for (l=1; l <= loadtj; l++)
        y[l] = x[l+lyoffset];


      /* === Calculate sum2_xy and sum_x arrays for loadtk x loadtj block === */
      for (lcol=1; lcol<=loadtk; lcol++)
      {
        lrow_o = (tj==tk) ? lcol:1;
        for (lrow=lrow_o; lrow <= loadtj; lrow++)
        {
          trow = tj+lrow-1;
          tcol = tk+lcol-1;
          if (bDebug)  printf("\n tk=%d, tj=%d ** lcol=%d, lrow_o=%d, lrow=%d ** tcol=%d, trow=%d ",
                              tk,tj,lcol,lrow_o,lrow,tcol,trow);
          for (i=0; i < nRows*nCols; i++)
            if (masktest(pMask[i]))
            {
              /* catch floating-point underflow errors w/ tolerance "TOL" */
              if ((fabs(x[lcol][i]) > TOL) && (fabs(y[lrow][i]) > TOL))
              {
                x.c.html#sum2_xy">sum2_xy[trow][tcol] += x[lcol][i]*y[lrow][i];
                if (trow==tcol)
                  x.c.html#sum_x">sum_x[tcol] += x[lcol][i];
              }
            }  /* END pMask if-block */
        } /* END lrow for-loop */
      } /* END lcol for-loop */

    } /* END tj for-loop */
  } /* END tk for-loop */


  /* ========== Calculate Mean Vector And Covariance Matrix =========== */
  for (tk=1; tk <= nTimeSteps; tk++)
    meant[tk] = sum_x[tk]/(float)MskTotCellcount;        /* mean */

  for (tk=1; tk <= nTimeSteps; tk++)
    for (tj=tk; tj <= nTimeSteps; tj++)
    {
      covtt[tj][tk] = sum2_xy[tj][tk]/(float)MskTotCellcount - meant[tj]*meant[tk];
      covtt[tk][tj] = covtt[tj][tk];
    }
}

/* ================================================================= */
/* =============== WRITE COVARIANCE MATRIX TO FILE ================= */
void Output(void)
{
  FILE *fout;

  fout = OUTPUTfile(pszOutputFile, ".cov");

  for (tj=1; tj <= nTimeSteps; tj++)
  {
    if (bDebug)  printf("  %f   ", meant[tj]);
    if (tj!=1)  fprintf(fout, "\n");
    for (tk=1; tk <= nTimeSteps; tk++)
    {
      if (tk!=1)  fprintf(fout, ",");
      /*    to output the correlation matrix, use:
       *    covtt[tj][tk] /= sqrt(covtt[tj][tj]*covtt[tk][tk])  */
      fprintf(fout, "%f", covtt[tj][tk]);
    }
  }

  fclose(fout);
}


/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free1d_d(meant,1,nTimeSteps);
  free1d_d(sum_x,1,nTimeSteps);

  free2d_f(x,1,2*LOADLIMIT,0,nRows*nCols-1);
  free2d_d(sum2_xy,1,nTimeSteps,1,nTimeSteps);
  free2d_d(covtt,1,nTimeSteps,1,nTimeSteps);
}