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

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

Fixed documentation and moved around code.

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