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

History | View | Annotate | Download (9.37 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 "../envir/envirbase.h" |
||

25 | #include <pthread.h> |
||

26 | |||

27 | ```
//
``` |
||

28 | ```
// Linux - <math.h> supplies a PI, MS does not...
``` |
||

29 | ```
//
``` |
||

30 | ```
#ifndef M_PI
``` |
||

31 | #define M_PI 3.1415926535897932384626433832795 |
||

32 | ```
#endif
``` |
||

33 | |||

34 | ```
#ifndef M_E
``` |
||

35 | #define M_E 2.7182818284590452353602874713527 |
||

36 | ```
#endif
``` |
||

37 | |||

38 | ```
cNumberGenerator::cNumberGenerator(unsigned nrRngs) :
``` |
||

39 | nrRngs(nrRngs) { |
||

40 | ```
if (simulation.isRunning()) {
``` |
||

41 | ```
throw cRuntimeError(
``` |
||

42 | ```
"Cannot (yet) create new instance of cNumberGenerator"
``` |
||

43 | ```
" at runtime due to potentially concurrent access to"
``` |
||

44 | ```
" the central seed generator.");
``` |
||

45 | } |
||

46 | |||

47 | ```
rngs = (cRNG**) new cMersenneTwister*[nrRngs];
``` |
||

48 | for (unsigned i = 0; i < nrRngs; i++) { |
||

49 | ```
rngs[i] = new cMersenneTwister();
``` |
||

50 | ```
EnvirBase* e = dynamic_cast<EnvirBase*> (&ev);
``` |
||

51 | unsigned long seed; |
||

52 | if (e != NULL) { |
||

53 | seed = e->getSeedGenerator()->intRand(); |
||

54 | ```
} else {
``` |
||

55 | ```
throw cRuntimeError(
``` |
||

56 | ```
"cNumberGenerator: Failed to initialize Number Generator, dynamic cast to EnvirBase failed.");
``` |
||

57 | } |
||

58 | |||

59 | rngs[i]->seed(seed); |
||

60 | } |
||

61 | } |
||

62 | |||

63 | |||

64 | cNumberGenerator::~cNumberGenerator() { |
||

65 | for (int i = 0; i < nrRngs; i++) |
||

66 | ```
delete rngs[i];
``` |
||

67 | |||

68 | } |
||

69 | |||

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

71 | ```
if (nrRngs < rngId)
``` |
||

72 | ```
throw cRuntimeError(
``` |
||

73 | ```
"cNumberGenerator only has %d < %d local generators.", nrRngs,
``` |
||

74 | rngId); |
||

75 | } |
||

76 | ```
//----------------------------------------------------------------------------
``` |
||

77 | ```
//
``` |
||

78 | ```
// C O N T I N U O U S
``` |
||

79 | ```
//
``` |
||

80 | ```
//----------------------------------------------------------------------------
``` |
||

81 | |||

82 | double cNumberGenerator::uniform(double a, double b, int rngId) { |
||

83 | checkBounds(rngId); |
||

84 | ```
return a + rngs[rngId]->doubleRand() * (b - a);
``` |
||

85 | |||

86 | } |
||

87 | |||

88 | double cNumberGenerator::exponential(double p, int rngId) { |
||

89 | checkBounds(rngId); |
||

90 | ```
return -p * log(rngs[rngId]->doubleRand());
``` |
||

91 | |||

92 | } |
||

93 | |||

94 | double cNumberGenerator::unit_normal(int rngId) { |
||

95 | checkBounds(rngId); |
||

96 | return sqrt(-2.0 * log(rngs[rngId]->doubleRand())) * cos(PI * 2 |
||

97 | * rngs[rngId]->doubleRand()); |
||

98 | |||

99 | } |
||

100 | |||

101 | double cNumberGenerator::normal(double m, double d, int rngId) { |
||

102 | checkBounds(rngId); |
||

103 | return m + d * sqrt(-2.0 * log(rngs[rngId]->doubleRand())) * cos(PI * 2 |
||

104 | * rngs[rngId]->doubleRand()); |
||

105 | |||

106 | } |
||

107 | |||

108 | double cNumberGenerator::truncnormal(double m, double d, int rngId) { |
||

109 | checkBounds(rngId); |
||

110 | ```
double res;
``` |
||

111 | ```
do {
``` |
||

112 | res = normal(m, d, rngId); |
||

113 | } while (res < 0); |
||

