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

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


////////////////////////////////////////////////////////////////////////////
/*!
  \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 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();
}



/////////////////////////////////////////////////////////////////
/*
**** Remarques sur complex< r_8 > ComplexGaussRan(double sig) ****

--- variables gaussiennes x,y independantes
x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2))
x,y independants --> pdf f(x,y) = f(x) f(y)
On a:
  <x>   = Integrate[x*f(x)]   = Mx
  <x^2> = Integrate[x^2*f(x)] = Mx^2 + Sx^2

--- On cherche la pdf g(r,t) du module et de la phase
  x = r cos(t) ,  y = r sin(t)
  r=sqrt(x^2+y^2 , t=atan2(y,x)
  (r,t) --> (x,y): le Jacobien = r

  g(r,t) = r f(x,y) = r f(x) f(y)
         = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2))

- Le cas general est complique
  (cf D.Pelat cours DEA "bruits et signaux" section 4.5)

- Cas ou "Mx = My = 0" et "Sx = Sy = S"
  c'est la pdf du module et de la phase d'un nombre complexe
     dont les parties reelles et imaginaires sont independantes
     et sont distribuees selon des gaussiennes de variance S^2
  g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2))
  La distribution de "r" est donc:
    g(r) = Integrate[g(r,t),{t,0,2Pi}]
         = r/S^2 exp(-r^2/(2 S^2))
  La distribution de "t" est donc:
    g(t) = Integrate[g(r,t),{r,0,Infinity}]
         = 1 / 2Pi  (distribution uniforme sur [0,2Pi[)
  Les variables aleatoires r,t sont independantes:
    g(r,t) = g(r) g(t)
On a:
  <r>   = Integrate[r*g(r)]   = sqrt(PI/2)*S
  <r^2> = Integrate[r^2*g(r)] = 2*S^2
  <r^3> = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3
  <r^4> = Integrate[r^4*g(r)] = 8*S^4

- Attention:
La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie:
              <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 S^2
Si on veut generer une variable complexe gaussienne telle que
     <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument

*/