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

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

////////////////////////////////////////////////////////////////////////////
//++
// Class        FunRan
// Lib  Outils++ 
// include      perandom.h
//
//      Tirage aleatoire sur un histogramme 1D.
//--

//++
FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin)
//
//      Createur. f is a probability density function (PDF).
// 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
// au centre du bin.
// Les valeurs retournees sont les valeurs du centre des bins.
// Si binw est la largeur du bin, les valeurs retournees
// vont de xmin+binw/2 a xmax-binw/2.
//--
: Histo(xMin,xMax,nBin)
{
 if(nBin<0)
   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");

 // On cree la fonction de distribution a partir de la PDF
 (*this)(0) = f(BinCenter(0));
 for(int_4 i=1;i<nBin;i++) (*this)(i) = (*this)(i-1) + f(BinCenter(i));

 if((*this)(nBin-1)<=0.)
   throw RangeCheckError("FunRan::FunRan not a distribution function last bin is <=0");

 for(int_4 i=0;i<nBin;i++) (*this)(i) /= (*this)(nBin-1);
}

//++
FunRan::FunRan(r_8 *tab, int_4 nBin)
//
//      Createur. tab is a probability density function
// The return random values will be between 0 and nBin-1
// See FunRan::FunRan(FunRan::Func...) for further comments.
//--
: Histo(-0.5,nBin-0.5,nBin)
{
 if(nBin<=0)
   throw RangeCheckError("FunRan::FunRan no bin in Histo");

 (*this)(0) = tab[0];
 for(int_4 i=1;i<nBin;i++) (*this)(i) = (*this)(i-1) + tab[i];

 if((*this)(nBin-1)<=0.)
   throw RangeCheckError("FunRan::FunRan not a distribution function last bin is <=0");

 for(int_4 i=0;i<nBin;i++) (*this)(i) /= (*this)(nBin-1);
}

//++
FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax)
//
//      Createur. tab is a probability density function
// 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.
//--
: Histo(xMin,xMax,nBin)
{
 if(nBin<=0)
   throw RangeCheckError("FunRan::FunRan no bin in Histo");

 (*this)(0) = tab[0];
 for(int_4 i=1;i<nBin;i++) (*this)(i) = (*this)(i-1) + tab[i];

 if((*this)(nBin-1)<=0.)
   throw RangeCheckError("FunRan::FunRan not a distribution function last bin is <=0");

 for(int_4 i=0;i<nBin;i++) (*this)(i) /= (*this)(nBin-1);
}

//++
//      Createur.
// If pdf=true, h is a probability density fonction.
// 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<=0)
   throw RangeCheckError("FunRan::FunRan(Histo) no bin in Histo");

 // On cree l'histo integre
 if(pdf)
   for(int_4 i=1;i<mBins;i++) (*this)(i) += (*this)(i-1);

 if((*this)(mBins-1)<=0.)
   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);
}

//++
int_4 FunRan::BinRandom()
//
//      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])
//--
{
 // 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;
}

//++
r_8 FunRan::Random()
//
//      Tirage avec retour abscisse du bin non interpole.
//--
{
 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);
}

                   

////////////////////////////////////////////////////////////////////////////
//++
// Class        FunRan2D
// Lib  Outils++ 
// include      perandom.h
//
//      Tirage aleatoire sur un histogramme 2D.
//--

//++
FunRan2D::FunRan2D(r_8 *tab, int_4 nBinX, int_4 nBinY)
//
//      Createur.
//--
{
  // 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;
}

//++
FunRan2D::FunRan2D(r_8 **tab, int_4 nBinX, int_4 nBinY)
//
//      Createur.
//--
{
  // 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;
}

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

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

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