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

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

fixes #133 removed all errors reported from Doxygen. Only one error left which says Index is not documented but I don't want it to be documented actually we use the Doxygens preprocessor to skip documenting that class, yet Doxygen complains that class is not documented huh. Only solution would be to move that class to its own file and not keep it together with SVM.

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