COPASI API  4.16.103
CRandom.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2014 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2002 - 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 #include <time.h>
16 #ifdef WIN32
17 # define _USE_MATH_DEFINES
18 # include <Windows.h>
19 #else
20 # include <unistd.h>
21 # include <sys/syscall.h>
22 #endif // WIN32
23 
24 #include <cmath>
25 #include <algorithm>
26 #include <string.h>
27 
28 #include "copasi.h"
29 #include "CRandom.h"
31 #include "utilities/CopasiTime.h"
32 
33 const std::string CRandom::TypeName[] =
34 {
35  "r250",
36  "Mersenne Twister",
37  "Mersenne Twister (HR)",
38  ""
39 };
40 
41 const char * CRandom::XMLType[] =
42 {
43  "r250",
44  "MersenneTwister",
45  "MersenneTwisterHR",
46  NULL
47 };
48 
50  unsigned C_INT32 seed)
51 {
52  CRandom * RandomGenerator = NULL;
53 
54  if (!seed)
55  seed = getSystemSeed();
56 
57  switch (type)
58  {
59  case r250:
60  RandomGenerator = new Cr250(seed);
61  RandomGenerator->mType = type;
62  break;
63 
64  case mt19937:
65  RandomGenerator = new Cmt19937(seed);
66  RandomGenerator->mType = type;
67  break;
68 
69  case mt19937HR:
70  RandomGenerator = new Cmt19937HR(seed);
71  RandomGenerator->mType = type;
72  break;
73 
74  default:
75  RandomGenerator = new Cmt19937(seed);
76  RandomGenerator->mType = type;
77  break;
78  }
79 
80  return RandomGenerator;
81 }
82 
84  mNumberU(0),
85  mNumberS(0),
86  mFloat(0.0),
87  mType(CRandom::unkown),
88  mModulus(1),
89  mModulusInv(1.0),
90  mModulusInv1(1.0)
91 {
92  varp.a0 = -0.5;
93  varp.a1 = 0.3333333;
94  varp.a2 = -0.2500068;
95  varp.a3 = 0.2000118;
96  varp.a4 = -0.1661269;
97  varp.a5 = 0.1421878;
98  varp.a6 = -0.1384794;
99  varp.a7 = 0.125006;
100  varp.muold = -1.0E37;
101  varp.muprev = -1.0E37;
102  varp.fact[0] = 1.0;
103  varp.fact[1] = 1.0;
104  varp.fact[2] = 2.0;
105  varp.fact[3] = 6.0;
106  varp.fact[4] = 24.0;
107  varp.fact[5] = 120.0;
108  varp.fact[6] = 720.0;
109  varp.fact[7] = 5040.0;
110  varp.fact[8] = 40320.0;
111  varp.fact[9] = 362880.0;
112 
113  vare.q[0] = 0.6931472;
114  vare.q[1] = 0.9333737;
115  vare.q[2] = 0.9888778;
116  vare.q[3] = 0.9984959;
117  vare.q[4] = 0.9998293;
118  vare.q[5] = 0.9999833;
119  vare.q[6] = 0.9999986;
120  vare.q[7] = 0.9999999;
121 
122  vare.q1 = vare.q;
123 }
124 
126 const CRandom::Type & CRandom::getType() const {return mType;}
127 const unsigned C_INT32 & CRandom::getModulus() const {return mModulus;}
128 
129 void CRandom::setModulus(const unsigned C_INT32 & modulus)
130 {
131  mModulus = modulus;
132  mModulusInv = 1.0 / mModulus;
133  mModulusInv1 = 1.0 / (mModulus + 1.0);
134 }
135 
137 {
138  unsigned C_INT32 ThreadId = 0;
139 
140 #ifdef WIN32
141  ThreadId = (unsigned C_INT32)(GetCurrentThreadId() & 0xffffffffUL);
142 #elif defined(SYS_thread_selfid)
143  ThreadId = (unsigned C_INT32)(::syscall(SYS_thread_selfid) & 0xffffffffUL);
144 #elif defined(SYS_gettid)
145  ThreadId = (unsigned C_INT32)(::syscall(SYS_gettid) & 0xffffffffUL);
146 #elif defined(SYS_getthrid)
147  ThreadId = (unsigned C_INT32)(syscall(SYS_getthrid) & 0xffffffffUL);
148 #endif
149 
150  // Invert Byte order so that we do not have accidental cancellations since both time and thread id are incremented.
151  ThreadId = (ThreadId & 0x000000ffUL) << 24 | (ThreadId & 0x0000ff00UL) << 8 |
152  (ThreadId & 0x00ff0000UL) >> 8 | (ThreadId & 0xff000000UL) >> 24;
153 
154  unsigned C_INT32 Time = (unsigned C_INT32)(CCopasiTimeVariable::getCurrentWallTime().getMicroSeconds() & 0xffffffffUL);
155 
156  // We use XOR so that we do not favor set or unset bits.
157  unsigned C_INT32 Seed = ThreadId ^ Time;
158 
159  return Seed;
160 }
161 
163 {
164  /* Every random number generator has to implement this. */
165  fatalError();
166  return;
167 }
168 
169 /**
170  * Get a random number in 0 <= n <= Modulus
171  * @return unsigned C_INT32 random
172  */
174 {
175  /* Every random number generator has to implement this. */
176  fatalError();
177  return mNumberU;
178 }
179 
180 /**
181  * Get a random number in 0 <= n <= (Modulus & 0x7ffffff)
182  * @return C_INT32 random
183  */
185 {
186  return getRandomU(std::min(mModulus, (unsigned C_INT32) 0x7ffffff));
187  // The method below only works for mModulus = 0xffffffff.
188  // return getRandomU() & 0x7ffffff;
189 }
190 
191 /**
192  * Get a random number in 0 <= n <= max
193  * Note: max must be smaller than Modulus (no check for performance)
194  * @param const unsigned C_INT32 & max
195  * @return unsigned C_INT32 random
196  */
197 unsigned C_INT32 CRandom::getRandomU(const unsigned C_INT32 & max)
198 {
199  unsigned C_INT32 Max = max + 1;
200  unsigned C_INT32 Limit = (mModulus / Max) * Max - 1;
201  unsigned C_INT32 NumberU;
202 
203  do
204  NumberU = getRandomU();
205 
206  while (NumberU >= Limit);
207 
208  return NumberU % Max;
209 }
210 
211 /**
212  * Get a random number in 0 <= n <= max
213  * Note: max must be smaller than Modulus (no check for performance)
214  * @param const C_INT32 & max
215  * @return C_INT32 random
216  */
218 {
219  unsigned C_INT32 Max = max + 1;
220  unsigned C_INT32 Limit = (mModulus / Max) * Max - 1;
221  unsigned C_INT32 NumberU;
222 
223  do
224  NumberU = getRandomU();
225 
226  while (NumberU >= Limit);
227 
228  return NumberU % Max;
229 }
230 
231 /**
232  * Produces a uniformly distributed random number in 0 <= x <= 1.
233  * @return C_FLOAT64 random
234  */
236 {
237  return getRandomU() * mModulusInv;
238 }
239 
240 /**
241  * Produces a uniformly distributed random number in 0 <= x < 1.
242  * Note: 0 < x <= 1 may be achieved by 1.0 - getRandomCO().
243  * @return C_FLOAT64 random
244  */
246 {
247  return getRandomU() * mModulusInv1;
248 }
249 
250 /**
251  * Produces a uniformly distributed random number in 0 < x < 1.
252  * @return C_FLOAT64 random
253  */
255 {
256  return (getRandomU() + .5) * mModulusInv1;
257 }
258 
259 /* a normal distribution with mean 0 and S.D. 1 */
260 
262 {
263  static bool HaveValue = true;
264  static C_FLOAT64 SavedValue;
265  C_FLOAT64 a, b, s;
266 
267  /* negate the flag */
268  HaveValue = !HaveValue;
269 
270  /* return the stored number (if one is there) */
271  if (HaveValue) return SavedValue;
272 
273  do
274  {
275  a = 2.0 * getRandomOO() - 1.0;
276  b = 2.0 * getRandomOO() - 1.0;
277  s = a * a + b * b;
278  }
279  while (s >= 1.0 || s == 0);
280 
281  s = sqrt(-2.0 * log(s) / s);
282 
283  // save one of the numbers for the next time
284  SavedValue = s * a;
285 
286  // and return the other
287  return s * b;
288 }
289 
290 /* a normal distribution with mean m and S.D. sd */
291 
293  const C_FLOAT64 & sd)
294 {return getRandomNormal01() * sd + mean;}
295 
296 /* a strictly positive normal distribution with mean m and S.D. sd */
297 
299  const C_FLOAT64 & sd)
300 {
301  C_FLOAT64 x;
302 
303  while ((x = getRandomNormal(mean, sd)) < 0);
304 
305  return x;
306 }
307 
308 /* a tentative normal distribution in logarithmic space */
309 
311  const C_FLOAT64 & sd)
312 {return mean * pow(10.0, getRandomNormal01() * sd);}
313 
315 {
316  if (mu == varp.muprev) goto S10;
317 
318  if (mu < 10.0) goto S120;
319 
320  varp.muprev = mu;
321  varp.s = sqrt(mu);
322  varp.d = 6.0 * mu * mu;
323  varp.ll = (long)(mu - 1.1484);
324 S10:
325  varp.g = mu + varp.s * getRandomNormal01();
326 
327  if (varp.g < 0.0) goto S20;
328 
329  varp.ignpoi = (long)(varp.g);
330 
331  if (varp.ignpoi >= varp.ll) return varp.ignpoi;
332 
333  varp.fk = (float)varp.ignpoi;
334  varp.difmuk = mu - varp.fk;
335  varp.u = getRandomCC();
336 
337  if (varp.d * varp.u >= varp.difmuk * varp.difmuk * varp.difmuk) return varp.ignpoi;
338 
339 S20:
340 
341  if (mu == varp.muold) goto S30;
342 
343  varp.muold = mu;
344  varp.omega = 0.3989423 / varp.s;
345  varp.b1 = 4.166667E-2 / mu;
346  varp.b2 = 0.3 * varp.b1 * varp.b1;
347  varp.c3 = 0.1428571 * varp.b1 * varp.b2;
348  varp.c2 = varp.b2 - 15.0 * varp.c3;
349  varp.c1 = varp.b1 - 6.0 * varp.b2 + 45.0 * varp.c3;
350  varp.c0 = 1.0 - varp.b1 + 3.0 * varp.b2 - 15.0 * varp.c3;
351  varp.c = 0.1069 / mu;
352 S30:
353 
354  if (varp.g < 0.0) goto S50;
355 
356  varp.kflag = 0;
357  goto S70;
358 S40:
359 
360  if (varp.fy - varp.u * varp.fy <= varp.py * exp(varp.px - varp.fx)) return varp.ignpoi;
361 
362 S50:
363  varp.e = getRandomExp();
364  varp.u = getRandomCC();
365  varp.u += (varp.u - 1.0);
366  varp. t = 1.8 + fsign(varp.e, varp.u);
367 
368  if (varp.t <= -0.6744) goto S50;
369 
370  varp.ignpoi = (long)(mu + varp.s * varp.t);
371  varp.fk = (float)varp.ignpoi;
372  varp.difmuk = mu - varp.fk;
373  varp.kflag = 1;
374  goto S70;
375 S60:
376 
377  if (varp.c * fabs(varp.u) > varp.py * exp(varp.px + varp.e) - varp.fy * exp(varp.fx + varp.e)) goto S50;
378 
379  return varp.ignpoi;
380 S70:
381 
382  if (varp.ignpoi >= 10) goto S80;
383 
384  varp.px = -mu;
385  varp.py = pow(mu, (double)varp.ignpoi) / *(varp.fact + varp.ignpoi);
386  goto S110;
387 S80:
388  varp.del = 8.333333E-2 / varp.fk;
389  varp.del -= (4.8 * varp.del * varp.del * varp.del);
390  varp.v = varp.difmuk / varp.fk;
391 
392  if (fabs(varp.v) <= 0.25) goto S90;
393 
394  varp.px = varp.fk * log(1.0 + varp.v) - varp.difmuk - varp.del;
395  goto S100;
396 S90:
397  varp.px = varp.fk * varp.v * varp.v * (((((((varp.a7 * varp.v + varp.a6) * varp.v + varp.a5) * varp.v\
398  + varp.a4) * varp.v + varp.a3) * varp.v + varp.a2) * varp.v + varp.a1) * varp.v + varp.a0) - varp.del;
399 S100:
400  varp.py = 0.3989423 / sqrt(varp.fk);
401 S110:
402  varp.x = (0.5 - varp.difmuk) / varp.s;
403  varp.xx = varp.x * varp.x;
404  varp.fx = -0.5 * varp.xx;
405  varp.fy = varp.omega * (((varp.c3 * varp.xx + varp.c2) * varp.xx + varp.c1) * varp.xx + varp.c0);
406 
407  if (varp.kflag <= 0) goto S40;
408 
409  goto S60;
410 S120:
411  varp.muprev = -1.0E37;
412 
413  if (mu == varp.muold) goto S130;
414 
415  if (mu >= 0.0) goto S125;
416 
417  // instead of quitting COPASI return a NaN
418  return std::numeric_limits<double>::quiet_NaN();
419 
420 S125:
421  varp.muold = mu;
422  varp.m = std::max(1L, (long)(mu));
423  varp.l = 0;
424  varp.p = exp(-mu);
425  varp.q = varp.p0 = varp.p;
426 S130:
427  varp.u = getRandomCC();
428  varp.ignpoi = 0;
429 
430  if (varp.u <= varp.p0) return varp.ignpoi;
431 
432  if (varp.l == 0) goto S150;
433 
434  varp.j = 1;
435 
436  if (varp.u > 0.458) varp.j = std::min(varp.l, varp.m);
437 
438  for (varp.k = varp.j; varp.k <= varp.l; varp.k++)
439  {
440  if (varp.u <= *(varp.pp + varp.k - 1)) goto S180;
441  }
442 
443  if (varp.l == 35) goto S130;
444 
445 S150:
446  varp.l += 1;
447 
448  for (varp.k = varp.l; varp.k <= 35; varp.k++)
449  {
450  varp.p = varp.p * mu / (float)varp.k;
451  varp.q += varp.p;
452  *(varp.pp + varp.k - 1) = varp.q;
453 
454  if (varp.u <= varp.q) goto S170;
455  }
456 
457  varp.l = 35;
458  goto S130;
459 S170:
460  varp.l = varp.k;
461 S180:
462  varp.ignpoi = varp.k;
463  return varp.ignpoi;
464 }
465 
467 {
468  if ((sign > 0.0f && num < 0.0f) || (sign < 0.0f && num > 0.0f))
469  return - num;
470  else return num;
471 }
472 
474 {
475  vare.a = 0.0;
476  vare.u = getRandomCC();
477  goto S30;
478 S20:
479  vare.a += *vare.q1;
480 S30:
481  vare.u += vare.u;
482 
483  if (vare.u < 1.0) goto S20;
484 
485  vare.u -= 1.0;
486 
487  if (vare.u > *vare.q1) goto S60;
488 
489  vare.sexpo = vare.a + vare.u;
490  return vare.sexpo;
491 S60:
492  vare.i = 1;
493  vare.ustar = getRandomCC();
494  vare.umin = vare.ustar;
495 S70:
496  vare.ustar = getRandomCC();
497 
498  if (vare.ustar < vare.umin) vare.umin = vare.ustar;
499 
500  vare.i += 1;
501 
502  if (vare.u > *(vare.q + vare.i - 1)) goto S70;
503 
504  vare.sexpo = vare.a + vare.umin**vare.q1;
505  return vare.sexpo;
506 }
507 
509 {
510  if (a < 1.0)
511  {
512  C_FLOAT64 v1, v2, v3;
513  C_FLOAT64 x, e;
514 
515  C_FLOAT64 v0 = M_E / (M_E + a);
516 
517  do
518  {
519  v1 = 1.0 - getRandomCO();
520  v2 = 1.0 - getRandomCO();
521  v3 = 1.0 - getRandomCO();
522 
523  if (v1 <= v0)
524  {
525  x = pow(v2, (1.0 / a));
526  e = v3 * pow(x, (a - 1.0));
527  }
528  else
529  {
530  x = 1 - log(v2);
531  e = v3 * exp(-x);
532  }
533  }
534  while (e > pow(x, (a - 1.0)) * exp(-x));
535 
536  return x;
537  }
538  else
539  {
540  C_FLOAT64 d, c, x, v, u;
541 
542  d = a - 1.0 / 3.0;
543  c = 1.0 / sqrt(9.0 * d);
544 
545  while (true)
546  {
547  do
548  {
549  x = getRandomNormal01();
550  v = 1.0 + c * x;
551  }
552  while (v <= 0.0);
553 
554  v = v * v * v;
555  u = getRandomOO();
556 
557  if (u < 1.0 - 0.0331 * (x * x) * (x * x))
558  {
559  break;
560  }
561 
562  if (log(u) < 0.5 * x * x + d * (1.0 - v + log(v)))
563  {
564  break;
565  }
566  }
567 
568  return (d * v);
569  }
570 }
571 
573 {
574  return scale * getRandomStdGamma(shape);
575 }
virtual ~CRandom()
Definition: CRandom.cpp:125
C_FLOAT64 q[8]
Definition: CRandom.h:49
unsigned int & syscall(int arg1)
Definition: Cr250.h:23
C_FLOAT64 xx
Definition: CRandom.h:43
C_FLOAT64 sexpo
Definition: CRandom.h:51
static unsigned C_INT32 getSystemSeed()
Definition: CRandom.cpp:136
C_FLOAT64 del
Definition: CRandom.h:43
virtual C_FLOAT64 getRandomNormal(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:292
C_FLOAT64 u
Definition: CRandom.h:51
virtual C_FLOAT64 getRandomOO()
Definition: CRandom.cpp:254
C_FLOAT64 a7
Definition: CRandom.h:36
virtual C_FLOAT64 getRandomExp()
Definition: CRandom.cpp:473
virtual C_FLOAT64 getRandomStdGamma(C_FLOAT64 shape)
Definition: CRandom.cpp:508
const CRandom::Type & getType() const
Definition: CRandom.cpp:126
C_FLOAT64 a
Definition: CRandom.h:51
#define fatalError()
virtual C_INT32 getRandomS()
Definition: CRandom.cpp:184
C_INT32 i
Definition: CRandom.h:50
C_FLOAT64 c3
Definition: CRandom.h:43
C_FLOAT64 pp[35]
Definition: CRandom.h:43
C_FLOAT64 a0
Definition: CRandom.h:29
#define C_UNUSED(p)
Definition: copasi.h:220
#define C_INT32
Definition: copasi.h:90
C_FLOAT64 muold
Definition: CRandom.h:37
C_FLOAT64 a6
Definition: CRandom.h:35
C_FLOAT64 px
Definition: CRandom.h:43
static const char * XMLType[]
Definition: CRandom.h:72
C_FLOAT64 umin
Definition: CRandom.h:51
C_FLOAT64 c2
Definition: CRandom.h:43
C_FLOAT64 a1
Definition: CRandom.h:30
virtual C_FLOAT64 getRandomNormal01()
Definition: CRandom.cpp:261
const unsigned C_INT32 & getModulus() const
Definition: CRandom.cpp:127
static const std::string TypeName[]
Definition: CRandom.h:67
C_FLOAT64 mModulusInv1
Definition: CRandom.h:108
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
C_FLOAT64 a2
Definition: CRandom.h:31
virtual C_FLOAT64 getRandomPoisson(const C_FLOAT64 &mean)
Definition: CRandom.cpp:314
C_FLOAT64 fact[10]
Definition: CRandom.h:39
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
virtual C_FLOAT64 getRandomNormalPositive(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:298
virtual unsigned C_INT32 getRandomU()
Definition: CRandom.cpp:173
C_FLOAT64 mModulusInv
Definition: CRandom.h:103
C_FLOAT64 p0
Definition: CRandom.h:43
C_FLOAT64 muprev
Definition: CRandom.h:38
C_FLOAT64 b2
Definition: CRandom.h:43
C_FLOAT64 ustar
Definition: CRandom.h:51
C_FLOAT64 a3
Definition: CRandom.h:32
C_FLOAT64 c0
Definition: CRandom.h:43
PoissonVars varp
Definition: CRandom.h:112
void setModulus(const unsigned C_INT32 &modulus)
Definition: CRandom.cpp:129
static CCopasiTimeVariable getCurrentWallTime()
Definition: CopasiTime.cpp:160
C_FLOAT64 difmuk
Definition: CRandom.h:43
C_FLOAT64 fy
Definition: CRandom.h:43
CRandom()
Definition: CRandom.cpp:83
C_FLOAT64 omega
Definition: CRandom.h:43
C_FLOAT64 * q1
Definition: CRandom.h:52
virtual C_FLOAT64 getRandomCO()
Definition: CRandom.cpp:245
C_FLOAT64 fsign(C_FLOAT64 num, C_FLOAT64 sign)
Definition: CRandom.cpp:466
virtual C_FLOAT64 getRandomNormalLog(const C_FLOAT64 &mean, const C_FLOAT64 &sd)
Definition: CRandom.cpp:310
C_FLOAT64 c1
Definition: CRandom.h:43
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 a4
Definition: CRandom.h:33
ExpVars vare
Definition: CRandom.h:113
C_INT64 getMicroSeconds(const bool &bounded=false) const
Definition: CopasiTime.cpp:112
C_FLOAT64 a5
Definition: CRandom.h:34
CRandom::Type mType
Definition: CRandom.h:93
C_FLOAT64 fx
Definition: CRandom.h:43
virtual C_FLOAT64 getRandomGamma(C_FLOAT64 shape, C_FLOAT64 scale)
Definition: CRandom.cpp:572
if(!yymsg) yymsg
unsigned C_INT32 mModulus
Definition: CRandom.h:98
unsigned C_INT32 mNumberU
Definition: CRandom.h:78
C_FLOAT64 py
Definition: CRandom.h:43
C_FLOAT64 b1
Definition: CRandom.h:43
C_FLOAT64 fk
Definition: CRandom.h:43
#define min(a, b)
Definition: f2c.h:175
virtual void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: CRandom.cpp:162
#define max(a, b)
Definition: f2c.h:176