source: Sophya/trunk/SophyaLib/HiStats/hist2err.cc@ 3597

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

namespace SOPHYA {

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

/********* Methode *********/
/*! Constructeur par defaut */
Histo2DErr::Histo2DErr(void)
: xmin_(1.), xmax_(-1.), dx_(0.) , ymin_(1.), ymax_(-1.), dy_(0.)
  , nx_(0), ny_(0), mMean(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)
: xmin_(1.), xmax_(-1.), dx_(0.) , ymin_(1.), ymax_(-1.), dy_(0.)
  , nx_(0), ny_(0), mMean(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_; 
 mMean = H.mMean;
}

/********* Methode *********/
/*! Destructeur */
Histo2DErr::~Histo2DErr(void)
{
 mMean = 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_;
  }
  mMean = 0;
}

/********* Methode *********/
/*!
  Remise a zero
*/
void Histo2DErr::Zero(void)
{
 if(nx_<=0 || ny_<=0) return;
 data_ = 0.;
 err2_ = 0.;
 ndata_ = 0.;
 mMean = 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
 WARNING: number of bins is kept, bin width is changed
*/
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 *********/
/*!
 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
 WARNING: bin widths are kept, numbers of bins are increased by 1
*/
void Histo2DErr::ReCenterBinW(void)
{
 CreateOrResize(xmin_-dx_/2.,xmax_+dx_/2.,nx_+1,ymin_-dy_/2.,ymax_+dy_/2.,ny_+1);
}

/********* Methode *********/
/*!
  Compute the mean histogram.
  Each bin content is divided by the number of entries in the bin.
  Each squared error is divided by the number of entries in the bin.
  The number of entries by bin is NOT set to 1
  (calling ToMean many time will change the histogram !)
*/
void Histo2DErr::ToMean(void)
{
 if(nx_<1 || ny_<1) return;
 mMean++;
 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 after ToMean action
*/
void Histo2DErr::FromMean(void)
{
 if(nx_<1 || ny_<1) return;
 mMean--;
 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 *********/
/*!
  Compute the mean histogram and replace the "error table" by the variance.
  This should be done if Add(x,w,w) has been used.
  The "value table" is divided by the number of entries to get the mean
  The "error table" is replace by the variance
  The number of entries by bin is NOT set to 1 
  (calling ToMean many time will change the histogram !)
  Mixing ToMean and ToVariance leads to unpredictable results
*/
void Histo2DErr::ToVariance(void)
{
 if(nx_<1 || ny_<1) return;
 mMean++;
 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) = err2_(i,j)/ndata_(i,j) - data_(i,j)*data_(i,j);
   }
 }
 return;
}

/********* Methode *********/
/*!
 Recompute back the original HistoErr after ToVariance action
  Mixing FromMean and FromVariance leads to unpredictable results
*/
void Histo2DErr::FromVariance(void)
{
 if(nx_<1 || ny_<1) return;
 mMean--;
 for(int_4 i=0;i<nx_;i++) {
   for(int_4 j=0;j<ny_;j++) {
     if(ndata_(i,j)<1.) continue;
     err2_(i,j) = ndata_(i,j)*(err2_(i,j) + data_(i,j)*data_(i,j));
     data_(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);
   }
 }
 mMean = hfrom.mMean;

}

/********* 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_;
  mMean = h.mMean;
  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(nmean="<<mMean<<")"<<endl
     <<"          nx="<<nx_<<" ["<<xmin_<<","<<xmax_<<"] dx="<<dx_<<endl
     <<"          ny="<<ny_<<" ["<<ymin_<<","<<ymax_<<"] dy="<<dy_<<endl;
}

/********* Methode *********/
/*!
  Write to an ASCII file
*/
int Histo2DErr::WriteASCII(string fname)
{
  FILE *file = fopen(fname.c_str(),"w");
  if(file==NULL) {
    cout<<"Histo2DErr::WriteASCII_Error: error opening "<<fname<<endl;
    return -1;
  }

  if(NBinX()<=0 || NBinY()<=0) {
    cout<<"Histo2DErr::WriteASCII_Error: wrong number of bins"<<endl;
    return -2;
  }

  fprintf(file,"%ld %.17e %.17e %.17e %ld %.17e %.17e %.17e %d\n"
         ,(long)NBinX(),XMin(),XMax(),WBinX()
         ,(long)NBinY(),YMin(),YMax(),WBinY()
         ,NMean());
  for(long i=0;i<NBinX();i++) for(long j=0;j<NBinY();j++) {
    // ligne = i*NY+j
    fprintf(file,"%ld %ld %.17e %.17e %.0f\n"
           ,i,j,(*this)(i,j),Error2(i,j),NEntBin(i,j));
  }

  fclose(file);
  return 0;
}

/*!
  Read from an ASCII file
*/
#define __LENLINE_Histo2DErr_ReadASCII__ 2048
int Histo2DErr::ReadASCII(string fname)
{
  FILE *file = fopen(fname.c_str(),"r");
  if(file==NULL) {
    cout<<"Histo2DErr::ReadASCII_Error: error opening "<<fname<<endl;
    return -1;
  }

  char line[__LENLINE_Histo2DErr_ReadASCII__];
  long n=0, nbinx=0, nbiny=0;

  while ( fgets(line,__LENLINE_Histo2DErr_ReadASCII__,file) != NULL ) {

    if(n==0) {

      r_8 xmin,xmax,wx, ymin,ymax,wy; long mnmean=1;
      sscanf(line,"%ld %lf %lf %lf %ld %lf %lf %lf %ld"
            ,&nbinx,&xmin,&xmax,&wx
            ,&nbiny,&ymin,&ymax,&wy
            ,&mnmean);
      if(nbinx<=0 || nbiny<=0) {
        cout<<"Histo2Err::ReadASCII_Error: wrong number of bins"<<endl;
        return -2;
      }
      CreateOrResize(xmin,xmax,nbinx,ymin,ymax,nbiny);
      SetMean(mnmean);

    } else {

      long i,j; r_8 v,e2,nb;
      sscanf(line,"%ld %ld %lf %lf %lf",&i,&j,&v,&e2,&nb);
      SetBin(i,j,v);
      SetErr2(i,j,e2);
      SetNentB(i,j,nb);

    }

    n++;
  }

  fclose(file);
  return 0;
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   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 ToMean/FromMean
is.Get(dobj->mMean);

// 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 ToMean/FromMean
os.Put(dobj->mMean);

// 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 ObjFileIO<Histo2DErr>;
#endif

} // FIN namespace SOPHYA