/* statsmaps.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          5/11/98
 *                *
 * COMPILE AS: cc -o statsmaps statsmaps.c functions.c util.c -lnetcdf  -lm
 * RUN AS: statsmaps <Name of output NC/binary file>  
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: statsmaps \n\
Program to extract descriptive statistics (mean, sd, min, max)\n\
per pixel, for different variables.\n\
\n\
Command-line switches specific to this program are:\n\
\n\
"};

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

void writeBinaryVar(char *filename, char *Type, float *array);
void writeNCvar(int ncID, char *parameter, char *Type, char *LongType, float *array);

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

#define  NUMPARAM 4


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Summary Statistics: Mean, SD, Min, Max");
  AssignString(&pszNChist, "stats file created by statsmaps");

  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,NUMPARAM-1,0,nVars-1,0,nRows*nCols-1);
  for (v=0; v < nVars;v++)
    for (i=0; i < nRows*nCols;i++)
    {
      out[MEAN][v][i] = 0.0;    /* mean     */
      out[SD][v][i]   = 0.0;    /* standard deviation */
      out[MIN][v][i]  = 0.0;    /* minimum      */
      out[MAX][v][i]  = 0.0;    /* maximum      */
    }

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

/* ================================================================= */
/* ======== READ DATA AND CALCULATE STATISTICAL PARAMETERS ========= */
void Process(void)
{
  int v, i, iout, col, row, n, i0, nyblock;   int bug=0;
  int tvar = 0;
  float fill_val, txyval;
  double stats_i[NUMPARAM], mean_i, sd_i;


  for (v = 0; v < nVars; v++)
  {
    printf("v=%d *** \n",v);
    for (nyblock = 0; nyblock < nYBlocks; nyblock++)
    {
      ReadTXYbl(TsXYblock, pstVars[v], 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);

          stats_i[SUM]  = 0.0;
          stats_i[SUM2] = 0.0;
          stats_i[MIN]  = 100000000.0;
          stats_i[MAX]  = -100000000.0;

          for (tvar = 0; tvar < nTimeSteps;tvar++)
          {
            txyval = TsXYblock[tvar][i-i0];
            /* catch floating-point underflow errors w/ tolerance "TOL" */
            if (txyval > TOL || txyval < -TOL)
            {
              stats_i[SUM]  += txyval;
              stats_i[SUM2] += txyval * txyval;
            }
            if (txyval < stats_i[MIN])
              stats_i[MIN] = txyval;
            if (txyval > stats_i[MAX])
              stats_i[MAX] = txyval;
          } /* END tvar FOR-LOOP (loop thru each time step) */

          mean_i = stats_i[SUM] / (double)nTimeSteps;
          sd_i   = sqrt(stats_i[SUM2]/(double)nTimeSteps - DSQR(mean_i));
          out[MEAN][v][iout] = (float) mean_i;
          out[SD][v][iout]   = (float) sd_i;
          out[MIN][v][iout]  = stats_i[MIN];
          out[MAX][v][iout]  = stats_i[MAX];

        } /* 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[MEAN][v][iout] = fill_val;
          out[SD][v][iout]   = fill_val;
          out[MIN][v][iout]  = fill_val;
          out[MAX][v][iout]  = fill_val;
        }
      } /* 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   */
  char FileName[NAMELEN];         /* Base File Name for binary files output */


  if (bOutBinary)    /* write output to binary files */
    for (v = 0; v < nVars; v++)
    {
      sprintf(FileName,"%s.%s", pszOutputFile, pstVars[v].param);
      writeBinaryVar(FileName, ".mean", out[MEAN][v]);
      writeBinaryVar(FileName, ".sd", out[SD][v]);
      writeBinaryVar(FileName, ".min", out[MIN][v]);
      writeBinaryVar(FileName, ".max", out[MAX][v]);
    }
  else
  {                              /* write output to a single netCDF file */
    SetupDefaultAxes(FALSE, FALSE);
    ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
    for (v = 0; v < nVars; v++)
    {
      writeNCvar(ncID, pstVars[v].param, "mean", "Means", out[MEAN][v]);
      writeNCvar(ncID, pstVars[v].param, "sd", "Standard Deviations", out[SD][v]);
      writeNCvar(ncID, pstVars[v].param, "min", "Minima", out[MIN][v]);
      writeNCvar(ncID, pstVars[v].param, "max", "Maxima", out[MAX][v]);
    }
    errnc = nc_close(ncID);
  }
}

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

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

/* ============================================================ */
void writeBinaryVar(char *filename, char *Type, float *array)
{
  FILE *fout;

  fout = OUTPUTfile(filename, Type);
  fwrite(array, sizeof(float), nRows*nCols, fout);
  fclose(fout);
}

/* ============================================================ */
void writeNCvar(int ncID, char *parameter, char *Type, char *LongType, float *array)
{
  int varID;                  /* variable ID    */
  char varName[NAMELEN];      /* Name of NC variable(s) to be written      */
  char varLongName[NAMELEN];  /* Long Name of NC variable(s) to be written */

  sprintf(varName,"%s_%s", parameter,Type);
  sprintf(varLongName,"%s %s", parameter,LongType);
  varID = DefNCvar(ncID, varName, EMU_FLOAT, "xy", varLongName, "");
  errnc = nc_put_var_float(ncID, varID, array);
}