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

Last change on this file since 865 was 865, checked in by Peter, 14 years ago

changing URL to http://trac.thep.lu.se/trac/yat

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