Project

General

Profile

Statistics
| Branch: | Revision:

root / src / sim / cnumgen.cc @ a3be1d55

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