114 | ```
return res;
``` |
||

115 | } |
||

116 | |||

117 | ```
/*
``` |
||

118 | ```
* internal, for alpha>1.
``` |
||

119 | ```
*
``` |
||

120 | ```
* From: "A Simple Method for Generating Gamma Variables", George Marsaglia and
``` |
||

121 | ```
* Wai Wan Tsang, ACM Transactions on Mathematical Software, Vol. 26, No. 3,
``` |
||

122 | ```
* September 2000. Available online.
``` |
||

123 | ```
*/
``` |
||

124 | double cNumberGenerator::gamma_Marsaglia2000(double a, int rngId) |
||

125 | { |
||

126 | checkBounds(rngId); |
||

127 | ```
ASSERT(a>1);
``` |
||

128 | |||

129 | ```
double d,c,x,v,u;
``` |
||

130 | d = a - 1.0/3.0; |
||

131 | c = 1.0/sqrt(9.0*d); |
||

132 | ```
for(;;)
``` |
||

133 | { |
||

134 | do {x = unit_normal(rngId); v = 1.0 + c*x;} while (v<=0); |
||

135 | v = v*v*v; u = rngs[rngId]->doubleRand(); |
||

136 | if (u < 1.0 - 0.0331*(x*x)*(x*x)) |
||

137 | ```
return d*v;
``` |
||

138 | if (log(u)<0.5*x*x + d*(1.0-v+log(v))) |
||

139 | ```
return d*v;
``` |
||

140 | } |
||

141 | } |
||

142 | |||

143 | |||

144 | ```
/*
``` |
||

145 | ```
* internal, for alpha<1.
``` |
||

146 | ```
*
``` |
||

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

148 | ```
* end of Section 6 of the Marsaglia2000 paper (see above).
``` |
||

149 | ```
*/
``` |
||

150 | double cNumberGenerator::gamma_MarsagliaTransf(double alpha, int rngId) |
||

151 | { |
||

152 | checkBounds(rngId); |
||

153 | ```
ASSERT(alpha<1);
``` |
||

154 | |||

155 | return gamma_Marsaglia2000(1+alpha,rngId) * pow(rngs[rngId]->doubleRand(), 1/alpha); |
||

156 | } |
||

157 | |||

158 | double cNumberGenerator::gamma_d(double alpha, double theta, int rngId) |
||

