/***************************************************************

      UNISTAT.C

      header file: VSTAT.H

      library: VSTAT_S.LIB

      These routines calculate the mean and variance
      for vector time series.

            Eric Firing
            88-02-16

      Modified 88-04-18 to include calculation of extrema.

****************************************************************/

#include "common.h"
#include "dbext.h"   /* min_val */
#include "vstat.h"
#include "ioserv.h"  /* check_error() */

static char mem_msg[] =  "allocate_unistat, out of memory";

void allocate_unistat(UNISTAT_TYPE *s, int n, char *name, int options)
{
   s->n = n;
   s->sum = (double *) calloc(s->n, sizeof(double));
   s->npts = (int *) calloc(s->n, sizeof(int));
   s->name = (char *) calloc(strlen(name)+1, sizeof(char));
   check_error( (s->sum == NULL
                  || s->npts == NULL || s->name == NULL), mem_msg);

   if (options & U_INDICES)    /* if index option is requested */
      options |= U_EXTREMA;    /*    then require extrema also */

   /* To print out the extrema divided by the standard deviations,
      we need to calculate the variances as well as the extrema.
   */
   if (options & U_NORM_EXTREMA)
       options = options | U_EXTREMA | U_VARIANCE;
   /* checking for invalid options should be added */
   s->options = options;

   if (s->options & U_VARIANCE)
   {
      s->sumsq = (double *) calloc(s->n, sizeof(double));
      check_error( (s->sumsq == NULL), mem_msg);
   }
   if (s->options & U_EXTREMA)
   {
      s->min = (float *) calloc(s->n, sizeof(float));
      s->max = (float *) calloc(s->n, sizeof(float));
      check_error( (s->min == NULL || s->max == NULL) , mem_msg);
   }
   if (s->options & U_INDICES)
   {
      s->imin = (int *) calloc(s->n, sizeof(int));
      s->imax = (int *) calloc(s->n, sizeof(int));
      check_error( (s->imin == NULL || s->imax == NULL) , mem_msg);
   }
   strcpy(s->name, name);
}                       /* allocate_unistat */

void free_unistat(UNISTAT_TYPE *s)
{
   free(s->sum);
   free(s->npts);
   free(s->name);
   if (s->options & U_VARIANCE)
      free(s->sumsq);
   if (s->options & U_EXTREMA)
   {
      free(s->min);
      free(s->max);
   }
   if (s->options & U_INDICES)
   {
      free(s->imin);
      free(s->imax);
   }
}                       /* free_unistat */

void zero_unistat(UNISTAT_TYPE *s)
{
   int i;

   s->n_increments = 0;
   for (i=0; i<s->n; i++)
   {
      (s->sum)[i] = 0.0;
      (s->npts)[i] = 0;
   }
   if (s->options & U_VARIANCE)
      for (i=0; i<s->n; i++)
         (s->sumsq)[i] = 0.0;
   if (s->options & U_EXTREMA)
      for (i=0; i<s->n; i++)
      {
         (s->min)[i] = BADFLOAT;
         (s->max)[i] = -BADFLOAT;
      }
   if (s->options & U_INDICES)
      for (i=0; i<s->n; i++)
      {
         (s->imin)[i] = 0;
         (s->imax)[i] = 0;
      }
}                       /* zero_unistat */

/* The following code may not be very efficient, but is
   written for readability first, so as to get it working.
*/

void update_unistat(UNISTAT_TYPE *s, float *array, int na)
{
   double f;
   int i;
   int is_min, is_max;

   na = min_val(na, s->n);
   for (i=0; i<na; i++)
   {
      if( (f=array[i]) < ADJ_BADFLOAT)
      {
         is_min = 0;
         is_max = 0;
         ((s->npts)[i])++;
         (s->sum)[i] += f;
         if (s->options & U_VARIANCE)
            (s->sumsq)[i] += (f*f);
         if (s->options & U_EXTREMA)
         {
            if (f < (s->min)[i])
            {
               (s->min)[i] = f;
               is_min = 1;
            }
            if (f > (s->max)[i])
            {
               (s->max)[i] = f;
               is_max = 1;
            }
            if (s->options & U_INDICES)
            {
               if (is_min)
                  (s->imin)[i] = s->n_increments;
               if (is_max)
                  (s->imax)[i] = s->n_increments;
            }
         }
      }
   }
   (s->n_increments)++ ;
}                       /* update_unistat */

