source: Sophya/trunk/SophyaLib/BaseTools/randinterf.cc@ 3992

Last change on this file since 3992 was 3838, checked in by ansari, 15 years ago

MAJ commentaires pour doxygen, passage methode Next() en pure virtual (=0), Reza 09/08/2010
File size: 30.2 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <math.h>
4#include <stdlib.h>
5#include <stdio.h>
6#include <sys/time.h>
7#include <time.h>
8#include <iostream>
9#include <typeinfo>
10
11#include "pexceptions.h"
12
13#include "randinterf.h"
14
15namespace SOPHYA {
16
17//-------------------------------------------------------------------------------
18// ------ Definition d'interface des classes de generateurs de nombres aleatoires
19/*!
20 \class RandomGeneratorInterface
21 \ingroup BaseTools
22 \brief Base class for random number generators
23
24 This class defines the interface for random number generator classes and
25 implements the generation of some specific distributions (Gaussian, Poisson ...)
26 through generation of random number with a flat distribution in the range [0,1[.
27
28 The sub classes inheriting from this class should implement the Next() method.
29
30 This base class manages also a global instance of a default generator.
31
32 \sa frand01 drand01 frandpm1 drandpm1
33 \sa Gaussian Poisson
34
35*/
36
37
38RandomGeneratorInterface* RandomGeneratorInterface::gl_rndgen_p = NULL;
39
40/*!
41 \brief: static method to set or change the intance of the global Random Generator object
42
43 This method should be called during initialization, before any call to global
44 functions for random number generation. The rgp object should be created using new.
45*/
46void RandomGeneratorInterface::SetGlobalRandGenP(RandomGeneratorInterface* rgp)
47{
48 if (rgp == NULL) return;
49 if (gl_rndgen_p) delete gl_rndgen_p;
50 gl_rndgen_p = rgp;
51 return;
52}
53
54RandomGeneratorInterface::RandomGeneratorInterface()
55{
56 SelectGaussianAlgo();
57 SelectPoissonAlgo();
58 SelectExponentialAlgo();
59}
60
61
62RandomGeneratorInterface::~RandomGeneratorInterface(void)
63{
64 // rien a faire
65}
66
67void RandomGeneratorInterface::ShowRandom()
68{
69 cout << " RandomGeneratorInterface::ShowRandom() typeid(this)=" << typeid(*this).name() << " @ "
70 << hex << (unsigned long)(this) << dec << endl;
71 return;
72}
73
74/////////////////////////////////////////////////////////////////////////
75/////////////////////////////////////////////////////////////////////////
76/////////////////////////////////////////////////////////////////////////
77
78/*
79r_8 RandomGeneratorInterface::Next()
80{
81 printf("RandomGeneratorInterface::Next(): undefined code !!!\n");
82 throw MathExc("RandomGeneratorInterface::Next(): undefined code !!!");
83}
84*/
85
86/////////////////////////////////////////////////////////////////////////
87/////////////////////////////////////////////////////////////////////////
88/////////////////////////////////////////////////////////////////////////
89void RandomGeneratorInterface::GenerateSeedVector(int nseed,vector<uint_2>& seed,int lp)
90// renvoie un vecteur de nseed+2 entiers 32 bits
91// [0 - 2] = codage sur 48 bits du nombre (melange) de microsec depuis l'origine
92// [3 -> 3+ngene-1] = entiers aleatoires (poor man generator)
93//
94// L'initialiseur est donne par un codage du nombre de millisecondes
95// ecoulees depuis le 0 heure le 1er Janvier 1970 UTC (cf gettimeofday).
96// Seuls les 48 bits de poids faible sont retenus.
97// Un melange des bits est ensuite effectue pour que les 3 nombres
98// (unsigned short) d'initialisation ne soient pas trop semblables.
99// Le nombre le plus grand que l'on peut mettre
100// dans un entier unsigned de N bits est: 2^N-1
101// 48 bits -> 2^48-1 = 281474976710655 musec = 3257.8j = 8.9y
102// -> meme initialisation tous les 8.9 ans a 1 microsec pres !
103{
104 if(lp>0) cout<<"RandomGeneratorInterface::GenerateSeedVector: nseed="<<nseed<<endl;
105
106 // ---
107 // --- les deux premiers mots remplis avec le temps
108 // ---
109 // On recupere le temps ecoule depuis l'origine code en sec+musec
110 struct timeval now;
111 gettimeofday(&now,0);
112 // Calcul du temps ecoule depuis l'origine en microsecondes
113 uint_8 tmicro70 = (uint_8)now.tv_sec*(uint_8)1000000 + (uint_8)now.tv_usec;
114 if(lp>1) cout<<"."<<now.tv_sec<<" sec + "<<now.tv_usec<<" musec = "<<tmicro70<<" musec"<<endl;
115 // Remplissage du tableau de 48 bits
116 uint_2 b[48]; uint_8 tdum = tmicro70;
117 for(int ip=0;ip<48;ip++) {b[ip] = tdum&1; tdum = (tdum>>1);}
118 if(lp>2) {
119 cout<<"..b= ";
120 for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
121 cout<<endl;
122 }
123 // Melange des bits qui varient vite (poids faible, microsec)
124 // avec ceux variant lentement (poids fort, sec)
125 for(int ip=0;ip<16;ip++) {
126 if(ip%3==1) swap(b[ip],b[32+ip]);
127 else if(ip%3==2) swap(b[ip],b[16-ip]);
128 }
129 if(lp>2) {
130 cout<<"..b= ";
131 for(int ip=47;ip>=0;ip--) {cout<<b[ip]; if(ip%32==0 || ip%16==0) cout<<" ";}
132 cout<<endl;
133 }
134 // Remplissage
135 seed.resize(0);
136 for(int i=0;i<3;i++) {
137 seed.push_back(0);
138 uint_2 w = 1;
139 for(int ip=0;ip<16;ip++) {seed[i] += w*b[i*16+ip]; w *= 2;}
140 }
141 if(lp>0) cout<<"seed(time): "<<seed[0]<<" "<<seed[1]<<" "<<seed[2]<<endl;
142
143 // ---
144 // --- generation des nombres aleatoires complementaires (poor man generator)
145 // ---
146 //----------------------------------------------------------------------------//
147 // Ran088: L'Ecuyer's 1996 three-component Tausworthe generator "taus88"
148 // Returns an integer random number uniformly distributed within [0,4294967295]
149 // The period length is approximately 2^88 (which is 3*10^26).
150 // This generator is very fast and passes all standard statistical tests.
151 // Reference:
152 // (1) P. L'Ecuyer, Maximally equidistributed combined Tausworthe generators,
153 // Mathematics of Computation, 65, 203-213 (1996), see Figure 4.
154 // (2) recommended in:
155 // P. L'Ecuyer, Random number generation, chapter 4 of the
156 // Handbook on Simulation, Ed. Jerry Banks, Wiley, 1997.
157 //----------------------------------------------------------------------------//
158 if(nseed<=0) return;
159 // initialize seeds using the given seed value taking care of
160 // the requirements. The constants below are arbitrary otherwise
161 uint_4 seed0 = uint_4(tmicro70&0xFFFFFFFFULL);
162 if(lp>2) cout<<"seed0(time): "<<seed0<<endl;
163 uint_4 state_s1, state_s2, state_s3;
164 state_s1 = 1243598713U ^ seed0; if (state_s1 < 2) state_s1 = 1243598713U;
165 state_s2 = 3093459404U ^ seed0; if (state_s2 < 8) state_s2 = 3093459404U;
166 state_s3 = 1821928721U ^ seed0; if (state_s3 < 16) state_s3 = 1821928721U;
167 int nfill = 0, ico=0;
168 while(nfill<nseed) {
169 uint_4 s1 = state_s1, s2 = state_s2, s3 = state_s3;
170 // generate a random 32 bit number
171 s1 = ((s1 & -2) << 12) ^ (((s1 << 13) ^ s1) >> 19);
172 s2 = ((s2 & -8) << 4) ^ (((s2 << 2) ^ s2) >> 25);
173 s3 = ((s3 & -16) << 17) ^ (((s3 << 3) ^ s3) >> 11);
174 state_s1 = s1; state_s2 = s2; state_s3 = s3;
175 // le nombre aleatoire sur 32 bits est: s1^s2^s3
176 if(ico<15) {ico++; continue;} // des tirages blancs
177 uint_2 s = uint_2( (s1^s2^s3)&0xFFFFU );
178 seed.push_back(s);
179 if(lp>0) cout<<"seed(t88): "<<seed[3+nfill]<<endl;
180 nfill++;
181 }
182
183}
184
185void RandomGeneratorInterface::AutoInit(int lp)
186{
187 printf("RandomGeneratorInterface::AutoInit(): undefined code !!!\n");
188 throw MathExc("RandomGeneratorInterface::AutoInit(): undefined code !!!");
189}
190
191/////////////////////////////////////////////////////////////////////////
192/////////////////////////////////////////////////////////////////////////
193/////////////////////////////////////////////////////////////////////////
194
195r_8 RandomGeneratorInterface::Gaussian()
196{
197 switch (usegaussian_) {
198 case C_Gaussian_BoxMuller :
199 return GaussianBoxMuller();
200 break;
201 case C_Gaussian_RandLibSNorm :
202 return GaussianSNorm();
203 break;
204 case C_Gaussian_PolarBoxMuller :
205 return GaussianPolarBoxMuller();
206 break;
207 case C_Gaussian_RatioUnif :
208 return GaussianRatioUnif();
209 break;
210 case C_Gaussian_LevaRatioUnif :
211 return GaussianLevaRatioUnif();
212 break;
213 default:
214 return GaussianBoxMuller();
215 break;
216 }
217}
218
219//--- Generation de nombre aleatoires suivant une distribution gaussienne
220r_8 RandomGeneratorInterface::GaussianBoxMuller()
221{
222 r_8 A=Next();
223 while (A==0.) A=Next();
224 return sqrt(-2.*log(A))*cos(2.*M_PI*Next());
225}
226
227//-------------------------------------------
228// Adapte de ranlib float snorm()
229// http://orion.math.iastate.edu/burkardt/c_src/ranlib/ranlib.c
230/*
231**********************************************************************
232 (STANDARD-) N O R M A L DISTRIBUTION
233**********************************************************************
234
235 FOR DETAILS SEE:
236
237 AHRENS, J.H. AND DIETER, U.
238 EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
239 SAMPLING FROM THE NORMAL DISTRIBUTION.
240 MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
241
242 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
243 (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
244
245 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
246 SUNIF. The argument IR thus goes away.
247
248**********************************************************************
249 THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
250 H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
251*/
252static double a_snorm[32] = {
253 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
254 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
255 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
256 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
257 1.862732,2.153875
258};
259static double d_snorm[31] = {
260 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
261 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
262 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
263 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
264};
265static float t_snorm[31] = {
266 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
267 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
268 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
269 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
270 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
271};
272static float h_snorm[31] = {
273 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
274 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
275 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
276 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
277 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
278};
279r_8 RandomGeneratorInterface::GaussianSNorm()
280{
281long i;
282double snorm,u,s,ustar,aa,w,y,tt;
283 u = Next();
284 s = 0.0;
285 if(u > 0.5) s = 1.0;
286 u += (u-s);
287 u = 32.0*u;
288 i = (long) (u);
289 if(i == 32) i = 31;
290 if(i == 0) goto S100;
291/*
292 START CENTER
293*/
294 ustar = u-(double)i;
295 aa = *(a_snorm+i-1);
296S40:
297 if(ustar <= *(t_snorm+i-1)) goto S60;
298 w = (ustar-*(t_snorm+i-1))**(h_snorm+i-1);
299S50:
300/*
301 EXIT (BOTH CASES)
302*/
303 y = aa+w;
304 snorm = y;
305 if(s == 1.0) snorm = -y;
306 return snorm;
307S60:
308/*
309 CENTER CONTINUED
310*/
311 u = Next();
312 w = u*(*(a_snorm+i)-aa);
313 tt = (0.5*w+aa)*w;
314 goto S80;
315S70:
316 tt = u;
317 ustar = Next();
318S80:
319 if(ustar > tt) goto S50;
320 u = Next();
321 if(ustar >= u) goto S70;
322 ustar = Next();
323 goto S40;
324S100:
325/*
326 START TAIL
327*/
328 i = 6;
329 aa = *(a_snorm+31);
330 goto S120;
331S110:
332 aa += *(d_snorm+i-1);
333 i += 1;
334S120:
335 u += u;
336 if(u < 1.0) goto S110;
337 u -= 1.0;
338S140:
339 w = u**(d_snorm+i-1);
340 tt = (0.5*w+aa)*w;
341 goto S160;
342S150:
343 tt = u;
344S160:
345 ustar = Next();
346 if(ustar > tt) goto S50;
347 u = Next();
348 if(ustar >= u) goto S150;
349 u = Next();
350 goto S140;
351}
352
353r_8 RandomGeneratorInterface::GaussianPolarBoxMuller()
354{
355double x1,x2,w;
356do {
357 x1 = 2.0 * Next() - 1.0;
358 x2 = 2.0 * Next() - 1.0;
359 w = x1 * x1 + x2 * x2;
360 } while ( w >= 1.0 || w==0. );
361return x1 * sqrt(-2.0*log(w)/w);
362}
363
364static double s2se_RatioUnif=sqrt(2./M_E) , epm135_RatioUnif=exp(-1.35) , ep1q_RatioUnif=exp(1./4.);
365r_8 RandomGeneratorInterface::GaussianRatioUnif()
366{
367double u,v,x;
368while(true) {
369 do {u = Next();} while ( u == 0. );
370 v = (2.0*Next()-1.0)*s2se_RatioUnif;
371 x = v/u;
372 if(x*x <= 5.0-4.0*ep1q_RatioUnif*u) break;
373 if(x*x<4.0*epm135_RatioUnif/u+1.4)
374 if(v*v<-4.0*u*u*log(u)) break;
375}
376return x;
377}
378
379r_8 RandomGeneratorInterface::GaussianLevaRatioUnif()
380{
381double u,v,x,y,q;
382 do {
383 u = 1.-Next(); // in ]0,1]
384 v = Next()-0.5; // in [-0.5, 0.5[
385 v *= 1.7156;
386 x = u - 0.449871;
387 y = ((v<0)?-v:v) + 0.386595;
388 q = x*x + y*(0.19600*y - 0.25472*x);
389 } while( q>=0.27597 && (q>0.27846 || v*v>-4.0*u*u*log(u)) );
390 return v/u;
391}
392
393r_8 RandomGeneratorInterface::GaussianTail(double s)
394{
395 /* Returns a gaussian random variable larger than a
396 * This implementation does one-sided upper-tailed deviates.
397 */
398
399 if (s < 1)
400 {
401 /* For small s, use a direct rejection method. The limit s < 1
402 can be adjusted to optimise the overall efficiency */
403 double x;
404 do
405 {
406 x = Gaussian();
407 }
408 while (x < s);
409 return x;
410 }
411 else
412 {
413 /* Use the "supertail" deviates from the last two steps
414 * of Marsaglia's rectangle-wedge-tail method, as described
415 * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11, p139,
416 * and the solution, p586.)
417 */
418 double u, v, x;
419 do
420 {
421 u = Next();
422 do
423 {
424 v = Next();
425 }
426 while (v == 0.0);
427 x = sqrt (s * s - 2 * log (v));
428 }
429 while (x * u > s);
430 return x;
431 }
432}
433
434/////////////////////////////////////////////////////////////////////////
435/////////////////////////////////////////////////////////////////////////
436/////////////////////////////////////////////////////////////////////////
437
438uint_8 RandomGeneratorInterface::Poisson(double mu, double mumax)
439{
440 switch (usepoisson_) {
441 case C_Poisson_Simple :
442 return PoissonSimple(mu,mumax);
443 break;
444 case C_Poisson_Ahrens :
445 return PoissonAhrens(mu);
446 break;
447 default:
448 return PoissonSimple(mu,mumax);
449 break;
450 }
451}
452
453
454//--- Generation de nombre aleatoires suivant une distribution de Poisson
455uint_8 RandomGeneratorInterface::PoissonSimple(double mu,double mumax)
456{
457 double pp,ppi,x;
458
459 if((mumax>0.)&&(mu>=mumax)) {
460 pp = sqrt(mu);
461 while( (x=pp*Gaussian()) < -mu );
462 return (uint_8)(mu+x+0.5);
463 }
464 else {
465 uint_8 n;
466 ppi = pp = exp(-mu);
467 x = Next();
468 n = 0;
469 while (x > ppi) {
470 n++;
471 pp = mu*pp/(double)n;
472 ppi += pp;
473 }
474 return n;
475 }
476 return 0; // pas necessaire ?
477}
478
479
480static double a0_poiahr = -0.5;
481static double a1_poiahr = 0.3333333;
482static double a2_poiahr = -0.2500068;
483static double a3_poiahr = 0.2000118;
484static double a4_poiahr = -0.1661269;
485static double a5_poiahr = 0.1421878;
486static double a6_poiahr = -0.1384794;
487static double a7_poiahr = 0.125006;
488static double fact_poiahr[10] = {
489 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0};
490uint_8 RandomGeneratorInterface::PoissonAhrens(double mu)
491/*
492**********************************************************************
493 long ignpoi(float mu)
494 GENerate POIsson random deviate
495 Function
496 Generates a single random deviate from a Poisson
497 distribution with mean AV.
498 Arguments
499 av --> The mean of the Poisson distribution from which
500 a random deviate is to be generated.
501 genexp <-- The random deviate.
502 Method
503 Renames KPOIS from TOMS as slightly modified by BWB to use RANF
504 instead of SUNIF.
505 For details see:
506 Ahrens, J.H. and Dieter, U.
507 Computer Generation of Poisson Deviates
508 From Modified Normal Distributions.
509 ACM Trans. Math. Software, 8, 2
510 (June 1982),163-179
511**********************************************************************
512**********************************************************************
513
514
515 P O I S S O N DISTRIBUTION
516
517
518**********************************************************************
519**********************************************************************
520
521 FOR DETAILS SEE:
522
523 AHRENS, J.H. AND DIETER, U.
524 COMPUTER GENERATION OF POISSON DEVIATES
525 FROM MODIFIED NORMAL DISTRIBUTIONS.
526 ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
527
528 (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
529
530**********************************************************************
531 INTEGER FUNCTION IGNPOI(IR,MU)
532 INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
533 MU=MEAN MU OF THE POISSON DISTRIBUTION
534 OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
535 MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
536 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
537 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
538 SEPARATION OF CASES A AND B
539*/
540{
541uint_8 ignpoi,j,k,kflag,l,m;
542double b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
543 t,u,v,x,xx,pp[35];
544
545 if(mu < 10.0) goto S120;
546/*
547 C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
548*/
549 s = sqrt(mu);
550 d = 6.0*mu*mu;
551/*
552 THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
553 PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
554 IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
555*/
556 l = (uint_8) (mu-1.1484);
557/*
558 STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
559*/
560 g = mu+s*Gaussian();
561 if(g < 0.0) goto S20;
562 ignpoi = (uint_8) (g);
563/*
564 STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
565*/
566 if(ignpoi >= l) return ignpoi;
567/*
568 STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
569*/
570 fk = (double)ignpoi;
571 difmuk = mu-fk;
572 u = Next();
573 if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
574S20:
575/*
576 STEP P. PREPARATIONS FOR STEPS Q AND H.
577 (RECALCULATIONS OF PARAMETERS IF NECESSARY)
578 .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
579 THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
580 APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
581 C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
582*/
583 omega = 0.3989423/s;
584 b1 = 4.166667E-2/mu;
585 b2 = 0.3*b1*b1;
586 c3 = 0.1428571*b1*b2;
587 c2 = b2-15.0*c3;
588 c1 = b1-6.0*b2+45.0*c3;
589 c0 = 1.0-b1+3.0*b2-15.0*c3;
590 c = 0.1069/mu;
591 if(g < 0.0) goto S50;
592/*
593 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
594*/
595 kflag = 0;
596 goto S70;
597S40:
598/*
599 STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
600*/
601 if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
602S50:
603/*
604 STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
605 DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
606 (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
607*/
608 e = Exponential();
609 u = Next();
610 u += (u-1.0);
611 //t = 1.8+fsign(e,u);
612 t = 1.8 + (((u>0. && e<0.) || (u<0. && e>0.))?-e:e);
613 if(t <= -0.6744) goto S50;
614 ignpoi = (uint_8) (mu+s*t);
615 fk = (double)ignpoi;
616 difmuk = mu-fk;
617/*
618 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
619*/
620 kflag = 1;
621 goto S70;
622S60:
623/*
624 STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
625*/
626 if(c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
627 return ignpoi;
628S70:
629/*
630 STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
631 CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
632*/
633 if(ignpoi >= 10) goto S80;
634 px = -mu;
635 py = pow(mu,(double)ignpoi)/ *(fact_poiahr+ignpoi);
636 goto S110;
637S80:
638/*
639 CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
640 A0-A7 FOR ACCURACY WHEN ADVISABLE
641 .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
642*/
643 del = 8.333333E-2/fk;
644 del -= (4.8*del*del*del);
645 v = difmuk/fk;
646 if(fabs(v) <= 0.25) goto S90;
647 px = fk*log(1.0+v)-difmuk-del;
648 goto S100;
649S90:
650 px = fk*v*v*(((((((a7_poiahr*v+a6_poiahr)*v+a5_poiahr)*v+a4_poiahr)*v+a3_poiahr)*v+a2_poiahr)*v+a1_poiahr)*v+a0_poiahr)-del;
651S100:
652 py = 0.3989423/sqrt(fk);
653S110:
654 x = (0.5-difmuk)/s;
655 xx = x*x;
656 fx = -0.5*xx;
657 fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
658 if(kflag <= 0) goto S40;
659 goto S60;
660S120:
661/*
662 C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
663*/
664// m = max(1L,(long) (mu));
665 m = (1ULL >= (uint_8)mu) ? 1ULL: (uint_8)mu;
666
667 l = 0;
668 p = exp(-mu);
669 q = p0 = p;
670S130:
671/*
672 STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
673*/
674 u = Next();
675 ignpoi = 0;
676 if(u <= p0) return ignpoi;
677/*
678 STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
679 PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
680 (0.458=PP(9) FOR MU=10)
681*/
682 if(l == 0) goto S150;
683 j = 1;
684//if(u > 0.458) j = min(l,m);
685 if(u > 0.458) j = ((l<=m)? l: m);
686 for(k=j; k<=l; k++) {
687 if(u <= *(pp+k-1)) goto S180;
688 }
689 if(l == 35) goto S130;
690S150:
691/*
692 STEP C. CREATION OF NEW POISSON PROBABILITIES P
693 AND THEIR CUMULATIVES Q=PP(K)
694*/
695 l += 1;
696 for(k=l; k<=35; k++) {
697 p = p*mu/(double)k;
698 q += p;
699 *(pp+k-1) = q;
700 if(u <= q) goto S170;
701 }
702 l = 35;
703 goto S130;
704S170:
705 l = k;
706S180:
707 ignpoi = k;
708 return ignpoi;
709}
710
711/////////////////////////////////////////////////////////////////////////
712/////////////////////////////////////////////////////////////////////////
713/////////////////////////////////////////////////////////////////////////
714
715r_8 RandomGeneratorInterface::Exponential()
716{
717 switch (useexpo_) {
718 case C_Exponential_Simple :
719 return ExpoSimple();
720 break;
721 case C_Exponential_Ahrens :
722 return ExpoAhrens();
723 break;
724 default:
725 return ExpoSimple();
726 break;
727 }
728}
729
730r_8 RandomGeneratorInterface::ExpoSimple(void)
731{
732 return -log(1.-Next());
733}
734
735
736static double q_expo[8] = {
737 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,1.0};
738r_8 RandomGeneratorInterface::ExpoAhrens(void)
739/*
740**********************************************************************
741**********************************************************************
742 (STANDARD-) E X P O N E N T I A L DISTRIBUTION
743**********************************************************************
744**********************************************************************
745
746 FOR DETAILS SEE:
747
748 AHRENS, J.H. AND DIETER, U.
749 COMPUTER METHODS FOR SAMPLING FROM THE
750 EXPONENTIAL AND NORMAL DISTRIBUTIONS.
751 COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
752
753 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
754 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
755
756 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
757 SUNIF. The argument IR thus goes away.
758
759**********************************************************************
760 Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
761 (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
762*/
763{
764long i;
765double sexpo,a,u,ustar,umin;
766double *q1 = q_expo;
767 a = 0.0;
768 while((u=Next())==0.);
769 goto S30;
770S20:
771 a += *q1;
772S30:
773 u += u;
774 if(u <= 1.0) goto S20;
775 u -= 1.0;
776 if(u > *q1) goto S60;
777 sexpo = a+u;
778 return sexpo;
779S60:
780 i = 1;
781 ustar = Next();
782 umin = ustar;
783S70:
784 ustar = Next();
785 if(ustar < umin) umin = ustar;
786 i += 1;
787 if(u > *(q_expo+i-1)) goto S70;
788 sexpo = a+umin**q1;
789 return sexpo;
790}
791
792/////////////////////////////////////////////////////////////////////////
793/////////////////////////////////////////////////////////////////////////
794/////////////////////////////////////////////////////////////////////////
795
796int RandomGeneratorInterface::Gaussian2DRho(double &x,double &y,double mx,double my,double sx,double sy,double ro)
797/*
798++
799| Tirage de 2 nombres aleatoires x et y distribues sur une gaussienne 2D
800| de centre (mx,my), de coefficient de correlation rho (ro) et telle que
801| les sigmas finals des variables x et y soient sx,sy (ce sont
802| les valeurs des distributions marginales des variables aleatoires x et y
803| c'est a dire les sigmas des projections x et y de l'histogramme 2D
804| de la gaussienne). Retourne 0 si ok.
805|
806| - La densite de probabilite (normalisee a 1) sur laquelle on tire est:
807| N*exp[-0.5*{[(dx/sx)^2-2*ro/(sx*sy)*dx*dy+(dy/sy)^2]/(1-ro^2)}]
808| avec dx = x-mx, dy = y-my et N = 1/[2Pi*sx*sy*sqrt(1-ro^2)]
809| - Dans ce cas la distribution marginale est (ex en X):
810| 1/(sqrt(2Pi)*sx) * exp[-0.5*{dx^2/sx^2}]
811| - La matrice des covariances C des variables x,y est:
812| | sx^2 ro*sx*sy |
813| | | et det(C) = (1-ro^2)*sx^2*sy^2
814| | ro*sx*sy sy^2 |
815| - La matrice inverse C^(-1) est:
816| | 1/sx^2 -ro/(sx*sy) |
817| | | * 1/(1-ro^2)
818| | -ro/(sx*sy) 1/sy^2 |
819|
820| - Remarque:
821| le sigma que l'on obtient quand on fait une coupe de la gaussienne 2D
822| en y=0 (ou x=0) est: SX0(y=0) = sx*sqrt(1-ro^2) different de sx
823| SY0(x=0) = sy*sqrt(1-ro^2) different de sy
824| La distribution qui correspond a des sigmas SX0,SY0
825| pour les coupes en y=0,x=0 de la gaussienne 2D serait:
826| N*exp[-0.5*{ (dx/SX0)^2-2*ro/(SX0*SY0)*dx*dy+(dy/SY0)^2 }]
827| avec N = sqrt(1-ro^2)/(2Pi*SX0*SY0) et les variances
828| des variables x,y sont toujours
829| sx=SX0/sqrt(1-ro^2), sy=SY0/sqrt(1-ro^2)
830--
831*/
832{
833double a,b,sa;
834
835if( ro <= -1. || ro >= 1. ) return 1;
836
837while( (b=Flat01()) == 0. );
838b = sqrt(-2.*log(b));
839a = 2.*M_PI * Flat01();
840sa = sin(a);
841
842x = mx + sx*b*(sqrt(1.-ro*ro)*cos(a)+ro*sa);
843y = my + sy*b*sa;
844
845return 0;
846}
847
848void RandomGeneratorInterface::Gaussian2DAng(double &x,double &y,double mx,double my,double sa,double sb,double teta)
849/*
850++
851| Tirage de 2 nombres aleatoires x et y distribues sur une gaussienne 2D
852| de centre (x=mx,y=my), de sigmas grand axe et petit axe (sa,sb)
853| et dont le grand axe fait un angle teta (radian) avec l'axe des x.
854|
855| - La densite de probabilite (normalisee a 1) sur laquelle on tire est:
856| N*exp[-0.5*{ (A/sa)**2+(C/sc)**2 }], N=1/(2Pi*sa*sc)
857| ou A et B sont les coordonnees selon le grand axe et le petit axe
858| et teta = angle(x,A), le resultat subit ensuite une rotation d'angle teta.
859| - La matrice des covariances C des variables A,B est:
860| | sa^2 0 |
861| | | et det(C) = (1-ro^2)*sa^2*sb^2
862| | 0 sb^2 |
863| - La distribution x,y resultante est:
864| N*exp[-0.5*{[(dx/sx)^2-2*ro/(sx*sy)*dx*dy+(dy/sy)^2]/(1-ro^2)}]
865| ou N est donne dans NormCo et sx,sy,ro sont calcules a partir
866| de sa,sc,teta (voir fonctions paramga ou gaparam). La matrice des
867| covariances des variables x,y est donnee dans la fonction NormCo.
868--
869*/
870{
871double c,s,X,Y;
872
873while( (s = Flat01()) == 0. );
874s = sqrt(-2.*log(s));
875c = 2.*M_PI * Flat01();
876
877X = sa*s*cos(c);
878Y = sb*s*sin(c);
879
880c = cos(teta); s = sin(teta);
881x = mx + c*X - s*Y;
882y = my + s*X + c*Y;
883}
884
885} /* namespace SOPHYA */
886
887
888
889/////////////////////////////////////////////////////////////////
890/*
891**** Remarques sur complex< r_8 > ComplexGaussian(double sig) ****
892
893--- variables gaussiennes x,y independantes
894x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
895y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2))
896x,y independants --> pdf f(x,y) = f(x) f(y)
897On a:
898 <x> = Integrate[x*f(x)] = Mx
899 <x^2> = Integrate[x^2*f(x)] = Mx^2 + Sx^2
900
901--- On cherche la pdf g(r,t) du module et de la phase
902 x = r cos(t) , y = r sin(t)
903 r=sqrt(x^2+y^2 , t=atan2(y,x)
904 (r,t) --> (x,y): le Jacobien = r
905
906 g(r,t) = r f(x,y) = r f(x) f(y)
907 = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2))
908
909- Le cas general est complique
910 (cf D.Pelat cours DEA "bruits et signaux" section 4.5)
911
912- Cas ou "Mx = My = 0" et "Sx = Sy = S"
913 c'est la pdf du module et de la phase d'un nombre complexe
914 dont les parties reelles et imaginaires sont independantes
915 et sont distribuees selon des gaussiennes de variance S^2
916 g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2))
917 La distribution de "r" est donc:
918 g(r) = Integrate[g(r,t),{t,0,2Pi}]
919 = r/S^2 exp(-r^2/(2 S^2))
920 La distribution de "t" est donc:
921 g(t) = Integrate[g(r,t),{r,0,Infinity}]
922 = 1 / 2Pi (distribution uniforme sur [0,2Pi[)
923 Les variables aleatoires r,t sont independantes:
924 g(r,t) = g(r) g(t)
925On a:
926 <r> = Integrate[r*g(r)] = sqrt(PI/2)*S
927 <r^2> = Integrate[r^2*g(r)] = 2*S^2
928 <r^3> = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3
929 <r^4> = Integrate[r^4*g(r)] = 8*S^4
930
931- Attention:
932La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie:
933 <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 S^2
934Si on veut generer une variable complexe gaussienne telle que
935 <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument
936
937*/
Note: See TracBrowser for help on using the repository browser.