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

Last change on this file since 2770 was 2770, checked in by Peter, 10 years ago

change back to RNG being a singleton, which now stores gsl_rng* in a thread_specific_ptr<gsl_rng>. refs #568

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