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

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

Added missing virtual destructors.

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