// This is histogram.h classes and methods to implement a simple
// histogram package.
// original code developed for AIX xlC 1995-6
// modified Dec 1996 for sun solaris 2.5 and g++
#ifndef HISTOGRAM_H
#define HISTOGRAM_H

//#include <iostream.h>
#include <iostream>

using namespace std;

#include <string.h>
#include <math.h>

#define NUMBER_HBINS 20
#define SP " "
#define TAB '\t'

class histogram
 {
  private:
   char Title[80];
   int bin[NUMBER_HBINS];
   int overflow;
   int underflow;
   float low_edge;
   float high_edge;
   double xsum;
   double x2sum;

  public:
   histogram( float lowEdge, float highEdge, char* inTitle);
   void initialize( float lowEdge, float highEdge, char* inTitle);
   void update( float value );
   void update( int value );
   void print();
   void printLong(); // both horizontal and vertical format ... mostly for testing
   int get_sumIn(){
	int sumIn=0; for (int i=0; i < NUMBER_HBINS; ++i) sumIn = sumIn + bin[i];return sumIn;}

   double mean() { return (double) xsum/get_sumIn(); }
   double stdev(){  return sqrt ( (x2sum -get_sumIn()*mean()*mean())/(get_sumIn()-1));}  
 };

 histogram::histogram(float lowEdge, float highEdge, char* inTitle )
 {
  cout << "In histogram contructor: inTitle = " << inTitle << endl;
  low_edge = lowEdge;
  high_edge = highEdge;
  underflow =0; overflow =0;
  for ( int i=0; i<NUMBER_HBINS; ++i) bin[i]=0;
  strcpy( Title, inTitle);
  xsum=0.0;
  x2sum=0.0;
 }

 void histogram::initialize( float lowEdge, float highEdge, char* inTitle )
 {
  low_edge = lowEdge;
  high_edge = highEdge;
  underflow =0; overflow =0;
  for ( int i=0; i<NUMBER_HBINS; ++i) bin[i]=0;
  strcpy( Title, inTitle);
  xsum=0.0;
  x2sum=0.0;
 }

