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

Last change on this file since 1334 was 1275, checked in by Jari Häkkinen, 13 years ago

Updating copyright statements.

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