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

      MULSTAT.C

      header file: VSTAT.H

      library: VSTAT_S.LIB

      These routines calculate means and covariance matrices
      for time series of 2 or more vector variables, up to a
      maximum that is now set at 10 but can easily be 
      increased.  For example, given a velocity profile time 
      series, one can calculate the mean u and v, the 
      variance of u and v, and their covariance, all as 
      functions of depth.  Given the vertical Fourier 
      transforms of the velocity profiles, one can calculate 
      the cross spectral matrix.  In this case, the real and 
      imaginary parts of the transforms are treated as 
      separate variables.

      For an example of how to call the routines, see the 
      test program at the end of this file.

            Eric Firing
            88-02-16

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

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

#if PROTOTYPE_ALLOWED
void allocate_multistat(MULTISTAT_TYPE *s, int n, int nvar, char *name,
  char *var_name[])
#else
void allocate_multistat(s, n, nvar, name, var_name)
   MULTISTAT_TYPE *s;
   int n, nvar;
   char *name;
   char *var_name[];
#endif
{
   static char err_msg[] = "allocate_multistat; out of memory";
   int i, j;

   s->n = n;
   s->nv = nvar;
   for (i=0; i < (s->nv); i++)
   {
      s->sum[i] = (double *) calloc(s->n, sizeof(double));
      check_error( (s->sum[i] == NULL), err_msg);
      s->var_name[i] = (char *) calloc(strlen(var_name[i])+1, sizeof(char));
      check_error( (s->var_name[i] == NULL), err_msg);
      strcpy(s->var_name[i], var_name[i]);
      for (j=i; j < (s->nv); j++)
      {
         s_sumpr(s,i,j) = (double *) calloc(s->n, sizeof(double));
         check_error( (s_sumpr(s,i,j) == NULL), err_msg);
         s_npts(s,i,j) = (int *) calloc(s->n, sizeof(int));
         check_error( (s_npts(s,i,j) == NULL), err_msg);
      }
   }
   s->name = (char *) calloc(strlen(name)+1, sizeof(char));
   check_error( (s->name == NULL), err_msg);
   strcpy(s->name, name);
   return;              /* normal exit */
}                       /* allocate_multistat */

#if PROTOTYPE_ALLOWED
void free_multistat(MULTISTAT_TYPE *s)
#else
void free_multistat(s)
MULTISTAT_TYPE *s;
#endif
{
   int i, j;

   for (i=0; i < (s->nv); i++)
   {
      free(s->sum[i]);
      free(s->var_name[i]);
      for (j=i; j < (s->nv); j++)
      {
         free(s_sumpr(s,i,j));
         free(s_npts(s,i,j));
      }
   }
   free(s->name);
}                       /* free_multistat */

#if PROTOTYPE_ALLOWED
void zero_double_array(double *d, int n)
#else
void zero_double_array(d, n)
   double *d;
   int n;
#endif
{
   int i;

   for (i=0; i<n; i++)
      *d++ = 0.0;
}                       /* zero_double_array */

#if PROTOTYPE_ALLOWED
void zero_int_array(int *ia, int n)
#else
void zero_int_array(ia, n)
   int *ia;
   int n;
#endif
{
   int i;

   for (i=0; i<n; i++)
      *ia++ = 0;
}                       /* zero_int_array */

#if PROTOTYPE_ALLOWED
void zero_multistat(MULTISTAT_TYPE *s)
#else
void zero_multistat(s)
   MULTISTAT_TYPE *s;
#endif
{
   int i, j;

   for (i=0; i < (s->nv); i++)
   {
      zero_double_array(s->sum[i], s->n);
      for (j=i; j < (s->nv); j++)
      {
         zero_double_array(s_sumpr(s,i,j), s->n);
         zero_int_array(s_npts(s,i,j), s->n);
      }
   }
}                       /* zero_multistat */


#if PROTOTYPE_ALLOWED
void update_multistat(MULTISTAT_TYPE *s, float *array[], int na)
#else
void update_multistat(s, array, na)
   MULTISTAT_TYPE *s;
   float *array[];
   int na;           /* number of points in each array */
