IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
random.cpp
Go to the documentation of this file.
1 
29 #include <itpp/base/random.h>
30 #include <itpp/base/itcompat.h>
32 
33 
34 namespace itpp
35 {
36 
37 // ----------------------------------------------------------------------
38 // Random_Generator (DSFMT_19937_RNG)
39 // ----------------------------------------------------------------------
40 
41 template <>
42 bool Random_Generator::initialized = false;
43 template <>
44 bool Random_Generator::bigendian = is_bigendian();
45 template <>
46 Random_Generator::w128_t Random_Generator::status[N + 1] = { };
47 template <>
48 int Random_Generator::idx = 0;
49 template <>
50 unsigned int Random_Generator::last_seed = 0U;
51 #if defined(__SSE2__)
52 template <>
53 __m128i Random_Generator::sse2_param_mask = _mm_set_epi32(0, 0, 0, 0);
54 #endif // __SSE2__
55 
56 
57 // Set the seed of the Global Random Number Generator
58 void RNG_reset(unsigned int seed)
59 {
60  Random_Generator RNG;
61  RNG.reset(seed);
62 }
63 
64 // Set the seed of the Global Random Number Generator to the same as last time
65 void RNG_reset()
66 {
67  Random_Generator RNG;
68  RNG.reset();
69 }
70 
71 // Set a random seed for the Global Random Number Generator
73 {
74  Random_Generator RNG;
75  RNG.randomize();
76 }
77 
78 // Save current full state of generator in memory
79 void RNG_get_state(ivec &state)
80 {
81  Random_Generator RNG;
82  state = RNG.get_state();
83 }
84 
85 // Resume the state saved in memory
86 void RNG_set_state(const ivec &state)
87 {
88  Random_Generator RNG;
89  RNG.set_state(state);
90 }
91 
93 // I_Uniform_RNG
95 
97 {
98  setup(min, max);
99 }
100 
102 {
103  if (min <= max) {
104  lo = min;
105  hi = max;
106  }
107  else {
108  lo = max;
109  hi = min;
110  }
111 }
112 
113 void I_Uniform_RNG::get_setup(int &min, int &max) const
114 {
115  min = lo;
116  max = hi;
117 }
118 
120 {
121  ivec vv(n);
122 
123  for (int i=0; i < n; i++)
124  vv(i) = sample();
125 
126  return vv;
127 }
128 
129 imat I_Uniform_RNG::operator()(int h, int w)
130 {
131  imat mm(h, w);
132  int i, j;
133 
134  for (i = 0; i < h; i++)
135  for (j = 0; j < w; j++)
136  mm(i, j) = sample();
137 
138  return mm;
139 }
140 
142 // Uniform_RNG
144 
146 {
147  setup(min, max);
148 }
149 
150 void Uniform_RNG::setup(double min, double max)
151 {
152  if (min <= max) {
153  lo_bound = min;
154  hi_bound = max;
155  }
156  else {
157  lo_bound = max;
158  hi_bound = min;
159  }
160 }
161 
162 void Uniform_RNG::get_setup(double &min, double &max) const
163 {
164  min = lo_bound;
165  max = hi_bound;
166 }
167 
169 // Exp_RNG
171 
173 {
174  setup(lambda);
175 }
176 
178 {
179  vec vv(n);
180 
181  for (int i=0; i < n; i++)
182  vv(i) = sample();
183 
184  return vv;
185 }
186 
188 {
189  mat mm(h, w);
190  int i, j;
191 
192  for (i = 0; i < h; i++)
193  for (j = 0; j < w; j++)
194  mm(i, j) = sample();
195 
196  return mm;
197 }
198 
200 // Gamma_RNG
202 
204 {
205  vec vv(n);
206  for (int i = 0; i < n; i++)
207  vv(i) = sample();
208  return vv;
209 }
210 
211 mat Gamma_RNG::operator()(int r, int c)
212 {
213  mat mm(r, c);
214  for (int i = 0; i < r * c; i++)
215  mm(i) = sample();
216  return mm;
217 }
218 
220 {
221  // A copy of rgamma code from the R package, adapted to IT++ by Vasek
222  // Smidl
223 
224  /* Constants : */
225  const static double sqrt32 = 5.656854;
226  const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
227 
228  const static double q1 = 0.04166669;
229  const static double q2 = 0.02083148;
230  const static double q3 = 0.00801191;
231  const static double q4 = 0.00144121;
232  const static double q5 = -7.388e-5;
233  const static double q6 = 2.4511e-4;
234  const static double q7 = 2.424e-4;
235 
236  const static double a1 = 0.3333333;
237  const static double a2 = -0.250003;
238  const static double a3 = 0.2000062;
239  const static double a4 = -0.1662921;
240  const static double a5 = 0.1423657;
241  const static double a6 = -0.1367177;
242  const static double a7 = 0.1233795;
243 
244  /* State variables [FIXME for threading!] :*/
245  static double aa = 0.;
246  static double aaa = 0.;
247  static double s, s2, d; /* no. 1 (step 1) */
248  static double q0, b, si, c;/* no. 2 (step 4) */
249 
250  double e, p, q, r, t, u, v, w, x, ret_val;
251  double a = alpha;
252  double scale = 1.0 / beta;
253 
254  it_error_if(!std::isfinite(a) || !std::isfinite(scale) || (a < 0.0)
255  || (scale <= 0.0), "Gamma_RNG wrong parameters");
256 
257  if (a < 1.) { /* GS algorithm for parameters a < 1 */
258  if (a == 0)
259  return 0.;
260  e = 1.0 + exp_m1 * a;
261  for (;;) { //VS repeat
262  p = e * RNG.genrand_open_open();
263  if (p >= 1.0) {
264  x = -std::log((e - p) / a);
265  if (-std::log(RNG.genrand_open_close()) >= (1.0 - a) * std::log(x))
266  break;
267  }
268  else {
269  x = std::exp(std::log(p) / a);
270  if (-std::log(RNG.genrand_open_close()) >= x)
271  break;
272  }
273  }
274  return scale * x;
275  }
276 
277  /* --- a >= 1 : GD algorithm --- */
278 
279  /* Step 1: Recalculations of s2, s, d if a has changed */
280  if (a != aa) {
281  aa = a;
282  s2 = a - 0.5;
283  s = std::sqrt(s2);
284  d = sqrt32 - s * 12.0;
285  }
286  /* Step 2: t = standard normal deviate, x = (s,1/2) -normal deviate. */
287  /* immediate acceptance (i) */
288  t = NRNG.sample();
289  x = s + 0.5 * t;
290  ret_val = x * x;
291  if (t >= 0.0)
292  return scale * ret_val;
293 
294  /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
295  u = RNG.genrand_close_open();
296  if ((d * u) <= (t * t * t))
297  return scale * ret_val;
298 
299  /* Step 4: recalculations of q0, b, si, c if necessary */
300  if (a != aaa) {
301  aaa = a;
302  r = 1.0 / a;
303  q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
304  + q2) * r + q1) * r;
305 
306  /* Approximation depending on size of parameter a */
307  /* The constants in the expressions for b, si and c */
308  /* were established by numerical experiments */
309  if (a <= 3.686) {
310  b = 0.463 + s + 0.178 * s2;
311  si = 1.235;
312  c = 0.195 / s - 0.079 + 0.16 * s;
313  }
314  else if (a <= 13.022) {
315  b = 1.654 + 0.0076 * s2;
316  si = 1.68 / s + 0.275;
317  c = 0.062 / s + 0.024;
318  }
319  else {
320  b = 1.77;
321  si = 0.75;
322  c = 0.1515 / s;
323  }
324  }
325 
326  /* Step 5: no quotient test if x not positive */
327  if (x > 0.0) {
328  /* Step 6: calculation of v and quotient q */
329  v = t / (s + s);
330  if (std::fabs(v) <= 0.25)
331  q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
332  + a3) * v + a2) * v + a1) * v;
333  else
334  q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
335 
336  /* Step 7: quotient acceptance (q) */
337  if (log(1.0 - u) <= q)
338  return scale * ret_val;
339  }
340 
341  for (;;) { //VS repeat
342  /* Step 8: e = standard exponential deviate
343  * u = 0,1 -uniform deviate
344  * t = (b,si)-double exponential (laplace) sample */
345  e = -std::log(RNG.genrand_open_close()); //see Exponential_RNG
346  u = RNG.genrand_open_close();
347  u = u + u - 1.0;
348  if (u < 0.0)
349  t = b - si * e;
350  else
351  t = b + si * e;
352  /* Step 9: rejection if t < tau(1) = -0.71874483771719 */
353  if (t >= -0.71874483771719) {
354  /* Step 10: calculation of v and quotient q */
355  v = t / (s + s);
356  if (std::fabs(v) <= 0.25)
357  q = q0 + 0.5 * t * t *
358  ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
359  + a2) * v + a1) * v;
360  else
361  q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
362  /* Step 11: hat acceptance (h) */
363  /* (if q not positive go to step 8) */
364  if (q > 0.0) {
365  // Try to use w = expm1(q); (Not supported on w32)
366  w = expm1(q);
367  /* ^^^^^ original code had approximation with rel.err < 2e-7 */
368  /* if t is rejected sample again at step 8 */
369  if ((c * std::fabs(u)) <= (w * std::exp(e - 0.5 * t * t)))
370  break;
371  }
372  }
373  } /* repeat .. until `t' is accepted */
374  x = s + 0.5 * t;
375  return scale * x * x;
376 }
377 
379 // Normal_RNG
381 
382 void Normal_RNG::get_setup(double &meanval, double &variance) const
383 {
384  meanval = mean;
385  variance = sigma * sigma;
386 }
387 
388 // (Ziggurat) tabulated values for the heigt of the Ziggurat levels
389 const double Normal_RNG::ytab[128] = {
390  1, 0.963598623011, 0.936280813353, 0.913041104253,
391  0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
392  0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
393  0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
394  0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
395  0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
396  0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
397  0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
398  0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
399  0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
400  0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
401  0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
402  0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
403  0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
404  0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
405  0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
406  0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
407  0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
408  0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
409  0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
410  0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
411  0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
412  0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
413  0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
414  0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
415  0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
416  0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
417  0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
418  0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
419  0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
420  0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
421  0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
422 };
423 
424 /*
425  * (Ziggurat) tabulated values for 2^24 times x[i]/x[i+1], used to accept
426  * for U*x[i+1]<=x[i] without any floating point operations
427  */
428 const unsigned int Normal_RNG::ktab[128] = {
429  0, 12590644, 14272653, 14988939,
430  15384584, 15635009, 15807561, 15933577,
431  16029594, 16105155, 16166147, 16216399,
432  16258508, 16294295, 16325078, 16351831,
433  16375291, 16396026, 16414479, 16431002,
434  16445880, 16459343, 16471578, 16482744,
435  16492970, 16502368, 16511031, 16519039,
436  16526459, 16533352, 16539769, 16545755,
437  16551348, 16556584, 16561493, 16566101,
438  16570433, 16574511, 16578353, 16581977,
439  16585398, 16588629, 16591685, 16594575,
440  16597311, 16599901, 16602354, 16604679,
441  16606881, 16608968, 16610945, 16612818,
442  16614592, 16616272, 16617861, 16619363,
443  16620782, 16622121, 16623383, 16624570,
444  16625685, 16626730, 16627708, 16628619,
445  16629465, 16630248, 16630969, 16631628,
446  16632228, 16632768, 16633248, 16633671,
447  16634034, 16634340, 16634586, 16634774,
448  16634903, 16634972, 16634980, 16634926,
449  16634810, 16634628, 16634381, 16634066,
450  16633680, 16633222, 16632688, 16632075,
451  16631380, 16630598, 16629726, 16628757,
452  16627686, 16626507, 16625212, 16623794,
453  16622243, 16620548, 16618698, 16616679,
454  16614476, 16612071, 16609444, 16606571,
455  16603425, 16599973, 16596178, 16591995,
456  16587369, 16582237, 16576520, 16570120,
457  16562917, 16554758, 16545450, 16534739,
458  16522287, 16507638, 16490152, 16468907,
459  16442518, 16408804, 16364095, 16301683,
460  16207738, 16047994, 15704248, 15472926
461 };
462 
463 // (Ziggurat) tabulated values of 2^{-24}*x[i]
464 const double Normal_RNG::wtab[128] = {
465  1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
466  3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
467  3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
468  4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
469  5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
470  5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
471  5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
472  6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
473  6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
474  6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
475  7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
476  7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
477  7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
478  8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
479  8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
480  8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
481  9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
482  9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
483  9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
484  1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
485  1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
486  1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
487  1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
488  1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
489  1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
490  1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
491  1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
492  1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
493  1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
494  1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
495  1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
496  1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
497 };
498 
499 // (Ziggurat) position of right-most step
500 const double Normal_RNG::PARAM_R = 3.44428647676;
501 
502 // Get a Normal distributed (0,1) sample
504 {
505  uint32_t u, sign, i, j;
506  double x, y;
507 
508  while (true) {
509  u = RNG.genrand_uint32();
510  sign = u & 0x80; // 1 bit for the sign
511  i = u & 0x7f; // 7 bits to choose the step
512  j = u >> 8; // 24 bits for the x-value
513 
514  x = j * wtab[i];
515 
516  if (j < ktab[i])
517  break;
518 
519  if (i < 127) {
520  y = ytab[i+1] + (ytab[i] - ytab[i+1]) * RNG.genrand_close_open();
521  }
522  else {
523  x = PARAM_R - std::log(1.0 - RNG.genrand_close_open()) / PARAM_R;
524  y = std::exp(-PARAM_R * (x - 0.5 * PARAM_R)) * RNG.genrand_close_open();
525  }
526 
527  if (y < std::exp(-0.5 * x * x))
528  break;
529  }
530  return sign ? x : -x;
531 }
532 
533 
535 // Laplace_RNG
537 
538 Laplace_RNG::Laplace_RNG(double meanval, double variance)
539 {
540  setup(meanval, variance);
541 }
542 
543 void Laplace_RNG::setup(double meanval, double variance)
544 {
545  mean = meanval;
546  var = variance;
547  sqrt_12var = std::sqrt(variance / 2.0);
548 }
549 
550 void Laplace_RNG::get_setup(double &meanval, double &variance) const
551 {
552  meanval = mean;
553  variance = var;
554 }
555 
556 
557 
559 {
560  vec vv(n);
561 
562  for (int i=0; i < n; i++)
563  vv(i) = sample();
564 
565  return vv;
566 }
567 
568 mat Laplace_RNG::operator()(int h, int w)
569 {
570  mat mm(h, w);
571  int i, j;
572 
573  for (i = 0; i < h; i++)
574  for (j = 0; j < w; j++)
575  mm(i, j) = sample();
576 
577  return mm;
578 }
579 
581 // AR1_Normal_RNG
583 
584 AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho)
585 {
586  setup(meanval, variance, rho);
587 }
588 
589 void AR1_Normal_RNG::setup(double meanval, double variance, double rho)
590 {
591  mean = meanval;
592  var = variance;
593  r = rho;
594  factr = -2.0 * var * (1.0 - rho * rho);
595  mem = 0.0;
596  odd = true;
597 }
598 
599 void AR1_Normal_RNG::get_setup(double &meanval, double &variance,
600  double &rho) const
601 {
602  meanval = mean;
603  variance = var;
604  rho = r;
605 }
606 
608 {
609  vec vv(n);
610 
611  for (int i=0; i < n; i++)
612  vv(i) = sample();
613 
614  return vv;
615 }
616 
617 mat AR1_Normal_RNG::operator()(int h, int w)
618 {
619  mat mm(h, w);
620  int i, j;
621 
622  for (i = 0; i < h; i++)
623  for (j = 0; j < w; j++)
624  mm(i, j) = sample();
625 
626  return mm;
627 }
628 
630 {
631  mem = 0.0;
632 }
633 
635 // Weibull_RNG
637 
638 Weibull_RNG::Weibull_RNG(double lambda, double beta)
639 {
640  setup(lambda, beta);
641 }
642 
643 void Weibull_RNG::setup(double lambda, double beta)
644 {
645  l = lambda;
646  b = beta;
647  mean = tgamma(1.0 + 1.0 / b) / l;
648  var = tgamma(1.0 + 2.0 / b) / (l * l) - mean;
649 }
650 
651 
653 {
654  vec vv(n);
655 
656  for (int i=0; i < n; i++)
657  vv(i) = sample();
658 
659  return vv;
660 }
661 
662 mat Weibull_RNG::operator()(int h, int w)
663 {
664  mat mm(h, w);
665  int i, j;
666 
667  for (i = 0; i < h; i++)
668  for (j = 0; j < w; j++)
669  mm(i, j) = sample();
670 
671  return mm;
672 }
673 
675 // Rayleigh_RNG
677 
679 {
680  setup(sigma);
681 }
682 
684 {
685  vec vv(n);
686 
687  for (int i=0; i < n; i++)
688  vv(i) = sample();
689 
690  return vv;
691 }
692 
693 mat Rayleigh_RNG::operator()(int h, int w)
694 {
695  mat mm(h, w);
696  int i, j;
697 
698  for (i = 0; i < h; i++)
699  for (j = 0; j < w; j++)
700  mm(i, j) = sample();
701 
702  return mm;
703 }
704 
706 // Rice_RNG
708 
709 Rice_RNG::Rice_RNG(double lambda, double beta)
710 {
711  setup(lambda, beta);
712 }
713 
715 {
716  vec vv(n);
717 
718  for (int i=0; i < n; i++)
719  vv(i) = sample();
720 
721  return vv;
722 }
723 
724 mat Rice_RNG::operator()(int h, int w)
725 {
726  mat mm(h, w);
727  int i, j;
728 
729  for (i = 0; i < h; i++)
730  for (j = 0; j < w; j++)
731  mm(i, j) = sample();
732 
733  return mm;
734 }
735 
736 } // namespace itpp
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:13 for IT++ by Doxygen 1.8.1.2