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

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

Addresses #170.

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