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

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

removed min parameter in uniform_int and added test to rnd_test

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