source: trunk/yat/random/random.h @ 676

Last change on this file since 676 was 675, checked in by Jari Häkkinen, 15 years ago

References #83. Changing project name to yat. Compilation will fail in this revision.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.3 KB
Line 
1#ifndef _theplu_random_
2#define _theplu_random_
3
4// $Id: random.h 675 2006-10-10 12:08:45Z jari $
5
6/*
7  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson
8
9  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "yat/statistics/Histogram.h"
28
29#include <gsl/gsl_rng.h>
30#include <gsl/gsl_randist.h>
31
32#include <string>
33
34namespace theplu {
35namespace random {
36
37  //forward declarion
38  class RNG_state;
39
40  ///
41  /// @brief Random Number Generator
42  ///
43  /// The RNG class is wrapper to the GSL random number generator
44  /// (rng). This class provides a single global instance of the rng,
45  /// and makes sure there is only one point of access to the
46  /// generator.
47  ///
48  /// There is information about how to change seeding and generators
49  /// at run time without recompilation using environment
50  /// variables. RNG of course support seeding at compile time if you
51  /// don't want to bother about environment variables and GSL.
52  ///
53  /// There are many different rng's available in GSL. Currently only
54  /// the default generator is implemented and no other one is
55  /// choosable through the class interface. This means that you have
56  /// to fall back to the use of environment variables as described in
57  /// the GSL documentation, or be bold and request support for other
58  /// rng's through the class interface.
59  ///
60  /// Not all GSL functionality is implemented, we'll add
61  /// functionality when needed and may do it when requested. Better
62  /// yet, supply us with code and we will probably add it to the code
63  /// (BUT remember to implement reasonable tests for your code and
64  /// follow the coding style.)
65  ///
66  /// This implementation may be thread safe (according to GSL
67  /// documentation), but should be checked to be so before trusting
68  /// thread safety.
69  ///
70  /// @see Design Patterns (the singleton and adapter pattern). GSL
71  /// documentation.
72  ///
73  class RNG
74  {
75  public:
76
77    virtual ~RNG(void);
78
79    ///
80    /// @brief Get an instance of the random number generator.
81    ///
82    /// Get an instance of the random number generator. If the random
83    /// number generator is not already created, the call will create
84    /// a new generator and use the default seed. The seed must be
85    /// changed with the seed or seed_from_devurandom member
86    /// functions.
87    ///
88    /// @return A pointer to the random number generator.
89    ///
90    /// @see seed and seed_from_devurandom
91    ///
92    static RNG* instance(void)
93      { if (!instance_) instance_=new RNG; return instance_; }
94
95    ///
96    /// @brief Returns the largest number that the random number
97    /// generator can return.
98    ///
99    inline u_long max(void) const { return gsl_rng_max(rng_); }
100
101    ///
102    /// @brief Returns the smallest number that the random number
103    /// generator can return.
104    ///
105    inline u_long min(void) const { return gsl_rng_min(rng_); }
106
107    ///
108    /// @brief Returns the name of the random number generator
109    ///
110    inline std::string name(void) const { return gsl_rng_name(rng_); }
111
112    ///
113    /// @return const pointer to underlying GSL random generator.
114    ///
115    inline const gsl_rng* rng(void) const { return rng_; }
116
117    ///
118    /// @brief Set the seed \a s for the rng.
119    ///
120    /// Set the seed \a s for the rng. If \a s is zero, a default
121    /// value from the rng's original implementation is used (cf. GSL
122    /// documentation).
123    ///
124    /// @see seed_from_devurandom
125    ///
126    inline void seed(u_long s) const { gsl_rng_set(rng_,s); }
127
128    ///
129    /// @brief Seed the rng using the /dev/urandom device.
130    ///
131    /// @return The seed acquired from /dev/urandom.
132    ///
133    u_long seed_from_devurandom(void);
134
135    ///
136    /// @brief set the state
137    ///
138    /// @return see gsl_rng_memcpy
139    ///
140    int set_state(const RNG_state&);
141
142  private:
143    RNG(void);
144
145    static RNG* instance_;
146    gsl_rng* rng_;
147  };
148
149
150  ///
151  /// @brief Class holding state of a random generator
152  ///
153  class RNG_state
154  {
155  public:
156    ///
157    /// @brief Constructor
158    ///
159    RNG_state(const RNG*);
160
161    ///
162    /// @brief Destructor
163    ///
164    ~RNG_state(void);
165
166    ///
167    /// @return const pointer to underlying GSL random generator.
168    ///
169    inline const gsl_rng* rng(void) const { return rng_; }
170
171
172  private:
173    gsl_rng* rng_;
174
175  };
176   
177
178  // --------------------- Discrete distribtuions ---------------------
179
180  ///
181  /// @brief Discrete random number distributions.
182  ///
183  /// Abstract base class for discrete random number
184  /// distributions. Given K discrete events with different
185  /// probabilities \f$ P[k] \f$, produce a random value k consistent
186  /// with its probability.
187  ///
188  class Discrete
189  {
190  public:
191    ///
192    /// @brief Constructor
193    ///
194    inline Discrete(void) { rng_=RNG::instance(); }
195
196    inline virtual ~Discrete(void) { }
197
198    ///
199    /// @brief Set the seed to \a s.
200    ///
201    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
202    /// default value from the rng's original implementation is used
203    /// (cf. GSL documentation).
204    ///
205    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
206    ///
207    inline void seed(u_long s) const { rng_->seed(s); }
208
209    ///
210    /// @brief Set the seed using the /dev/urandom device.
211    ///
212    /// @return The seed acquired from /dev/urandom.
213    ///
214    /// @see seed, RNG::seed_from_devurandom, RNG::seed
215    ///
216    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
217
218    ///
219    /// @return A random number.
220    ///
221    virtual u_long operator()(void) const = 0;
222   
223  protected:
224    /// GSL random gererator
225    RNG* rng_;
226  };
227
228  ///
229  /// @brief General
230  ///
231  class DiscreteGeneral : public Discrete
232  {
233  public:
234    ///
235    /// @brief Constructor
236    ///
237    /// @param hist is a Histogram defining the probability distribution
238    ///
239    DiscreteGeneral(const statistics::Histogram& hist);
240   
241    ///
242    /// @brief Destructor
243    ///
244    ~DiscreteGeneral(void);
245
246    ///
247    /// The generated number is an integer and proportinal to the
248    /// frequency in the corresponding histogram bin. In other words,
249    /// the probability that 0 is returned is proportinal to the size
250    /// of the first bin.
251    ///
252    /// @return A random number.
253    ///
254    inline u_long
255    operator()(void) const { return gsl_ran_discrete(rng_->rng(), gen_); }
256
257  private:
258     gsl_ran_discrete_t* gen_;
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  class DiscreteUniform : public Discrete
274  {
275  public:
276    ///
277    /// @brief Default constructor.
278    ///
279    DiscreteUniform(void) : range_(rng_->max()) {}
280
281    ///
282    /// @brief Constructor.
283    ///
284    /// The generator will generate integers from \f$ [0,n-1] \f$. If
285    /// \a n is larger than the maximum number the random number
286    /// generator can return, then (currently) \a n is adjusted
287    /// appropriately.
288    ///
289    /// @todo If a too large \a n is given an exception should be
290    /// thrown, i.e. the behaviour of this class will change. The case
291    /// when argument is 0 is not treated gracefully (underlying GSL
292    /// functionality will not return).
293    ///
294    DiscreteUniform(const u_long n) : range_(n)
295    { if ( range_>rng_->max() ) range_=rng_->max(); }
296
297    ///
298    /// This function returns a random integer from 0 to n-1
299    /// inclusive. All integers in the range [0,n-1] are equally
300    /// likely. n is set in constructor.
301    ///
302    inline u_long
303    operator()(void) const { return gsl_rng_uniform_int(rng_->rng(), range_); }
304
305    ///
306    /// This function returns a random integer from 0 to n-1
307    /// inclusive. All integers in the range [0,n-1] are equally
308    /// likely.
309    ///
310    inline u_long operator()(const u_long n) const
311    { return gsl_rng_uniform_int(rng_->rng(), n); }
312
313  private:
314    u_long range_;
315  };
316
317  ///
318  /// @brief Poisson Distribution
319  ///
320  /// Having a Poisson process (i.e. no memory), number of occurences
321  /// within a given time window is Poisson distributed. This
322  /// distribution is the limit of a Binomial distribution when number
323  /// of attempts is large, and the probability for one attempt to be
324  /// succesful is small (in such a way that the expected number of
325  /// succesful attempts is \f$ m \f$.
326  ///
327  /// Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
328  /// k  \f$ \n
329  /// Expectation value: \f$ m \f$ \n
330  /// Variance: \f$ m \f$
331  ///
332  class Poisson : public Discrete
333  {
334  public:
335    ///
336    /// @brief Constructor
337    ///
338    /// @param m is expectation value
339    inline Poisson(const double m=1) : m_(m) {}
340
341    ///
342    /// @return A Poisson distributed number.
343    ///
344    inline u_long
345    operator()(void) const { return gsl_ran_poisson(rng_->rng(), m_); }
346
347    ///
348    /// @return A Poisson distributed number with expectation value \a
349    /// m
350    ///
351    /// @note this operator ignores parameters set in Constructor
352    ///
353    inline u_long
354    operator()(const double m) const { return gsl_ran_poisson(rng_->rng(), m); }
355
356  private:
357    double m_;
358  };
359
360  // --------------------- Continuous distribtuions ---------------------
361
362  ///
363  /// @brief Continuous random number distributions.
364  ///
365  /// Abstract base class for continuous random number distributions.
366  ///
367  class Continuous
368  {
369  public:
370
371    ///
372    /// @brief Constructor
373    ///
374    inline Continuous(void) { rng_=RNG::instance(); }
375
376    inline virtual ~Continuous(void) { }
377
378    ///
379    /// @brief Set the seed to \a s.
380    ///
381    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
382    /// default value from the rng's original implementation is used
383    /// (cf. GSL documentation).
384    ///
385    /// @see seed_from_devurandom, RNG::seed_from_devurandom, RNG::seed
386    ///
387    inline void seed(u_long s) const { rng_->seed(s); }
388
389    ///
390    /// @brief Set the seed using the /dev/urandom device.
391    ///
392    /// @return The seed acquired from /dev/urandom.
393    ///
394    /// @see seed, RNG::seed_from_devurandom, RNG::seed
395    ///
396    u_long seed_from_devurandom(void) { return rng_->seed_from_devurandom(); }
397
398    ///
399    /// @return A random number
400    ///
401    virtual double operator()(void) const = 0;
402
403  protected:
404    /// pointer to GSL random generator
405    RNG* rng_;
406  };
407
408  ///
409  /// @brief Uniform distribution
410  ///
411  /// Class for generating a random number from a uniform distribution
412  /// in the range [0,1), i.e. zero is included but not 1.
413  ///
414  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
415  /// Expectation value: 0.5 \n
416  /// Variance: \f$ \frac{1}{12} \f$
417  ///
418  class ContinuousUniform : public Continuous
419  {
420  public:
421    inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());}
422  };
423
424  ///
425  /// Class to generate numbers from a histogram in a continuous manner.
426  ///
427  class ContinuousGeneral : public Continuous
428  {
429  public:
430    ///
431    /// @brief Constructor
432    ///
433    /// @param hist is a Histogram defining the probability distribution
434    ///
435    inline ContinuousGeneral(const statistics::Histogram& hist)
436      : discrete_(DiscreteGeneral(hist)), hist_(hist) {}
437
438    ///
439    /// @brief Destructor
440    ///
441    ~ContinuousGeneral(void);
442
443    ///
444    /// The number is generated in a two step process. First the bin
445    /// in the histogram is randomly selected (see
446    /// DiscreteGeneral). Then a number is generated uniformly from
447    /// the interval defined by the bin.
448    ///
449    /// @return A random number.
450    ///
451    inline double operator()(void) const 
452    { return hist_.observation_value(discrete_())+(u_()-0.5)*hist_.spacing(); }
453
454  private:
455    const DiscreteGeneral discrete_;
456    const statistics::Histogram& hist_;
457    ContinuousUniform u_;
458  };
459
460  ///
461  /// @brief Generator of random numbers from an exponential
462  /// distribution.
463  ///
464  /// The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
465  /// \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
466  /// variance \f$ m^2 \f$
467  ///
468  class Exponential : public Continuous
469  {
470  public:
471    ///
472    /// @brief Constructor
473    ///
474    /// @param m is the expectation value of the distribution.
475    ///
476    inline Exponential(const double m=1) : m_(m) {}
477
478    ///
479    /// @return A random number from exponential distribution.
480    ///
481    inline double
482    operator()(void) const { return gsl_ran_exponential(rng_->rng(), m_); }
483
484    ///
485    /// @return A random number from exponential distribution, with
486    /// expectation value \a m
487    ///
488    /// @note This operator ignores parameters given in constructor.
489    ///
490    inline double operator()(const double m) const
491    { return gsl_ran_exponential(rng_->rng(), m); }
492
493  private:
494    double m_;
495  };
496
497  ///
498  /// @brief Gaussian distribution
499  ///
500  /// Class for generating a random number from a Gaussian
501  /// distribution between zero and unity. Utilizes the Box-Muller
502  /// algorithm, which needs two calls to random generator.
503  ///
504  /// Distribution function \f$ f(x) =
505  /// \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
506  /// \f$ \n
507  /// Expectation value: \f$ \mu \f$ \n
508  /// Variance: \f$ \sigma^2 \f$
509  ///
510  class Gaussian : public Continuous
511  {
512  public:
513    ///
514    /// @brief Constructor
515    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
516    /// @param m is the expectation value \f$ \mu \f$ of the distribution
517    ///
518    inline Gaussian(const double s=1, const double m=0) : m_(m), s_(s) {}
519
520    ///
521    /// @return A random Gaussian number
522    ///
523    inline double
524    operator()(void) const { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }
525
526    ///
527    /// @return A random Gaussian number with standard deviation \a s
528    /// and expectation value 0.
529    /// @note this operator ignores parameters given in Constructor
530    ///
531    inline double
532    operator()(const double s) const { return gsl_ran_gaussian(rng_->rng(), s); }
533
534    ///
535    /// @return A random Gaussian number with standard deviation \a s
536    /// and expectation value \a m.
537    /// @note this operator ignores parameters given in Constructor
538    ///
539    inline double operator()(const double s, const double m) const
540    { return gsl_ran_gaussian(rng_->rng(), s)+m; }
541
542  private:
543    double m_;
544    double s_;
545  };
546
547
548}} // of namespace random and namespace theplu
549
550#endif
Note: See TracBrowser for help on using the repository browser.