root / src / sim / cnumgen.cc @ f857128b
History  View  Annotate  Download (9.17 KB)
1 
//==========================================================================


2 
// CNUMGEN.CC

3 
//

4 
// OMNeT++/OMNEST

5 
// Discrete System Simulation in C++

6 
//

7 
// Random variate generation

8 
//

9 
// @author Werner Sandmann (ws), Kay Michael Masslow (kmm)

10 
//

11 
// @date 11/26/2002 "doxygenification" (kmm)

12 
// @date 11/20/2002 some final comments (ws)

13 
// @date 10/22/2002 implemented various discrete distributions (kmm)

14 
//

15 
//==========================================================================

16  
17 
#include <float.h> 
18 
#include <math.h> 
19 
#include "distrib.h" 
20 
#include "regmacros.h" 
21 
#include "cexception.h" 
22 
#include "cnumgen.h" 
23 
#include "cmersennetwister.h" 
24 
#include <pthread.h> 
25  
26 
//

27 
// Linux  <math.h> supplies a PI, MS does not...

28 
//

29 
#ifndef M_PI

30 
#define M_PI 3.1415926535897932384626433832795 
31 
#endif

32  
33 
#ifndef M_E

34 
#define M_E 2.7182818284590452353602874713527 
35 
#endif

36  
37 
cNumberGenerator::cNumberGenerator() : nrRNGs(0) {

38 
if (simulation.isRunning()) {

39 
throw cRuntimeError(

40 
"Cannot (yet) create new instance of cNumberGenerator"

41 
" at runtime due to potentially concurrent access to"

42 
" the central seed generator.");

43 
} 
44 
} 
45  
46  
47 
cNumberGenerator::~cNumberGenerator() { 
48 
for (int i = 0; i < nrRNGs; i++) 
49 
delete rngs[i];

50  
51 
} 
52  
53 
void cNumberGenerator::checkBounds(int rngId) { 
54 
if (nrRNGs <= rngId)

55 
throw cRuntimeError(

56 
"cNumberGenerator only has %d < %d local generators.", nrRNGs,

57 
rngId); 
58 
} 
59  
60 
void cNumberGenerator::setupSeeds(unsigned nrRng, unsigned long * seeds) { 
61 
nrRNGs = nrRng; 
62 
rngs = (cRNG**) new cMersenneTwister*[nrRNGs];

63 
for (int i = 0; i < nrRNGs; i++) { 
64 
rngs[i] = new cMersenneTwister();

65 
rngs[i]>seed(seeds[i]); 
66 
} 
67 
} 
68 
//

69 
//

70 
// C O N T I N U O U S

71 
//

72 
//

73  
74 
double cNumberGenerator::uniform(double a, double b, int rngId) { 
75 
checkBounds(rngId); 
76 
return a + rngs[rngId]>doubleRand() * (b  a);

77  
78 
} 
79  
80 
double cNumberGenerator::exponential(double p, int rngId) { 
81 
checkBounds(rngId); 
82 
return p * log(rngs[rngId]>doubleRand());

83  
84 
} 
85  
86 
double cNumberGenerator::unit_normal(int rngId) { 
87 
checkBounds(rngId); 
88 
return sqrt(2.0 * log(rngs[rngId]>doubleRand())) * cos(PI * 2 
89 
* rngs[rngId]>doubleRand()); 
90  
91 
} 
92  
93 
double cNumberGenerator::normal(double m, double d, int rngId) { 
94 
checkBounds(rngId); 
95 
return m + d * sqrt(2.0 * log(rngs[rngId]>doubleRand())) * cos(PI * 2 
96 
* rngs[rngId]>doubleRand()); 
97  
98 
} 
99  
100 
double cNumberGenerator::truncnormal(double m, double d, int rngId) { 
101 
checkBounds(rngId); 
102 
double res;

103 
do {

104 
res = normal(m, d, rngId); 
105 
} while (res < 0); 
106 
return res;

107 
} 
108  
109 
/*

110 
* internal, for alpha>1.

111 
*

112 
* From: "A Simple Method for Generating Gamma Variables", George Marsaglia and

113 
* Wai Wan Tsang, ACM Transactions on Mathematical Software, Vol. 26, No. 3,

114 
* September 2000. Available online.

115 
*/

