Changeset 4018 in Sophya for trunk/SophyaLib/BaseTools

Timestamp:

Sep 21, 2011, 6:21:17 PM (14 years ago)

Author:

cmv

Message:

• introduction du tirage gaussien par la methode de la ziggurat
• ajout argument slim par defaut dans GaussianTail pour proteger d'une reccursion infinie genere par le tirage gaussien ziggurat dans le cas ou GaussianTail applique la methode par rejection (sdev<slim)

cmv, 21/09/2011

Location:

trunk/SophyaLib/BaseTools

Files:

• randinterf.cc (modified) (3 diffs)
• randinterf.h (modified) (3 diffs)
• srandgen.h (modified) (1 diff)

trunk/SophyaLib/BaseTools/randinterf.cc

@@ -211 +211 @@
       return GaussianLevaRatioUnif();
       break;
+    case C_Gaussian_Ziggurat128 :
+      return GaussianZiggurat128();
+      break;
     default:
       return GaussianBoxMuller();
@@ -391 +394 @@
 }
 
-r_8 RandomGeneratorInterface::GaussianTail(double s)
+//-------------------------------------------------------------------------------
+// Definition des tableaux pour tirage Gaussien methode Ziggurat N=128 bandes
+// G.Marsaglia and W.W.Tsang "The ziggurat method for generating random variables
+// D.Thomas et al. ACM Computing Surveys, Vol 39, No 4, Article 11, October 2007
+// http://en.wikipedia.org/wiki/Ziggurat_algorithm  (tres bien explique)
+// Calcul des tableaux avec programme "cmvziggurat.cc"
+//-------------------------------------------------------------------------------
+const int N_ZIGGURRAT_GAUSS128 = 128;
+static const double X_ZIGGURRAT_GAUSS128[N_ZIGGURRAT_GAUSS128+1] = {
+0.0000000000000000, 0.2723208647046734, 0.3628714310284243, 0.4265479863033096, 0.4774378372537916, 
+0.5206560387251481, 0.5586921783755209, 0.5929629424419807, 0.6243585973090908, 0.6534786387150446, 
+0.6807479186459064, 0.7064796113136101, 0.7309119106218833, 0.7542306644345121, 0.7765839878761502, 
+0.7980920606262765, 0.8188539066833194, 0.8389522142812090, 0.8584568431780525, 0.8774274290977171, 
+0.8959153525662399, 0.9139652510088031, 0.9316161966013551, 0.9489026254979132, 0.9658550793881319, 
+0.9825008035027615, 0.9988642334806447, 1.0149673952393006, 1.0308302360564565, 1.0464709007525812, 
+1.0619059636836206, 1.0771506248819389, 1.0922188768965548, 1.1071236475235364, 1.1218769225722551, 
+1.1364898520030764, 1.1509728421389769, 1.1653356361550478, 1.1795873846544616, 1.1937367078237728, 
+1.2077917504067583, 1.2217602305309634, 1.2356494832544818, 1.2494664995643345, 1.2632179614460288, 
+1.2769102735517004, 1.2905495919178738, 1.3041418501204223, 1.3176927832013436, 1.3312079496576772, 
+1.3446927517457137, 1.3581524543224235, 1.3715922024197329, 1.3850170377251492, 1.3984319141236070, 
+1.4118417124397606, 1.4252512545068619, 1.4386653166774619, 1.4520886428822168, 1.4655259573357950, 
+1.4789819769830983, 1.4924614237746157, 1.5059690368565506, 1.5195095847593711, 1.5330878776675565, 
+1.5467087798535037, 1.5603772223598409, 1.5740982160167500, 1.5878768648844011, 1.6017183802152775, 
+1.6156280950371333, 1.6296114794646788, 1.6436741568569830, 1.6578219209482079, 1.6720607540918526, 
+1.6863968467734867, 1.7008366185643014, 1.7153867407081171, 1.7300541605582440, 1.7448461281083769, 
+1.7597702248942324, 1.7748343955807697, 1.7900469825946195, 1.8054167642140493, 1.8209529965910054, 
+1.8366654602533845, 1.8525645117230873, 1.8686611409895424, 1.8849670357028696, 1.9014946531003181, 
+1.9182573008597323, 1.9352692282919006, 1.9525457295488893, 1.9701032608497135, 1.9879595741230611, 
+2.0061338699589673, 2.0246469733729340, 2.0435215366506698, 2.0627822745039639, 2.0824562379877247, 
+2.1025731351849992, 2.1231657086697902, 2.1442701823562618, 2.1659267937448412, 2.1881804320720208, 
+2.2110814088747279, 2.2346863955870573, 2.2590595738653296, 2.2842740596736570, 2.3104136836950024, 
+2.3375752413355309, 2.3658713701139877, 2.3954342780074676, 2.4264206455302118, 2.4590181774083506, 
+2.4934545220919508, 2.5300096723854670, 2.5690336259216395, 2.6109722484286135, 2.6564064112581929, 
+2.7061135731187225, 2.7611693723841539, 2.8231253505459666, 2.8943440070186708, 2.9786962526450171, 
+3.0832288582142140, 3.2230849845786187, 3.4426198558966523, 3.7130862467403638
+};
+static const double Y_ZIGGURRAT_GAUSS128[N_ZIGGURRAT_GAUSS128+1] = {
+1.0000000000000000, 0.9635996931557651, 0.9362826817083690, 0.9130436479920363, 0.8922816508023012, 
+0.8732430489268521, 0.8555006078850629, 0.8387836053106459, 0.8229072113952607, 0.8077382946961199, 
+0.7931770117838580, 0.7791460859417020, 0.7655841739092348, 0.7524415591857027, 0.7396772436833371, 
+0.7272569183545049, 0.7151515074204761, 0.7033360990258165, 0.6917891434460349, 0.6804918410064135, 
+0.6694276673577053, 0.6585820000586529, 0.6479418211185500, 0.6374954773431442, 0.6272324852578138, 
+0.6171433708265618, 0.6072195366326042, 0.5974531509518116, 0.5878370544418199, 0.5783646811267017, 
+0.5690299910747210, 0.5598274127106941, 0.5507517931210546, 0.5417983550317235, 0.5329626593899870, 
+0.5242405726789923, 0.5156282382498716, 0.5071220510813041, 0.4987186354765838, 0.4904148252893212, 
+0.4822076463348383, 0.4740943006982492, 0.4660721526945706, 0.4581387162728716, 0.4502916436869266, 
+0.4425287152802462, 0.4348478302546615, 0.4272469983095620, 0.4197243320540379, 0.4122780401070242, 
+0.4049064208114880, 0.3976078564980422, 0.3903808082413892, 0.3832238110598833, 0.3761354695144541, 
+0.3691144536682749, 0.3621594953730330, 0.3552693848515469, 0.3484429675498723, 0.3416791412350135, 
+0.3349768533169710, 0.3283350983761522, 0.3217529158792085, 0.3152293880681574, 0.3087636380092518, 
+0.3023548277894796, 0.2960021568498557, 0.2897048604458103, 0.2834622082260124, 0.2772735029218976, 
+0.2711380791410251, 0.2650553022581618, 0.2590245673987105, 0.2530452985097656, 0.2471169475146965, 
+0.2412389935477511, 0.2354109422657275, 0.2296323252343025, 0.2239026993871337, 0.2182216465563704, 
+0.2125887730737359, 0.2070037094418736, 0.2014661100762031, 0.1959756531181102, 0.1905320403209136, 
+0.1851349970107133, 0.1797842721249620, 0.1744796383324022, 0.1692208922389246, 0.1640078546849276, 
+0.1588403711409350, 0.1537183122095865, 0.1486415742436969, 0.1436100800919329, 0.1386237799858510, 
+0.1336826525846476, 0.1287867061971039, 0.1239359802039817, 0.1191305467087185, 0.1143705124498882, 
+0.1096560210158177, 0.1049872554103545, 0.1003644410295455, 0.0957878491225781, 0.0912578008276347, 
+0.0867746718955429, 0.0823388982429574, 0.0779509825146547, 0.0736115018847548, 0.0693211173941802, 
+0.0650805852136318, 0.0608907703485663, 0.0567526634815385, 0.0526674019035031, 0.0486362958602840, 
+0.0446608622008724, 0.0407428680747906, 0.0368843887869688, 0.0330878861465051, 0.0293563174402538, 
+0.0256932919361496, 0.0221033046161116, 0.0185921027371658, 0.0151672980106720, 0.0118394786579823, 
+0.0086244844129305, 0.0055489952208165, 0.0026696290839025, 0.0000000000000000
+};
+
+
+r_8 RandomGeneratorInterface::GaussianZiggurat128()
+//--------
+// On a "N = 128" bandes horizontales numerotees [0,N-1=127]
+// Les tableaux ont une taille "N + 1 = 129"
+// On tire un numero de bande dans [0,N-1=127]
+//--------
+// Pour choisir le signe sans avoir a retirer un aleatoire,
+// on utilise un digit du premier tirage qui n'est pas utilise
+// dans le choix du numero de bande "I":
+// U = [0,1[ , Ai=[0,9] (chiffres)
+// U = A1/10 + A2/100 + A3/1000 + A4/10000 + A5/100000 + ...
+// 128*U = A1*12.8 + A2*1.28 + A3*0.128 + A4*0.0128 + A5*0.00128
+// pour Ai le plus grand possible cad 9
+//   128*U = 115.2 + 11.52 + 1.152 + 0.1152 + 0.01152
+// On voit que pour Ai = 9
+//  1 terme,  A1       : I = int(128*U) = int(115.2) = 115
+//  2 termes, A1,2     : I = int(128*U) = int(115.2+11.52) = int(126.72) = 126
+//  3 termes, A1,2,3   : I = int(128*U) = int(115.2+11.52+1.152) = int(127.872) = 127
+//  4 termes, A1,2,3,4 : I = int(128*U) = int(115.2+11.52+1.152+0.1152) = int(127.9872) = 127
+//  ==> le digit A4 ne sert pas dans la determination de "I"
+//  On prend un digit de + pour avoir de la marge -> A5 cad le dernier digit de int(U*10^5)
+//--------
+{
+ while(1) {
+   double U;
+   // -- choix de l'intervalle et de l'abscisse "x"
+   int I = N_ZIGGURRAT_GAUSS128;
+   while(I>=N_ZIGGURRAT_GAUSS128) {
+     U = Next();
+     I = int(N_ZIGGURRAT_GAUSS128*U);
+   }
+
+   // -- choix du signe (cf blabla ci-dessus)
+   double s = ( (int(U*100000)&1) == 0 ) ? -1.: 1.;
+
+   // -- choix de l'abscisse "x" dans l'intervalle
+   double x = Next() * X_ZIGGURRAT_GAUSS128[I+1];
+
+   // -- x est dans l'interieur de la bande
+   if(x<X_ZIGGURRAT_GAUSS128[I]) return s * x;
+
+   // -- x n'est pas a l'interieur de la bande mais dans la partie a 2 possibilites
+
+   // l'intervalle est celui qui contient la queue a l'infini
+   // On s'assure que la partie "rejection sur la gaussienne" ne sera pas appelle
+   //   (cad que slim=1. < X[127]) pour eviter les recursions infinies (possibles?)
+   if(I==N_ZIGGURRAT_GAUSS128-1)  // cas de la bande de la queue I=127
+     return s * GaussianTail(X_ZIGGURRAT_GAUSS128[N_ZIGGURRAT_GAUSS128-1],1.);
+
+   // on tire "y" uniforme a l'interieur des ordonnees de la bande choisie
+   // et on regarde si on est en-dessous ou au-dessus de la pdf
+   double y = Y_ZIGGURRAT_GAUSS128[I+1]
+            + Next()*(Y_ZIGGURRAT_GAUSS128[I]-Y_ZIGGURRAT_GAUSS128[I+1]);
+   double pdf = exp(-0.5*x*x);
+   if(pdf>y) return s * x;
+
+   // echec, on est au-dessus de la pdf -> on re-essaye
+ }
+}
+
+r_8 RandomGeneratorInterface::GaussianTail(double s,double slim)
 {
   /* Returns a gaussian random variable larger than a
@@ -397 +526 @@
    */
 
-  if (s < 1)
+  if (s < slim)
     {
       /* For small s, use a direct rejection method. The limit s < 1

trunk/SophyaLib/BaseTools/randinterf.h

@@ -29 +29 @@
   C_Gaussian_PolarBoxMuller = 2,
   C_Gaussian_RatioUnif = 3,
-  C_Gaussian_LevaRatioUnif = 4
+  C_Gaussian_LevaRatioUnif = 4,
+  C_Gaussian_Ziggurat128 = 5
 };
 
@@ -88 +89 @@
   virtual r_8 GaussianRatioUnif(); 
   virtual r_8 GaussianLevaRatioUnif();
+  virtual r_8 GaussianZiggurat128(); 
   /*! \brief Return a random number following a gaussian distribution "sigma", (mean=0)*/ 
   inline r_8 Gaussian(double sigma) {return sigma*Gaussian();}
@@ -94 +96 @@
 
   /*! \brief Return a random number following a gaussian tail distribution for x>sdev */ 
-  virtual r_8 GaussianTail(double sdev);
+  virtual r_8 GaussianTail(double sdev,double slim=1.);
 
   // --- Le tirage sur une distribution de poisson

trunk/SophyaLib/BaseTools/srandgen.h

@@ -30 +30 @@
   {return RandomGeneratorInterface::GetGlobalRandGenP()->Gaussian(sigma,mu);}
 
+  inline double GaussianTailRand(double sdev,double slim=1.) 
+  {return RandomGeneratorInterface::GetGlobalRandGenP()->GaussianTail(sdev,slim);}
+
 inline complex< r_8 > ComplexGaussianRand()
   {return RandomGeneratorInterface::GetGlobalRandGenP()->ComplexGaussian();}