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 size_t	Discrete (const gsl_ran_discrete_t *g)
	Calling the GSL ran_discrete function. More...

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

   {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_bernoulli(r, p); 
 #else    
     double u = RAND::Uniform() ;
     
     if (u < p)
     {
       return 1 ;
     }
     else
     {
       return 0 ;
     }
 #endif
   }

References Uniform().

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

◆ Beta()

static double RAND::Beta	(	const double	a,
		const double	b
	)

inlinestatic

   {
     /*from Knuth*/
     double x1 = RAND::Gamma(a, 1.0);
     double x2 = RAND::Gamma(b, 1.0);
     
     return x1 / (x1 + x2);
   }

References Gamma().

Referenced by Binomial().

◆ Binomial()

static double RAND::Binomial	(	double	p,
		unsigned int	n
	)

inlinestatic

   {
     /*implements Knuth method*/
     unsigned int i, a, b, k = 0;
     
     while (n > 10)        /* This parameter is tunable */
     {
       double X;
       a = 1 + (n / 2);
       b = 1 + n - a;
       
       X = RAND::Beta((double) a, (double) b);
       
       if (X >= p)
       {
         n = a - 1;
         p /= X;
       }
       else
       {
         k += a;
         n = b - 1;
         p = (p - X) / (1 - X);
       }
     }
     
     for (i = 0; i < n; i++)
     {
       double u = RAND::Uniform();
       if (u < p)
         k++;
     }
     
     return k;    
     
   }

References Beta(), and Uniform().

Referenced by Binomial2(), TProtoNeutralGenes::get_num_mutations(), TProtoQuanti::get_num_mutations(), Multinomial(), MultinomialOnNormalizedValarray(), MultinomialOnNormalizedValarray_expandedOut(), and MultinomialOnNormalizedValarray_scrambleOut().

◆ Binomial2()

static unsigned int RAND::Binomial2	(	double	p,
		unsigned int	n
	)

inlinestatic

                                                                  {
     
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_binomial (r,  p,  n);
 #else
     return Binomial(p, n);
 #endif
   }

References Binomial().

◆ BivariateGaussian()

static void RAND::BivariateGaussian	(	double	sigma1,
		double	sigma2,
		double	rho,
		double *	out1,
		double *	out2
	)

inlinestatic

                                                                                                            {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     gsl_ran_bivariate_gaussian(r,sigma1,sigma2,rho,out1,out2);
 #else
     //gsl code:
     double u, v, r2, scale;
    //double *x = out, *y = (++out);
     do
     {
       /* choose x,y in uniform square (-1,-1) to (+1,+1) */
       
       u = -1 + 2 * Uniform ();
       v = -1 + 2 * Uniform ();
       
       /* see if it is in the unit circle */
       r2 = u * u + v * v;
     }
     while (r2 > 1.0 || r2 == 0);
     
     scale = sqrt (-2.0 * log (r2) / r2);
     
     *out1 = sigma1 * u * scale;
     *out2 = sigma2 * (rho * u + sqrt(1 - rho*rho) * v) * scale;
     
 #endif
   }

References Uniform().

Referenced by TProtoQuanti::getMutationEffectBivariateGaussian(), and TProtoQuanti::getMutationEffectBivariateGaussianLocSpec().

◆ Discrete()

static size_t RAND::Discrete ( const gsl_ran_discrete_t * g )

inlinestatic

Calling the GSL ran_discrete function.

   {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_discrete(r, g);
 #else
     /* GSL code here, modified to call our Uniform() */
     size_t c=0;
     double u,f;
     u = Uniform();
  
 #if KNUTH_CONVENTION
     c = (u*(g->K));
 #else
     u *= g->K;
     c = u;
     u -= c;
 #endif
  
     f = (g->F)[c];
     /* fprintf(stderr,"c,f,u: %d %.4f %f\n",c,f,u); */
     if (f == 1.0) return c;
  
     if (u < f) {
         return c;
     }
     else {
         return (g->A)[c];
     }
  
 #endif
   }

References Uniform().

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

◆ Exponential()

static double RAND::Exponential ( double mu )

inlinestatic

                                                {
     return -mu * log(RAND::Uniform());
   }

