//
// This file test the functionalities of the Block class.
// Interectly it also test the functionalities of the ColumnAttribs class.
//

// test harness
#include "TestManager.h"

#include "HepTuple/HepHBookNtuple.h"
#include "HepTuple/HepNtuple.h"
#include "HepTuple/Column.h"
#include "HepTuple/Block.h"

#include "HepTuple/HepHist1D.h"
#include "HepTuple/HepHist2D.h"
#include "HepTuple/HepHistProf.h"

ZM_USING_NAMESPACE( zmht )

#include "ZMutility/iostream"

// Due to a g++ bug we commented out the part of the test
// using useMethod.h.  The problem is that this part will
// not compile with the -g flag but will compile with
// -ggdb3

Int4 c;

Int4 func () {
  return c;
}

void readFunc (Int4 i) {
  c = i;
}

class Test {
  int _a;
public:
  Int4 a() {return _a;};
  void seta(Int4 f) { _a = f; };
};

#include <cmath>
#ifndef DBL_EPSILON
#define    DBL_EPSILON     2.2204460492503131E-16
#endif

#define IS_DBL_EQUAL(a,b) ( a == b \
			    || ( ::fabs( (a) - (b) ) < 100*DBL_EPSILON \
			    || ( ::fabs( (a) - (b) ) / (b) ) < 100*DBL_EPSILON ) )

#ifndef FLT_EPSILON
#define    FLT_EPSILON     2.2204460492503131E-6
#endif

#define IS_FLT_EQUAL(a,b) ( a == b \
			    || ( ::fabs( (a) - (b) ) < 100*FLT_EPSILON \
			    || ( ::fabs( (a) - (b) ) / (b) ) < 100*FLT_EPSILON ) )

#include "ZMutility/fstream"

USING( namespace std )

const Int4 rowNumber = 10;

const int RANDLIMIT(100);
std::string name1D ( "1d_histogram" );
int nBins = 100;
float low = -11.9;
float high = 75.2;
unsigned int seed1d = 11;

std::string name2D ( "2d_histogram" );
int nBinsX = 10;
int nBinsY = 50;
float lowX = -11.9;
float highX = 75.2;
float lowY = 780.89;
float highY = 907.32;
unsigned int seed2d = 22;

std::string nameProf ( "profile_histogram" );
int nBinsP = 78;
unsigned int seedp = 33;

void Create( HepFileManager* );
void Transfer( HepFileManager* from, HepFileManager* to, const string& where,
	       bool no_ntuple = false, bool no_prof = false );
HepNtuple& addNtuple( HepFileManager* );
HepHist1D& addHist1D( HepFileManager* );
HepHist2D& addHist2D( HepFileManager* );
HepHistProf& addHistProf( HepFileManager* );
bool Check( HepNtuple& );
bool Check( const HepHist1D& );
bool Check( const HepHist2D& );
bool Check( const HepHistProf& );
bool Compare( const HepHist1D& , const HepHist1D&);
bool Compare( const HepHist2D& , const HepHist2D&);
bool Compare( const HepHistProf& , const HepHistProf&);

const char objectTypes = 4;

int main (int argc, char** argv) {
  ZMxHepUnknownHepObj::setHandler ( ZMexIgnoreAlways() );


  HepFileManager* files[3];
  string managerType[3];
  bool unimplemented[3][objectTypes];
  Int4 nFiles = 0;
  Int4 i, j;

  for ( i = 0; i < 3 ; i++ )
    for ( j = 0; j < objectTypes; j++ )
      unimplemented[i][j] = false;

#ifdef USE_HBOOK
  files[nFiles] =  OPEN_HBOOK_FILE("transfer");
  managerType[nFiles] = "hbook";
  if ( files[nFiles] != NULL ) nFiles++;
#endif
#ifdef USE_HISTOSCOPE
  files[nFiles] =  OPEN_HISTO_FILE("transfer");
  managerType[nFiles] = "histoscope";
  unimplemented[nFiles][3] = true; // no histogram profile
  if ( files[nFiles] != NULL ) nFiles++;
#endif
#ifdef USE_ROOT
  files[nFiles] =  OPEN_ROOT_FILE("transfer");
  managerType[nFiles] = "root";
//  unimplemented[nFiles][0] = true; // no ntuple
  if ( files[nFiles] != NULL ) nFiles++;
#endif


  // build original content
  for ( i = 0; i < nFiles; i++ ) {
    Create( files[i] );
  }

  // transfer it to all the managers
  for ( i = 0; i < nFiles; i++ ) {
    Int4 j;
    for ( j = 0; j < nFiles; j++ ) {
      string new_dir = string("from_")+managerType[i];
      Transfer ( files[i], files[j], new_dir,
		 unimplemented[i][0] || unimplemented[j][0],
		 unimplemented[i][3] || unimplemented[j][3]);
    }
  }

   // List final content
  for ( i = 0; i < nFiles; i++ ) {
    files[i]->ls("","r");
  }
  for( i=0; i<nFiles; ++i ) {
    if( files[i] != NULL )  delete files[i];
  }
  return 0;
}


