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

Last change on this file since 2145 was 2145, checked in by Peter, 12 years ago

improve compile error message in random_shuffle. refs #489

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