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

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

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