GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GRan.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GRan.cpp - Random number generator class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2014 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GRan.cpp
23  * @brief Random number generator class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GRan.hpp"
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GException.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 
39 /* __ Macros _____________________________________________________________ */
40 
41 /* __ Coding definitions _________________________________________________ */
42 
43 /* __ Debug definitions __________________________________________________ */
44 
45 /* __ Constants __________________________________________________________ */
46 
47 
48 /*==========================================================================
49  = =
50  = Constructors/destructors =
51  = =
52  ==========================================================================*/
53 
54 /***********************************************************************//**
55  * @brief Void constructor
56  ***************************************************************************/
58 {
59  // Initialise private members
60  init_members();
61 
62  // Return
63  return;
64 }
65 
66 
67 /***********************************************************************//**
68  * @brief Seed constructor
69  *
70  * @param[in] seed Random number generator seed.
71  ***************************************************************************/
72 GRan::GRan(unsigned long long int seed)
73 {
74  // Initialise private
75  init_members(seed);
76 
77  // Return
78  return;
79 }
80 
81 
82 /***********************************************************************//**
83  * @brief Copy constructor
84  *
85  * @param[in] ran Random number generator.
86  ***************************************************************************/
87 GRan::GRan(const GRan& ran)
88 {
89  // Initialise private
90  init_members();
91 
92  // Copy members
93  copy_members(ran);
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief Destructor
102  ***************************************************************************/
104 {
105  // Free members
106  free_members();
107 
108  // Return
109  return;
110 }
111 
112 
113 /*==========================================================================
114  = =
115  = Operators =
116  = =
117  ==========================================================================*/
118 
119 /***********************************************************************//**
120  * @brief Assignment operator
121  *
122  * @param[in] ran Random number generator.
123  * @return Random number generator.
124  ***************************************************************************/
126 {
127  // Execute only if object is not identical
128  if (this != &ran) {
129 
130  // Free members
131  free_members();
132 
133  // Initialise private members for clean destruction
134  init_members();
135 
136  // Copy members
137  copy_members(ran);
138 
139  } // endif: object was not identical
140 
141  // Return
142  return *this;
143 }
144 
145 
146 /*==========================================================================
147  = =
148  = Public methods =
149  = =
150  ==========================================================================*/
151 
152 /***********************************************************************//**
153  * @brief Clear random number generator
154  *
155  * This method properly resets the object to an initial state using the
156  * default seed.
157  ***************************************************************************/
158 void GRan::clear(void)
159 {
160  // Free class members
161  free_members();
162 
163  // Initialise members
164  init_members();
165 
166  // Return
167  return;
168 }
169 
170 
171 /***********************************************************************//**
172  * @brief Clone random number generator
173  *
174  * @return Pointer to deep copy of random number generator.
175  **************************************************************************/
176 GRan* GRan::clone(void) const
177 {
178  return new GRan(*this);
179 }
180 
181 
182 /***********************************************************************//**
183  * @brief Set seed of random number generator
184  *
185  * @param[in] seed Random number generator seed.
186  *
187  * This method properly resets the object to an initial state using the
188  * specified seed.
189  ***************************************************************************/
190 void GRan::seed(unsigned long long int seed)
191 {
192  // Free class members
193  free_members();
194 
195  // Initialise members
196  init_members(seed);
197 
198  // Return
199  return;
200 }
201 
202 
203 /***********************************************************************//**
204  * @brief Return 32-bit random unsigned integer
205  ***************************************************************************/
206 unsigned long int GRan::int32(void)
207 {
208  // Return random value
209  return ((unsigned long int)int64());
210 }
211 
212 
213 /***********************************************************************//**
214  * @brief Return 64-bit random unsigned integer
215  ***************************************************************************/
216 unsigned long long int GRan::int64(void)
217 {
218  // Update u
219  m_value1 = m_value1 * 2862933555777941757LL + 7046029254386353087LL;
220 
221  // Shuffle v
222  m_value2 ^= m_value2 >> 17;
223  m_value2 ^= m_value2 << 31;
224  m_value2 ^= m_value2 >> 8;
225 
226  // Update w
227  m_value3 = 4294957665U * (m_value3 & 0xffffffff) + (m_value3 >> 32);
228 
229  // Compute x
230  unsigned long long int x = m_value1 ^ (m_value1 << 21);
231  x ^= x >> 35;
232  x ^= x << 4;
233 
234  // Return random value
235  return ((x + m_value2) ^ m_value3);
236 }
237 
238 
239 /***********************************************************************//**
240  * @brief Returns random double precision floating value in range 0 to 1
241  ***************************************************************************/
242 double GRan::uniform(void)
243 {
244  // Return random value
245  return (5.42101086242752217e-20 * int64());
246 }
247 
248 
249 /***********************************************************************//**
250  * @brief Returns normal deviates
251  *
252  * Method Used: Box-Muller transform, outlined here:
253  * http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
254  *
255  * Code from: http://www.design.caltech.edu/erik/Misc/Gaussian.html
256  ***************************************************************************/
257 double GRan::normal(void)
258 {
259  // Get random value
260  double x1;
261  double w;
262  do {
263  x1 = 2.0 * uniform() - 1.0;
264  double x2 = 2.0 * uniform() - 1.0;
265  w = x1 * x1 + x2 * x2;
266  } while (w >= 1.0);
267 
268  // Compute random value
269  double value = x1 * std::sqrt((-2.0 * std::log(w)) / w);
270 
271  // Return random value
272  return value;
273 }
274 
275 
276 /***********************************************************************//**
277  * @brief Returns exponential deviates
278  *
279  * @param[in] lambda Mean rate.
280  *
281  * Returns exponential deviates from the probability distribution
282  * \f[p(x) = \lambda \exp( -\lambda x )\f]
283  * where
284  * \f$\lambda\f$ is the parameter of the distribution.
285  * This method may be used to simulate the occurence time of an event, where
286  * \f$\lambda\f$ is the mean event rate. Convsersely, \f$1/\lambda\f$ is the
287  * mean waiting time between events.
288  *
289  * @todo Check that \f$\lambda>0\f$.
290  ***************************************************************************/
291 double GRan::exp(const double& lambda)
292 {
293  // Allocate argument
294  double x;
295 
296  // Get non-zero uniform deviate
297  do {
298  x = uniform();
299  } while (x == 0.0);
300 
301  // Return random value
302  return (-std::log(x)/lambda);
303 }
304 
305 
306 /***********************************************************************//**
307  * @brief Returns Poisson deviates
308  *
309  * @param[in] lambda Expectation value.
310  *
311  * Returns Poisson deviates for an expectation value.
312  *
313  * This method may be used to simulate the number of events in case that a
314  * given mean number of events is expected.
315  ***************************************************************************/
316 double GRan::poisson(const double& lambda)
317 {
318  // Declare result
319  double value;
320 
321  // Use direct method for small numbers ...
322  if (lambda < 12.0) {
323  if (lambda != m_old_lambda) {
324  m_old_lambda = lambda;
325  m_exp_lambda = std::exp(-lambda);
326  }
327  value = -1.0;
328  double tmp = 1.0;
329  do {
330  value += 1.0;
331  tmp *= uniform();
332  } while (tmp > m_exp_lambda);
333  } // endif: direct method used
334 
335  // ... otherwise use rejection method
336  else {
337  double tmp;
338  if (lambda != m_old_lambda) {
339  m_old_lambda = lambda;
340  m_sqrt_lambda = std::sqrt(2.0*lambda);
341  m_log_lambda = std::log(lambda);
342  m_exp_lambda = lambda * m_log_lambda - gammalib::gammln(lambda+1.0);
343  }
344  do {
345  double factor;
346  do {
347  factor = std::tan(gammalib::pi * uniform());
348  value = m_sqrt_lambda * factor + lambda;
349  } while (value < 0.0);
350  value = floor(value);
351  tmp = 0.9*(1.0+factor*factor) *
353  } while (uniform() > tmp);
354  }
355 
356  // Return random deviate
357  return value;
358 }
359 
360 
361 /***********************************************************************//**
362  * @brief Returns Chi2 deviates for 2 degrees of freedom
363  *
364  * Returns exponential deviates from the probability distribution
365  * \f[p(x) = \frac{1}{2\pi} x \exp( -\frac{1}{2} x^2 )\f]
366  * This method can be used to simulate the random radial offset of a measured
367  * source position from the true source position, assuming an azimuthally
368  * symmetric 2D Gaussian probability distribution.
369  ***************************************************************************/
370 double GRan::chisq2(void)
371 {
372  // Allocate argument
373  double x;
374 
375  // Get uniform deviate < 1
376  do {
377  x = uniform();
378  } while (x == 1.0);
379 
380  // Return random value
381  return (std::sqrt(-2.0*std::log(1.0-x)));
382 }
383 
384 
385 /***********************************************************************//**
386  * @brief Random sampling from a cumulative density function
387  *
388  * @param[in] cdf Array containing cumulative density function
389  ***************************************************************************/
390 int GRan::cdf(const std::vector<double>& cdf)
391 {
392  // Get uniform random number
393  double u = uniform();
394 
395  // Get pixel index according to random number. We use a bi-section
396  // method to find the corresponding pixel index
397  int low = 0;
398  int high = cdf.size();
399  while ((high - low) > 1) {
400  int mid = (low+high) / 2;
401  if (u < cdf[mid]) {
402  high = mid;
403  }
404  else if (cdf[mid] <= u) {
405  low = mid;
406  }
407  }
408 
409  // Return index
410  return low;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Random sampling from a cumulative density function
416  *
417  * @param[in] cdf Vector containing cumulative density function
418  *
419  * @todo Somehow merge with the above method to reduce code repetition
420  ***************************************************************************/
421 int GRan::cdf(const GVector& cdf)
422 {
423  // Get uniform random number
424  double u = uniform();
425 
426  // Get pixel index according to random number. We use a bi-section
427  // method to find the corresponding pixel index
428  int low = 0;
429  int high = cdf.size();
430  while ((high - low) > 1) {
431  int mid = (low+high) / 2;
432  if (u < cdf[mid]) {
433  high = mid;
434  }
435  else if (cdf[mid] <= u) {
436  low = mid;
437  }
438  }
439 
440  // Return index
441  return low;
442 }
443 
444 
445 /***********************************************************************//**
446  * @brief Print random number generator information
447  *
448  * @param[in] chatter Chattiness (defaults to NORMAL).
449  * @return String containing random number generator information.
450  ***************************************************************************/
451 std::string GRan::print(const GChatter& chatter) const
452 {
453  // Initialise result string
454  std::string result;
455 
456  // Continue only if chatter is not silent
457  if (chatter != SILENT) {
458 
459  // Append header
460  result.append("=== GRan ===");
461 
462  // Append information
463  result.append("\n"+gammalib::parformat("Seed")+gammalib::str(m_seed));
464  result.append("\n"+gammalib::parformat("Value 1")+gammalib::str(m_value1));
465  result.append("\n"+gammalib::parformat("Value 2")+gammalib::str(m_value2));
466  result.append("\n"+gammalib::parformat("Value 3")+gammalib::str(m_value3));
467 
468  } // endif: chatter was not silent
469 
470  // Return result
471  return result;
472 }
473 
474 
475 /*==========================================================================
476  = =
477  = Private methods =
478  = =
479  ==========================================================================*/
480 
481 /***********************************************************************//**
482  * @brief Initialise class members
483  ***************************************************************************/
484 void GRan::init_members(unsigned long long int seed)
485 {
486  // Initialise members
487  m_seed = seed;
488  m_value2 = 4101842887655102017LL;
489  m_value3 = 1;
490  m_old_lambda = -1.0;
491  m_sqrt_lambda = 0.0;
492  m_log_lambda = 0.0;
493  m_exp_lambda = 0.0;
494 
495  // Initialise values
497  int64();
498  m_value2 = m_value1;
499  int64();
500  m_value3 = m_value2;
501  int64();
502 
503  // Return
504  return;
505 }
506 
507 
508 /***********************************************************************//**
509  * @brief Copy class members
510  *
511  * @param[in] ran Random number generator.
512  ***************************************************************************/
513 void GRan::copy_members(const GRan& ran)
514 {
515  // Copy members
516  m_seed = ran.m_seed;
517  m_value1 = ran.m_value1;
518  m_value2 = ran.m_value2;
519  m_value3 = ran.m_value3;
524 
525  // Return
526  return;
527 }
528 
529 
530 /***********************************************************************//**
531  * @brief Delete class members
532  ***************************************************************************/
534 {
535  // Return
536  return;
537 }
double chisq2(void)
Returns Chi2 deviates for 2 degrees of freedom.
Definition: GRan.cpp:370
double m_exp_lambda
exp(-lambda)
Definition: GRan.hpp:89
Random number generator class definition.
std::string print(const GChatter &chatter=NORMAL) const
Print random number generator information.
Definition: GRan.cpp:451
const double pi
Definition: GMath.hpp:35
double m_sqrt_lambda
sqrt(2*lambda)
Definition: GRan.hpp:87
unsigned long long int int64(void)
Return 64-bit random unsigned integer.
Definition: GRan.cpp:216
void copy_members(const GRan &ran)
Copy class members.
Definition: GRan.cpp:513
GRan(void)
Void constructor.
Definition: GRan.cpp:57
int cdf(const std::vector< double > &cdf)
Random sampling from a cumulative density function.
Definition: GRan.cpp:390
unsigned long long int m_value2
Value 2.
Definition: GRan.hpp:82
Random number generator class.
Definition: GRan.hpp:44
Gammalib tools definition.
void free_members(void)
Delete class members.
Definition: GRan.cpp:533
double m_old_lambda
Old lambda value.
Definition: GRan.hpp:86
unsigned long long int m_value1
Value 1.
Definition: GRan.hpp:81
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1379
unsigned long int int32(void)
Return 32-bit random unsigned integer.
Definition: GRan.cpp:206
unsigned long long int m_value3
Value 3.
Definition: GRan.hpp:83
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
GRan * clone(void) const
Clone random number generator.
Definition: GRan.cpp:176
unsigned long long int seed(void) const
Return seed value.
Definition: GRan.hpp:111
GChatter
Definition: GTypemaps.hpp:33
double normal(void)
Returns normal deviates.
Definition: GRan.cpp:257
unsigned long long int m_seed
Random number generator seed.
Definition: GRan.hpp:80
void init_members(unsigned long long int seed=41L)
Initialise class members.
Definition: GRan.cpp:484
double gammln(const double &arg)
Computes logarithm of gamma function.
Definition: GMath.cpp:383
GRan & operator=(const GRan &ran)
Assignment operator.
Definition: GRan.cpp:125
virtual ~GRan(void)
Destructor.
Definition: GRan.cpp:103
double m_log_lambda
log(lambda)
Definition: GRan.hpp:88
double poisson(const double &lambda)
Returns Poisson deviates.
Definition: GRan.cpp:316
Exception handler interface definition.
void clear(void)
Clear random number generator.
Definition: GRan.cpp:158
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:178
Vector class.
Definition: GVector.hpp:46
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
Mathematical function definitions.
double exp(const double &lambda)
Returns exponential deviates.
Definition: GRan.cpp:291
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489