## root / src / sim / cnumgen.cc @ fbe00e73

History | View | Annotate | Download (9.17 KB)

1 | 01873262 | Georg Kunz | ```
//==========================================================================
``` |
---|---|---|---|

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 | f857128b | Simon Tenbusch | ```
cNumberGenerator::cNumberGenerator() : nrRNGs(0) {
``` |

38 | 01873262 | Georg Kunz | ```
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 | f857128b | Simon Tenbusch | for (int i = 0; i < nrRNGs; i++) |

49 | 01873262 | Georg Kunz | ```
delete rngs[i];
``` |

50 | |||

51 | } |
||

52 | |||

53 | void cNumberGenerator::checkBounds(int rngId) { |
||

54 | f857128b | Simon Tenbusch | ```
if (nrRNGs <= rngId)
``` |

55 | 01873262 | Georg Kunz | ```
throw cRuntimeError(
``` |

56 | f857128b | Simon Tenbusch | ```
"cNumberGenerator only has %d < %d local generators.", nrRNGs,
``` |

57 | 01873262 | Georg Kunz | rngId); |

58 | } |
||

59 | f857128b | Simon Tenbusch | |

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 | 01873262 | Georg Kunz | ```
//----------------------------------------------------------------------------
``` |

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.0-v+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(b-a+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 | } |