Empirical
Random.h
Go to the documentation of this file.
1 
11 #ifndef EMP_RANDOM_H
12 #define EMP_RANDOM_H
13 
14 #include <ctime>
15 #include <climits>
16 #include <cmath>
17 #include <iterator>
18 #include <unistd.h>
19 
20 #include "../base/assert.h"
21 #include "Range.h"
22 
23 namespace emp {
24 
26  class Random {
27  protected:
28  int seed;
29  int original_seed;
30  int inext;
31  int inextp;
32  int ma[56];
33 
34  // Members & functions for stat functions
35  double expRV; // Exponential Random Variable for the randNormal function
36 
37  // Constants ////////////////////////////////////////////////////////////////
38  // Statistical Approximation
39  static const int32_t _BINOMIAL_TO_NORMAL = 50; // if < n*p*(1-p)
40  static const int32_t _BINOMIAL_TO_POISSON = 1000; // if < n && !Normal approx Engine
41 
42  // Engine
43  static const int32_t _RAND_MBIG = 1000000000;
44  static const int32_t _RAND_MSEED = 161803398;
45 
46  // Internal functions
47 
48  // Setup, called on initialization and seed reset.
49  void init()
50  {
51  // Clear variables
52  for (int i = 0; i < 56; ++i) ma[i] = 0;
53 
54  int32_t mj = (_RAND_MSEED - seed) % _RAND_MBIG;
55  ma[55] = mj;
56  int32_t mk = 1;
57 
58  for (int32_t i = 1; i < 55; ++i) {
59  int32_t ii = (21 * i) % 55;
60  ma[ii] = mk;
61  mk = mj - mk;
62  if (mk < 0) mk += _RAND_MBIG;
63  mj = ma[ii];
64  }
65 
66  for (int32_t k = 0; k < 4; ++k) {
67  for (int32_t j = 1; j < 55; ++j) {
68  ma[j] -= ma[1 + (j + 30) % 55];
69  if (ma[j] < 0) ma[j] += _RAND_MBIG;
70  }
71  }
72 
73  inext = 0;
74  inextp = 31;
75 
76  // Setup variables used by Statistical Distribution functions
77  expRV = -log(Random::Get() / (double) _RAND_MBIG);
78  }
79 
80  // Basic Random number
81  // Returns a random number [0,_RAND_MBIG)
82  int32_t Get() {
83  if (++inext == 56) inext = 0;
84  if (++inextp == 56) inextp = 0;
85  int mj = ma[inext] - ma[inextp];
86  if (mj < 0) mj += _RAND_MBIG;
87  ma[inext] = mj;
88 
89  return mj;
90  }
91 
92  public:
99  Random(const int _seed = -1) : seed(0), original_seed(0), inext(0), inextp(0), expRV(0) {
100  for (int i = 0; i < 56; ++i) ma[i] = 0;
101  ResetSeed(_seed); // Calls init()
102  }
103 
104  ~Random() { ; }
105 
106 
110  inline int GetSeed() const { return seed; }
111 
115  inline int GetOriginalSeed() const { return original_seed; }
116 
124  inline void ResetSeed(const int _seed) {
125  original_seed = _seed;
126 
127  if (_seed <= 0) {
128  int seed_time = (int) time(NULL);
129  int seed_mem = (int) ((uint64_t) this);
130  seed = seed_time ^ seed_mem;
131  } else {
132  seed = _seed;
133  }
134 
135  if (seed < 0) seed *= -1;
136  seed %= _RAND_MSEED;
137 
138  init();
139  }
140 
141 
142  // Random Number Generation /////////////////////////////////////////////////
143 
149  inline double GetDouble() { return Get() / (double) _RAND_MBIG; }
150 
157  inline double GetDouble(const double max) {
158  // emp_assert(max <= (double) _RAND_MBIG, max, (double) _RAND_MBIG); // Precision will be too low past this point...
159  return GetDouble() * max;
160  }
161 
169  inline double GetDouble(const double min, const double max) {
170  emp_assert((max-min) <= (double) _RAND_MBIG, min, max); // Precision will be too low past this point...
171  return GetDouble() * (max - min) + min;
172  }
173 
180  inline double GetDouble(const Range<double> range) {
181  return GetDouble(range.GetLower(), range.GetUpper());
182  }
183 
190  template <typename T>
191  inline uint32_t GetUInt(const T max) {
192  emp_assert(max <= (T) _RAND_MBIG, max); // Precision will be too low past this point...
193  return static_cast<uint32_t>(GetDouble() * static_cast<double>(max));
194  }
195 
203  template <typename T>
204  inline uint64_t GetUInt64(const T max) {
205  if (max <= (T) _RAND_MBIG) return (uint64_t) GetUInt(max); // Don't need extra precision.
206  const double max2 = ((double) max) / (double) _RAND_MBIG;
207  emp_assert(max2 <= (T) _RAND_MBIG, max); // Precision will be too low past this point...
208 
209  return static_cast<uint64_t>(GetDouble() * static_cast<double>(max))
210  + static_cast<uint64_t>(GetDouble() * static_cast<double>(max2) * _RAND_MBIG);
211  }
212 
213 
221  template <typename T1, typename T2>
222  inline uint32_t GetUInt(const T1 min, const T2 max) {
223  return GetUInt<uint32_t>((uint32_t) max - (uint32_t) min) + (uint32_t) min;
224  }
225 
232  template <typename T>
233  inline uint32_t GetUInt(const Range<T> range) {
234  return GetUInt(range.GetLower(), range.GetUpper());
235  }
236 
244  inline int GetInt(const int max) { return static_cast<int>(GetUInt((uint32_t) max)); }
245  inline int GetInt(const int min, const int max) { return GetInt(max - min) + min; }
246  inline int GetInt(const Range<int> range) { return GetInt(range.GetLower(), range.GetUpper()); }
247 
248 
249  // Random Event Generation //////////////////////////////////////////////////
250 
253  inline bool P(const double p) {
254  emp_assert(p >= 0.0 && p <= 1.0, p);
255  return (Get() < (p * _RAND_MBIG));
256  }
257 
258 
259  // Statistical functions ////////////////////////////////////////////////////
260 
261  // Distributions //
262 
266  inline double GetRandNormal() {
267  // Draw from a Unit Normal Dist
268  // Using Rejection Method and saving of initial exponential random variable
269  double expRV2;
270  while (1) {
271  expRV2 = -log(GetDouble());
272  expRV -= (expRV2-1)*(expRV2-1)/2;
273  if (expRV > 0) break;
274  expRV = -log(GetDouble());
275  }
276  if (P(.5)) return expRV2;
277  return -expRV2;
278  }
279 
284  inline double GetRandNormal(const double mean, const double std) { return mean + GetRandNormal() * std; }
285 
289  inline uint32_t GetRandPoisson(const double n, double p) {
290  emp_assert(p >= 0.0 && p <= 1.0, p);
291  // Optimizes for speed and calculability using symetry of the distribution
292  if (p > .5) return (uint32_t)n - GetRandPoisson(n * (1 - p));
293  else return GetRandPoisson(n * p);
294  }
295 
301  inline uint32_t GetRandPoisson(const double mean) {
302  // Draw from a Poisson Dist with mean; if cannot calculate, return UINT_MAX.
303  // Uses Rejection Method
304  const double a = exp(-mean);
305  if (a <= 0) return UINT_MAX; // cannot calculate, so return UINT_MAX
306  uint32_t k = 0;
307  double u = GetDouble();
308  while (u >= a) {
309  u *= GetDouble();
310  ++k;
311  }
312  return k;
313  }
314 
321  inline uint32_t GetFullRandBinomial(const double n, const double p) { // Exact
322  emp_assert(p >= 0.0 && p <= 1.0, p);
323  // Actually try n Bernoulli events with probability p
324  uint32_t k = 0;
325  for (uint32_t i = 0; i < n; ++i) if (P(p)) k++;
326  return k;
327  }
328 
337  inline uint32_t GetRandBinomial(const double n, const double p) { // Approx
338  emp_assert(p >= 0.0 && p <= 1.0, p);
339  emp_assert(n >= 0.0, n);
340  // Approximate Binomial if appropriate
341  // if np(1-p) is large, use a Normal approx
342  if (n * p * (1 - p) >= _BINOMIAL_TO_NORMAL) {
343  return static_cast<uint32_t>(GetRandNormal(n * p, n * p * (1 - p)) + 0.5);
344  }
345  // elseif n is large, use a Poisson approx
346  if (n >= _BINOMIAL_TO_POISSON) {
347  uint32_t k = GetRandPoisson(n, p);
348  if (k < UINT_MAX) return k; // if approx worked
349  }
350  // otherwise, actually generate the randBinomial
351  return GetFullRandBinomial(n, p);
352  }
353  };
354 
355 
357  struct RandomStdAdaptor {
358  typedef int argument_type;
359  typedef int result_type;
360 
361  RandomStdAdaptor(Random& rng) : _rng(rng) { }
362  int operator()(int n) { return _rng.GetInt(n); }
363 
364  Random& _rng;
365  };
366 
367 
369  template <typename ForwardIterator, typename OutputIterator, typename RNG>
370  void sample_with_replacement(ForwardIterator first, ForwardIterator last, OutputIterator ofirst, OutputIterator olast, RNG rng) {
371  std::size_t range = std::distance(first, last);
372  while(ofirst != olast) {
373  *ofirst = *(first+rng(range));
374  ++ofirst;
375  }
376  }
377 
378 
379 } // END emp namespace
380 
381 #endif
int GetOriginalSeed() const
Definition: Random.h:115
int argument_type
Definition: Random.h:358
int GetInt(const Range< int > range)
Definition: Random.h:246
double GetRandNormal(const double mean, const double std)
Definition: Random.h:284
static const uint32_t _RAND_MSEED
Definition: ce_random.h:71
int GetInt(const int max)
Definition: Random.h:244
uint32_t GetUInt(const T max)
Definition: Random.h:191
Random(const int _seed=-1)
Definition: Random.h:99
A versatile and non-patterned pseudo-random-number generator (Mersenne Twister).
Definition: ce_random.h:52
Definition: BitVector.h:785
void init()
Definition: Random.h:49
constexpr uint32_t Get()
Definition: ce_random.h:109
double GetDouble()
Definition: Random.h:149
A range of values from a lower limit to and upper limit, of any provided type.
Definition: Range.h:23
int result_type
Definition: Random.h:359
T GetLower() const
Definition: Range.h:32
This is an adaptor to make Random behave like a proper STL random number generator.
Definition: ce_random.h:315
~Random()
Definition: Random.h:104
uint32_t GetFullRandBinomial(const double n, const double p)
Definition: Random.h:321
int operator()(int n)
Definition: Random.h:362
constexpr double GetRandNormal()
Definition: ce_random.h:226
constexpr uint32_t GetUInt(const uint32_t max)
Definition: ce_random.h:191
A simple way to track value ranges.
int seed
Current random number seed.
Definition: ce_random.h:55
constexpr void init()
Definition: ce_random.h:76
static const uint32_t _BINOMIAL_TO_POISSON
Definition: ce_random.h:67
static const uint32_t _BINOMIAL_TO_NORMAL
Definition: ce_random.h:66
RandomStdAdaptor(Random &rng)
Definition: Random.h:361
constexpr double GetDouble()
Definition: ce_random.h:166
uint32_t GetRandPoisson(const double mean)
Definition: Random.h:301
int GetInt(const int min, const int max)
Definition: Random.h:245
double GetDouble(const double max)
Definition: Random.h:157
int ma[56]
Internal state of RNG.
Definition: ce_random.h:59
double GetRandNormal()
Definition: Random.h:266
static const uint32_t _RAND_MBIG
Definition: ce_random.h:70
constexpr void ResetSeed(const int _seed)
Definition: ce_random.h:150
int GetSeed() const
Definition: Random.h:110
void ResetSeed(const int _seed)
Definition: Random.h:124
uint32_t GetRandBinomial(const double n, const double p)
Definition: Random.h:337
constexpr uint32_t GetRandPoisson(const double mean)
Definition: ce_random.h:251
uint32_t GetRandPoisson(const double n, double p)
Definition: Random.h:289
int32_t Get()
Definition: Random.h:82
uint32_t GetUInt(const Range< T > range)
Definition: Random.h:233
If we are in emscripten, make sure to include the header.
Definition: array.h:37
double GetDouble(const Range< double > range)
Definition: Random.h:180
int inext
First position in use in internal state.
Definition: ce_random.h:57
#define emp_assert(...)
Definition: assert.h:199
double expRV
Definition: ce_random.h:62
double GetDouble(const double min, const double max)
Definition: Random.h:169
uint64_t GetUInt64(const T max)
Definition: Random.h:204
T GetUpper() const
Definition: Range.h:33
constexpr bool P(const double _p)
Definition: ce_random.h:216
int original_seed
Orignal random number seed when object was first created.
Definition: ce_random.h:56
constexpr uint32_t GetFullRandBinomial(const double n, const double p)
Definition: ce_random.h:280
bool P(const double p)
Definition: Random.h:253
uint32_t GetUInt(const T1 min, const T2 max)
Definition: Random.h:222
void sample_with_replacement(ForwardIterator first, ForwardIterator last, OutputIterator ofirst, OutputIterator olast, RNG rng)
Draw a sample (with replacement) from an input range, copying to the output range.
Definition: ce_random.h:329
int inextp
Second position in use in internal state.
Definition: ce_random.h:58