void Create( HepFileManager* manager ) {
  cerr << "Creating " << manager->fileName() << endl;

  manager->mkdir("original");
  manager->cd("original");
  HepNtuple & ntuple  = addNtuple( manager );
  HepHist1D & h1 = addHist1D( manager );
  HepHist2D & h2 = addHist2D( manager );
  HepHistProf & hp = addHistProf( manager );

  if ( ntuple.isValid() ) Check ( ntuple );
  Check ( h1 );
  Check ( h2 );
  if ( hp.isValid() ) Check ( hp );
}


void Transfer( HepFileManager* from, HepFileManager* to, const string& where,
	       bool no_ntuple, bool no_prof) {
  cerr << "Transfer from " << from->fileName() << " to " << to->fileName() << endl;

  from->cd();
  from->cd("original"); // or is it "//"+from->fileName()+"/original"

  HepNtuple & origNtuple = from->retrieveNtuple("ntuple");
  HepHist1D & origHist1D = from->retrieveHist1D(name1D);
  HepHist2D & origHist2D = from->retrieveHist2D(name2D);
  HepHistProf & origHistProf = from->retrieveHistProf(nameProf);

  to->cd();
  to->mkdir(where);
  to->cd(where);

  if ( ! no_ntuple ) origNtuple.makeClone( to, "copyOfNtuple" );
  origHist1D.makeClone( to, "copyOfHist1D" );
  origHist2D.makeClone( to, "copyOfHist2D" );
  if ( ! no_prof ) origHistProf.makeClone( to, "copyOfHistProf" );

  if ( ! no_ntuple ) Check( to->retrieveNtuple("copyOfNtuple") );
  Compare( to->retrieveHist1D("copyOfHist1D"), origHist1D );
  Compare( to->retrieveHist2D("copyOfHist2D"), origHist2D );
  if ( ! no_prof ) Compare( to->retrieveHistProf("copyOfHistProf"), origHistProf );
}

// Use both by addNtuple and Check(HepNtuple&)
int dim1 = 10;
int dim2 = 10;

HepNtuple& addNtuple( HepFileManager* manager) {

  Int4 atInt1 = 55;
  Int4 atInt2 = 55;
  //  Float8 atDouble = 55.0;
  Test object;
  // Float8 arr[10][10];
  Float8 arr[100][10];

 int i;
  for ( i = 0; i < dim1 ; i++) {
    ((double*)arr)[i] = 1.0;
  }

  cerr << "Start ntuple creation \n";

  HepNtuple& hep = manager->ntuple("ntuple",1);

  hep.setColumnWise();

  hep.column("Col1",3.2,7.0);
  hep.column("Col2",8.2);
  hep.columnVia("colF",func,(Int4)0);
  c = 0;

#ifdef USEMETHOD_WORKS
  hep.columnVia("colM",newUseMethod(&object, object.a));
#endif

  hep.column("Col3",8.2,99.0);
  hep.column("Col4",8.2,99.0);

  hep.columnAt("colAt",&atInt1);
  hep.columnAt("colAt2",&atInt2);

  hep.columnArray("colArray",2,(Float8)8).dimension(dim1,dim2);

  hep.capture();
  hep.storeCapturedData();
  hep.clearData();

  cerr << hep.blockFormat(0) << endl;
  cerr << "Filling ntuple" << endl;

  for (i=0; i<rowNumber; i++) {
    //    cerr << i << ":";

    hep.capture("Col1",3.2+i);
    hep.capture("Col2",8.2+2*i);

    atInt1 = -3*i;
    c = i;
    object.seta(5*i);
    // cerr << "before loop " << dim1 << "," << dim2 << "\n";
    for (int j=0; j < dim1+1; j++) {
      for (int k=0; k< dim2+1; k++) {
	arr[j][k] = 8.2/(i+1)*j*k;
	// cerr << j << "," << k << ":" << arr[j][k] << ";";
      }
    }
    hep.capture("colArray",&(arr[0][0]));
    hep.capture();
    hep.storeCapturedData();
    hep.clearData();
  }

  return hep;
}


