source:trunk/lib/random/Random.h@367

Last change on this file since 367 was 367, checked in by Peter, 16 years ago

• Property svn:eol-style set to native
• Property svn:keywords set to Id
File size: 7.1 KB
Line
1 // $Id: Random.h 367 2005-08-05 11:55:16Z peter$
2
3#ifndef _theplu_random_
4#define _theplu_random_
5
6
7#include <gsl/gsl_rng.h>
8#include <gsl/gsl_randist.h>
9
10namespace theplu {
11namespace random {
12
13  ///
14  /// The RNG class provides a single global random number generator
16  ///
17  /// This is probably not thread safe.
18  ///
19  /// @see Design Patterns (the singleton and adapter pattern)
20  ///
21  class RNG
22  {
23  public:
24
25    virtual ~RNG(void);
26
27    static RNG* instance(int seed);
28
29    inline const gsl_rng* rng(void) const { return rng_; }
30
31  private:
32
33    RNG(void);
34
35    static RNG* instance_;
36    gsl_rng* rng_;
37  };
38
39
40
41  ///
42  /// @brief continuous random distributions.
43  ///
44  /// Abstract base class for continuous random distributions.
45  ///
46  class RandomContinuous
47  {
48  public:
49
50    ///
51    /// @brief Constructor
52    ///
53    inline RandomContinuous(void) { rng_=RNG::instance(89); }
54
55    ///
56    /// @return A random number
57    ///
58    virtual double operator()(void) const = 0;
59
60  protected:
61    RNG* rng_;
62  };
63
64
65  ///
66  /// @brief Uniform distribution
67  ///
68  /// Class for generating a random number from a uniform distribution
69  /// between zero and unity.
70  ///
71  /// Distribution function \f$f(x) = 1 \f$ for \f$0 \le x < 1 \f$ \n
72  /// Expectation value: 0.5 \n
73  /// Variance: \f$\frac{1}{12} \f$
74  ///
75  class RandomContinuousUniform : public RandomContinuous
76  {
77  public:
78
79    inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());}
80
81  };
82
83  ///
84  /// @brief Gaussian distribution
85  ///
86  /// Class for generating a random number from a Gaussian
87  /// distribution between zero and unity. Utilizes the Box-Muller
88  /// algorithm, which needs two calls to random generator.
89  ///
90  /// Distribution function \f$f(x) = 91 /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2}) 92 /// \f$ \n
93  /// Expectation value: \f$\mu \f$ \n
94  /// Variance: \f$\sigma^2 \f$
95  ///
96  class RandomGaussian : public RandomContinuous
97  {
98  public:
99    ///
100    /// @brief Constructor
101    /// @param s is the standard deviation \f$\sigma \f$ of distribution
102    /// m is the expectation value \f$\mu \f$ of the distribution
103    ///
104    inline RandomGaussian(const double s=1, const double m=0)
105      : m_(m), s_(s) {}
106
107    ///
108    /// @return A random Gaussian number
109    ///
110    inline double operator()(void) const
111    { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }
112
113    ///
114    /// @return A random Gaussian number with standard deviation \a s
115    /// and expectation value 0.
116    /// @note this operator ignores parameters given in Constructor
117    ///
118    inline double operator()(const double s) const
119    { return gsl_ran_gaussian(rng_->rng(), s); }
120
121    ///
122    /// @return A random Gaussian number with standard deviation \a s
123    /// and expectation value \a m.
124    /// @note this operator ignores parameters given in Constructor
125    ///
126    inline double operator()(const double s, const double m) const
127    { return gsl_ran_gaussian(rng_->rng(), s)+m; }
128
129  private:
130    double m_;
131    double s_;
132  };
133
134  ///
135  /// @brief Exponential distribution
136  ///
137  /// Class for generating a random number from a Exponential
138  /// distribution.
139  ///
140  /// Distribution function \f$f(x) = \frac{1}{m}\exp(-x/a) \f$ for
141  /// \f$x \f$ \n
142  /// Expectation value: \f$m \f$ \n
143  /// Variance: \f$m^2 \f$
144  ///
145  class RandomExponential : public RandomContinuous
146  {
147  public:
148    ///
149    /// @brief Constructor
150    /// @param m is the expectation value of the distribution.
151    ///
152    inline RandomExponential(const double m=1)
153      : m_(m){}
154
155    ///
156    /// @return A random number from exponential distribution
157    ///
158    inline double operator()(void) const
159    { return gsl_ran_exponential(rng_->rng(), m_); }
160
161    ///
162    /// @return A random number from exponential distribution, with
163    /// expectation value \a m
164    /// @note this operator ignores parameters given in Constructor
165    ///
166    inline double operator()(const double m) const
167    { return gsl_ran_exponential(rng_->rng(), m); }
168
169  private:
170    double m_;
171  };
172
173
174  ///
175  /// @brief discrete random distributions.
176  ///
177  /// Abstract Base Class for discrete random distributions. Given K
178  /// discrete events with different probabilities \f$P[k] \f$,
179  /// produce a random value k consistent with its probability.
180  ///
181  class RandomDiscrete
182  {
183  public:
184    ///
185    /// @brief Constructor
186    ///
187    inline RandomDiscrete(void) { rng_=RNG::instance(89); }
188
189    ///
190    /// @return A random number.
191    ///
192    virtual long operator()(void) const = 0;
193
194  protected:
195    RNG* rng_;
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 \f$[a,b) \f$, have equal probability to occur.
204  ///
205  /// Distribution function \f$p(k) = \frac{1}{b-a} \f$ for \f$a \le 206 /// k < b \f$ \n
207  /// Expectation value: \f$\frac{a+b}{2} \f$ \n
208  /// Variance: \f$\frac{1}{3}((b-a)^2-1) \f$
209  ///
210
211  class RandomDiscreteUniform : public RandomDiscrete
212  {
213  public:
214    ///
215    /// @brief Constructor.
216    ///
217    /// @param \a range is number of different integers that can be
218    /// generated \f$(b-a+1) \f$ \a min is the minimal integer that
219    /// can be generated \f$(a) \f$
220    ///
221    inline RandomDiscreteUniform(const u_long range, const long min=0)
222      : min_(min), range_(range) {}
223
224    ///
225    /// @return A random number between \a min_ and \a range_ - \a min_
226    ///
227    inline long operator()(void) const
228    { return gsl_rng_uniform_int(rng_->rng(), range_)+min_; }
229
230    ///
231    /// @return A random number between 0 and \a range.
232    ///
233    inline long operator()(const u_long range) const
234    { return gsl_rng_uniform_int(rng_->rng(), range); }
235
236    ///
237    /// @return A random number between \a min and \a range - \a min
238    /// @note this operator ignores parameters given in Constructor
239    ///
240    inline long operator()(const u_long range, const long min) const
241    { return gsl_rng_uniform_int(rng_->rng(), range)+min; }
242
243
244  private:
245    long min_;
246    u_long range_;
247  };
248
249
250  ///
251  /// @brief Poisson Distribution
252  ///
253  /// Having a Poisson process (no memory), number of occurences
254  /// within a given time window is Poisson distributed. This
255  /// distribution is the limit of a Binomial distribution when number
256  /// of attempts is large, and the probability for one attempt to be
257  /// succesful is small (in such a way that the expected number of
258  /// succesful attempts is \f$m \f$.
259  ///
260  /// Probability function \f$p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$0 \le 261 /// k \f$ \n
262  /// Expectation value: \f$m \f$ \n
263  /// Variance: \f$m \f$
264  ///
265  class RandomPoisson : public RandomDiscrete
266  {
267  public:
268    ///
269    /// @brief Constructor
270    ///
271    /// @param m is expectation value
272    inline RandomPoisson(const double m=1)
273      : m_(m){}
274
275    ///
276    /// @return a Poisson distributed number.
277    ///
278    inline double operator()(void) const
279    { return gsl_ran_poisson(rng_->rng(), m_); }
280
281    ///
282    /// @return a Poisson distributed number with expectation value \a
283    /// m
284    /// @note this operator ignores parameters set in Constructor
285    inline double operator()(const double m) const
286    { return gsl_ran_poisson(rng_->rng(), m); }
287
288  private:
289    double m_;
290  };
291
292  ///
293  /// @brief General
294  ///
295  class RandomDiscreteGeneral : public RandomDiscrete {
296  };
297
298
299}} // of namespace random and namespace theplu
300
301#endif
Note: See TracBrowser for help on using the repository browser.