Empirical
ce_random.h
Go to the documentation of this file.
1 // This file is part of Empirical, https://github.com/devosoft/Empirical
2 // Copyright (C) Michigan State University, 2016.
3 // Released under the MIT Software license; see doc/LICENSE
4 //
5 //
6 // A versatile and non-patterned pseudo-random-number generator.
7 //
8 // Status: DESIGN
9 //
10 // Constructor:
11 // Random(int _seed=-1)
12 // _seed is the random number seed that will produce a unique pseudo-random sequence.
13 // (a value of -1 indicates that the seed should be bassed off of a combination of time
14 // and the memory position of the random number generator, in case multiple generators
15 // start at the same time.)
16 //
17 // Other useful functions:
18 // double GetDouble()
19 // double GetDouble(double max)
20 // double GetDouble(double min, double max)
21 // Retrive a random double in the range [min, max). By default, min=0.0 and max=1.0.
22 //
23 // int GetInt(int max)
24 // int GetInt(int min, int max)
25 // uint32_t GetUInt(uint32_t max)
26 // uint32_t GetUInt(uint32_t min, uint32_t max)
27 // Retrive a random int or uint in the range [min, max). By default, min=0.
28 //
29 // bool P(double p)
30 // Tests a random value [0,1) against a given probability p, and returns true of false.
31 //
32 // double GetRandNormal(const double mean, const double std)
33 // uint32_t GetRandPoisson(const double n, double p)
34 // uint32_t GetRandPoisson(const double mean)
35 // uint32_t GetRandBinomial(const double n, const double p)
36 // Draw a value from the given distributions
37 //
38 
39 #ifndef EMP_CE_RANDOM_H
40 #define EMP_CE_RANDOM_H
41 
42 // #include <algorithm>
43 #include <ctime>
44 #include <climits>
45 #include <cmath>
46 #include <iterator>
47 #include <unistd.h>
48 
49 #include "../tools/math.h"
50 
51 namespace emp {
52  class Random {
53  protected:
54  // Internal members
55  int seed;
57  int inext;
58  int inextp;
59  int ma[56] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
60 
61  // Members & functions for stat functions
62  double expRV; // Exponential Random Variable for the randNormal function
63 
64  // Constants ////////////////////////////////////////////////////////////////
65  // Statistical Approximation
66  static const uint32_t _BINOMIAL_TO_NORMAL = 50; // if < n*p*(1-p)
67  static const uint32_t _BINOMIAL_TO_POISSON = 1000; // if < n && !Normal approx Engine
68 
69  // Engine
70  static const uint32_t _RAND_MBIG = 1000000000;
71  static const uint32_t _RAND_MSEED = 161803398;
72 
73  // Internal functions
74 
75  // Setup, called on initialization and seed reset.
76  constexpr void init()
77  {
78  // Clear variables
79  for (int i = 0; i < 56; ++i) ma[i] = 0;
80 
81  int mj = (_RAND_MSEED - seed) % _RAND_MBIG;
82  ma[55] = mj;
83  int mk = 1;
84 
85  for (int i = 1; i < 55; ++i) {
86  int ii = (21 * i) % 55;
87  ma[ii] = mk;
88  mk = mj - mk;
89  if (mk < 0) mk += _RAND_MBIG;
90  mj = ma[ii];
91  }
92 
93  for (int k = 0; k < 4; ++k) {
94  for (int j = 1; j < 55; ++j) {
95  ma[j] -= ma[1 + (j + 30) % 55];
96  if (ma[j] < 0) ma[j] += _RAND_MBIG;
97  }
98  }
99 
100  inext = 0;
101  inextp = 31;
102 
103  // Setup variables used by Statistical Distribution functions
104  expRV = -emp::Ln(Random::Get() / (double) _RAND_MBIG);
105  }
106 
107  // Basic Random number
108  // Returns a random number [0,_RAND_MBIG)
109  constexpr uint32_t Get() {
110  if (++inext == 56) inext = 0;
111  if (++inextp == 56) inextp = 0;
112  int mj = ma[inext] - ma[inextp];
113  if (mj < 0) mj += _RAND_MBIG;
114  ma[inext] = mj;
115 
116  return mj;
117  }
118 
119  public:
126  constexpr Random(const int _seed=-1) : seed(0), original_seed(0), inext(0), inextp(0), expRV(0) {
127  for (int i = 0; i < 56; ++i) ma[i] = 0;
128  ResetSeed(_seed); // Calls init()
129  }
130  ~Random() = default;
131 
132 
136  constexpr int GetSeed() const { return seed; }
137 
141  constexpr int GetOriginalSeed() const { return original_seed; }
142 
150  constexpr void ResetSeed(const int _seed) {
151  original_seed = _seed;
152  seed = _seed;
153  if (seed < 0) seed *= -1;
154  seed %= _RAND_MSEED;
155  init();
156  }
157 
158 
159  // Random Number Generation /////////////////////////////////////////////////
160 
166  constexpr double GetDouble() { return Get() / (double) _RAND_MBIG; }
167 
174  constexpr double GetDouble(const double max) { return GetDouble() * max; }
175 
183  constexpr double GetDouble(const double min, const double max) { return GetDouble() * (max - min) + min; }
184 
191  constexpr uint32_t GetUInt(const uint32_t max) { return static_cast<int>(GetDouble() * static_cast<double>(max)); }
192 
200  constexpr uint32_t GetUInt(const uint32_t min, const uint32_t max) { return GetUInt(max - min) + min; }
201 
209  constexpr int32_t GetInt(const int max) { return static_cast<int>(GetUInt(max)); }
210  constexpr int32_t GetInt(const int min, const int max) { return static_cast<int>(GetUInt(max - min)) + min; }
211 
212 
213  // Random Event Generation //////////////////////////////////////////////////
214 
215  // P(p) => if p < [0,1) random variable
216  constexpr bool P(const double _p) { return (Get() < (_p * _RAND_MBIG));}
217 
218 
219  // Statistical functions ////////////////////////////////////////////////////
220 
221  // Distributions //
222 
226  constexpr double GetRandNormal() {
227  // Draw from a Unit Normal Dist
228  // Using Rejection Method and saving of initial exponential random variable
229  double expRV2 = -emp::Ln(GetDouble());
230  while (1) {
231  expRV -= (expRV2-1)*(expRV2-1)/2;
232  if (expRV > 0) break;
233  expRV = -emp::Ln(GetDouble());
234  expRV2 = -emp::Ln(GetDouble());
235  }
236  if (P(.5)) return expRV2;
237  return -expRV2;
238  }
239 
244  constexpr double GetRandNormal(const double mean, const double std) { return mean + GetRandNormal() * std; }
245 
251  constexpr uint32_t GetRandPoisson(const double mean) {
252  // Draw from a Poisson Dist with mean; if cannot calculate, return UINT_MAX.
253  // Uses Rejection Method
254  const double a = emp::Exp(-mean);
255  if (a <= 0) return UINT_MAX; // cannot calculate, so return UINT_MAX
256  uint32_t k = 0;
257  double u = GetDouble();
258  while (u >= a) {
259  u *= GetDouble();
260  ++k;
261  }
262  return k;
263  }
264 
268  constexpr uint32_t GetRandPoisson(const double n, double p) {
269  // Optimizes for speed and calculability using symetry of the distribution
270  if (p > .5) return (uint32_t) n - GetRandPoisson(n * (1 - p));
271  else return GetRandPoisson(n * p);
272  }
273 
280  constexpr uint32_t GetFullRandBinomial(const double n, const double p) { // Exact
281  // Actually try n Bernoulli events with probability p
282  uint32_t k = 0;
283  for (uint32_t i = 0; i < n; ++i) if (P(p)) k++;
284  return k;
285  }
286 
295  constexpr uint32_t GetRandBinomial(const double n, const double p) { // Approx
296  // Approximate Binomial if appropriate
297  // if np(1-p) is large, use a Normal approx
298  if (n * p * (1 - p) >= _BINOMIAL_TO_NORMAL) {
299  return static_cast<uint32_t>(GetRandNormal(n * p, n * p * (1 - p)) + 0.5);
300  }
301  // elseif n is large, use a Poisson approx
302  if (n >= _BINOMIAL_TO_POISSON) {
303  uint32_t k = GetRandPoisson(n, p);
304  if (k < UINT_MAX) return k; // if approx worked
305  }
306  // otherwise, actually generate the randBinomial
307  return GetFullRandBinomial(n, p);
308  }
309  };
310 
311 
316  typedef int argument_type;
317  typedef int result_type;
318 
319  RandomStdAdaptor(Random& rng) : _rng(rng) { }
320  int operator()(int n) { return _rng.GetInt(n); }
321 
323  };
324 
325 
328  template <typename ForwardIterator, typename OutputIterator, typename RNG>
329  void sample_with_replacement(ForwardIterator first, ForwardIterator last, OutputIterator ofirst, OutputIterator olast, RNG rng) {
330  std::size_t range = std::distance(first, last);
331  while(ofirst != olast) {
332  *ofirst = *(first+rng(range));
333  ++ofirst;
334  }
335  }
336 
337 
338 } // END emp namespace
339 
340 #endif
constexpr int GetSeed() const
Definition: ce_random.h:136
int argument_type
Definition: ce_random.h:316
static const uint32_t _RAND_MSEED
Definition: ce_random.h:71
Random & _rng
Definition: ce_random.h:322
constexpr uint32_t GetUInt(const uint32_t min, const uint32_t max)
Definition: ce_random.h:200
A versatile and non-patterned pseudo-random-number generator (Mersenne Twister).
Definition: ce_random.h:52
Definition: BitVector.h:785
constexpr uint32_t Get()
Definition: ce_random.h:109
int result_type
Definition: ce_random.h:317
static constexpr double Ln(double x)
Compile-time natural log calculator.
Definition: math.h:128
This is an adaptor to make Random behave like a proper STL random number generator.
Definition: ce_random.h:315
int operator()(int n)
Definition: ce_random.h:320
constexpr double GetRandNormal()
Definition: ce_random.h:226
constexpr uint32_t GetUInt(const uint32_t max)
Definition: ce_random.h:191
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: ce_random.h:319
constexpr double GetDouble()
Definition: ce_random.h:166
constexpr int32_t GetInt(const int max)
Definition: ce_random.h:209
int ma[56]
Internal state of RNG.
Definition: ce_random.h:59
constexpr Random(const int _seed=-1)
Definition: ce_random.h:126
static const uint32_t _RAND_MBIG
Definition: ce_random.h:70
constexpr void ResetSeed(const int _seed)
Definition: ce_random.h:150
~Random()=default
constexpr uint32_t GetRandPoisson(const double n, double p)
Definition: ce_random.h:268
constexpr uint32_t GetRandPoisson(const double mean)
Definition: ce_random.h:251
constexpr double GetDouble(const double min, const double max)
Definition: ce_random.h:183
If we are in emscripten, make sure to include the header.
Definition: array.h:37
constexpr uint32_t GetRandBinomial(const double n, const double p)
Definition: ce_random.h:295
int inext
First position in use in internal state.
Definition: ce_random.h:57
static constexpr double Exp(double exp)
A fast method of calculating e^x.
Definition: math.h:181
double expRV
Definition: ce_random.h:62
constexpr double GetDouble(const double max)
Definition: ce_random.h:174
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
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
constexpr int32_t GetInt(const int min, const int max)
Definition: ce_random.h:210
constexpr int GetOriginalSeed() const
Definition: ce_random.h:141
constexpr double GetRandNormal(const double mean, const double std)
Definition: ce_random.h:244
int inextp
Second position in use in internal state.
Definition: ce_random.h:58