bool Check( HepNtuple& hep ) {
  cerr << "Check Ntuple ... ";

  // Float8 arr[10][10];
  Float8 arr[100][10];

  int i;
  for ( i = 0; i < dim1 ; i++) {
    ((double*)arr)[i] = 1.0;
  }

  // set the readColumn back destinations.

  Float8 col1 =0.0;
  hep.readColumn(1,"Col1",&col1);
  if ( ! IS_DBL_EQUAL( col1, 3.2) ) {
    cerr << "first col1 is wrong " << col1 << endl;
  }
  hep.readColumn(1,"Col2",&col1);
  if ( ! IS_DBL_EQUAL( col1, 8.2) ) {
    cerr << "first col2 is wrong " << col1 <<  endl;
  }
  hep.readColumn(1,"Col3",&col1);
  if ( ! IS_DBL_EQUAL( col1, 8.2) ) {
    cerr << "first col3 is wrong " << col1 <<  endl;
  }
  hep.readColumn(1,"Col4",&col1);
  if ( ! IS_DBL_EQUAL( col1, 8.2) ) {
    cerr << "first col4 is wrong " << col1 <<  endl;
  }
  Int4 creadColumn = -99;
  hep.readColumn(1,"colF",&creadColumn);
  if ( creadColumn != 0 ) {
    cerr << "first colF is wrong " << creadColumn <<endl;
  }
#ifdef USEMETHOD_WORKS
  hep.readColumn(1,"colM",&creadColumn);
  if (  creadColumn != 0 ) {
    cerr << "first colM is wrong " << creadColumn <<endl;
  }
#endif
  hep.readColumn(1,"colAt",&creadColumn);
  if ( creadColumn != 55 ) {
    cerr << "first colAt is wrong " << creadColumn <<endl;
  }
  hep.readColumn(1,"colAt2",&creadColumn);
  if ( creadColumn != 55 ) {
    cerr << "first colAt2 is wrong " << creadColumn <<endl;
  }
  hep.readColumnArray(1,"colArray",&(arr[0][0]));
  for ( i = 0; i < dim1*dim2 ; i++) {
    if ( ! IS_DBL_EQUAL( ((double*)arr)[i] , 8.0 ) ) {
      cerr << "first colArray of " << i << " is wrong " << ((double*)arr)[i]
	<< endl;
    }
  }

  Float8 float1,float2,float3,float4;
  Int4 int1,int2;

  hep.destinationAt("Col1",&float1);
  hep.destinationAt("Col2",&float2);
  hep.destinationVia("colF",readFunc);
#ifdef USEMETHOD_WORKS
  hep.destinationVia("colM",newSetMethod(&object, object.seta));
#endif
  hep.destinationAt("Col3",&float3);
  hep.destinationAt("Col4",&float4);
  hep.destinationAt("colAt",&int1);
  hep.destinationAt("colAt2",&int2);
  hep.destinationAt("colArray",&(arr[0][0]));


  // now readColumn back and compare !!
  for (i=0; i<rowNumber; i++) {
    // cerr << i << ":";

    hep.readRow(i+2); // the first row is "1"

    if ( ! IS_DBL_EQUAL( float1 ,3.2+i ) ) {
      cerr << "Col1 wrong " << float1 - (3.2+i) << endl;
    }
    if ( ! IS_DBL_EQUAL( float2 ,8.2+2*i ) ) {
      cerr << "Col2 wrong " << float2 - ( 8.2+2*i) << endl;
    }
    if ( ! IS_DBL_EQUAL( float3 ,99.0 ) ) {
      cerr << "Col3 wrong " << float3 << endl;
    }
    if ( ! IS_DBL_EQUAL( float4 ,99.0 ) ) {
      cerr << "Col4 wrong " << float4 << endl;
    }
    if ( int1 != -3*i ) {
      cerr << "colAt wrong " << int1 << endl;
    }
    if ( int2 != 55 ) {
      cerr << "colAt2 wrong " << int2 << endl;
    }
    if ( c!= i ) {
      cerr << "colF wrong " << c << endl;
    }
#ifdef USEMETHOD_WORKS
    if ( object.a() != 5*i ) {
      cerr << "colM wrong " << object.a() << endl;
    }
#endif
    for (int j=0; j < dim1; j++) {
      for (int k=0; k< dim2; k++) {
	if ( ! IS_DBL_EQUAL( arr[j][k], 8.2/(i+1)*j*k) ) {
	  cerr << "colArray wrong " << j << "," << k << ":" <<::fabs(arr[j][k] - (8.2/(i+1)*j*k)) /(8.2/(i+1)*j*k) << ";" << endl;
	}
      }
    }
  }
  hep.release();
  cerr << "\n";

  return 0;
}