#endif
{
   double *s_ptr;
   float *fi_ptr, *fj_ptr;
   int *n_ptr;
   int i, j, k;
   na = min_val(na, s->n);

   for (i=0; i < s->nv; i++)
   {
      s_ptr = s->sum[i];
      fi_ptr = array[i];
      /* accumulate the sum for each array element */
      for (k=0; k<na; k++, s_ptr++, fi_ptr++)        
         if (*fi_ptr <= ADJ_BADFLOAT)  
            *s_ptr += *fi_ptr;
      /* for each unique pair of variables: */
      for (j=i; j < (s->nv); j++)
      {
         /* initialize the pointers to the array beginnings */
         fi_ptr = array[i];
         fj_ptr = array[j];
         s_ptr  = s_sumpr(s,i,j);
         n_ptr  = s_npts(s,i,j);

         /* for each array element */
         for (k=0; k<na; k++, fi_ptr++, fj_ptr++, s_ptr++, n_ptr++)
         {
            if ((*fi_ptr <= ADJ_BADFLOAT) && (*fj_ptr <= ADJ_BADFLOAT))
            *s_ptr += (*fi_ptr) * (*fj_ptr);        /* sum of products */
            (*n_ptr)++;                            /* count */
         }
      }
   }
}                       /* update_multistat */

#if PROTOTYPE_ALLOWED
void calculate_multistat(MULTISTAT_TYPE *s)
#else
void calculate_multistat(s)
   MULTISTAT_TYPE *s;
#endif
{
   double *s_ptr, *mean_i_ptr, *mean_j_ptr;
   int *n_ptr;
   int i, j, k;

   for (i=0; i < s->nv; i++)
   {
      s_ptr = (s->sum)[i];
      n_ptr = s_npts(s, i, i);
      /* calculate the mean for each array element */
      for (k=0; k<s->n; k++, s_ptr++, n_ptr++)
         if (*n_ptr > 0)
            *s_ptr /= *n_ptr;
         else
            *s_ptr = BADFLOAT;
   }
   for (i=0; i < s->nv; i++)
   {
      /* for each unique pair of variables: */
      for (j=i; j < (s->nv); j++)
      {
         /* initialize the pointers to the array beginnings */
         s_ptr  = s_sumpr(s,i,j);
         n_ptr  = s_npts(s,i,j);
         mean_i_ptr = (s->sum)[i];    /* this is now the mean */
         mean_j_ptr = (s->sum)[j];    /* this is now the mean */

         /* for each array element, replace sum of products with covariance */
         for (k=0; k<s->n; k++, s_ptr++, n_ptr++, mean_i_ptr++, mean_j_ptr++)
         {
            if (*n_ptr > 0)
               *s_ptr = (*s_ptr)/(*n_ptr) - (*mean_i_ptr) * (*mean_j_ptr);
            else
               *s_ptr = BADFLOAT;
         }
      }
   }
}                       /* calculate_multistat */

/***************************************************************
      write_multistat()

      This is a fairly general output routine for the result 
      of a multistat calculation.  It does not currently 
      print the number of points used in each covariance 
      calculation.  Instead, it writes only the number of 
      points used in the variance and the mean.  In most 
      cases this should be adequate for standard output.  If 
      not, the option of writing n with each covariance can 
      be added.

      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) k+1 : d[k])
#define sqrt_(x) sqrt(max_val(x, 0.0))
#define cor_coef(s,i,j,k) s_sumpr(s,i,j)[k]/sqrt_(s_sumpr(s,i,i)[k]*s_sumpr(s,j,j)[k])
#define print_fixed(x) (x) >= ADJ_BADFLOAT ? fprintf(fpt, "      1E38") : fprintf(fpt, "%10.3f", (x))

#if PROTOTYPE_ALLOWED
void write_multistat(FILE *fpt, MULTISTAT_TYPE *s, double scale, float *d,
  char *d_name, int nd_prec, char cv_flag)
#else
void write_multistat(fpt, s, scale, d, d_name, nd_prec, cv_flag)
   FILE *fpt;
   MULTISTAT_TYPE *s;
   double scale;   /* e.g., set scale to 0.01 to convert m/s to cm/s */
   float *d;       /* if it is a null pointer, then the index will be used */
   char *d_name;   /* max 12 char; to be used in labelling the "depth" column */
   int  nd_prec;   /* precision (decimal places) for "depth" column */
   char cv_flag;   /* c for correlation, v for covariance */