References Uniform().

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

◆ free()

static void RAND::free ( )

inlinestatic

Memory de-allocation.

   {
  
 #if defined(HAS_SPRNG)
     // do nothing with pointer?
 #elif defined(HAS_GSL)
     gsl_rng_free(r);
 #endif
  
   }

Referenced by SimRunner::run().

◆ Gamma()

static double RAND::Gamma	(	double	a,
		double	b
	)

inlinestatic

                                                   {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_gamma(r, a, b);
 #else
     /*pasting GSL code in here*/
     
     /* assume a > 0 */
     
     if (a < 1)
     {
       double u = Uniform();//might give singularity at 0....
       return Gamma(1.0 + a, b) * pow (u, 1.0 / a);
     }
     
     {
       double x, v, u;
       double d = a - 1.0 / 3.0;
       double c = (1.0 / 3.0) / sqrt (d);
       
       while (1)
       {
         do
         {
           x = Gaussian(1.0); //should use the Ziggurat method?
           v = 1.0 + c * x;
         }
         while (v <= 0);
         
         v = v * v * v;
         u = Uniform();//might give singularity at 0....
         
         if (u < 1 - 0.0331 * x * x * x * x) 
           break;
         
         if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
           break;
       }
       
       return b * d * v;
     }
     
 #endif
   }

References Gaussian(), and Uniform().

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

◆ gammln()

static double RAND::gammln ( double xx )

inlinestatic

From the Numerical Recieps.

                                           {
     double x,y,tmp,ser=1.000000000190015;
     static double cof[6]={76.18009172947146,-86.50532032941677,
       24.01409824083091,-1.231739572450155,
       0.1208650973866179e-2,-0.5395239384953e-5};
     int j;
     y=x=xx;
     tmp=x+5.5;
     tmp -= (x+0.5)*log(tmp);
     for (j = 0; j < 6; ++j) ser += cof[j]/++y;
  
     return -tmp+log(2.5066282746310005*ser/x);
  
   }

Referenced by Poisson().

◆ Gaussian()

static double RAND::Gaussian ( double sigma )

inlinestatic

From the GSL.

   {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     
     return gsl_ran_gaussian_ziggurat (r, sigma);
     
 #else
 //    double x, y, r2;
 //    
 //    do
 //    {
 //      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
 //      
 //      x = -1 + 2 * Uniform ( );
 //      y = -1 + 2 * Uniform ( );
 //      
 //      /* see if it is in the unit circle */
 //      r2 = x * x + y * y;
 //    }
 //    while (r2 > 1.0 || r2 == 0);
 //    
 //    /* Box-Muller transform */
 //    return sigma * y * sqrt (-2.0 * log (r2) / r2);
     
   /*trying the ratio method from GSL*/
   /* see code in gsl/randist/gauss.c */  
     double u, v, x, y, Q;
     const double s = 0.449871;    /* Constants from Leva */
     const double t = -0.386595;
     const double a = 0.19600;
     const double b = 0.25472;
     const double r1 = 0.27597;
     const double r2 = 0.27846;
     
     do                            /* This loop is executed 1.369 times on average  */
     {
       /* Generate a point P = (u, v) uniform in a rectangle enclosing
        the K+M region v^2 <= - 4 u^2 log(u). */
       
       /* u in (0, 1] to avoid singularity at u = 0 */
       u = 1 - RAND::Uniform();
       
       /* v is in the asymmetric interval [-0.5, 0.5).  However v = -0.5
        is rejected in the last part of the while clause.  The
        resulting normal deviate is strictly symmetric about 0
        (provided that v is symmetric once v = -0.5 is excluded). */
       v = RAND::Uniform() - 0.5;
       
       /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
        much (for efficiency). */
       v *= 1.7156;
       
       /* Compute Leva's quadratic form Q */
       x = u - s;
       y = fabs (v) - t;
       Q = x * x + y * (a * y - b * x);
       
       /* Accept P if Q < r1 (Leva) */
       /* Reject P if Q > r2 (Leva) */
       /* Accept if v^2 <= -4 u^2 log(u) (K+M) */
       /* This final test is executed 0.012 times on average. */
     }
     while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));
     
     return sigma * (v / u);
 #endif
   }

