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

Last change on this file since 2268 was 2268, checked in by Peter, 12 years ago

make DiscreteUniform? adaptable. closes #629

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