HepHist1D& addHist1D( HepFileManager* mgr ) {

  HepHist1D& h = mgr->hist1D(name1D, nBins, low, high);
  cout << "Histogram " << name1D << " was created. " << endl;

  float inc = (high - low)/(float)nBins;

  ::srand(seed1d);

  // iterate through a loop that generates data
  // the number of loop iterations may not be correct
  int bob = 0; // temporary counter
  for (float i = low - inc; i <= high; i += inc/2) {
    // cout << "\tDatapoint: " << i << "     \t";
    for (int j = 0; j < ( ::rand() % RANDLIMIT); j++) {
      h.accumulate(i+inc/100);
      //      cout << ".";
      bob++;
    }
    // cerr << bob << endl; // " from " << r << endl;
    bob = 0;
  }

  return h;
}

#define ERROR( h, orig, what, val1, val2 )				\
cerr << "Error in " << h.manager()->fileName() << "::" << h.title() 	\
     << ",copy of " << orig.manager()->fileName() << ", for " << what  	\
     << ": " << val1 << " instead of " << val2 << endl;

#define COMP_HIST( h, orig, descript, funct ) 		\
if ( h.funct != orig.funct ) {				\
   ERROR ( h, orig, descript, h.funct, orig.funct );	\
   return false;					\
}


bool Compare(const HepHist1D& h1, const  HepHist1D& h2) {
  cerr << "Check Hist1D ... " << endl;

  if (h1.nBins() != h2.nBins()) {
    ERROR( h1, h2, "number of bins", h1.nBins(), h2.nBins() );
    return false;
  }
  bool anyFalse = false;
  for (int i = 0; i <= h1.nBins() + 1; i++) {
    //cout << "Bin " << i << ": ";
    if (IS_FLT_EQUAL(h1.bin(i),h2.bin(i))) {
      //cout << "Equal" << endl;
    }
    else {
      ERROR( h1, h2, "value of bin #" << i , h1.bin(i), h2.bin(i) );
      anyFalse = true;
    }
  }
  if (anyFalse) {
    cout << "The two histograms are *NOT* equal. " << endl;
    return false;
  }
  else {
    // cout << "The two histograms are equal. " << endl;
    return true;
  }
}

bool Check( const HepHist1D& h ) {
  cerr << "Check Hist1D ... " << endl;

  ::srand(seed1d);
  float low = h.min();
  float high = h.max();
  int nBins = h.nBins();

  float inc = (high - low)/(float)nBins;

  // iterate through a loop that generates data
  // the number of loop iterations may not be correct
  int bob = 0; // temporary counter
  float i;
  int bin = 0;
  bool next = false;
  for (i = low - inc; i <= high; i += inc/2) {
    for (int j = 0; j < ( ::rand() % RANDLIMIT); j++) {
      bob++;
    }
    if (next) {
      if ( bob != h.bin(  bin ) ) {
	cerr << "Error in bin #" << static_cast<int>( bin ) << " for " << h.title()
	     << " in " << h.manager()->fileName()
	     << " found " << h.bin(bin) << " instead of " << bob
	     << endl;
      }
      bin ++;
      bob = 0;
      next = false;
    } else {
      next = true;
    }
  }
  return true;
}

