/*  stax.c  (C) 1997 Mark A. Livingston */
/* Any non-commercial use of this code is permitted so long as */
/* the above copyright line appears in the file.  Please contact */
/* the author regarding permission for commercial use. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "stax.h"

static int numStat, numHist;
static Stat *stx;
static Hist *hst;

void
stxInit( int n )
{
    int i;

    stx = (Stat *) malloc( n * sizeof( Stat ) );
    if( stx == NULL ) {
      fprintf( stderr, "Unable to allocate statistical counters\n" );
      return;
    }

    for( i = 0; i < n; i++ ) {
      stx[i].s[ AVG ] = 0.0;
      stx[i].s[ STD ] = 0.0;
      stx[i].n = 0;
    }

    numStat = n;
    return;
}

void
stxUpdate( int i, double v )
{
    stx[i].s[ AVG ] += v;
    stx[i].s[ STD ] += SQR( v );
    if( stx[i].n > 0 ) {
      if( stx[i].s[ MIN ] > v )
        stx[i].s[ MIN ] = v;
      if( stx[i].s[ MAX ] < v )
        stx[i].s[ MAX ] = v;
    } else {
      stx[i].s[ MIN ] = v;
      stx[i].s[ MAX ] = v;
    }
    stx[i].n++;

    return;
}

void
stxPrint( int i, char *hdr, char *name )
{
    double std, avg, std_sq_numerator, nd;

    if( hdr != NULL )
      printf( "%d points %s\n", stx[i].n, hdr );
    if( stx[i].n > 0 ) {
      nd = (double) stx[i].n;
      if( hdr != NULL )
	printf( "Stat\t             AVG        STD        MIN        MAX\n" );
      if( stx[i].n > 1 ) {
	/* compute standard deviation from running sum of squares */
	std_sq_numerator = nd * stx[i].s[ STD ] - SQR( stx[i].s[ AVG ] );
	if( std_sq_numerator > 0.0 )
	  std = sqrt( std_sq_numerator / ( nd * ( nd - 1.0 ) ) );
	else
	  std = 0.0;
      } else
	std = 0.0;

      /* compute average from running sum */
      avg = stx[i].s[ AVG ] / nd;

      /* print the stats */
      if( name != NULL )
	printf( "%s %9.6lf  %9.6lf  %9.6lf  %9.6lf\n", name,
	        avg, std, stx[i].s[ MIN ], stx[i].s[ MAX ] );
      else
	printf( "         %9.6lf  %9.6lf  %9.6lf  %9.6lf\n",
	        avg, std, stx[i].s[ MIN ], stx[i].s[ MAX ] );
    }

    return;
}

void
stxClear( int i )
{
    /* min and max do not need reseting - the counter set to zero handles it */
    stx[i].s[ AVG ] = 0.0;
    stx[i].s[ STD ] = 0.0;
    stx[i].n = 0;
    return;
}

void
stxHistInit( int n, double *mins, double *maxs, double *widths )
{
    int i, j;

    hst = (Hist *) malloc( n * sizeof( Hist ) );
    if( hst == NULL ) {
      fprintf( stderr, "Unable to allocate histogram structures\n" );
      return;
    }

    for( i = 0; i < n; i++ ) {
      hst[i].min = mins[i];
      hst[i].max = maxs[i];
      hst[i].width = widths[i];
      hst[i].numSlots = ( maxs[i] - mins[i] ) / widths[i] + 2;
      /* assuming this works out to an integer, don't know what would happen */
      hst[i].slot = (int *) malloc( hst[i].numSlots * sizeof( int ) );
      for( j = 0; j < hst[i].numSlots; j++ )
	hst[i].slot[j] = 0;
    }

    return;
}

void
stxHistUpdate( int i, double v )
{
    int cell;

    if( v < hst[i].min ) {
      hst[i].slot[ H_MIN ]++;
    } else if( v > hst[i].max ) {
      hst[i].slot[ H_MAX ]++;
    } else {
      cell = H_BASE + (int) ( ( v - hst[i].min ) / hst[i].width );
      hst[i].slot[ cell ]++;
    }

    return;
}

void
stxHistPrint( int i )
{
    int j;
    double n1, n2;

    printf( "      <   %lf : %d\n", hst[i].min, hst[i].slot[ H_MIN ] );
    n1 = hst[i].min, n2 = hst[i].min + hst[i].width;
    for( j = 2; j < hst[i].numSlots;
	 j++, n1 += hst[i].width, n2 += hst[i].width )
      printf( "%lf, %lf : %d\n", n1, n2, hst[i].slot[j] );
    printf( "      >=  %lf : %d\n", hst[i].max, hst[i].slot[ H_MAX ] );

    return;
}