References Uniform().

Referenced by Gamma(), LCE_PhenotypeExpression::get_env_cue_noise(), LCE_Selection_base::getFitnessMultivariateGaussian_VE(), LCE_Selection_base::getFitnessUnivariateGaussian_VE(), LCE_Breed_base::getGaussianFecundity(), TProtoQuanti::getMutationEffectUnivariateGaussian(), TProtoQuanti::getMutationEffectUnivariateGaussianLocSpec(), TTDispersal::init_sequence(), TTQuanti_continuous_full_pleio::init_sequence(), TTQuanti_continuous_var_pleio::init_sequence(), TTQuanti_continuous_no_pleio::init_sequence(), TTQuanti_continuous_full_pleio_epistasis::init_sequence(), TTQuanti_continuous_no_pleio_epistasis::init_sequence(), LCE_Patch_Extinction::rand_gaussian(), ParamsParser::rnorm(), and TProtoQuanti::set_trait_value_VE().

◆ init()

static void RAND::init ( unsigned long seed )

inlinestatic

Initialize the random generator's seed.

                                         {
  
 #if defined(HAS_SPRNG)
     // initializing the SPRNG stream with default param and LFG (0) generator
  
     stream = SelectType(4); //MLFG Modified Lagged Fibonacci
  
     stream->init_sprng( _myenv->workerRank(), _myenv->workerCount()+1, seed, SPRNG_DEFAULT );
  
 #ifdef DEBUG_MPI
     message("--- initialized the SPRNG random generator on rank %i\n", _myenv->workerRank());
  
     stream->print_sprng();
  
     message("--- process %i random number: %.14f\n", _myenv->workerRank(), stream->sprng());
 #endif
  
 #elif defined(HAS_GSL)
     
     init_gsl( gsl_rng_mt19937, seed );
   
 #else
     
     Seed1 = seed;
     
 #endif
   }

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

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

◆ LogNormal()

static double RAND::LogNormal	(	double	zeta,
		double	sigma
	)

inlinestatic

                                                              {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_lognormal(r,zeta,sigma);
 #else
     //this is the GSL code:
     double u, v, r2, normal, z;
     
     do
     {
       /* choose x,y in uniform square (-1,-1) to (+1,+1) */
       
       u = -1 + 2 * Uniform();
       v = -1 + 2 * Uniform();
       
       /* see if it is in the unit circle */
       r2 = u * u + v * v;
     }
     while (r2 > 1.0 || r2 == 0);
     
     normal = u * sqrt (-2.0 * log (r2) / r2);
     
     z =  exp (sigma * normal + zeta);
     
     return z;
     
 #endif
   }

References Uniform().

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

◆ Multinomial()

static void RAND::Multinomial	(	size_t	K,
		unsigned int	N,
		const double	p[],
		unsigned int	n[]
	)

