source: trunk/lib/utility/Random.h @ 365

 // $Id: Random.h 365 2005-08-05 08:32:13Z peter $

#ifndef _theplu_random_
#define _theplu_random_


#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

namespace theplu {
namespace random {  

  ///
  /// The RNG class provides a single global random number generator
  /// instance with one point of access to the generator.
  ///
  /// This is probably not thread safe.
  ///
  /// @see Design Patterns (the singleton and adapter pattern)
  ///
  class RNG
  {
  public:

    virtual ~RNG(void);

    static RNG* instance(int seed);

    inline const gsl_rng* rng(void) const { return rng_; }

  private:

    RNG(void);

    static RNG* instance_;
    gsl_rng* rng_;
  };



  class RandomContinuous
  {
  public:

    inline RandomContinuous(void) { rng_=RNG::instance(89); }

    virtual double operator()(void) const = 0;

  protected:
    RNG* rng_;
  };



  class RandomContinuousUniform : public RandomContinuous
  {
  public:

    inline double operator()(void) const { return gsl_rng_uniform(rng_->rng());}

  };

  ///
  /// Class to generate Gaussian random numbers. 
  ///
  class RandomGaussian : public RandomContinuous
  {
  public:
    inline RandomGaussian(const double s=1, const double m=0) 
      : m_(m), s_(s) {}

    inline double operator()(void) const 
    { return gsl_ran_gaussian(rng_->rng(), s_)+m_; }

    inline double operator()(const double s) const 
    { return gsl_ran_gaussian(rng_->rng(), s); }

    inline double operator()(const double s, const double m) const 
    { return gsl_ran_gaussian(rng_->rng(), s)+m; }

  private:
    double m_;
    double s_;
  };

  ///
  /// @brief exponential
  ///
  class RandomExponential : public RandomContinuous
  {
  public:
    ///
    /// Constructor
    ///
    inline RandomExponential(const double m=1) 
      : m_(m){}

    ///
    ///
    ///
    inline double operator()(void) const 
    { return gsl_ran_exponential(rng_->rng(), m_); }

    ///
    ///
    ///
    inline double operator()(const double m) const 
    { return gsl_ran_exponential(rng_->rng(), m); }

  private:
    double m_;
  };


  ///
  /// Base class for discrete random distributions.
  ///
  class RandomDiscrete 
  {
  public:
    ///
    /// Constructor
    ///
    inline RandomDiscrete(void) { rng_=RNG::instance(89); }

    ///
    /// Common interface for all Discrete Random Classes
    ///
    virtual long operator()(void) const = 0;
    
  protected:
    RNG* rng_;
  };


  class RandomDiscreteUniform : public RandomDiscrete 
  {
  public:
    ///
    /// Constructor.
    ///
    inline RandomDiscreteUniform(const u_long range, const long min=0) 
      : min_(min), range_(range) {}

    ///
    /// @return A random number between \a min_ and \a range_ - \a min_
    ///
    inline long operator()(void) const 
    { return gsl_rng_uniform_int(rng_->rng(), range_)+min_; }

    ///
    /// @return A random number between 0 and \a range.
    ///
    inline long operator()(const u_long range) const 
    { return gsl_rng_uniform_int(rng_->rng(), range); }

    ///
    /// @return A random number between \a min and \a range - \a min
    ///
    inline long operator()(const u_long range, const long min) const 
    { return gsl_rng_uniform_int(rng_->rng(), range)+min; }


  private:
    long min_;
    u_long range_;
  };


  ///
  /// @brief Poisson
  ///
  class RandomPoisson : public RandomContinuous
  {
  public:
    ///
    /// Constructor
    ///
    inline RandomPoisson(const double m=1) 
      : m_(m){}

    ///
    ///
    ///
    inline double operator()(void) const 
    { return gsl_ran_poisson(rng_->rng(), m_); }

    ///
    ///
    ///
    inline double operator()(const double m) const 
    { return gsl_ran_poisson(rng_->rng(), m); }

  private:
    double m_;
  };

  ///
  /// @brief General
  ///
  class RandomDiscreteGeneral : public RandomDiscrete {
  };


}} // of namespace random and namespace theplu

#endif