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

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

refs #803; let std::random_shuffle do the concept check

  • 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 3518 2016-10-05 08:01:11Z 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 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  };
203
204
205  ///
206  /// @brief Class holding state of a random generator
207  ///
208  class RNG_state
209  {
210  public:
211    ///
212    /// @brief Constructor
213    ///
214    explicit RNG_state(const RNG*);
215
216    /**
217       Copy Constructor
218
219       \since Explicitely declared since yat 0.5
220     */
221    RNG_state(const RNG_state&);
222
223    ///
224    /// @brief Destructor
225    ///
226    ~RNG_state(void);
227
228    ///
229    /// @return const pointer to underlying GSL random generator.
230    ///
231    const gsl_rng* rng(void) const;
232
233    /**
234       Assignment operator
235
236       \since Explicitely declared since yat 0.5
237     */
238    RNG_state& operator=(const RNG_state&);
239
240  private:
241    gsl_rng* rng_;
242
243    void clone(const gsl_rng&);
244  };
245
246
247  // --------------------- Discrete distribtuions ---------------------
248
249  ///
250  /// @brief Discrete random number distributions.
251  ///
252  /// Abstract base class for discrete random number
253  /// distributions. Given K discrete events with different
254  /// probabilities \f$ P[k] \f$, produce a random value k consistent
255  /// with its probability.
256  ///
257  class Discrete
258  {
259  public:
260    /**
261       type returned by operator()
262
263       \since New in yat 0.10
264     */
265    typedef unsigned long int result_type;
266
267    ///
268    /// @brief Constructor
269    ///
270    Discrete(void);
271
272    ///
273    /// @brief The destructor
274    ///
275    virtual ~Discrete(void);
276
277    ///
278    /// @brief Set the seed to \a s.
279    ///
280    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
281    /// default value from the rng's original implementation is used
282    /// (cf. GSL documentation).
283    ///
284    /// \deprecated Provided for backward compatibility with the 0.7
285    /// API. Use RNG::instance()->seed(s) instead.
286    ///
287    void seed(unsigned long s) const YAT_DEPRECATE;
288
289    ///
290    /// @brief Set the seed using the /dev/urandom device.
291    ///
292    /// @return The seed acquired from /dev/urandom.
293    ///
294    /// \deprecated Provided for backward compatibility with the 0.7
295    /// API. Use RNG::instance()->seed_from_devurandom() instead.
296    ///
297    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
298
299    ///
300    /// @return A random number.
301    ///
302    virtual result_type operator()(void) const = 0;
303
304  protected:
305    /// GSL random gererator
306    RNG* rng_;
307  };
308
309  /**
310     \brief Binomial distribution
311
312     \see gsl_ran_binomial
313
314     \since New in yat 0.10
315   */
316  class Binomial : public Discrete
317  {
318  public:
319    /**
320       \brief Constructor
321
322       Create an object that generates random numbers from binomial
323       distribution, the number of successes in \n trials each with
324       success probability \a p.
325     */
326    Binomial(double p, unsigned int n);
327
328    /**
329       \return number from binomial distrubtion as parametrized in constructor.
330     */
331    unsigned long operator()(void) const;
332  private:
333    double p_;
334    unsigned int n_;
335  };
336
337  ///
338  /// @brief General
339  ///
340  class DiscreteGeneral : public Discrete
341  {
342  public:
343    ///
344    /// @brief Constructor
345    ///
346    /// @param hist is a Histogram defining the probability distribution
347    ///
348    explicit DiscreteGeneral(const statistics::Histogram& hist);
349
350    /**
351       \brief Copy constructor
352
353       \since Explicitely implemented in yat 0.5
354     */
355    DiscreteGeneral(const DiscreteGeneral&);
356
357    ///
358    /// @brief Destructor
359    ///
360    ~DiscreteGeneral(void);
361
362    ///
363    /// The generated number is an integer and proportional to the
364    /// frequency in the corresponding histogram bin. In other words,
365    /// the probability that 0 is returned is proportinal to the size
366    /// of the first bin.
367    ///
368    /// @return A random number.
369    ///
370    unsigned long operator()(void) const;
371
372    /**
373       \brief Assignment operator
374
375       \since Explicitely implemented in yat 0.5
376     */
377    DiscreteGeneral& operator=(const DiscreteGeneral&);
378
379  private:
380    void free(void);
381    void preproc(void);
382
383    gsl_ran_discrete_t* gen_;
384    std::vector<double> p_;
385  };
386
387  /**
388     @brief Discrete uniform distribution
389
390     Discrete uniform distribution also known as the "equally likely
391     outcomes" distribution. Each outcome, in this case an integer
392     from [0,n-1] , have equal probability to occur.
393
394     Distribution function \f$ p(k) = \frac{1}{n+1} \f$ for \f$ 0 \le
395     k < n \f$ \n
396     Expectation value: \f$ \frac{n-1}{2} \f$ \n
397     Variance: \f$ \frac{1}{12}(n-1)(n+1) \f$
398  */
399  class DiscreteUniform
400    : public Discrete,
401      public std::unary_function<unsigned long, unsigned long>
402  {
403  public:
404    /**
405       \brief Constructor.
406
407       The generator will generate integers within the range \f$
408       [0,n-1] \f$. If \a n is zero, then the whole range of the
409       underlying RNG will be used \f$ [min,max] \f$. Setting \a n to
410       zero is the preferred way to sample the whole range of the
411       underlying RNG, i.e. not setting \n to RNG.max.
412
413       \throw If \a n is larger than the maximum number the underlying
414       random number generator can return, then a GSL_error exception
415       is thrown.
416    */
417    explicit DiscreteUniform(unsigned long n=0);
418
419    /**
420       \brief Get a random number
421
422       The returned integer is either in the range [RNG.min,RNG.max]
423       or [0,n-1] depending on how the random number generator was
424       created.
425
426       \see DiscreteUniform(const unsigned long n=0)
427    */
428    unsigned long operator()(void) const;
429
430    /**
431       \brief Get a random integer in the range \f$ [0,n-1] \f$.
432
433       All integers in the range [0,n-1] are equally likely. This
434       function should be avoided for sampling the whole range of the
435       underlying RNG.
436
437       \throw GSL_error if \a n is larger than the range of the
438       underlying generator.
439    */
440    unsigned long operator()(unsigned long n) const;
441
442  private:
443    unsigned long range_;
444  };
445
446
447  /**
448     @brief Geomtric Distribution
449
450     The number of independent trials with probability \em p until the
451     first success.
452
453     Probability function \f$ p(k) = p (1-p)^(k-1) \f$
454
455     \since New in yat 0.14
456  */
457  class Geometric : public Discrete
458  {
459  public:
460    /**
461       \brief Constructor
462
463       \param p is probability for success in one trial
464    */
465    Geometric(double p);
466
467    /*
468      \return a number from Geomtric distribution
469    */
470    unsigned long operator()(void) const;
471
472    /**
473       \return a number from Geomtric distribution with success rate \a p
474
475       \note this operator ignores parameters set in Constructor
476    */
477    unsigned long operator()(double p) const;
478
479  private:
480    double p_;
481  };
482
483
484  /**
485     If we have \a n1 samples of type 1 and \a n2 samples of type 2
486     and draw \a t samples with replacement, number of drawn samples
487     of type 1 will follow the hyper geometric distribution.
488
489     \since New in yat 0.14
490   */
491  class HyperGeometric : public Discrete
492  {
493  public:
494    /**
495       \brief Defaul constructor
496     */
497    HyperGeometric(void);
498
499    /**
500       \brief Constructor
501       \param n1 number of samples of type 1
502       \param n2 number of samples of type 2
503       \param t number of samples to draw
504     */
505    HyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
506
507    /**
508       \return random number from hypergeometric distribution using
509       parameters set in constructor.
510     */
511    unsigned long int operator()(void) const;
512
513    /**
514       \return random number from hypergeometric distribution using
515       parameters passed.
516     */
517    unsigned long int operator()(unsigned int n1, unsigned int n2,
518                                 unsigned int t) const;
519  private:
520    unsigned int n1_;
521    unsigned int n2_;
522    unsigned int t_;
523  };
524
525  /**
526     We have \a n1 samples of type 1 and \a n2 samples of type 2.
527     Samples are drawn with replacement until \a t samles of type 2
528     are drawn. Then \a k, number of drawn samples of type 1, follows
529     negative hypergeometric distribution.
530
531     \since New in yat 0.14
532
533     \see HyperGeometric
534   */
535  class NegativeHyperGeometric : public Discrete
536  {
537  public:
538    /**
539       \brief Default constructor
540     */
541    NegativeHyperGeometric(void);
542
543    /**
544       \brief Constructor
545       \param n1 number of samples of type 1
546       \param n2 number of samples of type 2
547       \param t number of samples of type 2 to draw
548     */
549    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);
550
551    /**
552       \return random number from negative hypergeometric distribution
553       using parameters set in constructor.
554     */
555    unsigned long int operator()(void) const;
556
557    /**
558       \return random number from negative hypergeometric distribution
559       using parameters passed.
560     */
561    unsigned long int operator()(unsigned int n1, unsigned int n2,
562                                 unsigned int t) const;
563  private:
564    unsigned int n1_;
565    unsigned int n2_;
566    unsigned int t_;
567  };
568
569
570  /**
571     @brief Poisson Distribution
572
573     Having a Poisson process (i.e. no memory), number of occurences
574     within a given time window is Poisson distributed. This
575     distribution is the limit of a Binomial distribution when number
576     of attempts is large, and the probability for one attempt to be
577     succesful is small (in such a way that the expected number of
578     succesful attempts is \f$ m \f$.
579
580     Probability function \f$ p(k) = e^{-m}\frac{m^k}{k!} \f$ for \f$ 0 \le
581     k  \f$ \n
582     Expectation value: \f$ m \f$ \n
583     Variance: \f$ m \f$
584  */
585  class Poisson : public Discrete
586  {
587  public:
588    ///
589    /// @brief Constructor
590    ///
591    /// @param m is expectation value
592    ///
593    explicit Poisson(const double m=1);
594
595    ///
596    /// @return A Poisson distributed number.
597    ///
598    unsigned long operator()(void) const;
599
600    ///
601    /// @return A Poisson distributed number with expectation value
602    /// \a m
603    ///
604    /// @note this operator ignores parameters set in Constructor
605    ///
606    unsigned long operator()(const double m) const;
607
608  private:
609    double m_;
610  };
611
612  // --------------------- Continuous distribtuions ---------------------
613
614  ///
615  /// @brief Continuous random number distributions.
616  ///
617  /// Abstract base class for continuous random number distributions.
618  ///
619  class Continuous
620  {
621  public:
622    /**
623       type returned by operator()
624
625       \since New in yat 0.10
626     */
627    typedef double result_type;
628
629    ///
630    /// @brief Constructor
631    ///
632    Continuous(void);
633
634    ///
635    /// @brief The destructor
636    ///
637    virtual ~Continuous(void);
638
639    ///
640    /// @brief Set the seed to \a s.
641    ///
642    /// Set the seed to \a s in the underlying rng. If \a s is zero, a
643    /// default value from the rng's original implementation is used
644    /// (cf. GSL documentation).
645    ///
646    /// \deprecated Provided for backward compatibility with the 0.7
647    /// API. Use RNG::instance()->seed(s) instead.
648    ///
649    void seed(unsigned long s) const YAT_DEPRECATE;
650
651    ///
652    /// @brief Set the seed using the /dev/urandom device.
653    ///
654    /// @return The seed acquired from /dev/urandom.
655    ///
656    /// \deprecated Provided for backward compatibility with the 0.7
657    /// API. Use RNG::instance()->seed_from_devurandom() instead.
658    ///
659    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;
660
661    ///
662    /// @return A random number
663    ///
664    virtual result_type operator()(void) const = 0;
665
666  protected:
667    /// pointer to GSL random generator
668    RNG* rng_;
669  };
670
671  // ContinuousUniform is declared before ContinuousGeneral to avoid
672  // forward declaration
673  ///
674  /// @brief Uniform distribution
675  ///
676  /// Class for generating a random number from a uniform distribution
677  /// in the range [0,1), i.e. zero is included but not 1.
678  ///
679  /// Distribution function \f$ f(x) = 1 \f$ for \f$ 0 \le x < 1 \f$ \n
680  /// Expectation value: 0.5 \n
681  /// Variance: \f$ \frac{1}{12} \f$
682  ///
683  class ContinuousUniform : public Continuous
684  {
685  public:
686    double operator()(void) const;
687  };
688
689  ///
690  /// @brief Generates numbers from a histogram in a continuous manner.
691  ///
692  class ContinuousGeneral : public Continuous
693  {
694  public:
695    ///
696    /// @brief Constructor
697    ///
698    /// @param hist is a Histogram defining the probability distribution
699    ///
700    explicit ContinuousGeneral(const statistics::Histogram& hist);
701
702    ///
703    /// The number is generated in a two step process. First the bin
704    /// in the histogram is randomly selected (see
705    /// DiscreteGeneral). Then a number is generated uniformly from
706    /// the interval defined by the bin.
707    ///
708    /// @return A random number.
709    ///
710    double operator()(void) const;
711
712  private:
713    const DiscreteGeneral discrete_;
714    const statistics::Histogram hist_;
715    ContinuousUniform u_;
716  };
717
718  /**
719     \brief Generator of random numbers from an exponential
720     distribution.
721
722     The distribution function is \f$ f(x) = \frac{1}{m}\exp(-x/a)
723     \f$ for \f$ x \f$ with the expectation value \f$ m \f$ and
724     variance \f$ m^2 \f$
725  */
726  class Exponential : public Continuous
727  {
728  public:
729    ///
730    /// @brief Constructor
731    ///
732    /// @param m is the expectation value of the distribution.
733    ///
734    explicit Exponential(const double m=1);
735
736    ///
737    /// @return A random number from exponential distribution.
738    ///
739    double operator()(void) const;
740
741    ///
742    /// @return A random number from exponential distribution, with
743    /// expectation value \a m
744    ///
745    /// @note This operator ignores parameters given in constructor.
746    ///
747    double operator()(const double m) const;
748
749  private:
750    double m_;
751  };
752
753  /**
754     @brief Gaussian distribution
755
756     Class for generating a random number from a Gaussian distribution
757     between zero and unity. Utilizes the Box-Muller algorithm, which
758     needs two calls to random generator.
759
760     Distribution function \f$ f(x) =
761     \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})
762     \f$ \n
763     Expectation value: \f$ \mu \f$ \n
764     Variance: \f$ \sigma^2 \f$
765  */
766  class Gaussian : public Continuous
767  {
768  public:
769    ///
770    /// @brief Constructor
771    ///
772    /// @param s is the standard deviation \f$ \sigma \f$ of distribution
773    /// @param m is the expectation value \f$ \mu \f$ of the distribution
774    ///
775    explicit Gaussian(const double s=1, const double m=0);
776
777    ///
778    /// @return A random Gaussian number
779    ///
780    double operator()(void) const;
781
782    ///
783    /// @return A random Gaussian number with standard deviation \a s
784    /// and expectation value 0.
785    ///
786    /// @note this operator ignores parameters given in Constructor
787    ///
788    double operator()(const double s) const;
789
790    ///
791    /// @return A random Gaussian number with standard deviation \a s
792    /// and expectation value \a m.
793    ///
794    /// @note this operator ignores parameters given in Constructor
795    ///
796    double operator()(const double s, const double m) const;
797
798  private:
799    double m_;
800    double s_;
801  };
802
803  /**
804     \brief Convenience function to shuffle a range with singleton RNG.
805
806     Wrapper around std::random_shuffle using DiscreteUniform as
807     random generator and thereby using the underlying RNG class,
808     which is singleton.
809
810     Type Requirements:
811     - RandomAccessIterator is \random_access_iterator
812   */
813  template<typename RandomAccessIterator>
814  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
815  {
816    DiscreteUniform rnd;
817    std::random_shuffle(first, last, rnd);
818  }
819
820}}} // of namespace random, yat, and theplu
821
822#endif
Note: See TracBrowser for help on using the repository browser.