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

Last change on this file since 1614 was 1614, checked in by Peter, 13 years ago

implementing copying in RNG_state. Previously compiler generated versions were used which was not sound. Also disallowed assignment of singleton RNG. Assignment of a singleton class does not make sense.

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