116 
double cNumberGenerator::gamma_Marsaglia2000(double a, int rngId) 
117 
{ 
118 
checkBounds(rngId); 
119 
ASSERT(a>1);

120  
121 
double d,c,x,v,u;

122 
d = a  1.0/3.0; 
123 
c = 1.0/sqrt(9.0*d); 
124 
for(;;)

125 
{ 
126 
do {x = unit_normal(rngId); v = 1.0 + c*x;} while (v<=0); 
127 
v = v*v*v; u = rngs[rngId]>doubleRand(); 
128 
if (u < 1.0  0.0331*(x*x)*(x*x)) 
129 
return d*v;

130 
if (log(u)<0.5*x*x + d*(1.0v+log(v))) 
131 
return d*v;

132 
} 
133 
} 
134  
135  
136 
/*

137 
* internal, for alpha<1.

138 
*

139 
* We can derive the alpha<1 case from the alpha>1 case. See "Note" at the

140 
* end of Section 6 of the Marsaglia2000 paper (see above).

141 
*/

142 
double cNumberGenerator::gamma_MarsagliaTransf(double alpha, int rngId) 
143 
{ 
144 
checkBounds(rngId); 
145 
ASSERT(alpha<1);

146  
147 
return gamma_Marsaglia2000(1+alpha,rngId) * pow(rngs[rngId]>doubleRand(), 1/alpha); 
148 
} 
149  
150 
double cNumberGenerator::gamma_d(double alpha, double theta, int rngId) 
151 
{ 
152 
checkBounds(rngId); 
153 
if (alpha<=0  theta<=0) 
154 
throw cRuntimeError("gamma(): alpha and theta params must be positive " 
155 
"(alpha=%lg, theta=%lg)", alpha, theta);

156  
157 
if (fabs(alpha  1.0) <= DBL_EPSILON) 
158 
{ 
159 
return exponential(theta, rngId);

160 
} 
161 
else if (alpha < 1.0) 
162 
{ 
163 
return theta * gamma_MarsagliaTransf(alpha, rngId);

164 
} 
165 
else // if (alpha > 1.0) 
166 
{ 
167 
return theta * gamma_Marsaglia2000(alpha, rngId);

168 
} 
169 
} 
170  
171  
172 
double cNumberGenerator::beta(double alpha1, double alpha2, int rngId) 
173 
{ 
174 
checkBounds(rngId); 
175 
if (alpha1<=0  alpha2<=0) 
176 
throw cRuntimeError("beta(): alpha1 and alpha2 parameters must be positive " 
177 
"(alpha1=%lg, alpha2=%lg)", alpha1, alpha2);

178  
179 
double Y1 = gamma_d(alpha1, 1.0, rngId); 
180 
double Y2 = gamma_d(alpha2, 1.0, rngId); 
181  
182 
return Y1 / (Y1 + Y2);

183 
} 
184  
185  
186 
double cNumberGenerator::erlang_k(unsigned int k, double m, int rngId) 
187 
{ 
188 
checkBounds(rngId); 
189 
double U = 1.0; 
190 
for (unsigned int i = 0; i < k; i++) 
191 
U *= (1.0  rngs[rngId]>doubleRand()); 
192  
193 
return (m / (double) k) * log(U); 
194 
} 
195  
196  
197 
double cNumberGenerator::chi_square(unsigned int k, int rngId) 
198 
{ 
199 
checkBounds(rngId); 
200 
if (!(k % 2)) 
201 
return erlang_k(k >> 1, k, rngId); 
202 
else

203 
return gamma_d((double) k / 2.0, 2.0, rngId); 
204 
} 
205  
206  
207 
double cNumberGenerator::student_t(unsigned int i, int rngId) 
208 
{ 
209 
checkBounds(rngId); 
210 
double Z = normal(0, 1, rngId); 
211 
double W = sqrt(chi_square(i, rngId) / (double) i); 
212 
return Z / W;

213 
} 
214  
215  
216 
double cNumberGenerator::cauchy(double a, double b, int rngId) 
217 
{ 
218 
checkBounds(rngId); 
219 
if (b<=0) 
220 
throw cRuntimeError("cauchy(): parameters must be b>0 (a=%lg, b=%lg)", a, b); 
221  
222 
return a + b * tan(M_PI * rngs[rngId]>doubleRand());

223  
224 
} 
225  
226  
227 
double cNumberGenerator::triang(double a, double b, double c, int rngId) 
228 
{ 
229 
checkBounds(rngId); 
230 
if (b<a  c<b  a==c)

231 
throw cRuntimeError("triang(): parameters must be a<=b<=c, a<c (a=%lg, b=%lg, c=%lg)", a, b, c); 
232  
233 
double U, beta, T;

234  
235 
U = rngs[rngId]>doubleRand(); 
236 
beta = (b  a) / (c  a); 
237  
238 
if (U < beta)

239 
T = sqrt(beta * U); 
240 
else

241 
T = 1.0  sqrt((1.0  beta) * (1.0  U)); 
242  
243 
return a + (c  a) * T;

244 
} 
245  
246  
247 
// lognormal() is inline