inlinestatic

   {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_multinomial (r, K, N, p, n);
 #else
     // GSL code here:
     size_t k;
     double norm = 0.0;
     double sum_p = 0.0;
  
     unsigned int sum_n = 0;
  
     /* p[k] may contain non-negative weights that do not sum to 1.0.
      * Even a probability distribution will not exactly sum to 1.0
      * due to rounding errors.
      */
  
     for (k = 0; k < K; k++)
       {
         norm += p[k];
       }
  
     for (k = 0; k < K; k++)
       {
         if (p[k] > 0.0)
           {
             n[k] = Binomial (p[k] / (norm - sum_p), N - sum_n); //MOD
           }
         else
           {
             n[k] = 0;
           }
  
         sum_p += p[k];
         sum_n += n[k];
       }
  
 #endif
   }

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.

   {
     // GSL code here: FG 2020: modified for normalized arrays
     size_t k;
     double sum_p = 0.0;
  
     unsigned int sum_n = 0;
  
     for (k = 0; k < K; k++)
     {
       n[k] = RAND::Binomial (p[k] / (1.0 - sum_p), N - sum_n); //MOD FG
       sum_p += p[k];
       sum_n += n[k];
     }
   }

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.

   {
     // GSL code here: FG 2020: modified for normalized arrays
     size_t k;
     double sum_p = 0.0;
     unsigned int _n, sum_n = 0;
     unsigned int i = 0, pos = 0;
  
     for (k = 0; k < K; k++) {
       //draw number of events k with probability p[k]
       _n =  RAND::Binomial (p[k] / (1.0 - sum_p), N - sum_n);
  
       for(i = 0; i < _n; i++)
         n[pos++] = k;
  
       sum_p += p[k];
       sum_n += _n;
       assert(pos <= N);
     }
     assert(sum_n == N);
   }

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.

   {
     // GSL code here: FG 2020: modified for normalized arrays
     size_t k;
     double sum_p = 0.0;
     unsigned int _n, sum_n = 0;
     unsigned int i = 0, pos = 0;
  
     for (k = 0; k < K; k++) {
       //draw number of events k with probability p[k]
       _n =  RAND::Binomial (p[k] / (1.0 - sum_p), N - sum_n);
  
       for(i = 0; i < _n; i++)
         n[pos++] = k;
  
       sum_p += p[k];
       sum_n += _n;
       assert(pos <= N);
     }
     assert(sum_n == N);
  
     //scramble the output array
     RAND::ScrambleArrayUInt(N, n);
  
   }

References Binomial(), and ScrambleArrayUInt().

◆ 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.

   {
     size_t k;
     double sum = 0.0;
     double rand;
     unsigned int i = 0;
  
     for (k = 0; k < N; k++) {
  
       rand = std::min(0.99999995, RAND::Uniform()); // to avoid problems with rounding errors
       i = 0;
       sum = p[i];
  
       //zip through the proba array:
       while( sum < rand ) {
         i++;
         sum += p[i];
       }
  
       n[k] = i;
  
     }
     //scramble the output array
     RAND::ScrambleArrayUInt(N, n);
  
   }

References ScrambleArrayUInt(), and Uniform().

◆ Poisson()

static double RAND::Poisson ( double mean )

inlinestatic

From the Numerical Recieps.

                                              {
     
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_ran_poisson(r, mean);
 #else
     static double sq,alxm,g,oldm=(-1.0);
     double em,t,y;
  
     if (mean < 12.0)
     {
       if (mean != oldm){
         oldm=mean;
         g=exp(-mean);
       }
       em = -1;
       t=1.0;
       do {
         ++em;
         t *= Uniform();
       } while (t > g);
     }else
     {
       if (mean != oldm)
       {
         oldm=mean;
         sq=sqrt(2.0*mean);
         alxm=log(mean);
         g=mean*alxm-gammln(mean+1.0);
       }
       do {
         do {
           y=tan(M_PI*Uniform());
           em=sq*y+mean;
         } while (em < 0.0);
         em=floor(em);
         t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
       } while (Uniform() > t);
     }
     return em;
 #endif
   }

References gammln(), and Uniform().

Referenced by LCE_Breed_base::getPoissonFecundity(), TTDispersal::mutate(), TT_BDMI::mutate_diplo(), TT_BDMI::mutate_haplo(), 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().

◆ RandBool()

static bool RAND::RandBool ( )

inlinestatic

Returns a random boolean.

                                 {
     //generate a random int (or long)
 #ifdef HAS_SPRNG
     static int intrand = stream->isprng();
 #elif defined(HAS_GSL)
     static unsigned long int intrand =  gsl_rng_get(r);
 #else
     static int intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
 #endif
     //read up to the first 16 bits
     static unsigned int num = 0;
     //redraw an number after that limit
     if ( ++num > 16 ) {
       num = 0;
  
 #ifdef HAS_SPRNG
       intrand = stream->isprng();
 #elif defined(HAS_GSL)
       intrand =  gsl_rng_get(r);
 #else
       intrand = (int)(Uniform()*(double)std::numeric_limits<int>::max());
 #endif
  
     }
     return (intrand & (1 << num) );
  
   }

References Uniform().

◆ RandULong()

static unsigned long RAND::RandULong ( )

inlinestatic