159 | { |
||

160 | checkBounds(rngId); |
||

161 | if (alpha<=0 || theta<=0) |
||

162 | throw cRuntimeError("gamma(): alpha and theta params must be positive " |
||

163 | ```
"(alpha=%lg, theta=%lg)", alpha, theta);
``` |
||

164 | |||

165 | if (fabs(alpha - 1.0) <= DBL_EPSILON) |
||

166 | { |
||

167 | ```
return exponential(theta, rngId);
``` |
||

168 | } |
||

169 | else if (alpha < 1.0) |
||

170 | { |
||

171 | ```
return theta * gamma_MarsagliaTransf(alpha, rngId);
``` |
||

172 | } |
||

173 | else // if (alpha > 1.0) |
||

174 | { |
||

175 | ```
return theta * gamma_Marsaglia2000(alpha, rngId);
``` |
||

176 | } |
||

177 | } |
||

178 | |||

179 | |||

180 | double cNumberGenerator::beta(double alpha1, double alpha2, int rngId) |
||

181 | { |
||

182 | checkBounds(rngId); |
||

183 | if (alpha1<=0 || alpha2<=0) |
||

184 | throw cRuntimeError("beta(): alpha1 and alpha2 parameters must be positive " |
||

185 | ```
"(alpha1=%lg, alpha2=%lg)", alpha1, alpha2);
``` |
||

186 | |||

187 | double Y1 = gamma_d(alpha1, 1.0, rngId); |
||

188 | double Y2 = gamma_d(alpha2, 1.0, rngId); |
||

189 | |||

190 | ```
return Y1 / (Y1 + Y2);
``` |
||

191 | } |
||

192 | |||

193 | |||

194 | double cNumberGenerator::erlang_k(unsigned int k, double m, int rngId) |
||

195 | { |
||

196 | checkBounds(rngId); |
||

197 | double U = 1.0; |
||

198 | for (unsigned int i = 0; i < k; i++) |
||

199 | U *= (1.0 - rngs[rngId]->doubleRand()); |
||

200 | |||

201 | return -(m / (double) k) * log(U); |
||

202 | } |
||

203 | |||

204 | |||

205 | double cNumberGenerator::chi_square(unsigned int k, int rngId) |
||

206 | { |
||

207 | checkBounds(rngId); |
||

208 | if (!(k % 2)) |
||

209 | return erlang_k(k >> 1, k, rngId); |
||

210 | ```
else
``` |
||

211 | return gamma_d((double) k / 2.0, 2.0, rngId); |
||

212 | } |
||

213 | |||

214 | |||

215 | double cNumberGenerator::student_t(unsigned int i, int rngId) |
||

216 | { |
||

217 | checkBounds(rngId); |
||

218 | double Z = normal(0, 1, rngId); |
||

219 | double W = sqrt(chi_square(i, rngId) / (double) i); |
||

220 | ```
return Z / W;
``` |
||

221 | } |
||

222 | |||

223 | |||

224 | double cNumberGenerator::cauchy(double a, double b, int rngId) |
||

225 | { |
||

226 | checkBounds(rngId); |
||

227 | if (b<=0) |
||

228 | throw cRuntimeError("cauchy(): parameters must be b>0 (a=%lg, b=%lg)", a, b); |
||

229 | |||

230 | ```
return a + b * tan(M_PI * rngs[rngId]->doubleRand());
``` |
||

231 | |||

232 | } |
||

233 | |||

234 | |||

235 | double cNumberGenerator::triang(double a, double b, double c, int rngId) |
||

236 | { |
||

237 | checkBounds(rngId); |
||

238 | ```
if (b<a || c<b || a==c)
``` |
||

239 | throw cRuntimeError("triang(): parameters must be a<=b<=c, a<c (a=%lg, b=%lg, c=%lg)", a, b, c); |
||

240 | |||

241 | ```
double U, beta, T;
``` |
||

242 | |||

243 | U = rngs[rngId]->doubleRand(); |
||

244 | beta = (b - a) / (c - a); |
||

245 | |||

246 | ```
if (U < beta)
``` |
||

247 | T = sqrt(beta * U); |
||

248 | ```
else
``` |
||

249 | T = 1.0 - sqrt((1.0 - beta) * (1.0 - U)); |
||

250 | |||

251 | ```
return a + (c - a) * T;
``` |
||

252 | } |
||

253 | |||

254 | |||

255 | ```
// lognormal() is inline
``` |
||

256 | |||

257 | |||

258 | double cNumberGenerator::weibull(double a, double b, int rngId) |
||

259 | { |
||

260 | checkBounds(rngId); |
||

261 | if (a<=0 || b<=0) |
||

262 | throw cRuntimeError("weibull(): a,b parameters must be positive (a=%lg, b=%lg)", a, b); |
||

263 | |||

264 | return a * pow(-log(1.0 - rngs[rngId]->doubleRand()), 1.0 / b); |
||

265 | } |
||

266 | |||

267 | |||

268 | double cNumberGenerator::pareto_shifted(double a, double b, double c, int rngId) |
||

269 | { |
||

270 | checkBounds(rngId); |
||

271 | if (a==0) |
||

272 | throw cRuntimeError("pareto_shifted(): parameter a cannot be zero)"); |
||

273 | |||

274 | double u_pow = pow(1.0 - rngs[rngId]->doubleRand(), 1.0 / a); |
||

275 | |||

276 | ```
return (b - c * u_pow) / u_pow;
``` |
||

277 | } |
||

278 | |||

279 | ```
//----------------------------------------------------------------------------
``` |
||

280 | ```
//
``` |
||

281 | ```
// D I S C R E T E
``` |
||

282 | ```
//
``` |
||

283 | ```
//----------------------------------------------------------------------------
``` |
||

284 | |||

285 | |||

286 | int cNumberGenerator::intuniform(int a, int b, int rngId) |
||

287 | { |
||

288 | checkBounds(rngId); |
||

289 | return a + rngs[rngId]->intRand(b-a+1); |
||

290 | } |
||

291 | |||

292 | |||

293 | ```
// bernoulli() is inline
``` |
||

294 | |||

295 | |||

296 | int cNumberGenerator::binomial(int n, double p, int rngId) |
||

297 | { |
||

298 | checkBounds(rngId); |
||

299 | int X = 0; |
||

300 | ```
// sum up n bernoulli trials
``` |
||

301 | for (int i = 0; i < n; i++) |
||

302 | { |
||

303 | ```
double U = rngs[rngId]->doubleRand();
``` |
||

304 | ```
if (p > U)
``` |
||

305 | X++; |
||

306 | } |
||

307 | ```
return X;
``` |
||

308 | } |
||

309 | |||

310 | |||

311 | int cNumberGenerator::geometric(double p, int rngId) |
||

312 | { |
||

313 | checkBounds(rngId); |
||

314 | double a = 1.0 / (log(1.0 - p)); |
||

315 | return (int)floor(a * log(rngs[rngId]->doubleRand())); |
||

316 | } |
||

317 | |||

318 | |||

319 | int cNumberGenerator::negbinomial(int n, double p, int rngId) |
||

320 | { |
||

321 | checkBounds(rngId); |
||

322 | int X = 0; |
||

323 | for (int i = 0; i < n; i++) |
||

324 | { |
||

325 | X += geometric(p, rngId); |
||

326 | } |
||

327 | ```
return X;
``` |
||

328 | } |
||

329 | |||

330 | |||

331 | int cNumberGenerator::poisson(double lambda, int rngId) |
||

332 | { |
||

333 | checkBounds(rngId); |
||

334 | ```
int X;
``` |
||

335 | if (lambda > 30.0) |
||

336 | { |
||

337 | double a = M_PI * sqrt(lambda / 3.0); |
||

338 | ```
double b = a / lambda;
``` |
||

339 | double c = 0.767 - 3.36 / lambda; |
||

340 | ```
double d = log(c) - log(b) - lambda;
``` |
||

341 | ```
double U, V, Y;
``` |
||

342 | |||

343 | ```
do
``` |
||

344 | { |
||

345 | ```
do
``` |
||

346 | { |
||

347 | U = rngs[rngId]->doubleRand(); |
||

348 | Y = (a - log((1.0 - U) / U)) / b; |
||

349 | } |
||

350 | while (Y <= -0.5); |
||

351 | |||

352 | X = (int)floor(Y + 0.5); |
||

353 | V = rngs[rngId]->doubleRand(); |
||

354 | } |
||

355 | while (a - b * Y + log(V / 1.0 + pow(exp(a - b * Y), 2.0)) > d + X * log(lambda) - log(double(X))); |
||

356 | } |
||

357 | ```
else
``` |
||

358 | { |
||

359 | ```
double a = exp(-lambda);
``` |
||

360 | double p = 1.0; |
||

361 | ```
X = -1;
``` |
||

362 | |||

363 | ```
while (p > a)
``` |
||

364 | { |
||

365 | p *= rngs[rngId]->doubleRand(); |
||

366 | X++; |
||

367 | } |
||

368 | } |
||

369 | ```
return X;
``` |
||

370 | } |
||

371 | |||

372 | |||

373 | ```
// compatibility genk_ versions:
``` |
||

374 | |||

375 | double cNumberGenerator::genk_uniform(double rng, double a, double b) |
||

376 | { |
||

377 | return uniform(a, b, (int)rng); |
||

378 | } |
||

379 | |||

380 | double cNumberGenerator::genk_intuniform(double rng, double a, double b) |
||

381 | { |
||

382 | return intuniform((int)a, (int)b, (int)rng); |
||

383 | } |
||

384 | |||

385 | double cNumberGenerator::genk_exponential(double rng, double p) |
||

386 | { |
||

387 | return exponential(p, (int)rng); |
||

388 | } |
||

389 | |||

390 | double cNumberGenerator::genk_normal(double rng, double m, double d) |
||

391 | { |
||

392 | return normal(m, d, (int)rng); |
||

393 | } |
||

394 | |||

395 | double cNumberGenerator::genk_truncnormal(double rng, double m, double d) |
||

396 | { |
||

397 | return truncnormal(m, d, (int)rng); |
||

398 | } |
||

399 | |||

400 | |||

401 | |||

402 | long cNumberGenerator::intrand(long r) |
||

403 | { |
||

404 | return rngs[0]->intRand(r); |
||

405 | } |
||

406 | |||

407 | ```
double cNumberGenerator::dblrand()
``` |
||

408 | { |
||

409 | return rngs[0]->doubleRand(); |
||

410 | |||

411 | } |