Changeset 2840 in Sophya for trunk/SophyaLib/NTools/perandom.cc

Timestamp:

Nov 10, 2005, 5:53:42 PM (20 years ago)

Author:

cmv

Message:

Refonte totale de FunRan qui n'etait pas OK.
Arret de l'interpolation entre bin pour le tirage aleatoire
car ca pose un pb dans le premier et dernier bin.. A voir + tard
eventuellement. cmv 10/11/2005

File:

• trunk/SophyaLib/NTools/perandom.cc (modified) (5 diffs)

trunk/SophyaLib/NTools/perandom.cc

@@ -18 +18 @@
 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.
-//--
-: Histo(xMin, xMax, nBin)
-{
-  (*this)(0) = f(BinLowEdge(0));
-  for(int_4 i=1; i<nBin; i++)
-    (*this)(i) = (*this)(i-1) + f(BinLowEdge(i));
-    
-  for(int_4 j=0; j<nBin; j++)
-    (*this)(j) /= (*this)(nBin-1);
-}
-
-//++
-FunRan::FunRan(r_8 *tab, int_4 nBin)
-//
-//      Createur.
-//--
-: Histo(0, (r_8)(nBin), nBin)
-{
-  (*this)(0) = tab[0];
-  for(int_4 i=1; i<nBin; i++)
-    (*this)(i) = (*this)(i-1) + tab[i];
-
-  if ((*this)(nBin-1) == 0) {
-    cerr << "FunRan::FunRan : sum(prob) = 0" << endl;
-    exit(-1);
-  }
-
-  for(int_4 j=0; j<nBin; j++)
-    (*this)(j) /= (*this)(nBin-1);
-}
-
-FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax)
-: Histo(xMin, xMax, nBin)
-{
-  (*this)(0) = tab[0];
-  for(int_4 i=1; i<nBin; i++)
-    (*this)(i) = (*this)(i-1) + tab[i];
-
-  if ((*this)(nBin-1) == 0) {
-    cerr << "FunRan::FunRan : sum(prob) = 0" << endl;
-    exit(-1);
-  }
-
-  for(int_4 j=0; j<nBin; j++)
-    (*this)(j) /= (*this)(nBin-1);
-}
-
-FunRan::FunRan(Histo &h)
+// 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);
 }
 
@@ -74 +110 @@
 int_4 FunRan::BinRandom()
 //
-//      Tirage avec retour du numero de bin.
-//--
-{
-  r_8 z=drand01();
-  if (z <= 0) return 0;
-  if (z >= 1) return mBins-1;
-  
-  // recherche du premier bin plus grand que z
-  int_4 iBin = 0;
-  for (; iBin<mBins; iBin++)
-    if (z < (*this)(iBin)) break;
-
-  return iBin;
+//      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;
 }
 
@@ -92 +124 @@
 r_8 FunRan::Random()
 //
-//      Tirage avec retour abscisse du bin interpole.
-//--
-{
-  r_8 z=drand01();
-  if (z <= 0) return mMin;
-  if (z >= 1) return mMax;
-  // cas z <= tab[0] 
-  if (z <= (*this)(0)) {
-    r_8 t = mMin + binWidth/(*this)(0) * z;
-    return t;
-  }
-
-  // recherche du premier bin plus grand que z
-  int_4 iBin = 0;
-  for (; iBin<mBins; iBin++)
-    if (z < (*this)(iBin)) break;
-
-  // interpolation pour trouver la valeur du tirage aleatoire
-  r_8 t1 = (*this)(iBin-1);
-  r_8 x1 = BinLowEdge(iBin-1);
-  r_8 t2 = (*this)(iBin);
-  r_8 x2 = x1 + binWidth;
-  r_8 t = x1 + (x2-x1) / (t2-t1) * (z-t1);
-  if (t < mMin) t = mMin;
-  if (t > mMax) t = mMax;
-  return(t);
+//      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);
 }
 
@@ -200 +215 @@
 {
   x = ranX->BinRandom();
-  //  FAILNIL(ranY[x]);  Ne compile pas $CHECK$ Reza 22/04/99
   y = ranY[x]->BinRandom();
 }
@@ -213 +227 @@
   x = ranX->Random();
   int_4 i = int_4(ceil(x));
-  //  FAILNIL(ranY[i]);  Ne compile pas $CHECK$ Reza 22/04/99
   y = ranY[i]->Random();
 }