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 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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

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

  {
#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().

Here is the caller graph for this function:

RandBool()

static bool RAND::RandBool ( )

inlinestatic

Returns a random boolean.

One bit is consumed per call from a cached random word, which is redrawn once its bits are exhausted. The shift must be done on an unsigned long: with the former 1 << num the shift count wrapped at 32, so bits 32-63 returned the same values as bits 0-31 of the same draw and the sequence repeated every 32 calls. This must stay consistent with RandULong() below: while RandULong() filled only the lower 32 bits, that same wrap was what kept these reads on the random half of the word, so fixing either one alone biases the result (see the CHANGELOG for 2.4.2).

                                {
    //number of usable random bits per draw, and the cached draw itself
#ifdef HAS_SPRNG
    static const unsigned int NBITS = 31; //isprng() returns a non-negative int
    static unsigned long intrand = (unsigned long)stream->isprng();
#else
    static const unsigned int NBITS = (unsigned int)(CHAR_BIT * sizeof(unsigned long));
    static unsigned long intrand = RandULong();
#endif
    static unsigned int num = 0;
 
    //redraw once all bits of the cached word have been consumed
    if ( num == NBITS ) {
      num = 0;
 
#ifdef HAS_SPRNG
      intrand = (unsigned long)stream->isprng();
#else
      intrand = RandULong();
#endif
 
    }
    return ( intrand >> num++ ) & 1UL;
 
  }

References RandULong().

Referenced by LCE_Patch_Extinction::do_remove(), TProtoQuanti::getMutationEffectBivariateDiallelic(), LCE_Breed_base::getOffsprgSexFixed(), LCE_Breed_base::getOffsprgSexRandom(), TTDispersal::inherit(), TTNeutralGenes_bitstring::inherit(), TProtoBDMI::inherit_free(), TProtoDeletMutations_bitstring::inherit_free(), TProtoQuanti::inherit_free(), TProtoNeutralGenes::inherit_free(), TTNeutralGenes_bitstring::init_sequence(), TTQuanti_continuous_full_pleio::init_sequence(), TTQuanti_continuous_var_pleio::init_sequence(), TTQuanti_continuous_no_pleio::init_sequence(), TTQuanti_diallelic_no_pleio::init_sequence(), TTQuanti_diallelic_full_pleio::init_sequence(), TTQuanti_diallelic_var_pleio::init_sequence(), TTQuanti_diallelic_bitstring_no_pleio::init_sequence(), TTQuanti_diallelic_bitstring_full_pleio::init_sequence(), TTQuanti_diallelic_bitstring_var_pleio::init_sequence(), TTQuanti_continuous_full_pleio_epistasis::init_sequence(), TTQuanti_diallelic_full_pleio_epistasis::init_sequence(), TTQuanti_continuous_no_pleio_epistasis::init_sequence(), TTQuanti_diallelic_no_pleio_epistasis::init_sequence(), TTQuanti_diallelic_bitstring_no_pleio_epistasis::init_sequence(), TTQuanti_diallelic_bitstring_full_pleio_epistasis::init_sequence(), LCE_Disperse_EvolDisp::Migrate_SteppingStone1D(), TTDispersal::mutate(), TTNeutralGenes_bitstring::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(), 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(), 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.

The GSL integer generators (mt19937 used here) return only 32 random bits, so we compose two draws to fill all 64 bits. Without this the upper 32 bits are always 0, which biases the word-level mask recombination.

                                          {
#if defined(HAS_GSL) && !defined(HAS_SPRNG)
    return ((unsigned long)(gsl_rng_get(r) & 0xFFFFFFFFUL) << 32)
         |  (unsigned long)(gsl_rng_get(r) & 0xFFFFFFFFUL);
#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(), TProtoDeletMutations_bitstring::inherit_free(), and RandBool().

Here is the caller graph for this function:

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(), ParamsParser::sample(), 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

  {
    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().

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

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().

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.

                                  {
 
#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(), LCE_Breed_base::random_hermaphrodite(), LCE_Breed_base::RandomMating(), RandULong(), GeneticMap::recombine(), LCE_Resize::regulateAgeClassNoBackup(), LCE_Resize::regulateAgeClassWithBackup(), LCE_Regulation::regulatePatch(), ParamsParser::runif(), Sample(), ParamsParser::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:

• Uniform.h
• simulation.cc