/* update_unistat_d is identical to update_unistat except for
   the input array type.
*/
void update_unistat_d(UNISTAT_TYPE *s, double *array, int na)
{
   double f;
   int i;
   int is_min, is_max;

   na = min_val(na, s->n);
   for (i=0; i<na; i++)
   {
      if( (f=array[i]) < ADJ_BADFLOAT)
      {
         is_min = 0;
         is_max = 0;
         ((s->npts)[i])++;
         (s->sum)[i] += f;
         if (s->options & U_VARIANCE)
            (s->sumsq)[i] += (f*f);
         if (s->options & U_EXTREMA)
         {
            if (f < (s->min)[i])
            {
               (s->min)[i] = f;
               is_min = 1;
            }
            if (f > (s->max)[i])
            {
               (s->max)[i] = f;
               is_max = 1;
            }
            if (s->options & U_INDICES)
            {
               if (is_min)
                  (s->imin)[i] = s->n_increments;
               if (is_max)
                  (s->imax)[i] = s->n_increments;
            }
         }
      }
   }
   (s->n_increments)++ ;
}                       /* update_unistat_d */

void calculate_unistat(UNISTAT_TYPE *s)
{
   double sum, ss, mean, var, dnpts;
   int i;

   for (i=0; i<s->n; i++)
   {
      sum = (s->sum)[i];
      if (s->options & U_VARIANCE)
         ss = (s->sumsq)[i];

      if ((dnpts = (double) (s->npts)[i]) != 0)
      {
         mean = sum/dnpts;
         (s->sum)[i] = mean;
         if (s->options & U_VARIANCE)
         {
            var = ss/dnpts - mean*mean;
            (s->sumsq)[i] = max_val(var, 0.0);
         }
      }
      else /* 93/11/30 jr+ set to BADFLOATs if no points */
      {
         (s->sum)[i] = BADFLOAT;
         if (s->options & U_VARIANCE)
         {
            (s->sumsq)[i] = BADFLOAT;
         }
      }
   }
}                       /* calculate_unistat */


/***************************************************************
      write_unistat()

      This is a fairly general output routine for the result
      of a unistat calculation.

      The input variable d is a pointer to a float array
      used for labelling the rows of output.  If the
      statistics are calculated from velocity profiles, then
      d is the array of depths.  If fourier transformed
      profiles were used, then d is the array of vertical
      wavenumbers.  If d is a NULL pointer, then the array
      index, starting from 1, will be printed.

***************************************************************/

#define depth  (d==NULL ? (float) i+1 : d[i])

void write_unistat(FILE *fpt, UNISTAT_TYPE *s, double scale, float *d,
  char *d_name, int nd_prec)
{
   int i;
   double std_dev;

   check_error( (scale == 0.0), "write_unistat; scale is zero");
   fprintf(fpt, "\nStatistics for %s (scaled by %f)", s->name, scale);
   fprintf(fpt, "\n%12s    N       MEAN", d_name);
   if (s->options & U_VARIANCE)
      fprintf(fpt, "    STD DEV        SEM");
   if (s->options & U_EXTREMA)
      fprintf(fpt, "    minimum    maximum");
   if (s->options & U_NORM_EXTREMA)
      fprintf(fpt, " min/sd max/sd");
   if (s->options & U_INDICES)
      fprintf(fpt, " imin imax");
   for (i=0; i<s->n; i++)  if (depth <= ADJ_BADFLOAT)
   {
      if ((s->npts)[i]>0)
      {
         fprintf(fpt, "\n%10.*f %6d %10.3f", nd_prec,
                       depth, (s->npts)[i], (s->sum)[i]/scale);
         if (s->options & U_VARIANCE)
         {
            std_dev = sqrt((s->sumsq)[i]);
            fprintf(fpt, " %10.3f %10.3f",
                          std_dev/scale,
                          (std_dev/sqrt( (double) ((s->npts)[i]) ))/scale );
         }
         if (s->options & U_EXTREMA)
            fprintf(fpt, " %10.3f %10.3f", (s->min)[i]/scale, (s->max)[i]/scale);
         if (s->options & U_NORM_EXTREMA)
            fprintf(fpt, " %6.2f %6.2f",
            (std_dev != 0.0) ? ((s->min)[i]-(s->sum)[i])/std_dev : 0.0,
            (std_dev != 0.0) ? ((s->max)[i]-(s->sum)[i])/std_dev : 0.0);
         if (s->options & U_INDICES)
            fprintf(fpt, " %4d %4d", (s->imin)[i], (s->imax)[i]);
      }
   }
   fprintf(fpt, "\n");
}                       /* write_unistat */
#undef depth