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

Last change on this file since 831 was 831, checked in by Peter, 15 years ago

Refs #185.

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