source: trunk/lib/utility/random_singleton.h @ 364

// $Id: random_singleton.h 364 2005-08-05 08:25:21Z peter $

#ifndef _theplu_utility_random_singleton_ 
#define _theplu_utility_random_singleton_ 


// #include <cstdlib>
// #include <cstdio>
#include <string>

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

namespace theplu {
namespace gslapi {
  class vector;
}

namespace utility {  

  ///
  /// Class defining interface for different random stuff.
  ///
  class random_singleton
  {
  public:
    ///
    ///   Method used to instantiate an object from this class. The
    ///   seed is zero by default. If user wish to use time the user
    ///   should specify a negative seed. Note that if an object
    ///   already exists the method will return that object and hence
    ///   no new object is created.
    ///
    static random_singleton* instance( int seed = 0 );


    /// 
    /// A random unsigned long integer number is returned from a
    /// uniform distribution.
    ///
    /// @return If \a n is zero, or left out, a number between [ \a
    /// min, \a max ] is returned, otherwise a number between 0 and
    /// n-1 is returned.
    ///
    /// @see uniform_int_min() and uniform_int_max().
    ///
    inline u_long uniform_int(const u_long n=0) const
    { return ( n ? gsl_rng_uniform_int(r_,n) : gsl_rng_get( r_ )); }


    /// A random number (double precision) is returned from a uniform
    /// distribution between [ 0, 1 [.
    ///
    inline double uniform_double() const
    { return gsl_rng_uniform( r_ ); }

    /// 
    /// A random number (double precision) is returned from a
    /// gaussioan distribution between where \a sigma can be specified
    /// by caller. \a sigma (\f$\sigma\f$) is the standard deviation
    /// in a gaussian distribution:  
    ///
    /// \f$ f(x) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}\f$
    ///
    inline double gaussian( const double& sigma=1 ) const
    { return gsl_ran_gaussian( r_, sigma ); }


    /// 
    /// A random number (double precision) is returned from an
    /// exponential distribution. \a mu (\f$\mu\f$), the mean, can be
    /// specified by caller: \f$f(x) =
    /// \frac{1}{\mu}e^{-\frac{x}{\mu}}\f$.
    ///
    inline double exponential( const double& mu=1 ) const
    { return gsl_ran_exponential( r_, mu ); }


    ///
    /// @brief random number from Poisson distribution 
    ///
    /// Probability that returned number i \f$ k \f$ is \f$ P(k) =
    /// \frac{m^k}{k!}e^{-m} \f$
    ///
    inline u_int poisson(const double m) const
    { return gsl_ran_poisson(r_, m); }
    
    /// 
    /// @See uniform_int()!
    ///
    inline u_long uniform_min() const
    { return gsl_rng_min( r_ ); }

    /// 
    /// @See uniform_int()!
    ///
    inline u_long uniform_max() const
    { return gsl_rng_max( r_ ); }

    ///
    /// Return seed to user for documentation purposes.
    ///
    inline int seed() const
    { return gsl_rng_default_seed; }

    ///
    /// Set a new seed
    ///
    void seed( const int& seed );


    ///
    /// Set probabilities \a p for the general purpose distribution. 
    ///
    /// @see rnd_discrete
    ///
    void general_distribution_prob( const gslapi::vector& p);

    ///
    /// Get a random number (index) from the general purpose
    /// distribution
    ///
    /// @see general_distribution_prob
    ///
    inline size_t rnd_discrete()
    { return gsl_ran_discrete( r_, gen_ );}

    ///
    /// GSL has different types of random genenerators. This function
    /// returns the one used.
    ///
    std::string generator_type(void) const
    { return static_cast<std::string>( gsl_rng_name( r_ ) ); }

    ///
    /// Operator for STL
    ///
    ///inline double operator()(void) { return uniform_double();}


    ///
    ///   Destructor returning memory using GSL free functions
    ///
    ~random_singleton();

  private:
    ///
    /// pointer to an instance of the class.
    ///
    static random_singleton *singleobj_;

    ///
    /// GSL stuff
    ///
    gsl_rng* r_; 
    gsl_ran_discrete_t* gen_;  // general purpose distribution 

    ///
    /// Private constructor. Only way to create object is via
    /// instance().
    ///
    random_singleton();
  }; 
  
  struct my_uniform_rng
  {
  public:
    inline u_long operator()(u_long i) {
      return random_singleton::instance()->uniform_int(i);
    }
  };

}} // of namespace utility and namespace theplu

#endif