#include <stdio.h>
#include <stdlib.h>
#include "sopnamsp.h"
#include "math.h"
#include <iostream>
#include <string>

// Pour le test 
#include "histinit.h"
#include "histos.h"
#include "histos2.h"
#include "ntuple.h"


#include "fmath.h"
#include "srandgen.h"
#include "tarray.h"
#include "arrctcast.h"
#include "fioarr.h"
#include "timing.h"

#include <typeinfo> 
#include <datatype.h> 

/* ------------------------------------------------------------------
  Programme de test des fonctionalites de persistance (PPF) de SOPHYA 
  (C) LAL/IN2P3-CNRS, Univ. Paris Sud,   IRFU-CEA   1999-2010
  SOPHYA::PPersist .... 
  Utilisation:
  csh> tobjio --> mini help 
  csh> tobjio A 
  csh> tobjio B 
*/

static int ppftstA();
static int ppftstB();

/* -- Main program -- */
int main(int narg, char* arg[])
{

int rc = 0;
if (narg < 2) {
  cout << " Usage: tobjio A/B \n" 
       << "   A : Create tobjio.ppf file containing Histo/NTuple/Array \n" 
       << "   B : Check PPF I/O on TArray<T> \n" << endl;
  return 1;
}
string opt = arg[1]; 
try {
  SophyaInit();
  InitTim();
  cout << "Hash check, OMatrix : " << hex << PIOPersist::Hash("OMatrix") << dec << endl;
  if (opt=="A")  rc = ppftstA();
  else rc = ppftstB();
}
catch (PThrowable & exc) {
  cerr << " tobjio.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
       << " - Msg= " << exc.Msg() << endl;
  rc = 99;
 }
catch (std::exception & e) {
  cerr << " tobjio.cc: Catched std::exception what()=" << e.what() << endl;
  rc=98;
 }
catch (...) {
  cerr << " tobjio.cc: some other exception (...) was caught ! " << endl;
  rc = 97;
}
PrtTim("End tobjio " );
cout << " ---- Programme tobjio.cc-  FIN  (Rc=" << rc << ") --- " << endl;
return rc;

}




static sa_size_t NX = 15;
static sa_size_t NY = 20;

/* --Fonction-- */
template <class T>
int ckppf_arr(T x) 
{
  int rc=0;
  cout << " ---- ckppf_arr ----- T= " << typeid(x).name() << " , " 
       << DataTypeInfo<T>::getTypeName() << endl;
  TArray<T> a(10,12);
  a = RegularSequence();
  {
    POutPersist po("arrtobj.ppf");
    po << a;
  }
  TArray<T> b;
  {
    PInPersist pin("arrtobj.ppf");
    pin >> b;
  }
  TArray<T> d = a-b;
  T min, max;
  d.MinMax(min,max);
  cout << "  diff d=a-b, Min=" << min << " Max=" << max << endl;
  if ((min != T(0))||(max != T(0)))  rc=9;
  if (rc==0) 
    cout << " ===>>> PPF I/O TArray< " << DataTypeInfo<T>::getTypeName() << " > OK " << endl;
  else 
    cout << " ===>>> PPF I/O TArray< " << DataTypeInfo<T>::getTypeName() << " > ERROR !!!! " << endl;
  return rc;
}

