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

                            HISTO.C

   This is a set of routines for calculating simple statistics and
   histograms for each of several single variables.

                          Eric Firing
                           88-01-13

   07/12/90 - JR - changed 10 to 10.0 in call to pow()
                   removed cast to int in truncate_left macro so returns type x
   90-09-05 - JR - changed %f to %lf in printing scale and bin_lb
                   in write_histogram()
                 - restored cast to int in truncate_left macro--so will
                   truncate; forced caller to truncate_left to re-cast it
                   back to whatever

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

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

#define  h_bin(i,j)  (*(h->bin_array+i*(h->nbins+1)+j))
/* (the column dimension is h->nbins+1 because the last column
    is for the total count)
*/
#define includes(s,c)  (strchr(s,c) != NULL)

void allocate_histogram(HISTOGRAM_TYPE *h, int nrows, int nbins, char *name)
{
   h->nrows = nrows;
   h->nbins = nbins;
   h->bin_array = (unsigned int *) calloc((nrows * (nbins+1)), sizeof(int));
   h->name = (char *) calloc(strlen(name)+1, sizeof(char));
   if (h->bin_array == NULL || h->name == NULL)
   {
      printf("\n Not enough memory left for histogram of %s.\n", name);
      exit(-1);
   }
   strcpy(h->name, name);
}                       /* allocate_histogram */

void free_histogram(HISTOGRAM_TYPE *h)
{
   free(h->bin_array);
   free(h->name);
}                       /* free_histogram */

void zero_histogram(HISTOGRAM_TYPE *h)
{
   int i, n;
   n = h->nrows * (h->nbins + 1);
   for (i=0; i<n; i++)
      *(h->bin_array + i) = 0;
}                       /* zero_histogram */

void update_histogram(HISTOGRAM_TYPE *h, float *x, int n)
{
   int i,j;
   double fi;
   n = min_val(n, h->nrows);
   check_error( (h->bin_width == 0), "update_histogram; zero bin width");
   for (i=0; i<n; i++)
   {
      if ( (fi = x[i]) < ADJ_BADFLOAT )
      {
         j = (int)((fi-h->origin)/h->bin_width);
         j = max_val(j, 0);
         j = min_val(j, h->nbins-1);
         h_bin(i,j)++;
         h_bin(i,h->nbins)++; /* total count in last column */
      }
   }
}                       /* update_histogram */


#define truncate_left(x) (x)>0 ? (int) (x) : ( -(int)(x) -1 )
/* jr 7/12/90 - removed cast to int so will return whatever type x is */
/* jr 9/05/90 - restored cast to int */
#define depth  (d==NULL ? (float) i+1 : d[i])

void write_histogram(FILE *fp, HISTOGRAM_TYPE *h, int nd, float *d,
  char *d_name, int nd_prec, char *style)
{
   int i,j;
   int xi;
   double xd;
   double scale, exponent;
   double max_bin;    /* lower boundary of highest bin */
   double bin_lb;     /* bin lower boundary */
   int i_style = 0;   /* 0 = plain count, 1 = cumulative count,
                         2 = normalized,  3 = cum and norm  */

   if (includes(style, 'c')) i_style++;
   if (includes(style, 'n')) i_style += 2;

   max_bin = h->origin + h->bin_width * (h->nbins - 1);
   exponent = truncate_left(log10(max_bin));
   scale = pow(10.0, exponent);
   fprintf(fp,"\n\n%s  (bin lower boundaries: scale is %f)",h->name, scale);
   fprintf(fp, "\n%12s Total   ", d_name);
                        /* 6 extra spaces; shift to boundaries */
   for (j=1; j<h->nbins; j++)   /* skip first; no lower bdry for bin 0 */
   {
      bin_lb =  (h->origin + h->bin_width * j) / scale;
      if (!(j%5))  fprintf(fp, "  ");
      fprintf(fp, "%*.1f ", nd, bin_lb);
   }
   for (i=0; i<h->nrows; i++)  
   {
      xd = 0.0;
      xi = 0;
      if (!(i%10)) fprintf(fp, "\n");
      if (depth <= ADJ_BADFLOAT)
      {
         fprintf(fp, "\n%10.*f %5d", nd_prec, depth, h_bin(i, h->nbins));
         for (j=0; j<h->nbins; j++)
         {
            if (!(j%5)) fprintf(fp, "  ");
            switch (i_style)
            {
               case 0:          /* raw */
                  fprintf(fp,"%*d ", nd, h_bin(i,j));
                  break;
               case 1:          /* cumulative */
                  xi += h_bin(i,j);
                  fprintf(fp,"%*d ", nd, xi);
                  break;
               case 2:          /* normalized */
                  if (h_bin(i, h->nbins))    /* prevent division by zero */
                     fprintf(fp,"%*.*f ", nd, nd-2,
                           h_bin(i,j)/ (double)h_bin(i,h->nbins));
                  break;
               case 3:           /* normalized, cumulative */
                  if (h_bin(i, h->nbins))    /* prevent division by zero */
                  {
                     xd += h_bin(i,j)/ (double)h_bin(i,h->nbins);
                     fprintf(fp,"%*.*f ", nd, nd-2, xd);
                  }
                  break;
                  /* i_style is calculated in this routine, can only
                     be 0-3, so no need for default case.
                  */
            }   /* end case; printing each bin */
         }  /* end row printing */
      }  /* end of things to do if depth is valid */
   }   /* end of loop over the rows */
   fprintf(fp,"\n");
}                       /* write_histogram */