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

Last change on this file since 674 was 674, checked in by Peter, 15 years ago

fixes #152

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