/* --Fonction-- */
template <class T>
int ckppf_arrZ(T x) 
{
  int rc=0;
  cout << " ---- ckppf_arrZ ----- complex< T= " << typeid(x).name() << " , " 
       << DataTypeInfo<T>::getTypeName() << " > " << endl;
  TArray< complex<T> > a(10,12);
  TArray<T> rea = SDRealPart(a);
  TArray<T> ima = SDImagPart(a);
  
  rea = RegularSequence();
  ima = RegularSequence(0.5,2.);
  {
    POutPersist po("arrtobj.ppf");
    po << a;
  }
  TArray< complex<T> > b;
  {
    PInPersist pin("arrtobj.ppf");
    pin >> b;
  }
  TArray< complex<T> > d = a-b;
  T remin, remax;
  SDRealPart(d).MinMax(remin, remax);
  T immin, immax;
  SDImagPart(d).MinMax(immin, immax);

  cout << "  diff real[d=a-b], Min=" << remin << " Max=" << remax << endl;
  cout << "  diff imag[d=a-b], Min=" << immin << " Max=" << immax << endl;
  if ((remin != T(0))||(remax != T(0)) ||
      (immin != T(0))||(immax != T(0)) )  rc=9;
  if (rc==0) 
    cout << " ===>>> PPF I/O TArray< complex<" << DataTypeInfo<T>::getTypeName() << " > > OK " << endl;
  else 
    cout << " ===>>> PPF I/O TArray< complex< " << DataTypeInfo<T>::getTypeName() << " > > ERROR !!!! " << endl;
  return rc;
}

/* --Fonction-- */
int ppftstB()
{
  int rcg=0;
  int rc;
  int_4 i=1;
  int_8 l=1;
  r_4 f=1.;
  r_8 d=1.;
  rc=ckppf_arr(i);
  if (rc!=0) rcg+=10;
  rc=ckppf_arr(l);
  if (rc!=0) rcg+=10;
  rc=ckppf_arr(f);
  if (rc!=0) rcg+=10;
  rc=ckppf_arr(d);
  if (rc!=0) rcg+=10;
  rc=ckppf_arrZ(f);
  if (rc!=0) rcg+=10;
  rc=ckppf_arrZ(d);
  if (rc!=0) rcg+=10;
#ifdef SO_LDBLE128
  r_16 ld=1.;
  rc=ckppf_arr(ld);
  if (rc!=0) rcg+=10;
  rc=ckppf_arrZ(ld);
  if (rc!=0) rcg+=10;
#endif
  return rcg;
}


