MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
RNG.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "RNG.h"
27 #include "Helpers.h"
28 #include <limits>
29 
37 {
39  a_ = 1103515245;
40  c_ = 12345;
41  m_ = 1024 * 1024 * 1024;
43  p_ = 607;
44  q_ = 273;
47 
48  haveSavedBoxMuller_ = false;
49  savedBoxMuller_ = 0;
50 
51 }
52 
53 void RNG::setRandomSeed(unsigned long int new_seed)
54 {
57 }
58 
59 void RNG::read(std::istream& is)
60 {
61  std::string dummy;
62  unsigned int type;
63  is >> type;
64  type_ = static_cast<RNGType>(type);
65  is >> a_;
66  is >> c_;
67  is >> m_;
68  is >> p_;
69  is >> q_;
71  //note: the seeds for the LaggedFibonacciGenerator cannot be restarted currently.
73  //randomSeedLaggedFibonacciGenerator_.resize(p_);
74  //for (auto& v : randomSeedLaggedFibonacciGenerator_)
75  // is >> v;
76 }
77 
78 void RNG::write(std::ostream& os) const
79 {
80  os << " " << static_cast<unsigned int>(type_);
81  os << " " << a_;
82  os << " " << c_;
83  os << " " << m_;
84  os << " " << p_;
85  os << " " << q_;
87  //for (auto v : randomSeedLaggedFibonacciGenerator_)
88  // os << " " << v;
89 }
90 
91 void RNG::setLinearCongruentialGeneratorParmeters(unsigned const int a, unsigned const int c, unsigned const int m)
92 {
93  a_ = a;
94  c_ = c;
95  m_ = m;
96 }
97 
99 {
100 #ifdef MERCURY_USE_MPI
101  //First set a random seed on the root
102  if (PROCESSOR_ID == 0)
103  {
104  setRandomSeed(static_cast<unsigned long int>(time(nullptr)));
105  }
106 
107  //Communicate this to the rest of the processes
108  std::vector<int> values(7);
109  if (PROCESSOR_ID == 0)
110  {
111  values[0] = static_cast<unsigned int>(type_);
112  values[1] = a_;
113  values[2] = c_;
114  values[3] = m_;
115  values[4] = p_;
116  values[5] = q_;
118  }
119  MPIContainer::Instance().broadcast(values.data(),7,0);
120 
121  //Update the generators on the other processors
122  if (PROCESSOR_ID != 0)
123  {
124  type_ = static_cast<RNGType>(values[0]);
125  a_ = values[1];
126  c_ = values[2];
127  m_ = values[3];
128  p_ = values[4];
129  q_ = values[5];
131  }
133 #else
134  setRandomSeed(static_cast<unsigned long int>(time(nullptr)));
135 #endif
136 }
137 
139 {
140  type_ = type;
141 }
142 
144 {
145  return getRandomNumber(0, 1);
146 }
147 
149 {
150  logger.assert(min <= max, "getRandomNumber: min cannot be larger than max");
153  } else {
155  }
156 }
157 
158 /* \details This uses the Box--Muller transform
159  * https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
160  * The Box--Muller algorithm generates a pair of independently distributed
161  * normal random variables. One of these will be returned straight away and the
162  * other will be stored and returned the next time you call this function.
163  */
165 {
166  static const double epsilon = std::numeric_limits<Mdouble>::min();
167 
169  {
170  /* If we have already generated a normal variate, use it. */
171  haveSavedBoxMuller_ = false;
172  return savedBoxMuller_;
173  }
174  else
175  {
176  /* Otherwise, generate a pair of normal variates, return one of them,
177  * and save the other. */
178  Mdouble radius, theta;
179  do
180  {
181  radius = getRandomNumber(0, 1);
182  theta = getRandomNumber(0, 2 * constants::pi);
183  } while (radius <= epsilon);
184  // make sure that the radius generated is not too small
185  // (unlikely to happen, just a safety check)
186 
187  savedBoxMuller_ = sqrt(-2.0 * log(radius)) * sin(theta);
188  haveSavedBoxMuller_ = true;
189  return sqrt(-2.0 * log(radius)) * cos(theta);
190  }
191 }
192 
194 {
195  if (stdev == 0) {
196  logger(WARN,
197  "[RNG::getNormalVariate(Mdouble, Mdouble)] Zero stdev?");
198  return mean;
199  } else if (stdev < 0) {
200  logger(ERROR,
201  "[RNG::getNormalVariate(Mdouble, Mdouble)] Negative stdev is not allowed.");
202  exit(-1);
203  } else {
204  return getNormalVariate() * stdev + mean;
205  }
206 }
207 
212 unsigned int RNG::getPoissonVariate(Mdouble lambda)
213 {
214  if (lambda > 50)
215  {
216  logger(WARN, "[RNG::getPoissonVariate(Mdouble)] Knuth's algorithm for Poissons may be slow for lambda = %", lambda);
217  }
218  unsigned int k = 0;
219  Mdouble p = 1;
220  Mdouble u;
221  do
222  {
223  k++;
224  u = getRandomNumber(0, 1);
225  p *= u;
226  }
227  while (u > exp(-lambda));
228  return k-1;
229 }
230 
235 #pragma optimize( "", off )
237 {
238  //Update the random seed
240 
241  //Generate a random number in the required range
242 
243  Mdouble range = max - min;
244  Mdouble random_num = min + range * randomSeedLinearCongruentialGenerator_ / (static_cast<Mdouble>(m_) + 1.0);
245 
246  return random_num;
247 }
248 #pragma optimize( "", on )
249 
250 /**************************************
251  * This sets the seed for LFG using LCG
252  *
253  *
254  **************************************/
256 {
257  for (unsigned int i = 0; i < p_; i++)
258  {
260  }
261 }
262 
268 {
269 #pragma optimize( "", off )
271  static_cast<Mdouble>(1.0));
272  //Update the random seed
274  randomSeedLaggedFibonacciGenerator_.emplace_back(new_seed);
275 
276  //Generate a random number in the required range
277 
278  Mdouble random_num;
279 
280  Mdouble range = max - min;
281  random_num = min + range * new_seed;
282  return random_num;
283 #pragma optimize( "", on )
284 }
285 
292 {
293  //This are the fixed parameters that define the test
294  static unsigned int num_of_tests = 100000;
295  static Mdouble max_num = 100.0;
296  static unsigned int num_of_bins = 10;
297 
298  //This is the generated random_number
299  Mdouble rn;
300  //This is the bin the random number will lie in
301  unsigned int bin = 0;
302  //This is a vector of bins
303  std::vector<int> count;
304  count.resize(num_of_bins);
305 
306  //Initialisation of the bins
307  for (unsigned int i = 0; i < num_of_bins; i++)
308  {
309  count[bin] = 0;
310  }
311 
312  //Loop over a number of tests
313  for (unsigned int i = 0; i < num_of_tests; i++)
314  {
315  rn = getRandomNumber(0.0, max_num);
316  bin = static_cast<unsigned int>(std::floor(rn * num_of_bins / max_num));
317 
318  //Add one to the bin count
319  count[bin]++;
320 
321  }
322 
323  //Final post-process the result and report on the random number
324  Mdouble chi_cum = 0.0;
325  Mdouble expected = num_of_tests / num_of_bins;
326 
327  for (unsigned int i = 0; i < num_of_bins; i++)
328  {
329  chi_cum = chi_cum + (count[i] - expected) * (count[i] - expected) / expected;
330  std::cout << i << " : "
331  << count[i] << " : " << (count[i] - expected) * (count[i] - expected) / expected << std::endl;
332  }
333  //end for loop over computing the chi-squared value.
334  std::cout << "chi_cum " << chi_cum << std::endl;
335 
336  return mathsFunc::chi_squared_prob(chi_cum, num_of_bins);
337 }
338 
340 void RNG::setLaggedFibonacciGeneratorParameters(const unsigned int p, const unsigned int q)
341 {
342  //p must be greater than q so makes sure this is true. Not sure what happens if you set p=q, in the LFG alogrithm.
343  if (p < q)
344  {
345  p_ = q;
346  q_ = p;
347  }
348 
351 }
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:130
void setLinearCongruentialGeneratorParmeters(const unsigned int a, const unsigned int c, unsigned int m)
This functions set the parameters for the LCG random number generator. It goes multiplier, addition, mod.
Definition: RNG.cc:91
void read(std::istream &is)
Definition: RNG.cc:59
unsigned long int c_
Definition: RNG.h:164
Mdouble getRandomNumberFromLaggedFibonacciGenerator(Mdouble min, Mdouble max)
This is a Lagged Fibonacci Generator.
Definition: RNG.cc:267
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
void seedLaggedFibonacciGenerator()
This seed the LFG.
Definition: RNG.cc:255
double Mdouble
Definition: GeneralDefine.h:34
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
unsigned long int q_
Definition: RNG.h:169
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void setLaggedFibonacciGeneratorParameters(const unsigned int p, const unsigned int q)
This function sets the parameters for the LFG random number generator.
Definition: RNG.cc:340
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Definition: RNG.cc:143
unsigned int getPoissonVariate(Mdouble lambda)
Produces a random number according to a Poisson distribution.
Definition: RNG.cc:212
Mdouble log(Mdouble Power)
Mdouble chi_squared_prob(Mdouble x, unsigned int k)
This is the function which actually gives the probability back using a chi squared test...
void randomise()
sets the random variables such that they differ for each run
Definition: RNG.cc:98
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
void write(std::ostream &os) const
Definition: RNG.cc:78
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble savedBoxMuller_
A storage space for the so-far-unused variate from the pair generated by Box–Muller.
Definition: RNG.h:190
Mdouble getNormalVariate()
Produces a random number according to a normal distribution with mean 0 and standard deviation 1...
Definition: RNG.cc:164
Mdouble getRandomNumberFromLinearCongruentialGenerator(Mdouble min, Mdouble max)
This is a basic Linear Congruential Generator Random.
Definition: RNG.cc:236
#define PROCESSOR_ID
Definition: GeneralDefine.h:63
unsigned long int randomSeedLinearCongruentialGenerator_
This is the initial seed of the RNG.
Definition: RNG.h:154
unsigned long int p_
This are the parameters that control the LFG random generator.
Definition: RNG.h:169
void setRandomSeed(unsigned long int new_seed)
This is the seed for the random number generator (note the call to seed_LFG is only required really i...
Definition: RNG.cc:53
std::enable_if< std::is_scalar< T >::value, void >::type broadcast(T &t, int fromProcessor=0)
Broadcasts a scalar from the root to all other processors.
Definition: MpiContainer.h:429
Mdouble test()
This function tests the quality of random numbers, based on the chi-squared test. ...
Definition: RNG.cc:291
bool haveSavedBoxMuller_
A flag that keeps track of whether or not to generate a new pair of normal variates (using Box–Mulle...
Definition: RNG.h:185
std::vector< Mdouble > randomSeedLaggedFibonacciGenerator_
This is the seeds required for the LFG.
Definition: RNG.h:159
unsigned long int a_
This are the two parameters that control the LCG random generated.
Definition: RNG.h:164
void setRandomNumberGenerator(RNGType type)
Allows the user to set which random number generator is used.
Definition: RNG.cc:138
RNGType
Definition: RNG.h:38
RNGType type_
This is the type of random number generator.
Definition: RNG.h:174
RNG()
default constructor
Definition: RNG.cc:36
unsigned long int m_
Definition: RNG.h:164