#include "pdb.h"

int max_measurement = 800  ;
 int *st_count = NULL;
 double *st_mean ;
 double *st_sd ;
 double *st_top ;
 double *st_bottom ;
 double *st_sum ;
 double *st_sum_sq ;



#define LOOP_SOME_MEASUREMENTS(max,body) { int index ; for(index = 0 ; index < max; index++) body }
#define LOOP_MEASUREMENTS(body) LOOP_SOME_MEASUREMENTS(max_measurement,body) 



void init_measurement_vectors()
{
  
  st_count = (int *) malloc (sizeof(int) * max_measurement) ;
  
  st_mean = (double *) malloc (sizeof(double) * 6 * max_measurement) ;
  st_sd  = st_mean + max_measurement ;
  st_top = st_sd + max_measurement ;
  st_bottom = st_top + max_measurement ;
  st_sum = st_bottom + max_measurement ;
  st_sum_sq = st_sum + max_measurement ;
}

void add_measurement (int index, double value)
{
  st_count[index] ++ ;

  st_sum[index] += value ;
  if (value > st_top[index]) st_top[index]  = value;
  if (value < st_bottom[index]) st_bottom[index]  = value;
  st_sum_sq[index] += value * value ;
}

void init_measurement(int index) 
{
  if (st_count == NULL)
    init_measurement_vectors() ;
  st_count[index] = 0 ;
  st_sum[index] = 0.0 ;
  st_sum_sq[index] = 0.0 ;
  st_top[index] = -1e300; 
  st_bottom[index] = 1e300; 
}

void init_all_measurements()
{
  LOOP_MEASUREMENTS( init_measurement(index) ;) ;
}

void init_some_measurements(int max)
{
  if (max >= max_measurement)
    {
      if (st_count != NULL)
	{
	  free (st_count) ;
	  free (st_mean) ;
	  st_count = NULL ;
	}
      max_measurement = max ;
    }
  LOOP_SOME_MEASUREMENTS(max, init_measurement(index) ;) ;
}


void sum_measurement (int index)
{
  if (st_count[index])
    {
      int count = st_count[index] ;
      double mean = st_sum[index] /count ;
      st_mean[index] = mean;
      st_sd[index] = count > 1 ?sqrt( (st_sum_sq[index] - mean * mean * count) / (count - 1)) : 1e50 ;
    }
}

void sum_some_measurements(int max)
{
  LOOP_SOME_MEASUREMENTS(max, sum_measurement(index);) ;
}