/* --Fonction-- */
int ppftstA()
{
  cout << "Test Creation / PPersist Save NTuple, Histo, ..." << endl;

  POutPersist so("tobjio.ppf");

  NTuple *nt;
  float xnt[4];
  const char *ntn[4] = {"x","y","ex","ey"};
  int nent,i;
  nent = 1000;
  nt = new NTuple(2,ntn);  // Creation NTuple (AVEC new )
  for(i=0; i<nent; i++) {
    xnt[0] = i*50./nent; xnt[1] = 0.2*xnt[0]*xnt[0]+(NorRand()*xnt[0]/3.);
    nt->Fill(xnt);
  }

  string nom = "nt21";
  ObjFileIO<NTuple> oiont(nt);
  oiont.Write(so, nom);
  cout << "Writing " << nom << endl;

  string nx, ny, nz;
  string ex,ey,ez;
  nx = "x";  ny = "y";
  nent = 20;
  nt = new NTuple(4,ntn);  // Creation NTuple (AVEC new )
  for(i=0; i<nent; i++) {
    xnt[0] = i*50./nent; xnt[1] = 0.2*xnt[0]*xnt[0]+(NorRand()*xnt[0]/3.);
    xnt[2] = sqrt((double) xnt[0]); xnt[3] = xnt[0]*2;
    nt->Fill(xnt);
  }
  nom = "nt22";
  ObjFileIO<NTuple> oiont2(nt);
  oiont2.Write(so, nom);
  cout << "Writing " << nom << endl;

  {  //  Objet Histo-2D
  r_8 x,y,w;
  int i,j;
  Histo2D* h = new Histo2D(-10.,10.,25,-10.,10.,25);
  for(i=0; i<25; i++) for(j=0; j<25; j++) {
    h->BinCenter(i,j,x,y);
    w = 100.*exp(-0.5*(x*x/9. + y*y/9.));
    h->Add(x,y,w);
    w = 50.*exp(-0.5*((x+10.)*(x+10.)/9. + (y+10.)*(y+10.)/9.));
    h->Add(x,y,w);
    }
  nom = "h2d";
  ObjFileIO<Histo2D> oioh(h);
  oioh.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

  {
  Histo *h;
  float x;


  h = new Histo(0., 200., 100);  // Creation histo (AVEC new )
  for(i=0; i<100; i++)  // Remplissage d'histo
    { x = (2*i+1.);  h->Add(x, x*(200.-x)); }
  nom = "h1d";
  ObjFileIO<Histo> oioh(h);
  oioh.Write(so, nom);
  cout << "Writing " << nom << endl;
  
  }

  {   // Vecteurs 
  Vector* vpx, *vpy;
  vpx = new Vector(20);  vpy = new Vector(20);
  int i;   double x;
  for(i=0; i<20; i++) { x = 8.*i/20;  (*vpx)(i) = x;  (*vpy)(i) = 3*cos(x)+2*sin(x); }
  nom = "vx1";
  FIO_TArray<double>  oiov1(vpx);
  oiov1.Write(so, nom);
  cout << "Writing " << nom << endl;
  nom = "vy1";
  FIO_TArray<double>  oiov2(vpy);
  oiov2.Write(so, nom);
  cout << "Writing " << nom << endl;
  vpx = new Vector(30);  
  for(i=0; i<30; i++) { x = 8.*i/30+0.1;  (*vpx)(i) = sin(2.*x)/(2.*x); }
  nom = "vx2";
  FIO_TArray<double>  oiov3(vpx);
  oiov3.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

      
  {       // NTuple 3D
  float xnt[6];
  double ang;
  const char *ntn[6] = {"x","y","z","ex","ey","ez"};
  int nent,i;
      
  NTuple* nt = new NTuple(6,ntn);  // Creation NTuple (AVEC new )
  nent = 3000;
  for(i=0; i<nent; i++) {
  ang = 5.*6.28*(double)i/(double)nent;
  xnt[0] = cos(ang)*1.5+NorRand()*0.1; xnt[1] = sin(ang)*1.5+NorRand()*0.1;
  xnt[2] = ang/8.+NorRand()*0.05;   xnt[3] = xnt[4] = xnt[5] = 0.;
  nt->Fill(xnt);
  }

  nom = "nt31";
  ObjFileIO<NTuple> oiont(nt);
  oiont.Write(so, nom);
  cout << "Writing " << nom << endl;
  
  nt = new NTuple(6,ntn);  // Creation NTuple (AVEC new )
  nent = 100;
  for(i=0; i<nent; i++) {
    ang = 5.*6.28*(double)i/(double)nent;
    xnt[0] = cos(-ang)*1+NorRand()*0.0; xnt[1] = sin(-ang)*1.+NorRand()*0.0;
    xnt[2] = ang/8.+NorRand()*0.0;   xnt[3] = 0.1; xnt[4] = 0.15;  xnt[5] = 0.07;
    nt->Fill(xnt);
    }
  nom = "nt32";
  ObjFileIO<NTuple> oiont2(nt);
  oiont2.Write(so, nom);
  cout << "Writing " << nom << endl;
  }

  {    // Matrices 
  Matrix* mtx = new Matrix(50, 50);
  int i,j;
  double xp,yp,rp;
  for(i=0; i<50; i++)
     for(j=0; j<50; j++) { 
     xp = i-15;   yp = j-25;
     rp = (xp*xp)/16.+(yp*yp)/16.;
     (*mtx)(j,i) = 11.*exp(-rp);
     xp = i-35;   yp = j-15;
     rp = (xp*xp)/9.+(yp*yp)/9.;
     (*mtx)(j,i) += 3.*exp(-rp);
     xp = i-45;   yp = j-40;
     rp = (xp*xp)/10.+(yp*yp)/20.;
     (*mtx)(j,i) += 6.*exp(-rp);
     }
  nom = "mtx1";
  FIO_TArray<double> oiom(mtx);
  oiom.Write(so, nom);
  cout << "Writing " << nom << endl;
  }
  return 0;
}