random.h

#ifndef _theplu_yat_random_
#define _theplu_yat_random_

// $Id: random.h 3330 2014-10-14 08:03:25Z peter $

/*
  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "yat/statistics/Histogram.h"
#include "yat/utility/deprecate.h"

// always include this before <boost/excpetion_ptr.hpp> indirectly below
#include "yat/utility/boost_exception_ptr.h"

#include <boost/concept_check.hpp>
#include <boost/iterator/iterator_concepts.hpp>
#include <boost/thread.hpp>
#include <boost/thread/tss.hpp>

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

#include <algorithm>
#include <string>
#include <vector>

namespace theplu {
namespace yat {
namespace random {

  //forward declaration
  class RNG_state;

  class RNG
  {
  public:

    static RNG* instance(void);

    unsigned long max(void) const;

    unsigned long min(void) const;

    std::string name(void) const;

    const gsl_rng* rng(void) const;

    void seed(unsigned long s) const;

    unsigned long seed_from_devurandom(void);

    // return int for backward compatibility with yat 0.8
    int set_state(const RNG_state&);

  private:
    RNG(void);

    RNG(const RNG&);

    RNG& operator=(const RNG&);

    virtual ~RNG(void);
    void rng_alloc(void) const;

    static RNG* instance_;
    // holds one gsl_rng per thread. Access through rng(void) so a
    // gsl_rng is allocated if necessary.
    mutable boost::thread_specific_ptr<gsl_rng> rng_;
    mutable unsigned long seed_;
    // guard needs to be mutable because major mission for it is to protect seed_ against multi-access, and seed_ is mutable...
    mutable boost::mutex mutex_;
  };


  class RNG_state
  {
  public:
    explicit RNG_state(const RNG*);

    RNG_state(const RNG_state&);

    ~RNG_state(void);

    const gsl_rng* rng(void) const;

    RNG_state& operator=(const RNG_state&);

  private:
    gsl_rng* rng_;

    void clone(const gsl_rng&);
  };


  // --------------------- Discrete distribtuions ---------------------

  class Discrete
  {
  public:
    typedef unsigned long int result_type;

    Discrete(void);

    virtual ~Discrete(void);

    void seed(unsigned long s) const YAT_DEPRECATE;

    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;

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

  protected:
    RNG* rng_;
  };

  class Binomial : Discrete
  {
  public:
    Binomial(double p, unsigned int n);

    unsigned long operator()(void) const;
  private:
    double p_;
    unsigned int n_;
  };

  class DiscreteGeneral : public Discrete
  {
  public:
    explicit DiscreteGeneral(const statistics::Histogram& hist);

    DiscreteGeneral(const DiscreteGeneral&);

    ~DiscreteGeneral(void);

    unsigned long operator()(void) const;

    DiscreteGeneral& operator=(const DiscreteGeneral&);

  private:
    void free(void);
    void preproc(void);

    gsl_ran_discrete_t* gen_;
    std::vector<double> p_;
  };

  class DiscreteUniform
    : public Discrete,
      public std::unary_function<unsigned long, unsigned long>
  {
  public:
    explicit DiscreteUniform(unsigned long n=0);

    unsigned long operator()(void) const;

    unsigned long operator()(unsigned long n) const;

  private:
    unsigned long range_;
  };

  class Poisson : public Discrete
  {
  public:
    explicit Poisson(const double m=1);

    unsigned long operator()(void) const;

    unsigned long operator()(const double m) const;

  private:
    double m_;
  };

  // --------------------- Continuous distribtuions ---------------------

  class Continuous
  {
  public:
    typedef double result_type;

    Continuous(void);

    virtual ~Continuous(void);

    void seed(unsigned long s) const YAT_DEPRECATE;

    unsigned long seed_from_devurandom(void) YAT_DEPRECATE;

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

  protected:
    RNG* rng_;
  };

  // ContinuousUniform is declared before ContinuousGeneral to avoid
  // forward declaration
  class ContinuousUniform : public Continuous
  {
  public:
    double operator()(void) const;
  };

  class ContinuousGeneral : public Continuous
  {
  public:
    explicit ContinuousGeneral(const statistics::Histogram& hist);

    double operator()(void) const;

  private:
    const DiscreteGeneral discrete_;
    const statistics::Histogram hist_;
    ContinuousUniform u_;
  };

  class Exponential : public Continuous
  {
  public:
    explicit Exponential(const double m=1);

    double operator()(void) const;

    double operator()(const double m) const;

  private:
    double m_;
  };

  class Gaussian : public Continuous
  {
  public:
    explicit Gaussian(const double s=1, const double m=0);

    double operator()(void) const;

    double operator()(const double s) const;

    double operator()(const double s, const double m) const;

  private:
    double m_;
    double s_;
  };

  template<typename RandomAccessIterator>
  void random_shuffle(RandomAccessIterator first, RandomAccessIterator last)
  {
    typedef RandomAccessIterator rai;
    BOOST_CONCEPT_ASSERT((boost_concepts::SwappableIterator<rai>));
    BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<rai>));
    DiscreteUniform rnd;
    std::random_shuffle(first, last, rnd);
  }

}}} // of namespace random, yat, and theplu

#endif