HepHist2D& addHist2D( HepFileManager* mgr ) {

  HepHist2D& h = mgr->hist2D(name2D, nBinsX, lowX, highX,
			     nBinsY, lowY, highY);
  cout << "Histogram " << name2D << " was created. " << endl;

  float lowX, highX, lowY, highY;
  int rangeX, rangeY;

  // Gather data from histogram
  lowX = h.minX();
  highX = h.maxX();
  rangeX = h.nBinsX();
  lowY = h.minY();
  highY = h.maxY();
  rangeY = h.nBinsY();

  float incX = (highX - lowX)/(float)rangeX;
  float incY = (highY - lowY)/(float)rangeY;

  ::srand(seed2d);

  // iterate through a loop that generates data
  // the number of loop iterations may not be correct
  int bob = 0; // temporary counter
  float offsetX = incX; // change range of loops
  float offsetY = incY;
  for (float i = lowX - offsetX*.98; i <= highX + offsetX; i += incX) {
    for (float k = lowY - offsetY*.98; k <= highY + offsetY; k += incY) {
      for (int j = 0; j < ( ::rand() % RANDLIMIT); j++) {
	h.accumulate(i+incX/100, k+incY/100);
	bob++;
      }
      bob = 0;
    }
  }

  return h;
}

bool Compare(const HepHist2D& h1, const  HepHist2D& h2) {
  cerr << "Check Hist2D ... " << endl;

  COMP_HIST( h1, h2, "number of X bins", nBinsX() );
  COMP_HIST( h1, h2, "number of Y bins", nBinsY() );

  bool anyFalse = false;
  for (int i = 0; i <= h1.nBinsX() + 1; i++) {
    for (int j = 0; j <= h1.nBinsY() + 1; j++) {
      //cout << "Bin " << i << ": ";
      if (IS_FLT_EQUAL(h1.bin(i,j),h2.bin(i,j))) {
	//cout << "Equal" << endl;
      }
      else {
	ERROR( h1, h2, "value of bin #(" << i << "," << j << ")" , h1.bin(i,j), h2.bin(i,j) );
	//cout << "***Not Equal***" << endl;
	anyFalse = true;
      }
    }
  }
  if (anyFalse) {
    cout << "The two histograms are *NOT* equal. " << endl;
    return false;
  }
  else {
    // cout << "The two histograms are equal. " << endl;
    return true;
  }
}

bool Check( const HepHist2D& h ) {
  cerr << "Check Hist2D ... " << endl;

  float lowX, highX, lowY, highY;
  int rangeX, rangeY;

  // Gather data from histogram
  lowX = h.minX();
  highX = h.maxX();
  rangeX = h.nBinsX();
  lowY = h.minY();
  highY = h.maxY();
  rangeY = h.nBinsY();

  float incX = (highX - lowX)/(float)rangeX;
  float incY = (highY - lowY)/(float)rangeY;

  ::srand(seed2d);

  // iterate through a loop that generates data
  // the number of loop iterations may not be correct
  int bob = 0; // temporary counter
  int binX = 0;
  int binY = 0;
  float offsetX = incX; // change range of loops
  float offsetY = incY;
  for (float i = lowX - offsetX*.98; i <= highX + offsetX; i += incX) {
    for (float k = lowY - offsetY*.98; k <= highY + offsetY; k += incY) {
      for (int j = 0; j < ( ::rand() % RANDLIMIT); j++) {
	bob++;
      }
      if ( bob != h.bin(binX,binY) ) {
	cerr << "Error in bin #(" << binX << "," << binY << ") for " << h.title()
	     << " in " << h.manager()->fileName()
	     << " found " << h.bin(binX,binY) << " instead of " << bob
	     << endl;
      }
      binY ++;
      bob = 0;
    }
    binX++;
    binY=0;
  }
  return true;
}

