source: Sophya/trunk/SophyaLib/HiStats/histerr.cc@ 3307

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

namespace SOPHYA {

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

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

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

/********* Methode *********/
/*! Constructeur par copie */
HistoErr::HistoErr(const HistoErr& H)
: mMean(H.mMean)
{
 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)
{
 mMean = 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_;
  }
  mMean = 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
 WARNING: number of bins is kept, bin width is changed
*/
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 *********/
/*!
 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
 WARNING: bin width is kept, number of bins is increased by 1
*/
void HistoErr::ReCenterBinW(void)
{
 if(nx_<=1) return;
 CreateOrResize(xmin_-dx_/2.,xmax_+dx_/2.,nx_+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 HistoErr::ToMean(void)
{
 if(nx_<1) return;
 mMean++;
 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 after ToMean action
*/
void HistoErr::FromMean(void)
{
 if(nx_<1) return;
 mMean--;
 for(int_4 i=0;i<nx_;i++) {
   if(ndata_(i)<1.) continue;
   data_(i) *= ndata_(i);
   err2_(i) *= ndata_(i);
 }
 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 HistoErr::ToVariance(void)
{
 if(nx_<1) return;
 mMean++;
 for(int_4 i=0;i<nx_;i++) {
   if(ndata_(i)<1.) continue;
   data_(i) /= ndata_(i);
   err2_(i) = err2_(i)/ndata_(i) - data_(i)*data_(i);
 }
 return;
}

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

}

/********* Methode *********/
/*!
  Return the sum of bin value
*/
double HistoErr::Sum(void)
{
 double s = 0.;
 for(int_4 i=0;i<nx_;i++) {
   if(ndata_(i)<1.) continue;
   s += data_(i);
 }
 return s;
}

/********* Methode *********/
/*!
  Return the sum of the bin value squared
*/
double HistoErr::Sum2(void)
{
 double s = 0.;
 for(int_4 i=0;i<nx_;i++) {
   if(ndata_(i)<1.) continue;
   s += data_(i)*data_(i);
 }
 return s;
}

/********* Methode *********/
/*!
  Return the sum of the number of entries
*/
double HistoErr::SumN(void)
{
 double s = 0.;
 for(int_4 i=0;i<nx_;i++) {
   if(ndata_(i)<1.) continue;
   s += ndata_(i);
 }
 return s;
}

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

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

  if(NBins()<=0) {
    cout<<"HistoErr::WriteASCII_Error: nbin= "<<NBins()<<endl;
    return -2;
  }

  fprintf(file,"%ld %.17e %.17e %d\n",(long)NBins(),XMin(),XMax(),NMean());
  for(long i=0;i<NBins();i++) {
      fprintf(file,"%ld %.17e %.17e %.17e %.0f\n"
             ,i,BinCenter(i),(*this)(i),Error2(i),NEntBin(i));
  }

  fclose(file);
  return 0;
}

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

  char line[__LENLINE_HistoErr_ReadASCII__];
  long n=0, nbin=0;

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

    if(n==0) {

      r_8 xmin,xmax; long mnmean=1;
      sscanf(line,"%ld %lf %lf %ld",&nbin,&xmin,&xmax,&mnmean);
      if(nbin<=0) {
        cout<<"HistoErr::ReadASCII_Error: nbin= "<<nbin<<endl;
        return -2;
      }
      CreateOrResize(xmin,xmax,nbin);
      SetMean(mnmean);

    } else {

      long i; r_8 x,v,e2,nb;
      sscanf(line,"%ld %lf %lf %lf %lf",&i,&x,&v,&e2,&nb);
      SetBin(i,v);
      SetErr2(i,e2);
      SetNentB(i,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<HistoErr>::ReadSelf(PInPersist& is)
{
string strg;

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

// Lecture entete
is.GetStr(strg);

// Nombre d'appels a ToMean/FromMean
is.Get(dobj->mMean);

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

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

} // FIN namespace SOPHYA