OpenFPM_pdata  4.1.0 Project that contain the implementation of distributed structures
SimpleRNG.hpp
1
12 #ifndef SIMPLE_RNG_HPP
13 #define SIMPLE_RNG_HPP
14
15 #include <string>
16 #include <cmath>
17
18 class SimpleRNG
19 {
20 private:
21  unsigned int m_w;
22  unsigned int m_z;
23
24  // This is the heart of the generator.
25  // It uses George Marsaglia's MWC algorithm to produce an unsigned integer.
27  unsigned int GetUint()
28  {
29  m_z = 36969 * (m_z & 65535) + (m_z >> 16);
30  m_w = 18000 * (m_w & 65535) + (m_w >> 16);
31  return (m_z << 16) + m_w;
32  }
33
34 public:
35
36  SimpleRNG()
37  {
38  // These values are not magical, just the default values Marsaglia used.
39  // Any pair of unsigned integers should be fine.
40  m_w = 521288629;
41  m_z = 362436069;
42  }
43
44  // The random generator seed can be set three ways:
45  // 1) specifying two non-zero unsigned integers
46  // 2) specifying one non-zero unsigned integer and taking a default value for the second
47  // 3) setting the seed from the system time
48
49  void SetSeed(unsigned int u, unsigned int v)
50  {
51  if (u != 0) m_w = u;
52  if (v != 0) m_z = v;
53  }
54
55  void SetSeed(unsigned int u)
56  {
57  m_w = u;
58  }
59
60  void SetSeedFromSystemTime()
61  {
62  long x = clock();
63  SetSeed((unsigned int)(x >> 16), (unsigned int)(x % 4294967296));
64  }
65
66  // Produce a uniform random sample from the open interval (0, 1).
67  // The method will not return either end point.
68  double GetUniform()
69  {
70  // 0 <= u < 2^32
71  unsigned int u = GetUint();
72  // The magic number below is 1/(2^32 + 2).
73  // The result is strictly between 0 and 1.
74  return (u + 1.0) * 2.328306435454494e-10;
75  }
76
77  // Get normal (Gaussian) random sample with mean 0 and standard deviation 1
78  double GetNormal()
79  {
80  // Use Box-Muller algorithm
81  double u1 = GetUniform();
82  double u2 = GetUniform();
83  double r = std::sqrt( -2.0*std::log(u1) );
84  double theta = 2.0*3.14159265359*u2;
85  return r*std::sin(theta);
86  }
87
88  // Get normal (Gaussian) random sample with specified mean and standard deviation
89  double GetNormal(double mean, double standardDeviation)
90  {
91  if (standardDeviation <= 0.0)
92  {
93  std::cerr << __FILE__ << ":" << __LINE__ << " Shape must be positive. Received." << std::to_string(standardDeviation);
94  throw 100001;
95  }
96  return mean + standardDeviation*GetNormal();
97  }
98
99  // Get exponential random sample with mean 1
100  double GetExponential()
101  {
102  return -std::log( GetUniform() );
103  }
104
105  // Get exponential random sample with specified mean
106  double GetExponential(double mean)
107  {
108  if (mean <= 0.0)
109  {
110  std::cerr << __FILE__ << ":" << __LINE__ << "Mean must be positive. Received " << mean;
111  throw 100002;
112  }
113  return mean*GetExponential();
114  }
115
116  double GetGamma(double shape, double scale)
117  {
118  // Implementation based on "A Simple Method for Generating Gamma Variables"
119  // by George Marsaglia and Wai Wan Tsang. ACM Transactions on Mathematical Software
120  // Vol 26, No 3, September 2000, pages 363-372.
121
122  double d, c, x, xsquared, v, u;
123
124  if (shape >= 1.0)
125  {
126  d = shape - 1.0/3.0;
127  c = 1.0/std::sqrt(9.0*d);
128  for (;;)
129  {
130  do
131  {
132  x = GetNormal();
133  v = 1.0 + c*x;
134  }
135  while (v <= 0.0);
136  v = v*v*v;
137  u = GetUniform();
138  xsquared = x*x;
139  if (u < 1.0 -.0331*xsquared*xsquared || std::log(u) < 0.5*xsquared + d*(1.0 - v + std::log(v)))
140  return scale*d*v;
141  }
142  }
143  else if (shape <= 0.0)
144  {
145  std::cerr << __FILE__ << ":" << __LINE__ << " Shape must be positive. Received" << shape << "\n";
146  throw 100003;
147  }
148  else
149  {
150  double g = GetGamma(shape+1.0, 1.0);
151  double w = GetUniform();
152  return scale*g*std::pow(w, 1.0/shape);
153  }
154  }
155
156  double GetChiSquare(double degreesOfFreedom)
157  {
158  // A chi squared distribution with n degrees of freedom
159  // is a gamma distribution with shape n/2 and scale 2.
160  return GetGamma(0.5 * degreesOfFreedom, 2.0);
161  }
162
163  double GetInverseGamma(double shape, double scale)
164  {
165  // If X is gamma(shape, scale) then
166  // 1/Y is inverse gamma(shape, 1/scale)
167  return 1.0 / GetGamma(shape, 1.0 / scale);
168  }
169
170  double GetWeibull(double shape, double scale)
171  {
172  if (shape <= 0.0 || scale <= 0.0)
173  {
174  std::cerr << __FILE__ << ":" << __LINE__ << " Shape and scale parameters must be positive. Recieved " << shape << " and " << scale;
175  throw 100004;
176  }
177  return scale * std::pow(-std::log(GetUniform()), 1.0 / shape);
178  }
179
180  double GetCauchy(double median, double scale)
181  {
182  if (scale <= 0)
183  {
184  std::cerr << __FILE__ << ":" << __LINE__ << "Scale must be positive. Received " << scale << "\n";
185  throw 100005;
186  }
187
188  double p = GetUniform();
189
190  // Apply inverse of the Cauchy distribution function to a uniform
191  return median + scale*std::tan(3.14159265359*(p - 0.5));
192  }
193
194  double GetStudentT(double degreesOfFreedom)
195  {
196  if (degreesOfFreedom <= 0)
197  {
198  std::cerr << __FILE__ << ":" << __LINE__ << " Degrees of freedom must be positive. Received " << degreesOfFreedom;
199  throw 100006;
200  }
201
202  // See Seminumerical Algorithms by Knuth
203  double y1 = GetNormal();
204  double y2 = GetChiSquare(degreesOfFreedom);
205  return y1 / std::sqrt(y2 / degreesOfFreedom);
206  }
207
208  // The Laplace distribution is also known as the double exponential distribution.
209  double GetLaplace(double mean, double scale)
210  {
211  double u = GetUniform();
212  return (u < 0.5) ?
213  mean + scale*std::log(2.0*u) :
214  mean - scale*std::log(2*(1-u));
215  }
216
217  double GetLogNormal(double mu, double sigma)
218  {
219  return std::exp(GetNormal(mu, sigma));
220  }
221
222  double GetBeta(double a, double b)
223  {
224  if (a <= 0.0 || b <= 0.0)
225  {
226  std::cerr << __FILE__ << ":" << __LINE__ << "Beta parameters must be positive. Received " << a << " and " << b << "\n";
227  throw 100007;
228  }
229
230  // There are more efficient methods for generating beta samples.
231  // However such methods are a little more efficient and much more complicated.
232  // For an explanation of why the following method works, see
233  // http://www.johndcook.com/distribution_chart.html#gamma_beta
234
235  double u = GetGamma(a, 1.0);
236  double v = GetGamma(b, 1.0);
237  return u / (u + v);
238  }
239 };
240
241 #endif
SimpleRNG is a simple random number generator based on George Marsaglia's MWC (multiply with carry) g...
Definition: SimpleRNG.hpp:18