Statistics
| Branch: | Revision:

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.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
}