#endif
{
   int i, j, k;
   int cor = 0;        /* 1 for correlation output */
   double ssq;

   check_error( (scale == 0.0), "write_multistat; scale is zero");
   ssq = scale * scale;
   if (cv_flag == 'c')
      cor = 1;
   fprintf(fpt, "\nStatistics for %s (scaled by %f)", s->name, scale);
   for (i=0; i<s->nv; i++)
   {
      fprintf(fpt, "\n------- %s -------  %s\n", s->var_name[i],
                         (cor ? "Correlation (normalized cov)" :
                                 "Covariance"));
      fprintf(fpt, "%12s    N       MEAN ", d_name);
      for (j=i; j<s->nv; j++)
         fprintf(fpt, "%10s ", s->var_name[j]);
      for (k=0; k<s->n; k++)  if (depth <= ADJ_BADFLOAT)
      {
/*         if ((s->npts)[i]>0)  */
         if ( (s_npts(s, i, i)[k]) > 0 )
            fprintf(fpt, "\n%10.*f %6d %10.3f ", nd_prec,
                       depth, s_npts(s,i,i)[k], ((s->sum)[i])[k]/scale);
         for (j=i; j<s->nv; j++)
            print_fixed((cor ? cor_coef(s,i,j,k) :
                             (s_sumpr(s,i,j)[k])*ssq));
/*            fprintf(fpt, "%10.3f ", (cor ? cor_coef(s,i,j,k) :
                                          (s_sumpr(s,i,j)[k])*ssq)); */
      }
      fprintf(fpt, "\n");
   }
   fprintf(fpt, "\n");
}                       /* write multistat */
#undef depth
#undef cor_coef

#if 0
/***************************************************************

      MSTEST.C

      This is a quick and dirty test of the multistat routines.
      The command line should contain 3 frequencies in 
      cycles per record length, followed by the letter c or 
      v to indicate correlation or covariance output.

      e.g.    mstest 2 2.4 3 c

      Link with vector_s.lib.

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

#include <stdio.h>
#include <math.h>
#include "\turboc\source\vector.h"
#include "vstat.h"

#define ND 20   /* The number of depth elements per vector */


/* The array of input vectors for each call to 
   update_multistat:
*/
float data0[ND];
float data1[ND];
float data2[ND];

float depth[ND];

/* The array of vector names:
*/
char name0[] = "first";
char name1[] = "second";
char name2[] = "third";


MULTISTAT_TYPE ms;

#if PROTOTYPE_ALLOWED
void main(int argc, char **argv)
#else
void main(argc, argv)
int argc;
char **argv;
#endif
{
   int n = ND;
   double amp,
          cycles,
          degrees;
   char cv_flag;
   float *data[3];
   char  *names[3];
   int i;

   if (argc != 5)
   {
      printf("\n need 3 frequencies in cycles per record length \n");
      exit(-1);
   }
   cv_flag = *argv[4];
   for (i=0; i<ND; i++)
      depth[i] = i * 20.0;
   data[0] = data0;
   data[1] = data1;
   data[2] = data2;
   names[0] = name0;
   names[1] = name1;
   names[2] = name2;
   ms.n = ND;
   ms.nv = 3;
   allocate_multistat(&ms, "test ms", names);
   zero_multistat(&ms);

   for (i=0; i<36; i++)
   {
      degrees = i * 10.0;
      amp = 1.0;  cycles = atof(argv[1]);
      make_sinusoid(data0, n, amp, cycles, degrees);

      amp = 2.0;  cycles = atof(argv[2]);
      make_sinusoid(data1, n, amp, cycles, degrees);

      amp = 3.0;  cycles = atof(argv[3]);
      make_sinusoid(data2, n, amp, cycles, degrees);

      update_multistat(&ms, data, n);
   }

   calculate_multistat(&ms);
   printf("\n Multistat calculated.\n");

   write_multistat(stdout, &ms, 1.0, (float *) NULL, "Index", 0, cv_flag);
   printf("\n Multistat written with index. \n");

   write_multistat(stdout, &ms, 1.0, depth, "Depth", 1, cv_flag);
   printf("\n Multistat written with depth. \n");

   free_multistat(&ms);
   printf("\n Multistat freed.\n");
}
#endif