#include "sopnamsp.h"
#include "machdefs.h"
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "perrors.h"
#include "fioarr.h"
#include "hist2err.h"

/*!
  \class SOPHYA::Histo2DErr
  \ingroup HiStats
  Classe d'histogrammes 1D avec erreurs donnees par l'utilisateur
*/

/********* Methode *********/
/*! Constructeur par defaut */
Histo2DErr::Histo2DErr(void)
: xmin_(1.), xmax_(-1.), nx_(0), dx_(0.)
, ymin_(1.), ymax_(-1.), ny_(0), dy_(0.)
, mCorrel(0)
{
}

/********* Methode *********/
/*! Constructeur d'un histo */
Histo2DErr::Histo2DErr(r_8 xmin,r_8 xmax,int_4 nx,r_8 ymin,r_8 ymax,int_4 ny)
{
  CreateOrResize(xmin,xmax,nx,ymin,ymax,ny);
}

/********* Methode *********/
/*! Constructeur par copie */
Histo2DErr::Histo2DErr(const Histo2DErr& H)
: mCorrel(0)
{
 if(H.nx_<=0 || H.ny_<=0) return;
 CreateOrResize(H.xmin_,H.xmax_,H.nx_,H.ymin_,H.ymax_,H.ny_);
 data_  = H.data_; 
 err2_  = H.err2_; 
 ndata_ = H.ndata_; 
 mCorrel = H.mCorrel;
}

/********* Methode *********/
/*! Destructeur */
Histo2DErr::~Histo2DErr(void)
{
 mCorrel = 0;
}

/********* Methode *********/
/*! Gestion de l'allocation */
void Histo2DErr::CreateOrResize(r_8 xmin,r_8 xmax,int_4 nx,r_8 ymin,r_8 ymax,int_4 ny)
{
  xmin_ = xmin; xmax_ = xmax; nx_ = nx; dx_=0.;
  ymin_ = ymin; ymax_ = ymax; ny_ = ny; dy_=0.;
  if(nx_>0 && ny_>0) {
    data_.ReSize(nx_,ny_);  data_ = 0.;
    err2_.ReSize(nx_,ny_);  err2_ = 0.;
    ndata_.ReSize(nx_,ny_); ndata_ = 0.;
    dx_ = (xmax_-xmin_)/nx_;
    dy_ = (ymax_-ymin_)/ny_;
  }
  mCorrel = 0;
}

/********* Methode *********/
/*!
  Remise a zero
*/
void Histo2DErr::Zero(void)
{
 if(nx_<=0 || ny_<=0) return;
 data_ = 0.;
 err2_ = 0.;
 ndata_ = 0.;
}

/********* Methode *********/
/*!
 Recompute XMin (YMin) and XMax (YMax) so that
 the CENTER of the first bin is exactly XMin (YMin) and
 the CENTER of the last bin is exactly XMax (YMax).
 Remember that otherwise
 XMin (YMin) is the beginning of the first bin
 and XMax (YMax) is the end of the last bin
*/
void Histo2DErr::ReCenterBinX(void)
{
 if(nx_<=1) return;
 double dx = (xmax_-xmin_)/(nx_-1);
 xmin_ -= dx/2.;
 xmax_ += dx/2.;
 dx_ = (xmax_-xmin_)/nx_;
}

void Histo2DErr::ReCenterBinY(void)
{
 if(ny_<=1) return;
 double dy = (ymax_-ymin_)/(ny_-1);
 ymin_ -= dy/2.;
 ymax_ += dy/2.;
 dy_ = (ymax_-ymin_)/ny_;
}

void Histo2DErr::ReCenterBin(void)
{
 ReCenterBinX();
 ReCenterBinY();
}

/********* Methode *********/
/*!
  Compute the correlation histogram.
  Each bin content is divided by the number of entries in that bin.
  Each squared error is divided by the number of entries in that bin.
  The number of entries by bin is NOT set to 1 (calling ToCorrel
    many time will change the histogram !)
*/
void Histo2DErr::ToCorrel(void)
{
 if(nx_<1 || ny_<1) return;
 mCorrel++;
 for(int_4 i=0;i<nx_;i++) {
   for(int_4 j=0;j<ny_;j++) {
     if(ndata_(i,j)<1.) continue;
     data_(i,j) /= ndata_(i,j);
     err2_(i,j) /= ndata_(i,j);
   }
 }
 return;
}

