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

Last change on this file since 2804 was 2804, checked in by Peter, 9 years ago

RNG::set_state: prefer throw on error.

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