source: trunk/lib/random/random.h @ 562

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

Corrected copyright statement.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 13.9 KB
Line 
1// $Id: random.h 562 2006-03-15 09:19:51Z 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  /// @todo Get those exceptions in! Should classes be considered
72  /// const if underlying structures are changed such as GSL stuff?
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    inline const gsl_rng* rng(void) const { return rng_; }
114
115    ///
116    /// @brief Set the seed \a s for the rng.
117    ///
118    /// Set the seed \a s for the rng. If \a s is zero, a default
119    /// value from the rng's original implementation is used (cf. GSL
120    /// documentation).
121    ///
122    /// @see seed_from_devurandom
123    ///
124    inline void seed(u_long s) const { gsl_rng_set(rng_,s); }
125
126    ///
127    /// @brief Seed the rng using the /dev/urandom device.
128    ///
129    /// @return The seed acquired from /dev/urandom.
130    ///
131    u_long seed_from_devurandom(void);
132
133  private:
134    RNG(void);
135
136    static RNG* instance_;
137    gsl_rng* rng_;
138  };
139
140  // --------------------- Discrete distribtuions ---------------------
141
142  ///
143  /// @brief Discrete random number distributions.
144  ///
145  /// Abstract base class for discrete random number
146  /// distributions. Given K discrete events with different
147  /// probabilities \f$ P[k] \f$, produce a random value k consistent
148  /// with its probability.
149  ///
150  class Discrete
151  {
152  public:
153    ///
154    /// @brief Constructor
155    ///
156    inline Discrete(void) { rng_=RNG::instance(); }
157
158    inline virtual ~Discrete(void) { }
159
160    ///
161    /// @brief Set the seed to \a s.
162    ///
163    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
164    /// default value from the rng's original implementation is used
165    /// (cf. GSL documentation).
166    ///
167    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
168    ///
169    inline void seed(u_long s) const { rng_->seed(s); }
170
171    ///
172    /// @brief Set the seed using the /dev/urandom device.
173    ///
174    /// @return The seed acquired from /dev/urandom.
175    ///
176    /// @see seed, RNG::seed_from_devurandom, RNG::seed
177    ///
178    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
179
180    ///
181    /// @return A random number.
182    ///
183    virtual u_long operator()(void) const = 0;
184   
185  protected:
186    RNG* rng_;
187  };
188
189  ///
190  /// @brief General
191  ///
192  class DiscreteGeneral : public Discrete
193  {
194  public:
195    ///
196    /// @brief Constructor
197    ///
198    /// @param hist is a Histogram defining the probability distribution
199    ///
200    DiscreteGeneral(const statistics::Histogram& hist);
201   
202    ///
203    /// @brief Destructor
204    ///
205    ~DiscreteGeneral(void);
206
207    ///
208    /// The generated number is an integer and proportinal to the
209    /// frequency in the corresponding histogram bin. In other words,
210    /// the probability that 0 is returned is proportinal to the size
211    /// of the first bin.
212    ///
213    /// @return A random number.
214    ///
215    inline u_long
216    operator()(void) const { return gsl_ran_discrete(rng_->rng(), gen_); }
217
218  private:
219     gsl_ran_discrete_t* gen_;
220  };
221
222  ///
223  /// @brief Discrete uniform distribution
224  ///
225  /// Discrete uniform distribution also known as the "equally likely
226  /// outcomes" distribution. Each outcome, in this case an integer
227  /// from [0,n-1] , have equal probability to occur.
228  ///
229  /// Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
230  /// k < n \f$ \n
231  /// Expectation value: \f$ \frac{n-1}{2} \f$ \n
232  /// Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
233  ///
234  class DiscreteUniform : public Discrete
235  {
236  public:
237    ///
238    /// @brief Default constructor.
239    ///
240    DiscreteUniform(void) : range_(rng_->max()) {}
241
242    ///
243    /// @brief Constructor.
244    ///
245    /// The generator will generate integers from \f$ [0,n-1] \f$. If
246    /// \a n is larger than the maximum number the random number
247    /// generator can return, then (currently) \a n is adjusted
248    /// appropriately.
249    ///
250    /// @todo If a too large \a n is given an exception should be
251    /// thrown, i.e. the behaviour of this class will change. The case
252    /// when argument is 0 is not treated gracefully (underlying GSL
253    /// functionality will not return).
254    ///
255    DiscreteUniform(const u_long n) : range_(n)
256    { if ( range_>rng_->max() ) range_=rng_->max(); }
257
258    ///
259    /// This function returns a random integer from 0 to n-1
260    /// inclusive. All integers in the range [0,n-1] are equally
261    /// likely. n is set in constructor.
262    ///
263    inline u_long
264    operator()(void) const { return gsl_rng_uniform_int(rng_->rng(), range_); }
265
266    ///
267    /// This function returns a random integer from 0 to n-1
268    /// inclusive. All integers in the range [0,n-1] are equally
269    /// likely.
270    ///
271    inline u_long operator()(const u_long n) const
272    { return gsl_rng_uniform_int(rng_->rng(), n); }
273
274  private:
275    u_long range_;
276  };
277
278  ///
279  /// @brief Poisson Distribution
280  ///
281  /// Having a Poisson process (i.e. no memory), number of occurences
282  /// within a given time window is Poisson distributed. This
283  /// distribution is the limit of a Binomial distribution when number
284  /// of attempts is large, and the probability for one attempt to be
285  /// succesful is small (in such a way that the expected number of
286  /// succesful attempts is \f$ m \f$.
287  ///
288  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
289  /// k  \f$ \n
290  /// Expectation value: \f$ m \f$ \n
291  /// Variance: \f$ m \f$
292  ///
293  class Poisson : public Discrete
294  {
295  public:
296    ///
297    /// @brief Constructor
298    ///
299    /// @param m is expectation value
300    inline Poisson(const double m=1) : m_(m) {}
301
302    ///
303    /// @return A Poisson distributed number.
304    ///
305    inline u_long
306    operator()(void) const { return gsl_ran_poisson(rng_->rng(), m_); }
307
308    ///
309    /// @return A Poisson distributed number with expectation value \a
310    /// m
311    ///
312    /// @note this operator ignores parameters set in Constructor
313    ///
314    inline u_long
315    operator()(const double m) const { return gsl_ran_poisson(rng_->rng(), m); }
316
317  private:
318    double m_;
319  };
320
321  // --------------------- Continuous distribtuions ---------------------
322
323  ///
324  /// @brief Continuous random number distributions.
325  ///
326  /// Abstract base class for continuous random number distributions.
327  ///
328  class Continuous
329  {
330  public:
331
332    ///
333    /// @brief Constructor
334    ///
335    inline Continuous(void) { rng_=RNG::instance(); }
336
337    inline virtual ~Continuous(void) { }
338
339    ///
340    /// @brief Set the seed to \a s.
341    ///
342    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
343    /// default value from the rng's original implementation is used
344    /// (cf. GSL documentation).
345    ///
346    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
347    ///
348    inline void seed(u_long s) const { rng_->seed(s); }
349
350    ///
351    /// @brief Set the seed using the /dev/urandom device.
352    ///
353    /// @return The seed acquired from /dev/urandom.
354    ///
355    /// @see seed, RNG::seed_from_devurandom, RNG::seed
356    ///
357    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
358
359    ///
360    /// @return A random number
361    ///
362    virtual double operator()(void) const = 0;
363
364  protected:
365    RNG* rng_;
366  };
367
368  ///
369  /// @brief Uniform distribution
370  ///
371  /// Class for generating a random number from a uniform distribution
372  /// in the range [0,1), i.e. zero is included but not 1.
373  ///
374  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
375  /// Expectation value: 0.5 \n
376  /// Variance: \f$ \frac{1}{12} \f$
377  ///
378  class ContinuousUniform : public Continuous
379  {
380  public:
381    inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());}
382  };
383
384  ///
385  /// Class to generate numbers from a histogram in a continuous manner.
386  ///
387  class ContinuousGeneral : public Continuous
388  {
389  public:
390    ///
391    /// @brief Constructor
392    ///
393    /// @param hist is a Histogram defining the probability distribution
394    ///
395    inline ContinuousGeneral(const statistics::Histogram& hist)
396      : discrete_(DiscreteGeneral(hist)), hist_(hist) {}
397
398    ///
399    /// @brief Destructor
400    ///
401    ~ContinuousGeneral(void);
402
403    ///
404    /// The number is generated in a two step process. First the bin
405    /// in the histogram is randomly selected (see
406    /// DiscreteGeneral). Then a number is generated uniformly from
407    /// the interval defined by the bin.
408    ///
409    /// @return A random number.
410    ///
411    inline double operator()(void) const 
412    { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); }
413
414  private:
415    const DiscreteGeneral discrete_;
416    const statistics::Histogram& hist_;
417    ContinuousUniform u_;
418  };
419
420  ///
421  /// @brief Generator of random numbers from an exponential
422  /// distribution.
423  ///
424  /// The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
425  /// \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
426  /// variance \f$ m^2 \f$
427  ///
428  class Exponential : public Continuous
429  {
430  public:
431    ///
432    /// @brief Constructor
433    ///
434    /// @param \a m is the expectation value of the distribution.
435    ///
436    inline Exponential(const double m=1) : m_(m) {}
437
438    ///
439    /// @return A random number from exponential distribution.
440    ///
441    inline double
442    operator()(void) const { return gsl_ran_exponential(rng_->rng(), m_); }
443
444    ///
445    /// @return A random number from exponential distribution, with
446    /// expectation value \a m
447    ///
448    /// @note This operator ignores parameters given in constructor.
449    ///
450    inline double operator()(const double m) const
451    { return gsl_ran_exponential(rng_->rng(), m); }
452
453  private:
454    double m_;
455  };
456
457  ///
458  /// @brief Gaussian distribution
459  ///
460  /// Class for generating a random number from a Gaussian
461  /// distribution between zero and unity. Utilizes the Box-Muller
462  /// algorithm, which needs two calls to random generator.
463  ///
464  /// Distribution function \f$ f(x) =
465  /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
466  /// \f$ \n
467  /// Expectation value: \f$ \mu \f$ \n
468  /// Variance: \f$ \sigma^2 \f$
469  ///
470  class Gaussian : public Continuous
471  {
472  public:
473    ///
474    /// @brief Constructor
475    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
476    /// m is the expectation value \f$ \mu \f$ of the distribution
477    ///
478    inline Gaussian(const double s=1, const double m=0) : m_(m), s_(s) {}
479
480    ///
481    /// @return A random Gaussian number
482    ///
483    inline double
484    operator()(void) const { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }
485
486    ///
487    /// @return A random Gaussian number with standard deviation \a s
488    /// and expectation value 0.
489    /// @note this operator ignores parameters given in Constructor
490    ///
491    inline double
492    operator()(const double s) const { return gsl_ran_gaussian(rng_->rng(), s); }
493
494    ///
495    /// @return A random Gaussian number with standard deviation \a s
496    /// and expectation value \a m.
497    /// @note this operator ignores parameters given in Constructor
498    ///
499    inline double operator()(const double s, const double m) const
500    { return gsl_ran_gaussian(rng_->rng(), s)+m; }
501
502  private:
503    double m_;
504    double s_;
505  };
506
507
508}} // of namespace random and namespace theplu
509
510#endif
Note: See TracBrowser for help on using the repository browser.