IT++ Logo
random.cpp
Go to the documentation of this file.
1
29#include <itpp/base/random.h>
30#include <itpp/base/itcompat.h>
32
33
34namespace itpp
35{
36
37namespace random_details
38{
39
40//Thread-local context for thread-safe RNGs
41static ActiveDSFMT::Context thread_local_context;
42#pragma omp threadprivate(thread_local_context)
43//Thread-local context initialization flag
44static bool is_thread_local_context_initialized = false;
45#pragma omp threadprivate(is_thread_local_context_initialized)
46
48{
49 return thread_local_context;
50}
51
53{
54 return is_thread_local_context_initialized;
55}
56
58{
59 is_thread_local_context_initialized = true;
60}
61
68static unsigned int hash_time_to_seed(time_t t, clock_t c)
69{
70 static unsigned int differ = 0; // guarantee time-based seeds will change
71
72 unsigned int h1 = 0;
73 unsigned char *p = (unsigned char *) &t;
74 for(size_t i = 0; i < sizeof(t); ++i) {
75 h1 *= std::numeric_limits<unsigned char>::max() + 2U;
76 h1 += p[i];
77 }
78 unsigned int h2 = 0;
79 p = (unsigned char *) &c;
80 for(size_t j = 0; j < sizeof(c); ++j) {
81 h2 *= std::numeric_limits<unsigned char>::max() + 2U;
82 h2 += p[j];
83 }
84 return (h1 + differ++) ^ h2;
85}
86
87
88
89// ----------------------------------------------------------------------
90// ActiveDSFMT (DSFMT_19937_RNG)
91// ----------------------------------------------------------------------
92template <>
93const bool ActiveDSFMT::bigendian = is_bigendian();
94
95#if defined(__SSE2__)
96template <>
97const __m128i ActiveDSFMT::sse2_param_mask = _mm_set_epi32(ActiveDSFMT::MSK32_3, ActiveDSFMT::MSK32_4, ActiveDSFMT::MSK32_1, ActiveDSFMT::MSK32_2);
98#endif // __SSE2__
99}
100
101/*
102 *Global Seed Provider class definition.
103 *
104 * Provides unique seeds for thread-safe generators running in each thread.
105 */
106class GlobalSeedProvider
107{
108 static const unsigned int default_first_seed = 4257U;
109 typedef random_details::ActiveDSFMT DSFMT;
110public:
111 //constructor
112 GlobalSeedProvider(): _dsfmt(_c), _first_seed_given(false) {
113 _dsfmt.init_gen_rand(default_first_seed);
114 }
115 //set new seed
116 void set_seed(unsigned int s) {_dsfmt.init_gen_rand(s);}
117 //get preset seed
118 unsigned int get_seed() {return _c.last_seed;}
119 //reset provider state to previously set one
120 void reset() { if(_first_seed_given) _dsfmt.init_gen_rand(get_seed());}
121 //set previously saved state from ivec
122 void set_state(const ivec& st) {
123 int size = (DSFMT::N + 1) * 4;
124 it_assert(st.size() == size + 1, "GlobalSeedProvider::state(): "
125 "Invalid state initialization vector");
126 uint32_t *psfmt = &_c.status[0].u32[0];
127 for(int i = 0; i < size; ++i) {
128 psfmt[i] = static_cast<uint32_t>(st(i));
129 }
130 _c.idx = st(size);
131 _first_seed_given = true;
132 }
133 //get current provider state
134 ivec get_state() {
135 int size = (DSFMT::N + 1) * 4;
136 uint32_t *psfmt = &_c.status[0].u32[0];
137 ivec state(size + 1); // size + 1 to save idx variable in the same vec
138 for(int i = 0; i < size; ++i) {
139 state(i) = static_cast<int>(psfmt[i]);
140 }
141 state(size) = _c.idx;
142 return state;
143 }
144 //randomize current provider state with system time
145 void randomize() {_dsfmt.init_gen_rand(random_details::hash_time_to_seed(time(0), clock())); _first_seed_given = true;}
146 //generate new seed for random number generators
147 unsigned int generate() {
148 if(_first_seed_given)
149 return _dsfmt.genrand_uint32();
150 else {
151 //return default seed on first request.
152 //it is done in order not to breake the old-style itpp tests.
153 //Some tests rely on the default state equal to 4257U
154 _first_seed_given = true;
155 return default_first_seed;
156 }
157 }
158private:
159 //
160 DSFMT _dsfmt;
161 //
162 DSFMT::Context _c;
163 //
164 bool _first_seed_given;
165};
166
167
168//Global seed provider instance
169GlobalSeedProvider& global_seed_provider()
170{
171 static GlobalSeedProvider global_seed_provider_instance;
172 return global_seed_provider_instance;
173}
174
175
176void GlobalRNG_reset(unsigned int seed)
177{
178 #pragma omp critical
179 {
180 global_seed_provider().set_seed(seed);
181 }
182}
183
184
186{
187 #pragma omp critical
188 {
189 global_seed_provider().reset();
190 }
191}
192
194{
195 unsigned int s;
196 #pragma omp critical
197 {
198 s = global_seed_provider().generate();
199 }
200 return s;
201}
202
203
205{
206 #pragma omp critical
207 {
208 global_seed_provider().randomize();
209 }
210}
211
212
213void GlobalRNG_get_state(ivec &state)
214{
215 #pragma omp critical
216 {
217 state = global_seed_provider().get_state();
218 }
219}
220
221
222void GlobalRNG_set_state(const ivec &state)
223{
224 #pragma omp critical
225 {
226 global_seed_provider().set_state(state);
227 }
228}
229
230void RNG_reset(unsigned int seed)
231{
233 dsfmt.init_gen_rand(seed);
235}
236
237
239{
241
243 //already initialized. Reinit with last set seed;
244 dsfmt.init_gen_rand(random_details::lc_get().last_seed);
245 }
246 else {
247 //query global seed provider for new seed and init with it
250 }
251}
252
253
255{
257 dsfmt.init_gen_rand(random_details::hash_time_to_seed(time(0), clock()));
259}
260
261
262void RNG_get_state(ivec &state)
263{
264 int size = (random_details::ActiveDSFMT::N + 1) * 4;
265 uint32_t *psfmt = &random_details::lc_get().status[0].u32[0];
266 state.set_size(size + 1); // size + 1 to save idx variable in the same vec
267 for(int i = 0; i < size; ++i) {
268 state(i) = static_cast<int>(psfmt[i]);
269 }
271}
272
273
274void RNG_set_state(const ivec &state)
275{
276 int size = (random_details::ActiveDSFMT::N + 1) * 4;
277 it_assert(state.size() == size + 1, "RNG_set_state: "
278 "Invalid state initialization vector");
279 uint32_t *psfmt = &random_details::lc_get().status[0].u32[0];
280 for(int i = 0; i < size; ++i) {
281 psfmt[i] = static_cast<uint32_t>(state(i));
282 }
284}
285
287// I_Uniform_RNG
289
291{
292 setup(min, max);
293}
294
296{
297 if(min <= max) {
298 lo = min;
299 hi = max;
300 }
301 else {
302 lo = max;
303 hi = min;
304 }
305}
306
307void I_Uniform_RNG::get_setup(int &min, int &max) const
308{
309 min = lo;
310 max = hi;
311}
312
314{
315 ivec vv(n);
316
317 for(int i = 0; i < n; i++)
318 vv(i) = sample();
319
320 return vv;
321}
322
324{
325 imat mm(h, w);
326 int i, j;
327
328 for(i = 0; i < h; i++)
329 for(j = 0; j < w; j++)
330 mm(i, j) = sample();
331
332 return mm;
333}
334
336// Uniform_RNG
338
340{
341 setup(min, max);
342}
343
344void Uniform_RNG::setup(double min, double max)
345{
346 if(min <= max) {
347 lo_bound = min;
348 hi_bound = max;
349 }
350 else {
351 lo_bound = max;
352 hi_bound = min;
353 }
354}
355
356void Uniform_RNG::get_setup(double &min, double &max) const
357{
358 min = lo_bound;
359 max = hi_bound;
360}
361
363// Exp_RNG
365
367{
368 setup(lambda);
369}
370
372{
373 vec vv(n);
374
375 for(int i = 0; i < n; i++)
376 vv(i) = sample();
377
378 return vv;
379}
380
382{
383 mat mm(h, w);
384 int i, j;
385
386 for(i = 0; i < h; i++)
387 for(j = 0; j < w; j++)
388 mm(i, j) = sample();
389
390 return mm;
391}
392
394// Gamma_RNG
396void Gamma_RNG::init_state()
397{
398 const static double sqrt32 = 5.656854;
399
400 const static double q1 = 0.04166669;
401 const static double q2 = 0.02083148;
402 const static double q3 = 0.00801191;
403 const static double q4 = 0.00144121;
404 const static double q5 = -7.388e-5;
405 const static double q6 = 2.4511e-4;
406 const static double q7 = 2.424e-4;
407
408 double r = 1.0 / alpha;
409 _scale = 1.0 / beta;
410
411 it_error_if(!std::isfinite(alpha) || !std::isfinite(_scale) || (alpha < 0.0)
412 || (_scale <= 0.0), "Gamma_RNG::init_state() - wrong parameters");
413
414 _s2 = alpha - 0.5;
415 _s = std::sqrt(_s2);
416 _d = sqrt32 - _s * 12.0;
417
418 _q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
419 + q2) * r + q1) * r;
420
421 /* Approximation depending on size of parameter alpha */
422 /* The constants in the expressions for _b, _si and _c */
423 /* were established by numerical experiments */
424 if(alpha <= 3.686) {
425 _b = 0.463 + _s + 0.178 * _s2;
426 _si = 1.235;
427 _c = 0.195 / _s - 0.079 + 0.16 * _s;
428 }
429 else if(alpha <= 13.022) {
430 _b = 1.654 + 0.0076 * _s2;
431 _si = 1.68 / _s + 0.275;
432 _c = 0.062 / _s + 0.024;
433 }
434 else {
435 _b = 1.77;
436 _si = 0.75;
437 _c = 0.1515 / _s;
438 }
439
440}
442{
443 vec vv(n);
444 for(int i = 0; i < n; i++)
445 vv(i) = sample();
446 return vv;
447}
448
449mat Gamma_RNG::operator()(int r, int c)
450{
451 mat mm(r, c);
452 for(int i = 0; i < r * c; i++)
453 mm(i) = sample();
454 return mm;
455}
456
458{
459 // A copy of rgamma code from the R package, adapted to IT++ by Vasek
460 // Smidl
461
462 /* Constants : */
463 const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
464
465 const static double a1 = 0.3333333;
466 const static double a2 = -0.250003;
467 const static double a3 = 0.2000062;
468 const static double a4 = -0.1662921;
469 const static double a5 = 0.1423657;
470 const static double a6 = -0.1367177;
471 const static double a7 = 0.1233795;
472
473 double e, p, q, t, u, v, w, x, ret_val;
474 double a = alpha;
475 double scale = _scale;
476
477 if(a < 1.) { /* GS algorithm for parameters a < 1 */
478 if(a == 0)
479 return 0.;
480 e = 1.0 + exp_m1 * a;
481 for(;;) { //VS repeat
482 p = e * RNG.genrand_open_open();
483 if(p >= 1.0) {
484 x = -std::log((e - p) / a);
485 if(-std::log(RNG.genrand_open_close()) >= (1.0 - a) * std::log(x))
486 break;
487 }
488 else {
489 x = std::exp(std::log(p) / a);
490 if(-std::log(RNG.genrand_open_close()) >= x)
491 break;
492 }
493 }
494 return scale * x;
495 }
496
497 /* --- a >= 1 : GD algorithm --- */
498
499 /* Step 1: t = standard normal deviate, x = (s,1/2) -normal deviate. */
500 /* immediate acceptance (i) */
501 t = NRNG.sample();
502 x = _s + 0.5 * t;
503 ret_val = x * x;
504 if(t >= 0.0)
505 return scale * ret_val;
506
507 /* Step 2: u = 0,1 - uniform sample. squeeze acceptance (s) */
508 u = RNG.genrand_close_open();
509 if((_d * u) <= (t * t * t))
510 return scale * ret_val;
511
512
513 /* Step 3: no quotient test if x not positive */
514 if(x > 0.0) {
515 /* Step 4: calculation of v and quotient q */
516 v = t / (_s + _s);
517 if(std::fabs(v) <= 0.25)
518 q = _q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
519 + a3) * v + a2) * v + a1) * v;
520 else
521 q = _q0 - _s * t + 0.25 * t * t + (_s2 + _s2) * log(1.0 + v);
522
523 /* Step 5: quotient acceptance (q) */
524 if(log(1.0 - u) <= q)
525 return scale * ret_val;
526 }
527
528 for(;;) { //VS repeat
529 /* Step 6: e = standard exponential deviate
530 * u = 0,1 -uniform deviate
531 * t = (b,si)-double exponential (laplace) sample */
532 e = -std::log(RNG.genrand_open_close()); //see Exponential_RNG
533 u = RNG.genrand_open_close();
534 u = u + u - 1.0;
535 if(u < 0.0)
536 t = _b - _si * e;
537 else
538 t = _b + _si * e;
539 /* Step 7: rejection if t < tau(1) = -0.71874483771719 */
540 if(t >= -0.71874483771719) {
541 /* Step 8: calculation of v and quotient q */
542 v = t / (_s + _s);
543 if(std::fabs(v) <= 0.25)
544 q = _q0 + 0.5 * t * t *
545 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
546 + a2) * v + a1) * v;
547 else
548 q = _q0 - _s * t + 0.25 * t * t + (_s2 + _s2) * log(1.0 + v);
549 /* Step 9: hat acceptance (h) */
550 /* (if q not positive go to step 6) */
551 if(q > 0.0) {
552 // Try to use w = expm1(q); (Not supported on w32)
553 w = expm1(q);
554 /* ^^^^^ original code had approximation with rel.err < 2e-7 */
555 /* if t is rejected sample again at step 6 */
556 if((_c * std::fabs(u)) <= (w * std::exp(e - 0.5 * t * t)))
557 break;
558 }
559 }
560 } /* repeat .. until `t' is accepted */
561 x = _s + 0.5 * t;
562 return scale * x * x;
563}
564
566// Normal_RNG
568
569void Normal_RNG::get_setup(double &meanval, double &variance) const
570{
571 meanval = mean;
572 variance = sigma * sigma;
573}
574
575// (Ziggurat) tabulated values for the heigt of the Ziggurat levels
576const double Normal_RNG::ytab[128] = {
577 1, 0.963598623011, 0.936280813353, 0.913041104253,
578 0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
579 0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
580 0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
581 0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
582 0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
583 0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
584 0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
585 0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
586 0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
587 0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
588 0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
589 0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
590 0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
591 0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
592 0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
593 0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
594 0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
595 0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
596 0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
597 0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
598 0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
599 0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
600 0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
601 0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
602 0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
603 0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
604 0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
605 0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
606 0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
607 0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
608 0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
609};
610
611/*
612 * (Ziggurat) tabulated values for 2^24 times x[i]/x[i+1], used to accept
613 * for U*x[i+1]<=x[i] without any floating point operations
614 */
615const unsigned int Normal_RNG::ktab[128] = {
616 0, 12590644, 14272653, 14988939,
617 15384584, 15635009, 15807561, 15933577,
618 16029594, 16105155, 16166147, 16216399,
619 16258508, 16294295, 16325078, 16351831,
620 16375291, 16396026, 16414479, 16431002,
621 16445880, 16459343, 16471578, 16482744,
622 16492970, 16502368, 16511031, 16519039,
623 16526459, 16533352, 16539769, 16545755,
624 16551348, 16556584, 16561493, 16566101,
625 16570433, 16574511, 16578353, 16581977,
626 16585398, 16588629, 16591685, 16594575,
627 16597311, 16599901, 16602354, 16604679,
628 16606881, 16608968, 16610945, 16612818,
629 16614592, 16616272, 16617861, 16619363,
630 16620782, 16622121, 16623383, 16624570,
631 16625685, 16626730, 16627708, 16628619,
632 16629465, 16630248, 16630969, 16631628,
633 16632228, 16632768, 16633248, 16633671,
634 16634034, 16634340, 16634586, 16634774,
635 16634903, 16634972, 16634980, 16634926,
636 16634810, 16634628, 16634381, 16634066,
637 16633680, 16633222, 16632688, 16632075,
638 16631380, 16630598, 16629726, 16628757,
639 16627686, 16626507, 16625212, 16623794,
640 16622243, 16620548, 16618698, 16616679,
641 16614476, 16612071, 16609444, 16606571,
642 16603425, 16599973, 16596178, 16591995,
643 16587369, 16582237, 16576520, 16570120,
644 16562917, 16554758, 16545450, 16534739,
645 16522287, 16507638, 16490152, 16468907,
646 16442518, 16408804, 16364095, 16301683,
647 16207738, 16047994, 15704248, 15472926
648};
649
650// (Ziggurat) tabulated values of 2^{-24}*x[i]
651const double Normal_RNG::wtab[128] = {
652 1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
653 3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
654 3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
655 4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
656 5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
657 5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
658 5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
659 6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
660 6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
661 6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
662 7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
663 7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
664 7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
665 8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
666 8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
667 8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
668 9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
669 9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
670 9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
671 1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
672 1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
673 1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
674 1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
675 1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
676 1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
677 1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
678 1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
679 1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
680 1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
681 1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
682 1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
683 1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
684};
685
686// (Ziggurat) position of right-most step
687const double Normal_RNG::PARAM_R = 3.44428647676;
688
689// Get a Normal distributed (0,1) sample
691{
692 uint32_t u, sign, i, j;
693 double x, y;
694
695 while(true) {
696 u = RNG.genrand_uint32();
697 sign = u & 0x80; // 1 bit for the sign
698 i = u & 0x7f; // 7 bits to choose the step
699 j = u >> 8; // 24 bits for the x-value
700
701 x = j * wtab[i];
702
703 if(j < ktab[i])
704 break;
705
706 if(i < 127) {
707 y = ytab[i + 1] + (ytab[i] - ytab[i + 1]) * RNG.genrand_close_open();
708 }
709 else {
710 x = PARAM_R - std::log(1.0 - RNG.genrand_close_open()) / PARAM_R;
711 y = std::exp(-PARAM_R * (x - 0.5 * PARAM_R)) * RNG.genrand_close_open();
712 }
713
714 if(y < std::exp(-0.5 * x * x))
715 break;
716 }
717 return sign ? x : -x;
718}
719
720
722// Laplace_RNG
724
725Laplace_RNG::Laplace_RNG(double meanval, double variance)
726{
727 setup(meanval, variance);
728}
729
730void Laplace_RNG::setup(double meanval, double variance)
731{
732 mean = meanval;
733 var = variance;
734 sqrt_12var = std::sqrt(variance / 2.0);
735}
736
737void Laplace_RNG::get_setup(double &meanval, double &variance) const
738{
739 meanval = mean;
740 variance = var;
741}
742
743
744
746{
747 vec vv(n);
748
749 for(int i = 0; i < n; i++)
750 vv(i) = sample();
751
752 return vv;
753}
754
755mat Laplace_RNG::operator()(int h, int w)
756{
757 mat mm(h, w);
758 int i, j;
759
760 for(i = 0; i < h; i++)
761 for(j = 0; j < w; j++)
762 mm(i, j) = sample();
763
764 return mm;
765}
766
768// AR1_Normal_RNG
770
771AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho)
772{
773 setup(meanval, variance, rho);
774}
775
776void AR1_Normal_RNG::setup(double meanval, double variance, double rho)
777{
778 mean = meanval;
779 var = variance;
780 r = rho;
781 factr = -2.0 * var * (1.0 - rho * rho);
782 mem = 0.0;
783 odd = true;
784}
785
786void AR1_Normal_RNG::get_setup(double &meanval, double &variance,
787 double &rho) const
788{
789 meanval = mean;
790 variance = var;
791 rho = r;
792}
793
795{
796 vec vv(n);
797
798 for(int i = 0; i < n; i++)
799 vv(i) = sample();
800
801 return vv;
802}
803
805{
806 mat mm(h, w);
807 int i, j;
808
809 for(i = 0; i < h; i++)
810 for(j = 0; j < w; j++)
811 mm(i, j) = sample();
812
813 return mm;
814}
815
817{
818 mem = 0.0;
819}
820
822// Weibull_RNG
824
825Weibull_RNG::Weibull_RNG(double lambda, double beta)
826{
827 setup(lambda, beta);
828}
829
830void Weibull_RNG::setup(double lambda, double beta)
831{
832 l = lambda;
833 b = beta;
834 mean = tgamma(1.0 + 1.0 / b) / l;
835 var = tgamma(1.0 + 2.0 / b) / (l * l) - mean;
836}
837
838
840{
841 vec vv(n);
842
843 for(int i = 0; i < n; i++)
844 vv(i) = sample();
845
846 return vv;
847}
848
849mat Weibull_RNG::operator()(int h, int w)
850{
851 mat mm(h, w);
852 int i, j;
853
854 for(i = 0; i < h; i++)
855 for(j = 0; j < w; j++)
856 mm(i, j) = sample();
857
858 return mm;
859}
860
862// Rayleigh_RNG
864
866{
867 setup(sigma);
868}
869
871{
872 vec vv(n);
873
874 for(int i = 0; i < n; i++)
875 vv(i) = sample();
876
877 return vv;
878}
879
881{
882 mat mm(h, w);
883 int i, j;
884
885 for(i = 0; i < h; i++)
886 for(j = 0; j < w; j++)
887 mm(i, j) = sample();
888
889 return mm;
890}
891
893// Rice_RNG
895
896Rice_RNG::Rice_RNG(double lambda, double beta)
897{
898 setup(lambda, beta);
899}
900
902{
903 vec vv(n);
904
905 for(int i = 0; i < n; i++)
906 vv(i) = sample();
907
908 return vv;
909}
910
911mat Rice_RNG::operator()(int h, int w)
912{
913 mat mm(h, w);
914 int i, j;
915
916 for(i = 0; i < h; i++)
917 for(j = 0; j < w; j++)
918 mm(i, j) = sample();
919
920 return mm;
921}
922
923} // namespace itpp
double operator()()
Get a single random sample.
Definition: random.h:665
void reset()
Set memory contents to zero.
Definition: random.cpp:816
AR1_Normal_RNG(double meanval=0.0, double variance=1.0, double rho=0.0)
Constructor. Set mean, variance, and correlation.
Definition: random.cpp:771
void get_setup(double &meanval, double &variance, double &rho) const
Get mean, variance and correlation.
Definition: random.cpp:786
void setup(double meanval, double variance, double rho)
Set mean, variance, and correlation.
Definition: random.cpp:776
Exponential_RNG(double lambda=1.0)
constructor. Set lambda.
Definition: random.cpp:366
double operator()()
Get one sample.
Definition: random.h:415
void setup(double lambda)
Set lambda.
Definition: random.h:411
double operator()()
Get one sample.
Definition: random.h:518
double sample()
Get a sample.
Definition: random.cpp:457
void setup(int min, int max)
set min and max values
Definition: random.cpp:295
void get_setup(int &min, int &max) const
get the parameters
Definition: random.cpp:307
I_Uniform_RNG(int min=0, int max=1)
constructor. Sets min and max values.
Definition: random.cpp:290
int sample()
Return a single value from this random generator.
Definition: random.h:338
int operator()()
Get one sample.
Definition: random.h:332
void get_setup(double &meanval, double &variance) const
Get mean and variance.
Definition: random.cpp:737
Laplace_RNG(double meanval=0.0, double variance=1.0)
Constructor. Set mean and variance.
Definition: random.cpp:725
void setup(double meanval, double variance)
Set mean and variance.
Definition: random.cpp:730
double operator()()
Get one sample.
Definition: random.h:555
double sample()
Returns a single sample.
Definition: random.h:561
double sample()
Get a Normal distributed (0,1) sample.
Definition: random.cpp:690
void get_setup(double &meanval, double &variance) const
Get mean and variance.
Definition: random.cpp:569
double genrand_open_open()
Generate uniform (0, 1) double pseudorandom number.
Definition: random.h:254
double genrand_close_open()
Generate uniform [0, 1) double pseudorandom number.
Definition: random.h:234
uint32_t genrand_uint32()
Generate uniform [0, UINT_MAX) integer pseudorandom number.
Definition: random.h:213
double genrand_open_close()
Generate uniform (0, 1] double pseudorandom number.
Definition: random.h:244
Rayleigh_RNG(double sigma=1.0)
Constructor. Set sigma.
Definition: random.cpp:865
double operator()()
Get one sample.
Definition: random.h:743
void setup(double sigma)
Set sigma.
Definition: random.h:739
double operator()()
Get one sample.
Definition: random.h:773
Rice_RNG(double sigma=1.0, double v=1.0)
Constructor. Set sigma, and v (if v = 0, Rice -> Rayleigh).
Definition: random.cpp:896
void setup(double sigma, double v)
Set sigma, and v (if v = 0, Rice -> Rayleigh).
Definition: random.h:769
Uniform_RNG(double min=0, double max=1.0)
Constructor. Set min, max and seed.
Definition: random.cpp:339
void setup(double min, double max)
set min and max
Definition: random.cpp:344
void get_setup(double &min, double &max) const
get parameters
Definition: random.cpp:356
double operator()()
Get one sample.
Definition: random.h:715
Weibull_RNG(double lambda=1.0, double beta=1.0)
Constructor. Set lambda and beta.
Definition: random.cpp:825
void setup(double lambda, double beta)
Set lambda, and beta.
Definition: random.cpp:830
C++ implementation of dSFMT random number generator.
Definition: random_dsfmt.h:99
void init_gen_rand(unsigned int seed)
Initialise the generator with a new seed.
Definition: random_dsfmt.h:154
Elementary mathematical functions - header file.
#define it_error_if(t, s)
Abort if t is true.
Definition: itassert.h:117
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94
vec log(const vec &x)
The natural logarithm of the elements.
Definition: log_exp.h:241
int size(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:55
bool is_bigendian()
Returns true if machine endianness is BIG_ENDIAN.
Definition: misc.cpp:50
double sign(double x)
Signum function.
Definition: elem_math.h:81
T min(const Vec< T > &in)
Minimum value of vector.
Definition: min_max.h:125
T max(const Vec< T > &v)
Maximum value of vector.
Definition: min_max.h:45
void RNG_randomize()
Set a random seed for all Random Number Generators in the current thread.
Definition: random.cpp:254
void lc_mark_initialized()
Function to mark thread-local context as initialized.
Definition: random.cpp:57
void GlobalRNG_get_state(ivec &state)
Save current full state of global seed provider in memory.
Definition: random.cpp:213
void RNG_set_state(const ivec &state)
Resume Random Number generation in the current thread with previously stored context.
Definition: random.cpp:274
unsigned int GlobalRNG_get_local_seed()
Get new seed to initialize thread-local generators.
Definition: random.cpp:193
void RNG_reset()
Reset the seed to the previously set value for all Random Number Generators in the current thread.
Definition: random.cpp:238
bool lc_is_initialized()
Function to check if thread-local context is initialized.
Definition: random.cpp:52
void GlobalRNG_reset()
Reset the internal seed of the Global Seed Provider to the previously set value.
Definition: random.cpp:185
void RNG_get_state(ivec &state)
Save Random Number generation context used in the current thread.
Definition: random.cpp:262
DSFMT_19937_RNG ActiveDSFMT
Active Generator for random (stochastic) sources.
Definition: random_dsfmt.h:395
ActiveDSFMT::Context & lc_get()
Function to access thread-local context for random numbers generation.
Definition: random.cpp:47
void GlobalRNG_randomize()
Set a random seed for the Global Seed Provider seed.
Definition: random.cpp:204
void GlobalRNG_set_state(const ivec &state)
Resume the global seed provider state saved in memory.
Definition: random.cpp:222
double variance(const cvec &v)
The variance of the elements in the vector. Normalized with N-1 to be unbiased.
Definition: misc_stat.cpp:159
IT++ compatibility types and functions.
itpp namespace
Definition: itmex.h:37
static unsigned int hash_time_to_seed(time_t t, clock_t c)
Get an unsigned int from time variables t and c.
Definition: random.cpp:68
Definition of classes for random number generators.
DSFMT context structure.
Definition: random_dsfmt.h:122
int idx
State array indexing.
Definition: random_dsfmt.h:138
w128_t status[N+1]
128-bit internal state array
Definition: random_dsfmt.h:136
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.6