Return a random unsigned long, from uniform distribution.

                                           {
 #if defined(HAS_GSL) && !defined(HAS_SPRNG)
     return gsl_rng_get(r);
 #else
     unsigned long rnd, limit = 0x10000000; //=2^28 this gives ~7% redraws
     do{
       rnd = (unsigned long)(Uniform()*ULONG_MAX);
     }while(rnd < limit);
  
     return rnd;
 #endif
   }

References Uniform().

Referenced by TTNeutralGenes_bitstring::inherit().

◆ 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

from	starting number of the series
to	last number of the series
num	number of elements to draw within [from, to)
result	container to hold the resulting randomized sequence of integers

   {
     assert(from < to);
  
     unsigned int seq_length = to - from; //don't include last (to)
  
     assert(num <= seq_length);
  
     // we build a sequence of indexes to choose from
     int *seq = new int [seq_length];
  
     seq[0] = from;
  
     for(unsigned int i = 1; seq[i-1] < to && i < seq_length; ++i)
       seq[i] = seq[i-1] + 1;
  
     if(!replace) { //without replacement
  
       unsigned int size = seq_length, pos, last = seq_length - 1;
  
       for (unsigned int i = 0; i < num; i++) {
         pos = RAND::Uniform(size);
         result[i] = seq[pos];
         seq[pos] = seq[last];
         //                seq[last] = result[i]; useless operation
         size--; last--;
       }
  
     } else { //with replacement
       for (unsigned int i = 0; i < num; i++)
         result[i] = seq[ RAND::Uniform(seq_length) ];
     }
  
     delete [] seq;
   }

References Uniform().

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

◆ SampleSeq()

static void RAND::SampleSeq	(	int	from,
		int	to,
		int	by,
		unsigned int	num,
		int *	result,
		bool	replace = `false`
	)

inlinestatic

   {
     assert(from < to && by < to && by > 0);
  
     unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
  
     assert(num <= seq_length);
  
     int *seq = new int [seq_length];
  
     seq[0] = from;
  
     for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
  
     if(!replace) { //without replacement
  
       unsigned int size = seq_length, pos, last = seq_length - 1;
  
       for (unsigned int i = 0; i < num; i++) {
         pos = RAND::Uniform(size);
         result[i] = seq[pos];
         seq[pos] = seq[last];
         //        seq[last] = result[i]; useless operation
         size--; last--;
       }
  
     } else { //with replacement
       for (unsigned int i = 0; i < num; i++) result[i] = seq[ RAND::Uniform(seq_length) ];
     }
  
     delete [] seq;
   }

References Uniform().

◆ SampleSeqWithReciprocal()

static void RAND::SampleSeqWithReciprocal	(	int	from,
		int	to,
		int	by,
		unsigned int	num1,
		int *	result1,
		unsigned int	num2,
		int *	result2
	)

inlinestatic

   {
     assert(from < to && by < to && by > 0);
  
     unsigned int seq_length = (int)((to - from) / by); //don't include last (to)
  
     assert(num1 + num2 == seq_length);
  
     int *seq = new int [seq_length];
  
     seq[0] = from;
  
     for(unsigned int i = 1; seq[i-1] + by < to && i < seq_length; ++i) seq[i] = seq[i-1] + by;
  
     unsigned int size = seq_length, pos, last = seq_length - 1;
  
     for (unsigned int i = 0; i < num1; i++) {
       pos = RAND::Uniform(size);
       result1[i] = seq[pos];
       seq[pos] = seq[last];
       seq[last] = result1[i];
       size--; last--;
     }
  
     for (unsigned int i = 0; i < num2 && i < size; i++)
       result2[i] = seq[i];
  
     delete [] seq;
   }

References Uniform().

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

◆ ScrambleArrayUInt()

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

inlinestatic

Randomize the elements within an array.

Parameters

length	the length of the array
array	the array to scramble

   {
     unsigned int size = length, pos, last = length - 1, num = length -1, el;
  
     // with stop after 2 elements are left
     // this algo allows swapping an element with itself, which is fine
     for (unsigned int i = 0; i < num; i++) {
       pos = RAND::Uniform(size);
       el = array[last];
       array[last] = array[pos];
       array[pos] = el;
       size--; last--;
     }
  
      assert(size == 1); //we stopped before swapping the last (first) element with itself
   }

