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

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

Minor documentation change.

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