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

Last change on this file since 3579 was 3579, checked in by Peter, 5 years ago

merge release 0.14 into trunk

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