HepHistProf& addHistProf( HepFileManager* mgr ) {

  HepHistProf& h = mgr->histProf(nameProf, nBinsX, lowX, highX,
				 lowY, highY);
  if ( ! h.isValid() ) return h;
  cout << "Histogram " << nameProf << " was created. " << endl;

  float lowX, highX, lowY, highY;
  int rangeX;

  // Gather data from histogram
  lowX = h.minX();
  highX = h.maxX();
  rangeX = h.nBinsX();
  lowY = h.minY();
  highY = h.maxY();

  float incX = (highX - lowX)/(float)rangeX;
  float incY = (highY - lowY)/(float)rangeX;

  ::srand(seedp);

  // iterate through a loop that generates data
  // the number of loop iterations may not be correct
  int bob = 0; // temporary counter
  float offsetX = incX; // change range of loops
  float offsetY = incY;
#ifdef BUG_IN_CDFSGA_CERN_VERSION_97a
  int n = 0;
#endif
  for (float i = lowX - offsetX*.98; i <= highX + offsetX; i += incX) {
    for (float k = lowY - offsetY*.98; k <= highY + offsetY; k += incY) {
      //      std::cerr << "accu " << i+incX/100 << "," << k+incY/100;
      for (int j = 0; j < ( ::rand() % RANDLIMIT); j++) {
	h.accumulate(i+incX/100, k+incY/100);
	bob++;
      }
      //      std::cerr << " : " << bob << std::endl;
      bob = 0;
    }
#ifdef BUG_IN_CDFSGA_CERN_VERSION_97a
    cerr << "After i==" << i << std::endl;
    cerr << "\tbin# " << n-1 << " : " << h.bin(n-1) << "," << h.weight(n-1) << std::endl;
   cerr << "\tbin# " << n << " : " << h.bin(n) << "," << h.weight(n) << std::endl;
    cerr << "\tbin# " << n+1 << " : " << h.bin(n+1) << "," << h.weight(n+1) << std::endl;
    n++;
#endif
  }

  return h;
}

bool Compare(const HepHistProf& h1, const  HepHistProf& h2) {
  cerr << "Check HistProf ... " << endl;

  COMP_HIST( h1, h2, "number of X bins", nBinsX() );

  bool anyFalse = false;

  for (int i = 1; i <= h1.nBinsX(); i++) {
    if ( IS_FLT_EQUAL(h1.bin(i) ,h2.bin(i)) ) {
      //cout << "Equal" << endl;
    }
    else {
      ERROR( h1, h2, "value of bin #" << i  , h1.bin(i), h2.bin(i) );
      //cout << "***Not Equal***" << endl;
      anyFalse = true;
    }
  }
  if (anyFalse) {
    cout << "The two histograms are *NOT* equal. " << endl;
    return false;
  }
  else {
    // cout << "The two histograms are equal. " << endl;
    return true;
  }
}

bool Check( const HepHistProf& h ) {
  cerr << "Check HistProf ... " << endl;

  float lowX, highX, lowY, highY;
  int rangeX, rangeY;

  // Gather data from histogram
  lowX = h.minX();
  highX = h.maxX();
  rangeX = h.nBinsX();
  lowY = h.minY();
  highY = h.maxY();
  rangeY = rangeX;

  float incX = (highX - lowX)/(float)rangeX;
  float incY = (highY - lowY)/(float)rangeY;

  ::srand(seedp);

  // iterate through a loop that generates data
  // the number of loop iterations may not be correct
  int bob = 0; // temporary counter
  int binX = 0;
  float offsetX = incX; // change range of loops
  float offsetY = incY;
  float sum = 0.0;
  for (float i = lowX - offsetX*.98; i <= highX + offsetX; i += incX) {
    for (float k = lowY - offsetY*.98; k <= highY + offsetY; k += incY) {
      for (int j = 0; j < ( ::rand() % RANDLIMIT); j++) {
	if ( k>lowY && k<highY) {
	  sum += k;
	  bob++;
	}
      }
    }
    if (binX != 0 && binX != (rangeX+1)) {
      if ( bob != h.weight(binX) ) {
	cerr << "Error in bin #" << binX << " for " << h.title()
	     << " in " << h.manager()->fileName()
	     << " found " << h.weight(binX) << " instead of " << bob
	     << " bin number is : " << binX << " " << h.bin(binX)
	     << " " << sum/bob
	     << endl;
      }
      bob = 0;
      sum = 0.0;
    } else {
      sum = 0.0;
      bob = 0;
    }
    binX++;
  }
  return true;
}