source: branches/0.9-stable/yat/random/random.h @ 2833

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

fixes #720

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