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

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

fixes #762. workaround for broken boost header (in v1.41)

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