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

#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 *********/
/*! 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.
*/
FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin)
: 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);
}

/********* Methode *********/
/*! 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.
*/
FunRan::FunRan(r_8 *tab, int_4 nBin)
: 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);
}

/********* Methode *********/
/*! 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.
*/
FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax)
: 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);
}

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

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

                   

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



/////////////////////////////////////////////////////////////////
/*
--- Remarque sur complex< r_8 > ComplexGaussRan(double sig)
x = r cos(t) tire gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
y = r sin(t) tire 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 cherche la pdf g(r,t) du module et de la phase
  (r=sqrt(x^2+y^2,t=atan2(y,x)) --> (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"
  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) 
- Attention:
La variable complexe "c=x+iy" ainsi definie verifie:
              <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 sig^2
Si on veut generer une variable complexe gaussienne telle que
     <c c*> = s^2 alors il faut sig = s/sqrt(2) comme argument
*/