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

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

/********* Methode *********/
/*! Constructeur par defaut */
HistoErr::HistoErr(void)
: xmin_(1.), xmax_(-1.), nx_(0), dx_(0.)
, mCorrel(0)
{
}

/********* Methode *********/
/*! Constructeur d'un histo */
HistoErr::HistoErr(r_8 xmin,r_8 xmax,int_4 nx)
: mCorrel(0)
{
  CreateOrResize(xmin,xmax,nx);
}

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

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

/********* Methode *********/
/*! Gestion de l'allocation */
void HistoErr::CreateOrResize(r_8 xmin,r_8 xmax,int_4 nx)
{
  xmin_ = xmin; xmax_ = xmax; nx_ = nx; dx_=0.;
  if(nx_>0) {
    data_.ReSize(nx_);  data_ = 0.;
    err2_.ReSize(nx_);  err2_ = 0.;
    ndata_.ReSize(nx_); ndata_ = 0.;
    dx_ = (xmax_-xmin_)/nx_;
  }
  mCorrel = 0;
}

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

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

/********* 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 HistoErr::ToCorrel(void)
{
 if(nx_<1) return;
 mCorrel++;
 for(int_4 i=0;i<nx_;i++) {
   if(ndata_(i)<1.) continue;
   data_(i) /= ndata_(i);
   err2_(i) /= ndata_(i);
 }
 return;
}

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

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

 Zero();

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

}

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

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

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

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

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

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

// Lecture entete
is.GetStr(strg);

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

// Lecture des parametres HistoErr
is.Get(dobj->xmin_);
is.Get(dobj->xmax_);
is.Get(dobj->nx_);
is.Get(dobj->dx_);

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

return;
}

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

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

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

// Ecriture des parametres HistoErr
os.Put(dobj->xmin_);
os.Put(dobj->xmax_);
os.Put(dobj->nx_);
os.Put(dobj->dx_);

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

return;
}

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

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