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

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

Fixes #8 and #179, addresses #2.

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