/********* Methode *********/
/*!
 Recompute back the original Histo2DErr before ToCorrel action
*/
void Histo2DErr::FromCorrel(void)
{
 if(nx_<1 || ny_<1) return;
 mCorrel--;
 for(int_4 i=0;i<nx_;i++) {
   for(int_4 j=0;j<ny_;j++) {
     if(ndata_(i,j)<1.) continue;
     data_(i,j) *= ndata_(i,j);
     err2_(i,j) *= ndata_(i,j);
   }
 }
 return;
}

/********* Methode *********/
/*!
 Fill the histogram with an other histogram
*/
void Histo2DErr::FillFrHErr(Histo2DErr& hfrom)
{
 if(nx_<=0 || ny_<=0) return;
 if(hfrom.nx_<=0 || hfrom.ny_<=0) return;

 Zero();

 for(int_4 i=0;i<hfrom.nx_;i++) {
   for(int_4 j=0;j<hfrom.ny_;j++) {
     r_8 x,y;  hfrom.BinCenter(i,j,x,y);
     int ii,jj; FindBin(x,y,ii,jj);
     if(jj<0 || jj>=ny_ || ii<0 || ii>=nx_) continue;
     data_(ii,jj)  += hfrom.data_(ii,jj);
     err2_(ii,jj)  += hfrom.err2_(ii,jj);
     ndata_(ii,jj) += hfrom.ndata_(ii,jj);
   }
 }
 mCorrel = hfrom.mCorrel;

}

/********* Methode *********/
/*!
  Operateur egal Histo2DErr = Histo2DErr
*/
Histo2DErr& Histo2DErr::operator = (const Histo2DErr& h)
{
  if(this==&h) return *this;
  CreateOrResize(h.xmin_,h.xmax_,h.nx_,h.ymin_,h.ymax_,h.ny_);
  data_  = h.data_;
  err2_  = h.err2_;
  ndata_ = h.ndata_;
  mCorrel = h.mCorrel;
  return *this;
}

/********* Methode *********/
/*!
  Operateur de multiplication par une constante
*/
Histo2DErr& Histo2DErr::operator *= (r_8 b)
{
r_8 b2 = b*b;
for(int_4 i=0;i<nx_;i++) {
  for(int_4 j=0;j<ny_;j++) {
    data_(i,j) *= b;
    err2_(i,j) *= b2;
  }
}
return *this;
}

/********* Methode *********/
/*!
  Print info
*/
void Histo2DErr::Show(ostream & os) const
{
  os <<"Histo2DErr(ncorrel="<<mCorrel<<")"<<endl
     <<"          nx="<<nx_<<" ["<<xmin_<<","<<xmax_<<"] dx="<<dx_<<endl
     <<"          ny="<<ny_<<" ["<<ymin_<<","<<ymax_<<"] dy="<<dy_<<endl;
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   Les objets delegues pour la gestion de persistance 
// --------------------------------------------------------
///////////////////////////////////////////////////////////

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<Histo2DErr>::ReadSelf(PInPersist& is)
{
string strg;

if(dobj==NULL) dobj = new Histo2DErr;

// Lecture entete
is.GetStr(strg);

// Nombre d'appels a ToCorrel/FromCorrel
is.Get(dobj->mCorrel);

// Lecture des parametres Histo2DErr
is.Get(dobj->xmin_);
is.Get(dobj->xmax_);
is.Get(dobj->nx_);
is.Get(dobj->dx_);
is.Get(dobj->ymin_);
is.Get(dobj->ymax_);
is.Get(dobj->ny_);
is.Get(dobj->dy_);

// Lecture des donnees
if(dobj->nx_>0 && dobj->ny_>0) {
  is >> dobj->data_;
  is >> dobj->err2_;
  is >> dobj->ndata_;
}

return;
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<Histo2DErr>::WriteSelf(POutPersist& os) const
{
if(dobj == NULL)   return;
string strg;

// Ecriture entete
strg = "Hist2DErr";
os.PutStr(strg);

// Nombre d'appels a ToCorrel/FromCorrel
os.Put(dobj->mCorrel);

// Ecriture des parametres Histo2DErr
os.Put(dobj->xmin_);
os.Put(dobj->xmax_);
os.Put(dobj->nx_);
os.Put(dobj->dx_);
os.Put(dobj->ymin_);
os.Put(dobj->ymax_);
os.Put(dobj->ny_);
os.Put(dobj->dy_);

// Ecriture des donnees
if(dobj->nx_>0 && dobj->ny_>0) {
  os << dobj->data_;
  os << dobj->err2_;
  os << dobj->ndata_;
}

return;
}

#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<Histo2DErr>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class SOPHYA::ObjFileIO<Histo2DErr>;
#endif