source: trunk/c++_tools/random/random.h @ 614

Last change on this file since 614 was 614, checked in by Jari Häkkinen, 15 years ago

Minor fixes in doxygen documentation.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.7 KB
Line 
1// $Id: random.h 614 2006-08-31 05:08:05Z jari $
2
3/*
4  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson
5
6  This file is part of the thep c++ tools library,
7                                http://lev.thep.lu.se/trac/c++_tools
8
9  The c++ tools library is free software; you can redistribute it
10  and/or modify it under the terms of the GNU General Public License
11  as published by the Free Software Foundation; either version 2 of
12  the License, or (at your option) any later version.
13
14  The c++ tools library is distributed in the hope that it will be
15  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
16  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  General Public License for more details.
18
19  You should have received a copy of the GNU General Public License
20  along with this program; if not, write to the Free Software
21  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
22  02111-1307, USA.
23*/
24
25#ifndef _theplu_random_
26#define _theplu_random_
27
28#include <c++_tools/statistics/Histogram.h>
29
30#include <gsl/gsl_rng.h>
31#include <gsl/gsl_randist.h>
32
33#include <string>
34
35namespace theplu {
36namespace random {
37
38  ///
39  /// @brief Random Number Generator
40  ///
41  /// The RNG class is wrapper to the GSL random number generator
42  /// (rng). This class provides a single global instance of the rng,
43  /// and makes sure there is only one point of access to the
44  /// generator.
45  ///
46  /// There is information about how to change seeding and generators
47  /// at run time without recompilation using environment
48  /// variables. RNG of course support seeding at compile time if you
49  /// don't want to bother about environment variables and GSL.
50  ///
51  /// There are many different rng's available in GSL. Currently only
52  /// the default generator is implemented and no other one is
53  /// choosable through the class interface. This means that you have
54  /// to fall back to the use of environment variables as described in
55  /// the GSL documentation, or be bold and request support for other
56  /// rng's through the class interface.
57  ///
58  /// Not all GSL functionality is implemented, we'll add
59  /// functionality when needed and may do it when requested. Better
60  /// yet, supply us with code and we will probably add it to the code
61  /// (BUT remember to implement reasonable tests for your code and
62  /// follow the coding style.)
63  ///
64  /// This implementation may be thread safe (according to GSL
65  /// documentation), but should be checked to be so before trusting
66  /// thread safety.
67  ///
68  /// @see Design Patterns (the singleton and adapter pattern). GSL
69  /// documentation.
70  ///
71  class RNG
72  {
73  public:
74
75    virtual ~RNG(void);
76
77    ///
78    /// @brief Get an instance of the random number generator.
79    ///
80    /// Get an instance of the random number generator. If the random
81    /// number generator is not already created, the call will create
82    /// a new generator and use the default seed. The seed must be
83    /// changed with the seed or seed_from_devurandom member
84    /// functions.
85    ///
86    /// @return A pointer to the random number generator.
87    ///
88    /// @see seed and seed_from_devurandom
89    ///
90    static RNG* instance(void)
91      { if (!instance_) instance_=new RNG; return instance_; }
92
93    ///
94    /// @brief Returns the largest number that the random number
95    /// generator can return.
96    ///
97    inline u_long max(void) const { return gsl_rng_max(rng_); }
98
99    ///
100    /// @brief Returns the smallest number that the random number
101    /// generator can return.
102    ///
103    inline u_long min(void) const { return gsl_rng_min(rng_); }
104
105    ///
106    /// @brief Returns the name of the random number generator
107    ///
108    inline std::string name(void) const { return gsl_rng_name(rng_); }
109
110    inline const gsl_rng* rng(void) const { return rng_; }
111
112    ///
113    /// @brief Set the seed \a s for the rng.
114    ///
115    /// Set the seed \a s for the rng. If \a s is zero, a default
116    /// value from the rng's original implementation is used (cf. GSL
117    /// documentation).
118    ///
119    /// @see seed_from_devurandom
120    ///
121    inline void seed(u_long s) const { gsl_rng_set(rng_,s); }
122
123    ///
124    /// @brief Seed the rng using the /dev/urandom device.
125    ///
126    /// @return The seed acquired from /dev/urandom.
127    ///
128    u_long seed_from_devurandom(void);
129
130  private:
131    RNG(void);
132
133    static RNG* instance_;
134    gsl_rng* rng_;
135  };
136
137  // --------------------- Discrete distribtuions ---------------------
138
139  ///
140  /// @brief Discrete random number distributions.
141  ///
142  /// Abstract base class for discrete random number
143  /// distributions. Given K discrete events with different
144  /// probabilities \f$ P[k] \f$, produce a random value k consistent
145  /// with its probability.
146  ///
147  class Discrete
148  {
149  public:
150    ///
151    /// @brief Constructor
152    ///
153    inline Discrete(void) { rng_=RNG::instance(); }
154
155    inline virtual ~Discrete(void) { }
156
157    ///
158    /// @brief Set the seed to \a s.
159    ///
160    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
161    /// default value from the rng's original implementation is used
162    /// (cf. GSL documentation).
163    ///
164    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
165    ///
166    inline void seed(u_long s) const { rng_->seed(s); }
167
168    ///
169    /// @brief Set the seed using the /dev/urandom device.
170    ///
171    /// @return The seed acquired from /dev/urandom.
172    ///
173    /// @see seed, RNG::seed_from_devurandom, RNG::seed
174    ///
175    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
176
177    ///
178    /// @return A random number.
179    ///
180    virtual u_long operator()(void) const = 0;
181   
182  protected:
183    RNG* rng_;
184  };
185
186  ///
187  /// @brief General
188  ///
189  class DiscreteGeneral : public Discrete
190  {
191  public:
192    ///
193    /// @brief Constructor
194    ///
195    /// @param hist is a Histogram defining the probability distribution
196    ///
197    DiscreteGeneral(const statistics::Histogram& hist);
198   
199    ///
200    /// @brief Destructor
201    ///
202    ~DiscreteGeneral(void);
203
204    ///
205    /// The generated number is an integer and proportinal to the
206    /// frequency in the corresponding histogram bin. In other words,
207    /// the probability that 0 is returned is proportinal to the size
208    /// of the first bin.
209    ///
210    /// @return A random number.
211    ///
212    inline u_long
213    operator()(void) const { return gsl_ran_discrete(rng_->rng(), gen_); }
214
215  private:
216     gsl_ran_discrete_t* gen_;
217  };
218
219  ///
220  /// @brief Discrete uniform distribution
221  ///
222  /// Discrete uniform distribution also known as the "equally likely
223  /// outcomes" distribution. Each outcome, in this case an integer
224  /// from [0,n-1] , have equal probability to occur.
225  ///
226  /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
227  /// k < n \f$ \n
228  /// Expectation value: \f$ \frac{n-1}{2} \f$ \n
229  /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
230  ///
231  class DiscreteUniform : public Discrete
232  {
233  public:
234    ///
235    /// @brief Default constructor.
236    ///
237    DiscreteUniform(void) : range_(rng_->max()) {}
238
239    ///
240    /// @brief Constructor.
241    ///
242    /// The generator will generate integers from \f$ [0,n-1] \f$. If
243    /// \a n is larger than the maximum number the random number
244    /// generator can return, then (currently) \a n is adjusted
245    /// appropriately.
246    ///
247    /// @todo If a too large \a n is given an exception should be
248    /// thrown, i.e. the behaviour of this class will change. The case
249    /// when argument is 0 is not treated gracefully (underlying GSL
250    /// functionality will not return).
251    ///
252    DiscreteUniform(const u_long n) : range_(n)
253    { if ( range_>rng_->max() ) range_=rng_->max(); }
254
255    ///
256    /// This function returns a random integer from 0 to n-1
257    /// inclusive. All integers in the range [0,n-1] are equally
258    /// likely. n is set in constructor.
259    ///
260    inline u_long
261    operator()(void) const { return gsl_rng_uniform_int(rng_->rng(), range_); }
262
263    ///
264    /// This function returns a random integer from 0 to n-1
265    /// inclusive. All integers in the range [0,n-1] are equally
266    /// likely.
267    ///
268    inline u_long operator()(const u_long n) const
269    { return gsl_rng_uniform_int(rng_->rng(), n); }
270
271  private:
272    u_long range_;
273  };
274
275  ///
276  /// @brief Poisson Distribution
277  ///
278  /// Having a Poisson process (i.e. no memory), number of occurences
279  /// within a given time window is Poisson distributed. This
280  /// distribution is the limit of a Binomial distribution when number
281  /// of attempts is large, and the probability for one attempt to be
282  /// succesful is small (in such a way that the expected number of
283  /// succesful attempts is \f$ m \f$.
284  ///
285  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
286  /// k  \f$ \n
287  /// Expectation value: \f$ m \f$ \n
288  /// Variance: \f$ m \f$
289  ///
290  class Poisson : public Discrete
291  {
292  public:
293    ///
294    /// @brief Constructor
295    ///
296    /// @param m is expectation value
297    inline Poisson(const double m=1) : m_(m) {}
298
299    ///
300    /// @return A Poisson distributed number.
301    ///
302    inline u_long
303    operator()(void) const { return gsl_ran_poisson(rng_->rng(), m_); }
304
305    ///
306    /// @return A Poisson distributed number with expectation value \a
307    /// m
308    ///
309    /// @note this operator ignores parameters set in Constructor
310    ///
311    inline u_long
312    operator()(const double m) const { return gsl_ran_poisson(rng_->rng(), m); }
313
314  private:
315    double m_;
316  };
317
318  // --------------------- Continuous distribtuions ---------------------
319
320  ///
321  /// @brief Continuous random number distributions.
322  ///
323  /// Abstract base class for continuous random number distributions.
324  ///
325  class Continuous
326  {
327  public:
328
329    ///
330    /// @brief Constructor
331    ///
332    inline Continuous(void) { rng_=RNG::instance(); }
333
334    inline virtual ~Continuous(void) { }
335
336    ///
337    /// @brief Set the seed to \a s.
338    ///
339    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
340    /// default value from the rng's original implementation is used
341    /// (cf. GSL documentation).
342    ///
343    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
344    ///
345    inline void seed(u_long s) const { rng_->seed(s); }
346
347    ///
348    /// @brief Set the seed using the /dev/urandom device.
349    ///
350    /// @return The seed acquired from /dev/urandom.
351    ///
352    /// @see seed, RNG::seed_from_devurandom, RNG::seed
353    ///
354    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
355
356    ///
357    /// @return A random number
358    ///
359    virtual double operator()(void) const = 0;
360
361  protected:
362    RNG* rng_;
363  };
364
365  ///
366  /// @brief Uniform distribution
367  ///
368  /// Class for generating a random number from a uniform distribution
369  /// in the range [0,1), i.e. zero is included but not 1.
370  ///
371  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
372  /// Expectation value: 0.5 \n
373  /// Variance: \f$ \frac{1}{12} \f$
374  ///
375  class ContinuousUniform : public Continuous
376  {
377  public:
378    inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());}
379  };
380
381  ///
382  /// Class to generate numbers from a histogram in a continuous manner.
383  ///
384  class ContinuousGeneral : public Continuous
385  {
386  public:
387    ///
388    /// @brief Constructor
389    ///
390    /// @param hist is a Histogram defining the probability distribution
391    ///
392    inline ContinuousGeneral(const statistics::Histogram& hist)
393      : discrete_(DiscreteGeneral(hist)), hist_(hist) {}
394
395    ///
396    /// @brief Destructor
397    ///
398    ~ContinuousGeneral(void);
399
400    ///
401    /// The number is generated in a two step process. First the bin
402    /// in the histogram is randomly selected (see
403    /// DiscreteGeneral). Then a number is generated uniformly from
404    /// the interval defined by the bin.
405    ///
406    /// @return A random number.
407    ///
408    inline double operator()(void) const 
409    { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); }
410
411  private:
412    const DiscreteGeneral discrete_;
413    const statistics::Histogram& hist_;
414    ContinuousUniform u_;
415  };
416
417  ///
418  /// @brief Generator of random numbers from an exponential
419  /// distribution.
420  ///
421  /// The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
422  /// \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
423  /// variance \f$ m^2 \f$
424  ///
425  class Exponential : public Continuous
426  {
427  public:
428    ///
429    /// @brief Constructor
430    ///
431    /// @param \a m is the expectation value of the distribution.
432    ///
433    inline Exponential(const double m=1) : m_(m) {}
434
435    ///
436    /// @return A random number from exponential distribution.
437    ///
438    inline double
439    operator()(void) const { return gsl_ran_exponential(rng_->rng(), m_); }
440
441    ///
442    /// @return A random number from exponential distribution, with
443    /// expectation value \a m
444    ///
445    /// @note This operator ignores parameters given in constructor.
446    ///
447    inline double operator()(const double m) const
448    { return gsl_ran_exponential(rng_->rng(), m); }
449
450  private:
451    double m_;
452  };
453
454  ///
455  /// @brief Gaussian distribution
456  ///
457  /// Class for generating a random number from a Gaussian
458  /// distribution between zero and unity. Utilizes the Box-Muller
459  /// algorithm, which needs two calls to random generator.
460  ///
461  /// Distribution function \f$ f(x) =
462  /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
463  /// \f$ \n
464  /// Expectation value: \f$ \mu \f$ \n
465  /// Variance: \f$ \sigma^2 \f$
466  ///
467  class Gaussian : public Continuous
468  {
469  public:
470    ///
471    /// @brief Constructor
472    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
473    /// m is the expectation value \f$ \mu \f$ of the distribution
474    ///
475    inline Gaussian(const double s=1, const double m=0) : m_(m), s_(s) {}
476
477    ///
478    /// @return A random Gaussian number
479    ///
480    inline double
481    operator()(void) const { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }
482
483    ///
484    /// @return A random Gaussian number with standard deviation \a s
485    /// and expectation value 0.
486    /// @note this operator ignores parameters given in Constructor
487    ///
488    inline double
489    operator()(const double s) const { return gsl_ran_gaussian(rng_->rng(), s); }
490
491    ///
492    /// @return A random Gaussian number with standard deviation \a s
493    /// and expectation value \a m.
494    /// @note this operator ignores parameters given in Constructor
495    ///
496    inline double operator()(const double s, const double m) const
497    { return gsl_ran_gaussian(rng_->rng(), s)+m; }
498
499  private:
500    double m_;
501    double s_;
502  };
503
504
505}} // of namespace random and namespace theplu
506
507#endif
Note: See TracBrowser for help on using the repository browser.