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

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

Fixed bug in DiscreteUniform?.
Added seeding support for classes derived from Discrete and Continuous
virtual classes.

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