source: Sophya/trunk/SophyaLib/NTools/perandom.cc@ 4033

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include "pexceptions.h"
#include "srandgen.h"
#include "perandom.h"


////////////////////////////////////////////////////////////////////////////
/*!
  \class SOPHYA::FunRan
  \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
*/
FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf)
  : Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = f(BinCenter(i));
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator from a ClassFunc
See FunRan::FunRan(FunRan::Func...) for further comments.
*/
FunRan::FunRan(ClassFunc& f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf)
  : Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = f(BinCenter(i));
 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 FunRan::FunRan(FunRan::Func...) for further comments.
*/
FunRan::FunRan(r_8 *tab, int_4 nBin, bool pdf)
: Histo(-0.5,nBin-0.5,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab[i];
 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::FunRan(FunRan::Func...) for further comments.
*/
FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf)
: Histo(xMin,xMax,nBin)
{
 if(nBin<=1)
   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");
 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab[i];
 create_DF(pdf);
}

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

/********* Methode *********/
/*! Creator from a TVector
*/
FunRan::FunRan(TVector<r_8>& tab, r_8 xMin, r_8 xMax, bool pdf)
: Histo(xMin,xMax,tab.Size())
{
 if(tab.Size()<=1)
   throw RangeCheckError("TsFunRan::TsFunRan less than 2 bins requested");
 for(int_4 i=0;i<tab.Size();i++) (*this)(i) = tab(i);
 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 FunRan::FunRan(FunRan::Func...) for further comments.
*/
FunRan::FunRan(Histo &h, bool pdf)
  : Histo(h)
{
 if(mBins<=1)
   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");
 create_DF(pdf);
}

/********* Methode *********/
/*! Creator by copy */
FunRan::FunRan(const FunRan& fh)
  : Histo(fh)
{
}

/********* Methode *********/
/*! Creator by default */
FunRan::FunRan(void)
{
}

/********* Methode *********/
void FunRan::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<<"FunRan::FunRan(Histo) not a distribution function last bin is <=0"<<endl;
    throw RangeCheckError("FunRan::FunRan(Histo) not a distribution function last bin is <=0");
  }
 for(int_4 i=0;i<mBins;i++) (*this)(i) /= (*this)(mBins-1);
}

/********* 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 FunRan::BinRandom()
{
 // recherche du premier bin plus grand ou egal a z
 r_8 z=drand01();
 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 FunRan::Random()
{
 r_8 z=drand01();
 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 FunRan::RandomInterp(void)
{
 r_8 z=drand01();
 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;
 
}

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

/********* Methode *********/
/*! Creator for random from a table */
FunRan2D::FunRan2D(r_8 *tab, int_4 nBinX, int_4 nBinY)
{
  // Tirage en X, somme sur les Y.
   r_8* tabX = new r_8[nBinX];
   for (int_4 i=0; i<nBinX; i++) {
     tabX[i] = 0;
     for (int_4 j=0; j<nBinY; j++) {
       tabX[i] += tab[i*nBinY +j];
     }
   }
   ranX = new FunRan(tabX, nBinX);
   delete[] tabX;
   
   ranY = new FunRan* [nBinX];
   
   for (int_4 k=0; k<nBinX; k++)
      ranY[k] = new FunRan(tab + nBinY*k, nBinY);
   
   nx = nBinX;
}

/********* Methode *********/
/*! Creator for random from a table */
FunRan2D::FunRan2D(r_8 **tab, int_4 nBinX, int_4 nBinY)
{
  // Tirage en X, somme sur les Y.
   r_8* tabX = new r_8[nBinX];
   for (int_4 i=0; i<nBinX; i++) {
     tabX[i] = 0;
     for (int_4 j=0; j<nBinY; j++) {
       tabX[i] += tab[i][j];
     }
   }
   ranX = new FunRan(tabX, nBinX);
   
   ranY = new FunRan* [nBinX];
   
   for (int_4 k=0; k<nBinX; k++)
    if (tabX[k] != 0) 
      ranY[k] = new FunRan(tab[k], nBinY);
    else ranY[k] = NULL;
   
   delete[] tabX;
   nx = nBinX;
}

/********* Methode *********/
/*! Destructor */
FunRan2D::~FunRan2D()
{
  for (int_4 i=nx-1; i>=0; i--)
    delete ranY[i];
    
  delete[] ranY;
  
  delete ranX;
}

/********* Methode *********/
/*! Tirage avec retour du numeros de bin. */
void FunRan2D::BinRandom(int_4& x, int_4& y)
{
  x = ranX->BinRandom();
  y = ranY[x]->BinRandom();
}

/********* Methode *********/
/*! Tirage avec retour abscisse et ordonnee du bin interpole. */
void FunRan2D::Random(r_8& x, r_8& y)
{
  x = ranX->Random();
  int_4 i = int_4(ceil(x));
  y = ranY[i]->Random();
}