Nemo  2.3.46
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 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 void Sample (int from, int to, unsigned int num, int *result, bool replace)
 
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 ( )
private

Member Function Documentation

static double RAND::Bernoulli ( double  p)
inlinestatic

References Uniform().

Referenced by TTNeutralGenes::init_sequence().

407  {
408 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
409  return gsl_ran_bernoulli(r, p);
410 #else
411  double u = RAND::Uniform() ;
412 
413  if (u < p)
414  {
415  return 1 ;
416  }
417  else
418  {
419  return 0 ;
420  }
421 #endif
422  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static void RAND::BivariateGaussian ( double  sigma1,
double  sigma2,
double  rho,
double *  out1,
double *  out2 
)
inlinestatic

References Uniform().

Referenced by TProtoQuanti::getMutationEffectBivariateGaussian().

307  {
308 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
309  gsl_ran_bivariate_gaussian(r,sigma1,sigma2,rho,out1,out2);
310 #else
311  //gsl code:
312  double u, v, r2, scale;
313  //double *x = out, *y = (++out);
314  do
315  {
316  /* choose x,y in uniform square (-1,-1) to (+1,+1) */
317 
318  u = -1 + 2 * Uniform ();
319  v = -1 + 2 * Uniform ();
320 
321  /* see if it is in the unit circle */
322  r2 = u * u + v * v;
323  }
324  while (r2 > 1.0 || r2 == 0);
325 
326  scale = sqrt (-2.0 * log (r2) / r2);
327 
328  *out1 = sigma1 * u * scale;
329  *out2 = sigma2 * (rho * u + sqrt(1 - rho*rho) * v) * scale;
330 
331 #endif
332  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static double RAND::Exponential ( double  mu)
inlinestatic

References Uniform().

Referenced by TProtoDeletMutations_bitstring::set_effects_exp().

424  {
425  return -mu * log(RAND::Uniform());
426  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static double RAND::Gamma ( double  a,
double  b 
)
inlinestatic

References Gaussian(), and Uniform().

Referenced by TProtoDeletMutations_bitstring::set_effects_gamma().

362  {
363 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
364  return gsl_ran_gamma(r, a, b);
365 #else
366  /*pasting GSL code in here*/
367 
368  /* assume a > 0 */
369 
370  if (a < 1)
371  {
372  double u = Uniform();//might give singularity at 0....
373  return Gamma(1.0 + a, b) * pow (u, 1.0 / a);
374  }
375 
376  {
377  double x, v, u;
378  double d = a - 1.0 / 3.0;
379  double c = (1.0 / 3.0) / sqrt (d);
380 
381  while (1)
382  {
383  do
384  {
385  x = Gaussian(1.0); //should use the Ziggurat method?
386  v = 1.0 + c * x;
387  }
388  while (v <= 0);
389 
390  v = v * v * v;
391  u = Uniform();//might give singularity at 0....
392 
393  if (u < 1 - 0.0331 * x * x * x * x)
394  break;
395 
396  if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
397  break;
398  }
399 
400  return b * d * v;
401  }
402 
403 #endif
404  }
static double Gaussian(double sigma)
Definition: Uniform.h:238
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static double Gamma(double a, double b)
Definition: Uniform.h:362
static double RAND::gammln ( double  xx)
inlinestatic

From the Numerical Recieps.

Referenced by Poisson().

181  {
182  double x,y,tmp,ser=1.000000000190015;
183  static double cof[6]={76.18009172947146,-86.50532032941677,
184  24.01409824083091,-1.231739572450155,
185  0.1208650973866179e-2,-0.5395239384953e-5};
186  int j;
187  y=x=xx;
188  tmp=x+5.5;
189  tmp -= (x+0.5)*log(tmp);
190  for (j = 0; j < 6; ++j) ser += cof[j]/++y;
191 
192  return -tmp+log(2.5066282746310005*ser/x);
193 
194  }
static double RAND::Gaussian ( double  sigma)
inlinestatic

From the GSL.

References Uniform().

Referenced by Gamma(), LCE_Selection_base::getFitnessMultivariateGaussian_VE(), LCE_Selection_base::getFitnessUnivariateGaussian_VE(), LCE_Breed_base::getGaussianFecundity(), TProtoQuanti::getMutationEffectUnivariateGaussian(), TTQuanti::init_sequence(), LCE_Patch_Extinction::rand_gaussian(), and TTQuanti::set_value().

239  {
240 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
241 
242  return gsl_ran_gaussian_ratio_method (r, sigma);
243 
244 #else
245 
246 // double x, y, r2;
247 //
248 // do
249 // {
250 // /* choose x,y in uniform square (-1,-1) to (+1,+1) */
251 //
252 // x = -1 + 2 * Uniform ( );
253 // y = -1 + 2 * Uniform ( );
254 //
255 // /* see if it is in the unit circle */
256 // r2 = x * x + y * y;
257 // }
258 // while (r2 > 1.0 || r2 == 0);
259 //
260 // /* Box-Muller transform */
261 // return sigma * y * sqrt (-2.0 * log (r2) / r2);
262 
263  /*trying the ratio method from GSL*/
264  /* see code in gsl/randist/gauss.c */
265  double u, v, x, y, Q;
266  const double s = 0.449871; /* Constants from Leva */
267  const double t = -0.386595;
268  const double a = 0.19600;
269  const double b = 0.25472;
270  const double r1 = 0.27597;
271  const double r2 = 0.27846;
272 
273  do /* This loop is executed 1.369 times on average */
274  {
275  /* Generate a point P = (u, v) uniform in a rectangle enclosing
276  the K+M region v^2 <= - 4 u^2 log(u). */
277 
278  /* u in (0, 1] to avoid singularity at u = 0 */
279  u = 1 - RAND::Uniform();
280 
281  /* v is in the asymmetric interval [-0.5, 0.5). However v = -0.5
282  is rejected in the last part of the while clause. The
283  resulting normal deviate is strictly symmetric about 0
284  (provided that v is symmetric once v = -0.5 is excluded). */
285  v = RAND::Uniform() - 0.5;
286 
287  /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
288  much (for efficiency). */
289  v *= 1.7156;
290 
291  /* Compute Leva's quadratic form Q */
292  x = u - s;
293  y = fabs (v) - t;
294  Q = x * x + y * (a * y - b * x);
295 
296  /* Accept P if Q < r1 (Leva) */
297  /* Reject P if Q > r2 (Leva) */
298  /* Accept if v^2 <= -4 u^2 log(u) (K+M) */
299  /* This final test is executed 0.012 times on average. */
300  }
301  while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));
302 
303  return sigma * (v / u);
304 #endif
305  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static void RAND::init ( unsigned long  seed)
inlinestatic

Initialize the random generator's seed.

Referenced by SimRunner::init_random_seed().

80  {
81 
82 #if defined(HAS_SPRNG)
83 
84  init_sprng( SPRNG_LFG, seed, SPRNG_DEFAULT );
85 
86 #elif defined(HAS_GSL)
87 
88  init_gsl( gsl_rng_mt19937, seed );
89 
90 #else
91 
92  Seed1 = seed;
93 
94 #endif
95  }
static long Seed1
Definition: Uniform.h:76
static double RAND::LogNormal ( double  zeta,
double  sigma 
)
inlinestatic

References Uniform().

Referenced by LCE_Patch_Extinction::rand_lognormal(), and TProtoDeletMutations_bitstring::set_effects_lognorm().

334  {
335 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
336  return gsl_ran_lognormal(r,zeta,sigma);
337 #else
338  //this is the GSL code:
339  double u, v, r2, normal, z;
340 
341  do
342  {
343  /* choose x,y in uniform square (-1,-1) to (+1,+1) */
344 
345  u = -1 + 2 * Uniform();
346  v = -1 + 2 * Uniform();
347 
348  /* see if it is in the unit circle */
349  r2 = u * u + v * v;
350  }
351  while (r2 > 1.0 || r2 == 0);
352 
353  normal = u * sqrt (-2.0 * log (r2) / r2);
354 
355  z = exp (sigma * normal + zeta);
356 
357  return z;
358 
359 #endif
360  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static double RAND::Poisson ( double  mean)
inlinestatic

From the Numerical Recieps.

References gammln(), and Uniform().

Referenced by LCE_Breed_base::getPoissonFecundity(), TT_BDMI::mutate_diplo(), TT_BDMI::mutate_haplo(), TTQuanti::mutate_HC(), TTQuanti::mutate_noHC(), TTDeletMutations_bitstring::mutate_noredraw(), TTDeletMutations_bitstring::mutate_redraw(), LCE_Patch_Extinction::rand_poisson(), GeneticMap::recombine(), and LCE_Breed_Disperse::stochasticLogisticGrowth().

196  {
197 
198 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
199  return gsl_ran_poisson(r, mean);
200 #else
201  static double sq,alxm,g,oldm=(-1.0);
202  double em,t,y;
203 
204  if (mean < 12.0)
205  {
206  if (mean != oldm){
207  oldm=mean;
208  g=exp(-mean);
209  }
210  em = -1;
211  t=1.0;
212  do {
213  ++em;
214  t *= Uniform();
215  } while (t > g);
216  }else
217  {
218  if (mean != oldm)
219  {
220  oldm=mean;
221  sq=sqrt(2.0*mean);
222  alxm=log(mean);
223  g=mean*alxm-gammln(mean+1.0);
224  }
225  do {
226  do {
227  y=tan(M_PI*Uniform());
228  em=sq*y+mean;
229  } while (em < 0.0);
230  em=floor(em);
231  t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
232  } while (Uniform() > t);
233  }
234  return em;
235 #endif
236  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static double gammln(double xx)
From the Numerical Recieps.
Definition: Uniform.h:181
static bool RAND::RandBool ( )
inlinestatic

Returns a random boolean.

References Uniform().

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

139  {
140  //generate a random int (or long)
141 #ifdef HAS_SPRNG
142  static int intrand = isprng();
143 #elif defined(HAS_GSL)
144  static unsigned long int intrand = gsl_rng_get(r);
145 #else
146  static int intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
147 #endif
148  //read up to the first 16 bits
149  static unsigned int num = 0;
150  //redraw an number after that limit
151  if ( ++num > 16 ) {
152  num = 0;
153 
154 #ifdef HAS_SPRNG
155  intrand = isprng();
156 #elif defined(HAS_GSL)
157  intrand = gsl_rng_get(r);
158 #else
159  intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
160 #endif
161 
162  }
163  return (intrand & (1 << num) );
164 
165  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static unsigned long RAND::RandULong ( )
inlinestatic

Return a random unsigned long, from uniform distribution.

References Uniform().

168  {
169 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
170  return gsl_rng_get(r);
171 #else
172  unsigned long rnd, limit = 0x10000000; //=2^28 this gives ~7% redraws
173  do{
174  rnd = (unsigned long)(Uniform()*ULONG_MAX);
175  }while(rnd < limit);
176 
177  return rnd;
178 #endif
179  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static void RAND::Sample ( int  from,
int  to,
unsigned int  num,
int result,
bool  replace 
)
inlinestatic

References Uniform().

Referenced by MPFileHandler::createAndPrintSample().

513  {
514  assert(from < to);
515 
516  unsigned int seq_length = to - from; //don't include last (to)
517 
518  assert(num <= seq_length);
519 
520  int *seq = new int [seq_length];
521 
522  seq[0] = from;
523 
524  for(unsigned int i = 1; seq[i-1] < to && i < seq_length; ++i) seq[i] = seq[i-1] + 1;
525 
526  if(!replace) { //without replacement
527 
528  unsigned int size = seq_length, pos, last = seq_length - 1;
529 
530  for (unsigned int i = 0; i < num; i++) {
531  pos = RAND::Uniform(size);
532  result[i] = seq[pos];
533  seq[pos] = seq[last];
534 // seq[last] = result[i]; useless operation
535  size--; last--;
536  }
537 
538  } else { //with replacement
539  for (unsigned int i = 0; i < num; i++) result[i] = seq[ RAND::Uniform(seq_length) ];
540  }
541 
542  delete [] seq;
543  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static void RAND::SampleSeq ( int  from,
int  to,
int  by,
unsigned int  num,
int result,
bool  replace = false 
)
inlinestatic

References Uniform().

546  {
547  assert(from < to && by < to && by > 0);
548 
549  unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
550 
551  assert(num <= seq_length);
552 
553  int *seq = new int [seq_length];
554 
555  seq[0] = from;
556 
557  for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
558 
559  if(!replace) { //without replacement
560 
561  unsigned int size = seq_length, pos, last = seq_length - 1;
562 
563  for (unsigned int i = 0; i < num; i++) {
564  pos = RAND::Uniform(size);
565  result[i] = seq[pos];
566  seq[pos] = seq[last];
567 // seq[last] = result[i]; useless operation
568  size--; last--;
569  }
570 
571  } else { //with replacement
572  for (unsigned int i = 0; i < num; i++) result[i] = seq[ RAND::Uniform(seq_length) ];
573  }
574 
575  delete [] seq;
576  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
static void RAND::SampleSeqWithReciprocal ( int  from,
int  to,
int  by,
unsigned int  num1,
int result1,
unsigned int  num2,
int result2 
)
inlinestatic

References Uniform().

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

579  {
580  assert(from < to && by < to && by > 0);
581 
582  unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
583 
584  assert(num1 + num2 == seq_length);
585 
586  int *seq = new int [seq_length];
587 
588  seq[0] = from;
589 
590  for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
591 
592  unsigned int size = seq_length, pos, last = seq_length - 1;
593 
594  for (unsigned int i = 0; i < num1; i++) {
595  pos = RAND::Uniform(size);
596  result1[i] = seq[pos];
597  seq[pos] = seq[last];
598  seq[last] = result1[i];
599  size--; last--;
600  }
601 
602  for (unsigned int i = 0; i < num2 && i < size; i++)
603  result2[i] = seq[i];
604 
605  delete [] seq;
606  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101
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.

References Seed2.

Referenced by Bernoulli(), 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(), TTQuanti::init_sequence(), TTNeutralGenes::init_sequence(), TTDeletMutations_bitstring::init_sequence(), TT_BDMI::init_sequence(), LogNormal(), LCE_Breed_Selection::makeOffspringWithSelection(), LCE_Breed_Disperse::mate_selfing(), LCE_Disperse_EvolDisp::Migrate_Island(), LCE_Disperse_EvolDisp::Migrate_Island_Propagule(), TTWolbachia::mutate(), TTNeutralGenes::mutate_2all(), TT_BDMI::mutate_diplo(), TT_BDMI::mutate_haplo(), TTQuanti::mutate_HC(), TTNeutralGenes::mutate_KAM(), TTQuanti::mutate_noHC(), TTDeletMutations_bitstring::mutate_noredraw(), 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(), Sample(), LCE_Cross::sampleAmongPop(), SampleSeq(), SampleSeqWithReciprocal(), LCE_Cross::sampleWithinPop(), LCE_Disperse_base::setPropaguleTargets(), TTProtoWithMap::setRecombinationMapRandom(), setSpatialMatrix(), LCE_Breed_Disperse::stochasticFecundityGrowth(), Uniform(), and LCE_Breed_Wolbachia::wolbachia_model_1().

101  {
102 
103 #ifdef HAS_SPRNG
104  return sprng();
105 #elif defined(HAS_GSL)
106  return gsl_rng_uniform(r);
107 #else
108  register long z, w;
109 
110  do{
111  w = Seed1 / 53668;
112 
113  Seed1 = 40014 * (Seed1 - w * 53668) - w * 12211;
114 
115  if (Seed1 < 0) Seed1 += 2147483563;
116 
117  w = (Seed2 / 52774);
118 
119  Seed2 = 40692 * (Seed2 - w * 52774) - w * 3791;
120 
121  if (Seed2 < 0) Seed2 += 2147483399;
122 
123  z = Seed1 - Seed2;
124 
125  if (z < 1) z += 2147483562;
126 
127  }while (!((z * 4.656613e-10) < 1.0));
128 
129  return (z * 4.656613e-10);
130 #endif
131  }
static long Seed2
Definition: Uniform.h:76
static long Seed1
Definition: Uniform.h:76
static unsigned int RAND::Uniform ( unsigned int  max)
inlinestatic

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

References Uniform().

134  {
135  return (unsigned int)(Uniform() * max);
136  }
static double Uniform()
Generates a random number from [0.0, 1.0[ uniformly distributed.
Definition: Uniform.h:101

Member Data Documentation

long RAND::Seed1 = 0
static
long RAND::Seed2 = 98280582
static

Referenced by Uniform().


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

Generated for Nemo v2.3.0 by  doxygen 1.8.8 --
Catalogued on GSR