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

Last change on this file since 3591 was 3591, checked in by Peter, 6 years ago

fixes #877

Generalize YAT_CXX_RVALUE into YAT_CXX_TRY_CXX11 and use that to
implement new macro YAT_CXX_ATOMIC. Use this macro in configure.ac and
use std::atomic<> conditoionally in RNG singleton.

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