248  
249  
250 
double cNumberGenerator::weibull(double a, double b, int rngId) 
251 
{ 
252 
checkBounds(rngId); 
253 
if (a<=0  b<=0) 
254 
throw cRuntimeError("weibull(): a,b parameters must be positive (a=%lg, b=%lg)", a, b); 
255  
256 
return a * pow(log(1.0  rngs[rngId]>doubleRand()), 1.0 / b); 
257 
} 
258  
259  
260 
double cNumberGenerator::pareto_shifted(double a, double b, double c, int rngId) 
261 
{ 
262 
checkBounds(rngId); 
263 
if (a==0) 
264 
throw cRuntimeError("pareto_shifted(): parameter a cannot be zero)"); 
265  
266 
double u_pow = pow(1.0  rngs[rngId]>doubleRand(), 1.0 / a); 
267  
268 
return (b  c * u_pow) / u_pow;

269 
} 
270  
271 
//

272 
//

273 
// D I S C R E T E

274 
//

275 
//

276  
277  
278 
int cNumberGenerator::intuniform(int a, int b, int rngId) 
279 
{ 
280 
checkBounds(rngId); 
281 
return a + rngs[rngId]>intRand(ba+1); 
282 
} 
283  
284  
285 
// bernoulli() is inline

286  
287  
288 
int cNumberGenerator::binomial(int n, double p, int rngId) 
289 
{ 
290 
checkBounds(rngId); 
291 
int X = 0; 
292 
// sum up n bernoulli trials

293 
for (int i = 0; i < n; i++) 
294 
{ 
295 
double U = rngs[rngId]>doubleRand();

296 
if (p > U)

297 
X++; 
298 
} 
299 
return X;

300 
} 
301  
302  
303 
int cNumberGenerator::geometric(double p, int rngId) 
304 
{ 
305 
checkBounds(rngId); 
306 
double a = 1.0 / (log(1.0  p)); 
307 
return (int)floor(a * log(rngs[rngId]>doubleRand())); 
308 
} 
309  
310  
311 
int cNumberGenerator::negbinomial(int n, double p, int rngId) 
312 
{ 
313 
checkBounds(rngId); 
314 
int X = 0; 
315 
for (int i = 0; i < n; i++) 
316 
{ 
317 
X += geometric(p, rngId); 
318 
} 
319 
return X;

320 
} 
321  
322  
323 
int cNumberGenerator::poisson(double lambda, int rngId) 
324 
{ 
325 
checkBounds(rngId); 
326 
int X;

327 
if (lambda > 30.0) 
328 
{ 
329 
double a = M_PI * sqrt(lambda / 3.0); 
330 
double b = a / lambda;

331 
double c = 0.767  3.36 / lambda; 
332 
double d = log(c)  log(b)  lambda;

333 
double U, V, Y;

334  
335 
do

336 
{ 
337 
do

338 
{ 
339 
U = rngs[rngId]>doubleRand(); 
340 
Y = (a  log((1.0  U) / U)) / b; 
341 
} 
342 
while (Y <= 0.5); 
343  
344 
X = (int)floor(Y + 0.5); 
345 
V = rngs[rngId]>doubleRand(); 
346 
} 
347 
while (a  b * Y + log(V / 1.0 + pow(exp(a  b * Y), 2.0)) > d + X * log(lambda)  log(double(X))); 
348 
} 
349 
else

350 
{ 
351 
double a = exp(lambda);

352 
double p = 1.0; 
353 
X = 1;

354  
355 
while (p > a)

356 
{ 
357 
p *= rngs[rngId]>doubleRand(); 
358 
X++; 
359 
} 
360 
} 
361 
return X;

362 
} 
363  
364  
365 
// compatibility genk_ versions:

366  
367 
double cNumberGenerator::genk_uniform(double rng, double a, double b) 
368 
{ 
369 
return uniform(a, b, (int)rng); 
370 
} 
371  
372 
double cNumberGenerator::genk_intuniform(double rng, double a, double b) 
373 
{ 
374 
return intuniform((int)a, (int)b, (int)rng); 
375 
} 
376  
377 
double cNumberGenerator::genk_exponential(double rng, double p) 
378 
{ 
379 
return exponential(p, (int)rng); 
380 
} 
381  
382 
double cNumberGenerator::genk_normal(double rng, double m, double d) 
383 
{ 
384 
return normal(m, d, (int)rng); 
385 
} 
386  
387 
double cNumberGenerator::genk_truncnormal(double rng, double m, double d) 
388 
{ 
389 
return truncnormal(m, d, (int)rng); 
390 
} 
391  
392  
393  
394 
long cNumberGenerator::intrand(long r) 
395 
{ 
396 
return rngs[0]>intRand(r); 
397 
} 
398  
399 
double cNumberGenerator::dblrand()

400 
{ 
401 
return rngs[0]>doubleRand(); 
402  
403 
} 