Nemo  2.3.56
Simulate forward-in-time genetic evolution in a spatially explicit, individual-based stochastic simulator
RAND Class Reference

Random number generation class, uses various types of random generators depending on the implementation. More...

#include <Uniform.h>

+ Collaboration diagram for RAND:

Static Public Member Functions

static void init (unsigned long seed)
 Initialize the random generator's seed. More...
 
static void free ()
 Memory de-allocation. More...
 
static double Uniform ()
 Generates a random number from [0.0, 1.0[ uniformly distributed. More...
 
static unsigned int Uniform (unsigned int max)
 Returns a uniformly distributed random number from [0.0, max[. More...
 
static bool RandBool ()
 Returns a random boolean. More...
 
static unsigned long RandULong ()
 Return a random unsigned long, from uniform distribution. More...
 
static double gammln (double xx)
 From the Numerical Recieps. More...
 
static double Poisson (double mean)
 From the Numerical Recieps. More...
 
static double Gaussian (double sigma)
 
static void BivariateGaussian (double sigma1, double sigma2, double rho, double *out1, double *out2)
 
static double LogNormal (double zeta, double sigma)
 
static double Gamma (double a, double b)
 
static double Bernoulli (double p)
 
static double Exponential (double mu)
 
static double Binomial (double p, unsigned int n)
 
static double Beta (const double a, const double b)
 
static unsigned int Binomial2 (double p, unsigned int n)
 
static void Multinomial (size_t K, unsigned int N, const double p[], unsigned int n[])
 
static void MultinomialOnNormalizedValarray (size_t K, unsigned int N, const std::valarray< double > &p, unsigned int n[])
 Multinomial draw assuming the probabilities sum to 1.0 and are all > 0. More...
 
static void MultinomialOnNormalizedValarray_expandedOut (size_t K, unsigned int N, const std::valarray< double > &p, unsigned int n[])
 Multinomial draw assuming the probabilities sum to 1.0 and are all > 0, the output is an array of size N. More...
 
static void MultinomialOnNormalizedValarray_scrambleOut (size_t K, unsigned int N, const std::valarray< double > &p, unsigned int n[])
 Multinomial draw assuming the probabilities sum to 1.0 and are all > 0, the output is an array of size N. More...
 
static void MultinomialOnNormalizedValarrayZipper_scrambleOut (size_t K, unsigned int N, const std::valarray< double > &p, unsigned int n[])
 Multinomial draw assuming the probabilities sum to 1.0 and are all > 0, the output is an array of size N. More...
 
static void ScrambleArrayUInt (const int length, unsigned int *array)
 Randomize the elements within an array. More...
 
static void Sample (const int from, const int to, const unsigned int num, int *result, bool replace)
 Creates a sample of integers within range [from, to), with or without replacement. More...
 
static void SampleSeq (int from, int to, int by, unsigned int num, int *result, bool replace=false)
 
static void SampleSeqWithReciprocal (int from, int to, int by, unsigned int num1, int *result1, unsigned int num2, int *result2)
 

Static Public Attributes

static long Seed1 = 0
 
static long Seed2 = 98280582
 

Private Member Functions

 RAND ()
 

Detailed Description

Random number generation class, uses various types of random generators depending on the implementation.

Constructor & Destructor Documentation

◆ RAND()

RAND::RAND ( )
private

Member Function Documentation

◆ Bernoulli()

static double RAND::Bernoulli ( double  p)
inlinestatic
431 {
432#if defined(HAS_GSL) && !defined(HAS_SPRNG)
433 return gsl_ran_bernoulli(r, p);
434#else
435 double u = RAND::Uniform() ;
436
437 if (u < p)
438 {
439 return 1 ;
440 }
441 else
442 {
443 return 0 ;
444 }
445#endif
446 }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:125

References Uniform().

Referenced by TTNeutralGenes::init_sequence(), and ParamsParser::rbernoul().

+ Here is the caller graph for this function:

◆ Beta()

static double RAND::Beta ( const double  a,
const double  b 
)
inlinestatic
528 {
529 /*from Knuth*/
530 double x1 = RAND::Gamma(a, 1.0);
531 double x2 = RAND::Gamma(b, 1.0);
532
533 return x1 / (x1 + x2);
534 }
static double Gamma(double a, double b)
Definition: Uniform.h:386

References Gamma().

Referenced by Binomial().

+ Here is the caller graph for this function:

◆ Binomial()

static double RAND::Binomial ( double  p,
unsigned int  n 
)
inlinestatic
491 {
492 /*implements Knuth method*/
493 unsigned int i, a, b, k = 0;
494
495 while (n > 10) /* This parameter is tunable */
496 {
497 double X;
498 a = 1 + (n / 2);
499 b = 1 + n - a;
500
501 X = RAND::Beta((double) a, (double) b);
502
503 if (X >= p)
504 {
505 n = a - 1;
506 p /= X;
507 }
508 else
509 {
510 k += a;
511 n = b - 1;
512 p = (p - X) / (1 - X);
513 }
514 }
515
516 for (i = 0; i < n; i++)
517 {
518 double u = RAND::Uniform();
519 if (u < p)
520 k++;
521 }
522
523 return k;
524
525 }
static double Beta(const double a, const double b)
Definition: Uniform.h:527

References Beta(), and Uniform().

Referenced by Binomial2(), Multinomial(), MultinomialOnNormalizedValarray(), MultinomialOnNormalizedValarray_expandedOut(), MultinomialOnNormalizedValarray_scrambleOut(), TTNeutralGenes::mutate_2all(), TTNeutralGenes::mutate_KAM(), and TTNeutralGenes::mutate_SSM().

+ Here is the caller graph for this function:

◆ Binomial2()

static unsigned int RAND::Binomial2 ( double  p,
unsigned int  n 
)
inlinestatic
536 {
537
538#if defined(HAS_GSL) && !defined(HAS_SPRNG)
539 return gsl_ran_binomial (r, p, n);
540#else
541 return Binomial(p, n);
542#endif
543 }
static double Binomial(double p, unsigned int n)
Definition: Uniform.h:490

References Binomial().

◆ BivariateGaussian()

static void RAND::BivariateGaussian ( double  sigma1,
double  sigma2,
double  rho,
double *  out1,
double *  out2 
)
inlinestatic
331 {
332#if defined(HAS_GSL) && !defined(HAS_SPRNG)
333 gsl_ran_bivariate_gaussian(r,sigma1,sigma2,rho,out1,out2);
334#else
335 //gsl code:
336 double u, v, r2, scale;
337 //double *x = out, *y = (++out);
338 do
339 {
340 /* choose x,y in uniform square (-1,-1) to (+1,+1) */
341
342 u = -1 + 2 * Uniform ();
343 v = -1 + 2 * Uniform ();
344
345 /* see if it is in the unit circle */
346 r2 = u * u + v * v;
347 }
348 while (r2 > 1.0 || r2 == 0);
349
350 scale = sqrt (-2.0 * log (r2) / r2);
351
352 *out1 = sigma1 * u * scale;
353 *out2 = sigma2 * (rho * u + sqrt(1 - rho*rho) * v) * scale;
354
355#endif
356 }

References Uniform().

Referenced by TProtoQuanti::getMutationEffectBivariateGaussian().

+ Here is the caller graph for this function:

◆ Exponential()

static double RAND::Exponential ( double  mu)
inlinestatic
448 {
449 return -mu * log(RAND::Uniform());
450 }

References Uniform().

Referenced by ParamsParser::rexp(), and TProtoDeletMutations_bitstring::set_effects_exp().

+ Here is the caller graph for this function:

◆ free()

static void RAND::free ( )
inlinestatic

Memory de-allocation.

110 {
111
112#if defined(HAS_SPRNG)
113 // do nothing with pointer?
114#elif defined(HAS_GSL)
115 gsl_rng_free(r);
116#endif
117
118 }

Referenced by SimRunner::run().

+ Here is the caller graph for this function:

◆ Gamma()

static double RAND::Gamma ( double  a,
double  b 
)
inlinestatic
386 {
387#if defined(HAS_GSL) && !defined(HAS_SPRNG)
388 return gsl_ran_gamma(r, a, b);
389#else
390 /*pasting GSL code in here*/
391
392 /* assume a > 0 */
393
394 if (a < 1)
395 {
396 double u = Uniform();//might give singularity at 0....
397 return Gamma(1.0 + a, b) * pow (u, 1.0 / a);
398 }
399
400 {
401 double x, v, u;
402 double d = a - 1.0 / 3.0;
403 double c = (1.0 / 3.0) / sqrt (d);
404
405 while (1)
406 {
407 do
408 {
409 x = Gaussian(1.0); //should use the Ziggurat method?
410 v = 1.0 + c * x;
411 }
412 while (v <= 0);
413
414 v = v * v * v;
415 u = Uniform();//might give singularity at 0....
416
417 if (u < 1 - 0.0331 * x * x * x * x)
418 break;
419
420 if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
421 break;
422 }
423
424 return b * d * v;
425 }
426
427#endif
428 }
static double Gaussian(double sigma)
Definition: Uniform.h:262

References Gamma(), Gaussian(), and Uniform().

Referenced by Beta(), Gamma(), ParamsParser::rgamma(), and TProtoDeletMutations_bitstring::set_effects_gamma().

+ Here is the caller graph for this function:

◆ gammln()

static double RAND::gammln ( double  xx)
inlinestatic

From the Numerical Recieps.

205 {
206 double x,y,tmp,ser=1.000000000190015;
207 static double cof[6]={76.18009172947146,-86.50532032941677,
208 24.01409824083091,-1.231739572450155,
209 0.1208650973866179e-2,-0.5395239384953e-5};
210 int j;
211 y=x=xx;
212 tmp=x+5.5;
213 tmp -= (x+0.5)*log(tmp);
214 for (j = 0; j < 6; ++j) ser += cof[j]/++y;
215
216 return -tmp+log(2.5066282746310005*ser/x);
217
218 }

Referenced by Poisson().

+ Here is the caller graph for this function:

◆ Gaussian()

static double RAND::Gaussian ( double  sigma)
inlinestatic

From the GSL.

263 {
264#if defined(HAS_GSL) && !defined(HAS_SPRNG)
265
266 return gsl_ran_gaussian_ratio_method (r, sigma);
267
268#else
270// double x, y, r2;
271//
272// do
273// {
274// /* choose x,y in uniform square (-1,-1) to (+1,+1) */
275//
276// x = -1 + 2 * Uniform ( );
277// y = -1 + 2 * Uniform ( );
278//
279// /* see if it is in the unit circle */
280// r2 = x * x + y * y;
281// }
282// while (r2 > 1.0 || r2 == 0);
283//
284// /* Box-Muller transform */
285// return sigma * y * sqrt (-2.0 * log (r2) / r2);
286
287 /*trying the ratio method from GSL*/
288 /* see code in gsl/randist/gauss.c */
289 double u, v, x, y, Q;
290 const double s = 0.449871; /* Constants from Leva */
291 const double t = -0.386595;
292 const double a = 0.19600;
293 const double b = 0.25472;
294 const double r1 = 0.27597;
295 const double r2 = 0.27846;
296
297 do /* This loop is executed 1.369 times on average */
298 {
299 /* Generate a point P = (u, v) uniform in a rectangle enclosing
300 the K+M region v^2 <= - 4 u^2 log(u). */
301
302 /* u in (0, 1] to avoid singularity at u = 0 */
303 u = 1 - RAND::Uniform();
304
305 /* v is in the asymmetric interval [-0.5, 0.5). However v = -0.5
306 is rejected in the last part of the while clause. The
307 resulting normal deviate is strictly symmetric about 0
308 (provided that v is symmetric once v = -0.5 is excluded). */
309 v = RAND::Uniform() - 0.5;
310
311 /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
312 much (for efficiency). */
313 v *= 1.7156;
314
315 /* Compute Leva's quadratic form Q */
316 x = u - s;
317 y = fabs (v) - t;
318 Q = x * x + y * (a * y - b * x);
319
320 /* Accept P if Q < r1 (Leva) */
321 /* Reject P if Q > r2 (Leva) */
322 /* Accept if v^2 <= -4 u^2 log(u) (K+M) */
323 /* This final test is executed 0.012 times on average. */
324 }
325 while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));
326
327 return sigma * (v / u);
328#endif
329 }

References Uniform().

Referenced by Gamma(), LCE_Selection_base::getFitnessMultivariateGaussian_VE(), LCE_Selection_base::getFitnessUnivariateGaussian_VE(), LCE_Breed_base::getGaussianFecundity(), TProtoQuanti::getMutationEffectUnivariateGaussian(), TTDispersal::init_sequence(), TTQuanti::init_sequence(), LCE_Patch_Extinction::rand_gaussian(), ParamsParser::rnorm(), TProtoQuanti::set_trait_value_VE(), and TProtoQuanti::setDominanceParameters().

+ Here is the caller graph for this function:

◆ init()

static void RAND::init ( unsigned long  seed)
inlinestatic

Initialize the random generator's seed.

83 {
84
85#if defined(HAS_SPRNG)
86 // initializing the SPRNG stream with default param and LFG (0) generator
87
88 stream = SelectType(4); //MLFG Modified Lagged Fibonacci
89
90 stream->init_sprng( _myenv->workerRank(), _myenv->workerCount()+1, seed, SPRNG_DEFAULT );
91
92// message("--- initialized the SPRNG random generator on rank %i\n", _myenv->workerRank());
93//
94// stream->print_sprng();
95//
96// message("--- process %i random number: %.14f\n", _myenv->workerRank(), stream->sprng());
97
98#elif defined(HAS_GSL)
99
100 init_gsl( gsl_rng_mt19937, seed );
101
102#else
103
104 Seed1 = seed;
105
106#endif
107 }
MPIenv * _myenv
Definition: MPImanager.cc:36
unsigned int workerCount() const
Definition: MPImanager.h:128
unsigned int workerRank() const
Definition: MPImanager.h:129
static long Seed1
Definition: Uniform.h:79

References _myenv, Seed1, MPIenv::workerCount(), and MPIenv::workerRank().

Referenced by SimRunner::init_random_seed(), and main().

+ Here is the caller graph for this function:

◆ LogNormal()

static double RAND::LogNormal ( double  zeta,
double  sigma 
)
inlinestatic
358 {
359#if defined(HAS_GSL) && !defined(HAS_SPRNG)
360 return gsl_ran_lognormal(r,zeta,sigma);
361#else
362 //this is the GSL code:
363 double u, v, r2, normal, z;
364
365 do
366 {
367 /* choose x,y in uniform square (-1,-1) to (+1,+1) */
368
369 u = -1 + 2 * Uniform();
370 v = -1 + 2 * Uniform();
371
372 /* see if it is in the unit circle */
373 r2 = u * u + v * v;
374 }
375 while (r2 > 1.0 || r2 == 0);
376
377 normal = u * sqrt (-2.0 * log (r2) / r2);
378
379 z = exp (sigma * normal + zeta);
380
381 return z;
382
383#endif
384 }

References Uniform().

Referenced by LCE_Breed_base::getLogNormalFecundity(), LCE_Patch_Extinction::rand_lognormal(), ParamsParser::rlognorm(), and TProtoDeletMutations_bitstring::set_effects_lognorm().

+ Here is the caller graph for this function:

◆ Multinomial()

static void RAND::Multinomial ( size_t  K,
unsigned int  N,
const double  p[],
unsigned int  n[] 
)
inlinestatic
546 {
547#if defined(HAS_GSL) && !defined(HAS_SPRNG)
548 return gsl_ran_multinomial (r, K, N, p, n);
549#else
550 // GSL code here:
551 size_t k;
552 double norm = 0.0;
553 double sum_p = 0.0;
554
555 unsigned int sum_n = 0;
556
557 /* p[k] may contain non-negative weights that do not sum to 1.0.
558 * Even a probability distribution will not exactly sum to 1.0
559 * due to rounding errors.
560 */
561
562 for (k = 0; k < K; k++)
563 {
564 norm += p[k];
565 }
566
567 for (k = 0; k < K; k++)
568 {
569 if (p[k] > 0.0)
570 {
571 n[k] = Binomial (p[k] / (norm - sum_p), N - sum_n); //MOD
572 }
573 else
574 {
575 n[k] = 0;
576 }
577
578 sum_p += p[k];
579 sum_n += n[k];
580 }
581
582#endif
583 }

References Binomial().

◆ MultinomialOnNormalizedValarray()

static void RAND::MultinomialOnNormalizedValarray ( size_t  K,
unsigned int  N,
const std::valarray< double > &  p,
unsigned int  n[] 
)
inlinestatic

Multinomial draw assuming the probabilities sum to 1.0 and are all > 0.

587 {
588 // GSL code here: FG 2020: modified for normalized arrays
589 size_t k;
590 double sum_p = 0.0;
591
592 unsigned int sum_n = 0;
593
594 for (k = 0; k < K; k++)
595 {
596 n[k] = RAND::Binomial (p[k] / (1.0 - sum_p), N - sum_n); //MOD FG
597 sum_p += p[k];
598 sum_n += n[k];
599 }
600 }

References Binomial().

◆ MultinomialOnNormalizedValarray_expandedOut()

static void RAND::MultinomialOnNormalizedValarray_expandedOut ( size_t  K,
unsigned int  N,
const std::valarray< double > &  p,
unsigned int  n[] 
)
inlinestatic

Multinomial draw assuming the probabilities sum to 1.0 and are all > 0, the output is an array of size N.

606 {
607 // GSL code here: FG 2020: modified for normalized arrays
608 size_t k;
609 double sum_p = 0.0;
610 unsigned int _n, sum_n = 0;
611 unsigned int i = 0, pos = 0;
612
613 for (k = 0; k < K; k++) {
614 //draw number of events k with probability p[k]
615 _n = RAND::Binomial (p[k] / (1.0 - sum_p), N - sum_n);
616
617 for(i = 0; i < _n; i++)
618 n[pos++] = k;
619
620 sum_p += p[k];
621 sum_n += _n;
622 assert(pos <= N);
623 }
624 assert(sum_n == N);
625 }

References Binomial().

◆ MultinomialOnNormalizedValarray_scrambleOut()

static void RAND::MultinomialOnNormalizedValarray_scrambleOut ( size_t  K,
unsigned int  N,
const std::valarray< double > &  p,
unsigned int  n[] 
)
inlinestatic

Multinomial draw assuming the probabilities sum to 1.0 and are all > 0, the output is an array of size N.

631 {
632 // GSL code here: FG 2020: modified for normalized arrays
633 size_t k;
634 double sum_p = 0.0;
635 unsigned int _n, sum_n = 0;
636 unsigned int i = 0, pos = 0;
637
638 for (k = 0; k < K; k++) {
639 //draw number of events k with probability p[k]
640 _n = RAND::Binomial (p[k] / (1.0 - sum_p), N - sum_n);
641
642 for(i = 0; i < _n; i++)
643 n[pos++] = k;
644
645 sum_p += p[k];
646 sum_n += _n;
647 assert(pos <= N);
648 }
649 assert(sum_n == N);
650
651 //scramble the output array
653
654 }
static void ScrambleArrayUInt(const int length, unsigned int *array)
Randomize the elements within an array.
Definition: Uniform.h:690

References Binomial(), and ScrambleArrayUInt().

Referenced by LCE_Breed_Selection::do_breed_selection_WrightFisher_1sex(), and LCE_Breed_Selection::do_breed_selection_WrightFisher_2sex().

+ Here is the caller graph for this function:

◆ MultinomialOnNormalizedValarrayZipper_scrambleOut()

static void RAND::MultinomialOnNormalizedValarrayZipper_scrambleOut ( size_t  K,
unsigned int  N,
const std::valarray< double > &  p,
unsigned int  n[] 
)
inlinestatic

Multinomial draw assuming the probabilities sum to 1.0 and are all > 0, the output is an array of size N.

660 {
661 size_t k;
662 double sum = 0.0;
663 double rand;
664 unsigned int i = 0;
665
666 for (k = 0; k < N; k++) {
667
668 rand = std::min(0.99999995, RAND::Uniform()); // to avoid problems with rounding errors
669 i = 0;
670 sum = p[i];
671
672 //zip through the proba array:
673 while( sum < rand ) {
674 i++;
675 sum += p[i];
676 }
677
678 n[k] = i;
679
680 }
681 //scramble the output array
683
684 }

References ScrambleArrayUInt(), and Uniform().

◆ Poisson()

static double RAND::Poisson ( double  mean)
inlinestatic

From the Numerical Recieps.

220 {
221
222#if defined(HAS_GSL) && !defined(HAS_SPRNG)
223 return gsl_ran_poisson(r, mean);
224#else
225 static double sq,alxm,g,oldm=(-1.0);
226 double em,t,y;
227
228 if (mean < 12.0)
229 {
230 if (mean != oldm){
231 oldm=mean;
232 g=exp(-mean);
233 }
234 em = -1;
235 t=1.0;
236 do {
237 ++em;
238 t *= Uniform();
239 } while (t > g);
240 }else
241 {
242 if (mean != oldm)
243 {
244 oldm=mean;
245 sq=sqrt(2.0*mean);
246 alxm=log(mean);
247 g=mean*alxm-gammln(mean+1.0);
248 }
249 do {
250 do {
251 y=tan(M_PI*Uniform());
252 em=sq*y+mean;
253 } while (em < 0.0);
254 em=floor(em);
255 t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
256 } while (Uniform() > t);
257 }
258 return em;
259#endif
260 }
static double gammln(double xx)
From the Numerical Recieps.
Definition: Uniform.h:205

References gammln(), and Uniform().

Referenced by LCE_Breed_base::getPoissonFecundity(), TTDispersal::mutate(), TProtoQuanti::mutate_diallelic_HC(), TT_BDMI::mutate_diplo(), TT_BDMI::mutate_haplo(), TProtoQuanti::mutate_HC(), TProtoQuanti::mutate_noHC(), TTDeletMutations_bitstring::mutate_noredraw(), TTDeletMutations_bitstring::mutate_noredraw_noBackMutation(), TTDeletMutations_bitstring::mutate_redraw(), LCE_Patch_Extinction::rand_poisson(), GeneticMap::recombine(), ParamsParser::rpoiss(), and LCE_Breed_Disperse::stochasticLogisticGrowth().

+ Here is the caller graph for this function:

◆ RandBool()

static bool RAND::RandBool ( )
inlinestatic

Returns a random boolean.

163 {
164 //generate a random int (or long)
165#ifdef HAS_SPRNG
166 static int intrand = stream->isprng();
167#elif defined(HAS_GSL)
168 static unsigned long int intrand = gsl_rng_get(r);
169#else
170 static int intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
171#endif
172 //read up to the first 16 bits
173 static unsigned int num = 0;
174 //redraw an number after that limit
175 if ( ++num > 16 ) {
176 num = 0;
177
178#ifdef HAS_SPRNG
179 intrand = stream->isprng();
180#elif defined(HAS_GSL)
181 intrand = gsl_rng_get(r);
182#else
183 intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
184#endif
185
186 }
187 return (intrand & (1 << num) );
188
189 }

References Uniform().

Referenced by LCE_Patch_Extinction::do_remove(), TProtoQuanti::getMutationEffectBivariateDiallelic(), TProtoQuanti::getMutationEffectUnivariateDiallelic(), LCE_Breed_base::getOffsprgSexFixed(), LCE_Breed_base::getOffsprgSexRandom(), TTDispersal::inherit(), TProtoDeletMutations_bitstring::inherit_free(), TProtoQuanti::inherit_free(), TProtoNeutralGenes::inherit_free(), LCE_Disperse_EvolDisp::Migrate_SteppingStone1D(), TTDispersal::mutate(), TTNeutralGenes::mutate_2all(), TProtoQuanti::mutate_diallelic_HC(), TT_BDMI::mutate_diplo(), TProtoQuanti::mutate_HC(), TTNeutralGenes::mutate_KAM(), TProtoQuanti::mutate_noHC(), TTDeletMutations_bitstring::mutate_noredraw(), TTDeletMutations_bitstring::mutate_noredraw_noBackMutation(), TTDeletMutations_bitstring::mutate_redraw(), TTNeutralGenes::mutate_SSM(), GeneticMap::recombine(), LCE_Breed_Wolbachia::wolbachia_model_1(), and LCE_Breed_Wolbachia::wolbachia_model_2().

◆ RandULong()

static unsigned long RAND::RandULong ( )
inlinestatic

Return a random unsigned long, from uniform distribution.

192 {
193#if defined(HAS_GSL) && !defined(HAS_SPRNG)
194 return gsl_rng_get(r);
195#else
196 unsigned long rnd, limit = 0x10000000; //=2^28 this gives ~7% redraws
197 do{
198 rnd = (unsigned long)(Uniform()*ULONG_MAX);
199 }while(rnd < limit);
200
201 return rnd;
202#endif
203 }

References Uniform().

◆ Sample()

static void RAND::Sample ( const int  from,
const int  to,
const unsigned int  num,
int result,
bool  replace 
)
inlinestatic

Creates a sample of integers within range [from, to), with or without replacement.

Can be used to scramble an array (without replacement). The random sequence of integers is placed in the 'results' array.

Parameters
fromstarting number of the series
tolast number of the series
numnumber of elements to draw within [from, to)
resultcontainer to hold the resulting randomized sequence of integers
715 {
716 assert(from < to);
717
718 unsigned int seq_length = to - from; //don't include last (to)
719
720 assert(num <= seq_length);
721
722 // we build a sequence of indexes to choose from
723 int *seq = new int [seq_length];
724
725 seq[0] = from;
726
727 for(unsigned int i = 1; seq[i-1] < to && i < seq_length; ++i)
728 seq[i] = seq[i-1] + 1;
729
730 if(!replace) { //without replacement
731
732 unsigned int size = seq_length, pos, last = seq_length - 1;
733
734 for (unsigned int i = 0; i < num; i++) {
735 pos = RAND::Uniform(size);
736 result[i] = seq[pos];
737 seq[pos] = seq[last];
738 // seq[last] = result[i]; useless operation
739 size--; last--;
740 }
741
742 } else { //with replacement
743 for (unsigned int i = 0; i < num; i++)
744 result[i] = seq[ RAND::Uniform(seq_length) ];
745 }
746
747 delete [] seq;
748 }

References Uniform().

Referenced by MPFileHandler::createAndPrintSample(), and FileServices::subSamplePatch().

+ Here is the caller graph for this function:

◆ SampleSeq()

static void RAND::SampleSeq ( int  from,
int  to,
int  by,
unsigned int  num,
int result,
bool  replace = false 
)
inlinestatic
751 {
752 assert(from < to && by < to && by > 0);
753
754 unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
755
756 assert(num <= seq_length);
757
758 int *seq = new int [seq_length];
759
760 seq[0] = from;
761
762 for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
763
764 if(!replace) { //without replacement
765
766 unsigned int size = seq_length, pos, last = seq_length - 1;
767
768 for (unsigned int i = 0; i < num; i++) {
769 pos = RAND::Uniform(size);
770 result[i] = seq[pos];
771 seq[pos] = seq[last];
772 // seq[last] = result[i]; useless operation
773 size--; last--;
774 }
775
776 } else { //with replacement
777 for (unsigned int i = 0; i < num; i++) result[i] = seq[ RAND::Uniform(seq_length) ];
778 }
779
780 delete [] seq;
781 }

References Uniform().

◆ SampleSeqWithReciprocal()

static void RAND::SampleSeqWithReciprocal ( int  from,
int  to,
int  by,
unsigned int  num1,
int result1,
unsigned int  num2,
int result2 
)
inlinestatic
784 {
785 assert(from < to && by < to && by > 0);
786
787 unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
788
789 assert(num1 + num2 == seq_length);
790
791 int *seq = new int [seq_length];
792
793 seq[0] = from;
794
795 for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
796
797 unsigned int size = seq_length, pos, last = seq_length - 1;
798
799 for (unsigned int i = 0; i < num1; i++) {
800 pos = RAND::Uniform(size);
801 result1[i] = seq[pos];
802 seq[pos] = seq[last];
803 seq[last] = result1[i];
804 size--; last--;
805 }
806
807 for (unsigned int i = 0; i < num2 && i < size; i++)
808 result2[i] = seq[i];
809
810 delete [] seq;
811 }

References Uniform().

Referenced by LCE_NtrlInit::init_allele_freq(), and LCE_QuantiInit::init_allele_freq().

+ Here is the caller graph for this function:

◆ ScrambleArrayUInt()

static void RAND::ScrambleArrayUInt ( const int  length,
unsigned int array 
)
inlinestatic

Randomize the elements within an array.

Parameters
lengththe length of the array
arraythe array to scramble
691 {
692 unsigned int size = length, pos, last = length - 1, num = length -1, el;
693
694 // with stop after 2 elements are left
695 // this algo allows swapping an element with itself, which is fine
696 for (unsigned int i = 0; i < num; i++) {
697 pos = RAND::Uniform(size);
698 el = array[last];
699 array[last] = array[pos];
700 array[pos] = el;
701 size--; last--;
702 }
703
704 assert(size == 1); //we stopped before swapping the last (first) element with itself
705 }

References Uniform().

Referenced by MultinomialOnNormalizedValarray_scrambleOut(), and MultinomialOnNormalizedValarrayZipper_scrambleOut().

+ Here is the caller graph for this function:

◆ Uniform() [1/2]

static double RAND::Uniform ( )
inlinestatic

Generates a random number from [0.0, 1.0[ uniformly distributed.

If SPRNG or GSL libraries are not used, implement a random generator from: L'Ecuyer, 1988, "Efficient and Portable Combined Random Number Generators", Communication of the ACM, 31(6):742-774.

125 {
126
127#ifdef HAS_SPRNG
128 return stream->sprng();
129#elif defined(HAS_GSL)
130 return gsl_rng_uniform(r);
131#else
132 register long z, w;
133
134 do{
135 w = Seed1 / 53668;
136
137 Seed1 = 40014 * (Seed1 - w * 53668) - w * 12211;
138
139 if (Seed1 < 0) Seed1 += 2147483563;
140
141 w = (Seed2 / 52774);
142
143 Seed2 = 40692 * (Seed2 - w * 52774) - w * 3791;
144
145 if (Seed2 < 0) Seed2 += 2147483399;
146
147 z = Seed1 - Seed2;
148
149 if (z < 1) z += 2147483562;
150
151 }while (!((z * 4.656613e-10) < 1.0));
152
153 return (z * 4.656613e-10);
154#endif
155 }
static long Seed2
Definition: Uniform.h:79

References Seed1, and Seed2.

Referenced by Bernoulli(), Binomial(), BivariateGaussian(), LCE_Breed_base::checkPolygyny(), LCE_Cross::create_individual_ancestors(), LCE_Patch_Extinction::do_remove(), LCE_Selection_base::doViabilitySelection(), LCE_Disperse_EvolDisp::evoldisp(), LCE_Aging::execute(), LCE_Patch_Extinction::execute(), Exponential(), BinaryDataLoader::extractPop(), Metapop::fillPopulationFromSource(), LCE_Disperse_EvolDisp::fixdisp(), LCE_Breed_base::fullPolyginy_manyMales(), Gamma(), Gaussian(), LCE_Breed_Disperse::get_parent(), LCE_Disperse_base::getMigrationPatchBackward(), LCE_Disperse_base::getMigrationPatchForward(), TProtoQuanti::getMutationEffectBivariateDiallelic(), TT_BDMI::init_sequence(), TTDeletMutations_bitstring::init_sequence(), TTDispersal::init_sequence(), TTNeutralGenes::init_sequence(), TTQuanti::init_sequence(), LogNormal(), LCE_Breed_Selection::makeOffspringWithSelection(), LCE_Breed_Disperse::mate_selfing(), LCE_Disperse_EvolDisp::Migrate_Island(), LCE_Disperse_EvolDisp::Migrate_Island_Propagule(), LCE_Disperse_EvolDisp::Migrate_Lattice(), LCE_Disperse_ConstDisp::MigratePatchByNumber(), MultinomialOnNormalizedValarrayZipper_scrambleOut(), TTDispersal::mutate(), TTWolbachia::mutate(), TTNeutralGenes::mutate_2all(), TProtoQuanti::mutate_diallelic_HC(), TT_BDMI::mutate_diplo(), TT_BDMI::mutate_haplo(), TProtoQuanti::mutate_HC(), TTNeutralGenes::mutate_KAM(), TProtoQuanti::mutate_noHC(), TTDeletMutations_bitstring::mutate_noredraw(), TTDeletMutations_bitstring::mutate_noredraw_noBackMutation(), TTDeletMutations_bitstring::mutate_redraw(), TTNeutralGenes::mutate_SSM(), LCE_Breed_base::partialMonoginy(), LCE_Breed_base::partialPolyginy(), LCE_Breed_base::partialPolyginy_manyMales(), LCE_Breed_base::partialSelfing(), Poisson(), LCE_Patch_Extinction::rand_exp(), LCE_Patch_Extinction::rand_uniform(), RandBool(), LCE_Breed_base::random_hermaphrodite(), LCE_Breed_base::RandomMating(), RandULong(), GeneticMap::recombine(), LCE_Resize::regulateAgeClassNoBackup(), LCE_Resize::regulateAgeClassWithBackup(), LCE_Regulation::regulatePatch(), ParamsParser::runif(), Sample(), LCE_Cross::sampleAmongPop(), SampleSeq(), SampleSeqWithReciprocal(), LCE_Cross::sampleWithinPop(), ScrambleArrayUInt(), LCE_Disperse_base::setPropaguleTargets(), TTProtoWithMap::setRecombinationMapRandom(), setSpatialMatrix(), LCE_Breed_Disperse::stochasticFecundityGrowth(), Uniform(), LCE_Breed_Wolbachia::wolbachia_model_1(), LCE_Breed_base::WrightFisherPopulation(), and LCE_Breed_Quanti::WrightFisherPopulation().

◆ Uniform() [2/2]

static unsigned int RAND::Uniform ( unsigned int  max)
inlinestatic

Returns a uniformly distributed random number from [0.0, max[.

158 {
159 return (unsigned int)(Uniform() * max);
160 }

References Uniform().

Member Data Documentation

◆ Seed1

long RAND::Seed1 = 0
static

Referenced by init(), and Uniform().

◆ Seed2

long RAND::Seed2 = 98280582
static

Referenced by Uniform().


The documentation for this class was generated from the following files:

Generated for Nemo v2.3.56 by  doxygen 1.9.0 -- Nemo is hosted on  Download Nemo

Locations of visitors to this page
Catalogued on GSR