source: Sophya/trunk/SophyaLib/NTools/tsfunran.cc@ 3600

#include "sopnamsp.h"
#include "machdefs.h"
#include "pexceptions.h"
#include "tsfunran.h"
#include <iostream>


////////////////////////////////////////////////////////////////////////////
/*!
  \class SOPHYA::TsFunRan
  \ingroup NTools
  Classe for generating random variables from 1D function
*/


/********* Methode *********/
/*! Creator from a function.
\verbatim
 - if pdf==true: f est une densite de probabilite (PDF)
                 non necessairement normalisee
   if pdf==false: f est  une fonction de distribution (DF).
                  non necessairement normalisee
 - Le tirage aleatoire est fait sur un histogramme
   Histo(xMin,xMax,nBin) (voir convention dans Histo).
 - Chaque bin de l'histogramme contient la valeur de la PDF
   (ou de la DF) au centre du bin: h(i)=f(BinCenter(i))
 - Les valeurs retournees sont les valeurs du centre des bins
   pour le tirage non interpole et toutes les valeurs
   entre [xmin,xmax] pour le tirage interpole
 - La pdf doit etre interpretee comme etant nulle
   pour des x<=xmin et x>=xmax
 - Dans le bin "i" entre [x1,x2[ et de centre x0, h(i)=pdf(x0).
   Pour le tirage interpole, la DF est approximee par un segment
   et pdf(x0) est l'exces de proba entre x1 et x2:
   bin 0 entre [xmin,BinHighEdge(0)[ :
               la pdf va de 0 a pdf(BinCenter(0))
   bin 1 entre [BinLowEdge(1),BinHighEdge(1)[:
               la pdf va de pdf(BinCenter(0)) a pdf(BinCenter(1))
   ...
   bin n-1 entre [BinLowEdge(n-1),xmax[:
               la pdf va de pdf(BinCenter(n-2)) a pdf(BinCenter(n-1))
\endverbatim
*/
TsFunRan::TsFunRan(TsFunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf)
  : Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = f(BinCenter(i));
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from a ClassFunc
See TsFunRan::TsFunRan(TsFunRan::Func...) for further comments.
*/
TsFunRan::TsFunRan(ClassFunc& f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf)
  : Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = f(BinCenter(i));
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from an table.
If pdf=true, tab is a probability density fonction (not necessary normalised).
If pdf=false, tab is a distribution function (not necessarly normalized to 1).
The return random values will be between 0 and nBin-1.
See TsFunRan::TsFunRan(TsFunRan::Func...) for further comments.
*/
TsFunRan::TsFunRan(r_8 *tab, int_4 nBin, bool pdf)
: Histo(-0.5,nBin-0.5,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab[i];
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from an table.
If pdf=true, tab is a probability density fonction (not necessary normalised).
If pdf=false, tab is a distribution function (not necessarly normalized to 1).
The content of tab is identified has the content of
an Histogram define by Histo(xMin,xMax,nBin).
See FunRan::TsFunRan(TsFunRan::Func...) for further comments.
*/
TsFunRan::TsFunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf)
: Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab[i];
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from a TVector
*/
TsFunRan::TsFunRan(TVector<r_8>& tab, int_4 nBin, bool pdf)
: Histo(-0.5,nBin-0.5,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab(i);
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from a TVector
*/
TsFunRan::TsFunRan(TVector<r_8>& tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf)
: Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab(i);
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from an histogram
If pdf=true, h is a probability density fonction (not necessary normalised).
If pdf=false, h is a distribution function (not necessarly normalized to 1).
See TsFunRan::TsFunRan(TsFunRan::Func...) for further comments.
*/
TsFunRan::TsFunRan(Histo &h, bool pdf)
  : Histo(h)
{
 if(mBins<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 rg_ = new RandomGenerator(1,false);
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator by copy */
TsFunRan::TsFunRan(const TsFunRan& fh)
  : Histo(fh)
{
  rg_ = new RandomGenerator(*fh.rg_);
}

/********* Methode *********/
/*! Creator by default */
TsFunRan::TsFunRan(void)
{
  rg_ = new RandomGenerator(1,false);
}

/********* Methode *********/
TsFunRan::~TsFunRan(void)
{
  if(rg_!=NULL) delete rg_;
}

/********* Methode *********/
void TsFunRan::create_DF(bool pdf)
// Creation (si necessaire) et normalisation de la DF
{
  // On fabrique la FD si necessaire
  if(pdf)
    for(int_4 i=1;i<mBins;i++) (*this)(i) += (*this)(i-1);

  // On normalise la FD
  if((*this)(mBins-1)<=0.) {
    cout<<"TsFunRan::TsFunRan(Histo) not a distribution function last bin is <=0"<<endl;
    throw RangeCheckError("TsFunRan::TsFunRan(Histo) not a distribution function last bin is <=0");
  }
 for(int_4 i=0;i<mBins;i++) (*this)(i) /= (*this)(mBins-1);
}

/********* Methode *********/
void TsFunRan::SetRandomGenerator(RandomGenerator& rg)
{
  if(rg_!=NULL) delete rg_;
  rg_ = new RandomGenerator(rg);
}

/********* Methode *********/
/*! Tirage avec retour du numero de bin entre 0 et mBins-1.
It returns the first bin whose content is greater or equal
to the random uniform number (in [0,1[)
*/
int_4 TsFunRan::BinRandom()
{
 // recherche du premier bin plus grand ou egal a z
 r_8 z=rg_->Flat01();
 for(int_4 i=0;i<mBins;i++) if(z<(*this)(i)) return i;
 return mBins-1;
}

/********* Methode *********/
/*! Tirage avec retour abscisse du bin non interpole. */
r_8 TsFunRan::Random()
{
 r_8 z=rg_->Flat01();
 int ibin = mBins-1;
 for(int_4 i=0;i<mBins;i++) if(z<(*this)(i)) {ibin=i; break;}
 return BinCenter(ibin);
}

/********* Methode *********/
/*! Tirage avec retour abscisse du bin interpole. */
r_8 TsFunRan::RandomInterp(void)
{
 r_8 z=rg_->Flat01();
 int ibin = mBins-1;
 r_8 z1=0., z2;
 for(int_4 i=0;i<mBins;i++) {
   z2 = (*this)(i);
   if(z<z2) {ibin=i; break;}
   z1 = z2;
 }

 // l'algorithme garanti que "z2-z1 != 0" et "z1<z2"
 return BinLowEdge(ibin) + (z-z1)/(z2-z1)*binWidth;
 
}