References Uniform().

Referenced by MultinomialOnNormalizedValarray_scrambleOut(), and MultinomialOnNormalizedValarrayZipper_scrambleOut().

◆ 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.

                                   {
  
 #ifdef HAS_SPRNG
     return stream->sprng();
 #elif defined(HAS_GSL)
     return gsl_rng_uniform(r);
 #else
     register long z, w;
  
     do{
       w = Seed1 / 53668;
  
       Seed1 = 40014 * (Seed1 - w * 53668) - w * 12211;
  
       if (Seed1 < 0) Seed1 += 2147483563;
  
       w = (Seed2 / 52774);
  
       Seed2 = 40692 * (Seed2 - w * 52774) - w * 3791;
  
       if (Seed2 < 0) Seed2 += 2147483399;
  
       z = Seed1 - Seed2;
  
       if (z < 1) z += 2147483562;
  
     }while (!((z * 4.656613e-10) < 1.0));
  
     return (z * 4.656613e-10);
 #endif
   }

References Seed1, and Seed2.

Referenced by Bernoulli(), Binomial(), BivariateGaussian(), LCE_Breed_base::checkPolygyny(), LCE_Cross::create_individual_ancestors(), Discrete(), LCE_Patch_Extinction::do_remove(), LCE_Selection_base::doViabilitySelection(), LCE_Disperse_EvolDisp::evoldisp(), LCE_BreedAssortativeMating::execute(), 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::getMigrationIndex(), LCE_Disperse_base::getMigrationPatchBackward(), LCE_Disperse_base::getMigrationPatchForward(), TProtoQuanti::getMutationEffectBivariateDiallelic(), LCE_QuantiInit::init_allele_freq(), TT_BDMI::init_sequence(), TTDeletMutations_bitstring::init_sequence(), TTDispersal::init_sequence(), TTNeutralGenes_byte::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(), TTNeutralGenes_bitstring::mutate(), TTWolbachia::mutate(), TTNeutralGenes_byte::mutate_2all(), TProtoQuanti::mutate_diallelic_no_pleio(), TProtoQuanti::mutate_diallelic_pleio(), TProtoQuanti::mutate_diallelic_var_pleio(), TT_BDMI::mutate_diplo(), TProtoQuanti::mutate_full_pleio(), TT_BDMI::mutate_haplo(), TProtoQuanti::mutate_inplace_full_pleio(), TProtoQuanti::mutate_inplace_no_pleio(), TProtoQuanti::mutate_inplace_var_pleio(), TTNeutralGenes_byte::mutate_KAM(), TProtoQuanti::mutate_no_pleio(), TTDeletMutations_bitstring::mutate_noredraw(), TTDeletMutations_bitstring::mutate_noredraw_noBackMutation(), TTDeletMutations_bitstring::mutate_redraw(), TTNeutralGenes_byte::mutate_SSM(), TProtoQuanti::mutate_var_pleio(), 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_BreedAssortativeMating::ScrambleContainer(), 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[.

   {
     return (unsigned int)(Uniform() * max);
   }

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:

Static Public Member Functions

Static Public Attributes

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ RAND()

Member Function Documentation

◆ Bernoulli()

◆ Beta()

◆ Binomial()

◆ Binomial2()

◆ BivariateGaussian()

◆ Discrete()

◆ Exponential()

◆ free()

◆ Gamma()

◆ gammln()

◆ Gaussian()

◆ init()

◆ LogNormal()

◆ Multinomial()

◆ MultinomialOnNormalizedValarray()

◆ MultinomialOnNormalizedValarray_expandedOut()

◆ MultinomialOnNormalizedValarray_scrambleOut()

◆ MultinomialOnNormalizedValarrayZipper_scrambleOut()

◆ Poisson()

◆ RandBool()

◆ RandULong()

◆ Sample()

◆ SampleSeq()

◆ SampleSeqWithReciprocal()

◆ ScrambleArrayUInt()

◆ Uniform() [1/2]

◆ Uniform() [2/2]

Member Data Documentation

◆ Seed1

◆ Seed2