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

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

Changed filenames lib/random/Random.* to lib/random/random.*

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.7 KB
Line 
1// $Id: random.h 375 2005-08-07 21:51:58Z 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.