// a bin function for update based on khist.c
/* calculate a bin number for an array histogram[n+2]	*/
/* thus the histogram is 0 | 1, 2,   ,n | n+1		*/
/*                     under             over		*/
int histogramBin(float x, float a, float b, int n)
    {
        int bin; float y;
        if (x<a) bin=0;
        if (x>=b) bin = n+1;
        if (x>=a && x < b)
            {
                y = (n*(x-a))/(b-a);
                bin = (int) y + 1;
            }
        return bin;
     }

 void histogram::update( float value )
 {
  int binIndex = histogramBin( value, low_edge, high_edge, NUMBER_HBINS );
  if ( binIndex == 0 ) underflow++;
  if ( binIndex == NUMBER_HBINS + 1 ) overflow++;
  if ( binIndex > 0 && binIndex < NUMBER_HBINS + 1 ) {
	bin[binIndex-1]++;
	xsum = xsum + (double) value;
	x2sum = x2sum + (double) value * (double) value;
  }
  //cout << "In histogram::update " << value << " binIndex=" << binIndex << endl;
 }

 void histogram::update( int ivalue )
 {
  float value = ivalue;
  int binIndex = histogramBin( value, low_edge, high_edge, NUMBER_HBINS );
  if ( binIndex == 0 ) underflow++;
  if ( binIndex == NUMBER_HBINS + 1 ) overflow++;
  if ( binIndex > 0 && binIndex < NUMBER_HBINS + 1 ) {
	bin[binIndex-1]++;
	xsum = xsum + (double) ivalue;
	x2sum = x2sum + (double)  ivalue * (double) ivalue ;
  }
  //cout << "In histogram::update " << value << '\t' << binIndex << endl;
 }

 void histogram::printLong()
 {
  int i,j, dbin;
  float xmean, xstddev;
  cout << endl;
  cout << "Histogram Title: " << Title << endl;
  cout << "Underflow = " << underflow 
	<< " Overflow = " << overflow;
  int sumIn, sumUnderInOver;
  sumIn = 0; sumUnderInOver = 0;
  for ( i=0; i < NUMBER_HBINS; ++i) sumIn = sumIn + bin[i];
  sumUnderInOver = underflow + sumIn + overflow;
  cout << " Number in Histogram = " << sumIn
	<< " Total = " << sumUnderInOver;
  cout << endl;
  cout << NUMBER_HBINS << " bins with ";
  cout << "Low Edge value = " <<  low_edge 
	<< " Bin width = " << (high_edge - low_edge)/NUMBER_HBINS
	<< " High Edge value = " << high_edge;
  cout << endl;

  xmean = xsum/sumIn;
  //cout << " sumIn = " << sumIn << endl;
  xstddev = sqrt ( (x2sum -sumIn*xmean*xmean)/(sumIn-1)  );
  //cout << "x2sum = " << x2sum << endl;
  //cout << " sqrt ( "<< x2sum/(sumIn-1) - xmean*2<< " )"<<endl;
  cout << "Mean based on values = " << xmean << endl;
  cout << "StdDev based on values = " << xstddev << endl;

  // print in horizantal style...
  for ( i=0; i < NUMBER_HBINS; i++) cout << bin[i] << " ";
  cout << endl;
  // use vertical style for testing....
	 int cumulativeSum = 0;
	 float localEdge = 0.0;
         float localBinCenter = 0.0;
	 float binWidth = (high_edge - low_edge)/NUMBER_HBINS;
         float wgtMean = 0.0, wgtMean2 = 0.0, wgtSum = 0.0, wgtStdDev = 0.0;
  cout << "binNumb edgeVal Content cumSum " << endl;
  int maxValue = -1;
  for ( i=1; i<NUMBER_HBINS; i++) if (bin[i] >= maxValue ) maxValue=bin[i];
  
 if ( sumIn > 0 ) {
  for (int i=0; i < NUMBER_HBINS; i++)
	{
	 cumulativeSum = cumulativeSum + bin[i];
         localEdge = low_edge + i*binWidth;
         localBinCenter = localEdge + binWidth/2.0;
         wgtSum = wgtSum + bin[i];
         wgtMean = wgtMean + localBinCenter*bin[i];
         //wgtMean2 = wgtMean2 + localBinCenter*bin[i]*localBinCenter*bin[i];
	 wgtMean2 = wgtMean2 + localBinCenter*localBinCenter*bin[i];
	 dbin=histogramBin( (float) bin[i], 0.0, 1.2*maxValue, 50);
	 cout << i << TAB << localEdge << TAB << bin[i] << TAB << cumulativeSum << TAB;
		for ( j=1; j<dbin; j++) cout << "*";	
		cout << endl;
	}
 } // end of if ( sumIn > 9) 
  float temp = 0.0;
  if ( wgtSum > 0.0 ) wgtMean = wgtMean/wgtSum;
  else wgtMean = 0.0;
  if ( wgtSum > (float) NUMBER_HBINS ) 
	{
	 temp = wgtMean2/wgtSum - wgtMean*wgtMean;
         temp = (wgtSum/(wgtSum - NUMBER_HBINS))*temp;
//         temp = (wgtSum/(wgtSum - 1))*temp;
	}
  else temp = 0.0;
  if ( temp > 0.0 ) wgtStdDev = sqrt ( temp );
  else wgtStdDev = 0.0;
  cout << "Based on binned values mean = " << wgtMean 
       << " standard deviation = " << wgtStdDev << endl;
 }

 void histogram::print()
 {
  float xmean, xstddev;
  cout << "---------------------------- Histogram PrintOut ---------------------" << endl;
  cout << "Histogram Title: " << Title << endl;
  cout << "Underflow = " << underflow 
	<< " Overflow = " << overflow;
  int sumIn, sumUnderInOver;
  sumIn = 0; sumUnderInOver = 0;
  for ( int i=0; i < NUMBER_HBINS; ++i) sumIn = sumIn + bin[i];
  sumUnderInOver = underflow + sumIn + overflow;
  cout << " Number in Histogram = " << sumIn
	<< " Total = " << sumUnderInOver;
  cout << endl;
  cout << NUMBER_HBINS << " bins with ";
  cout << "Low Edge value = " <<  low_edge 
	<< " Bin width = " << (high_edge - low_edge)/NUMBER_HBINS
	<< " High Edge value = " << high_edge;
  cout << endl;
  // print in horizantal style...
  for (int i=0; i < NUMBER_HBINS; i++) cout << bin[i] << " ";
  cout << endl;
	 float localEdge = 0.0;
         float localBinCenter = 0.0;
	 float binWidth = (high_edge - low_edge)/NUMBER_HBINS;
         float wgtMean = 0.0, wgtMean2 = 0.0, wgtSum = 0.0, wgtStdDev = 0.0;
  for (int i=0; i < NUMBER_HBINS; i++)
	{
         localEdge = low_edge + i*binWidth;
         localBinCenter = localEdge + binWidth/2.0;
         wgtSum = wgtSum + bin[i];
         wgtMean = wgtMean + localBinCenter*bin[i];
         wgtMean2 = wgtMean2 + localBinCenter*localBinCenter*bin[i];
	}
  float temp = 0.0;
  if ( wgtSum > 0.0 ) wgtMean = wgtMean/wgtSum;
  else wgtMean = 0.0;
  if ( wgtSum > (float) NUMBER_HBINS ) 
	{
	 temp = wgtMean2/wgtSum - wgtMean*wgtMean;
         temp = (wgtSum/(wgtSum - NUMBER_HBINS))*temp;
	}
  else temp = 0.0;
  if ( temp > 0.0 ) wgtStdDev = sqrt ( temp );
  else wgtStdDev = 0.0;
  cout << "Based on binned values mean = " << wgtMean 
       << " standard deviation = " << wgtStdDev << endl;
  xmean = xsum/sumIn;
  //cout << " sumIn = " << sumIn << endl;
  xstddev = sqrt ( (x2sum -sumIn*xmean*xmean)/(sumIn-1)  );
  //cout << "x2sum = " << x2sum << endl;
  //cout << " sqrt ( "<< x2sum/(sumIn-1) - xmean*2<< " )"<<endl;
  cout << "Mean based on values = " << xmean << " ";
  cout << "StdDev based on values = " << xstddev << endl;
  cout << "---------------------------- Histogram PrintOut ---------------------" << endl;
 }




int khist( float x, float a, float b, int nbins ) {
     int index;
     float y;
     y = ( (nbins-0)*x + b-a - (nbins-0)*a )/(b-a);
     index = (int) y;
     if ( x < a ) index = 0;
     if ( x > b ) index = nbins + 1;
     //printf("IN khist: x = %f y= %f index = %d \n", x,y,index);
     return index;
} // end of khist


double flat( double a, double b ) {
	double x01;
	x01 = (double) rand() / RAND_MAX;
	return  (( b - a)*x01 + a);
} // end of flat

double normal( double mean, double s ) {
	double x, x01;
	x=0.0;
	for ( int i=0; i<12; i++) {
		x01 = (double) rand() / RAND_MAX;
		x = x01 + x;
	}
	x = s*(x-6.0)+mean;
	return x;
}
		

#endif //end of HISTOGRAM_H