source: trunk/lib/random/Random.h @ 374

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

Added documentation.
Added missing Makefile.am.
Implemented seeding of the rng. Implemented seed acquring from /dev/urandom.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.7 KB
Line 
1// $Id: Random.h 374 2005-08-05 22:30:57Z 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    /// @param \a range is number of different integers that can be
279    /// generated \f$ (b-a+1) \f$ \a min is the minimal integer that
280    /// can be generated \f$ (a) \f$
281    ///
282    inline DiscreteUniform(const u_long range, const long min=0) 
283      : min_(min), range_(range) {}
284
285    ///
286    /// @return A random number between \a min_ and \a range_ - \a min_
287    ///
288    inline long operator()(void) const 
289    { return gsl_rng_uniform_int(rng_->rng(), range_)+min_; }
290
291    ///
292    /// @return A random number between 0 and \a range.
293    ///
294    inline long operator()(const u_long range) const 
295    { return gsl_rng_uniform_int(rng_->rng(), range); }
296
297    ///
298    /// @return A random number between \a min and \a range - \a min
299    /// @note this operator ignores parameters given in Constructor
300    ///
301    inline long operator()(const u_long range, const long min) const 
302    { return gsl_rng_uniform_int(rng_->rng(), range)+min; }
303
304
305  private:
306    long min_;
307    u_long range_;
308  };
309
310
311  ///
312  /// @brief Poisson Distribution
313  ///
314  /// Having a Poisson process (no memory), number of occurences
315  /// within a given time window is Poisson distributed. This
316  /// distribution is the limit of a Binomial distribution when number
317  /// of attempts is large, and the probability for one attempt to be
318  /// succesful is small (in such a way that the expected number of
319  /// succesful attempts is \f$ m \f$.
320  ///
321  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
322  /// k  \f$ \n
323  /// Expectation value: \f$ m \f$ \n
324  /// Variance: \f$ m \f$
325  ///
326  class Poisson : public Discrete
327  {
328  public:
329    ///
330    /// @brief Constructor
331    ///
332    /// @param m is expectation value
333    inline Poisson(const double m=1) 
334      : m_(m){}
335
336    ///
337    /// @return A Poisson distributed number.
338    ///
339    inline long operator()(void) const 
340    { return gsl_ran_poisson(rng_->rng(), m_); }
341
342    ///
343    /// @return A Poisson distributed number with expectation value \a
344    /// m
345    /// @note this operator ignores parameters set in Constructor
346    inline u_long operator()(const double m) const 
347    { return gsl_ran_poisson(rng_->rng(), m); }
348
349  private:
350    double m_;
351  };
352
353  ///
354  /// @brief General
355  ///
356  class DiscreteGeneral : public Discrete
357  {
358  public:
359    ///
360    /// @brief Constructor
361    ///
362    /// @param hist is a Histogram defining the probability distribution
363    ///
364    DiscreteGeneral(const statistics::Histogram& hist) ;
365   
366    ///
367    /// @brief Destructor
368    ///
369    virtual ~DiscreteGeneral(void);
370
371    ///
372    /// The generated number is an integer and proportinal to the
373    /// frequency in the corresponding histogram bin. In other words,
374    /// the probability that 0 is returned is proportinal to the size
375    /// of the first bin.
376    ///
377    /// @return A random number.
378    ///
379    inline long operator()(void) const 
380    { return gsl_ran_discrete(rng_->rng(), gen_); }
381
382  private:
383     gsl_ran_discrete_t* gen_;
384  };
385
386
387  ///
388  /// Class to generate numbers from a histogram in a continuous manner.
389  ///
390  class ContinuousGeneral : public Continuous
391  {
392  public:
393    ///
394    /// @brief Constructor
395    ///
396    /// @param hist is a Histogram defining the probability distribution
397    ///
398    inline ContinuousGeneral(const statistics::Histogram& hist) 
399      : discrete_(DiscreteGeneral(hist)), hist_(hist) {}
400   
401    ///
402    /// @brief Destructor
403    ///
404    virtual ~ContinuousGeneral(void);
405
406    ///
407    /// The number is generated in a two step process. First the bin
408    /// in the histogram is randomly selected (see
409    /// DiscreteGeneral). Then a number is generated uniformly from
410    /// the interval defined by the bin.
411    ///
412    /// @return A random number.
413    ///
414    inline double operator()(void) const 
415    { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); }
416
417  private:
418    const DiscreteGeneral discrete_;
419    const statistics::Histogram& hist_;
420    ContinuousUniform u_;
421  };
422
423
424}} // of namespace random and namespace theplu
425
426#endif
Note: See TracBrowser for help on using the repository browser.