Index: trunk/SophyaLib/NTools/perandom.cc
===================================================================
--- trunk/SophyaLib/NTools/perandom.cc	(revision 3097)
+++ trunk/SophyaLib/NTools/perandom.cc	(revision 3109)
@@ -16,89 +16,111 @@
 
 /********* 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.
+/*! 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)
+FunRan::FunRan(r_8 *tab, int_4 nBin, bool pdf)
 : 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
+ 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)
+FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf)
 : 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(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<=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);
-
+  : Histo(h)
+{
+ if(mBins<=1)
+   throw RangeCheckError("FunRan::FunRan less than 2 bins requested");
+ create_DF(pdf);
+}
+
+/********* 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.)
    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);
 }
@@ -107,5 +129,5 @@
 /*! 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])
+to the random uniform number (in [0,1[)
 */
 int_4 FunRan::BinRandom()
@@ -123,11 +145,25 @@
  r_8 z=drand01();
  int ibin = mBins-1;
- for(int_4 i=0;i<mBins;i++)
-   if(z<(*this)(i)) {ibin=i; break;}
-
+ 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;
+ 
+}
 
 ////////////////////////////////////////////////////////////////////////////
@@ -219,15 +255,29 @@
 /////////////////////////////////////////////////////////////////
 /*
---- 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))
+**** 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 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
+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:
@@ -238,9 +288,16 @@
          = 1 / 2Pi  (distribution uniforme sur [0,2Pi[)
   Les variables aleatoires r,t sont independantes:
-    g(r,t) = g(r) g(t) 
+    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" ainsi definie verifie:
-              <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 sig^2
+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 sig = s/sqrt(2) comme argument
-*/
+     <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument
+
+*/
