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

Last change on this file since 389 was 389, checked in by Peter, 17 years ago

moved kernel to regression namespace and tried to fix some dox issues

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