random.h

#ifndef _theplu_yat_random_
#define _theplu_yat_random_

// $Id: random.h 3959 2020-07-30 12:02:27Z peter $

/*
  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2020 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/utility/config_public.h"

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

#include <boost/concept_check.hpp>

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

#include <algorithm>
#include <atomic>
#include <memory>
#include <mutex>
#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 std::atomic<RNG*> instance_;
    // holds one gsl_rng per thread. Access through rng(void) so a
    // gsl_rng is allocated if necessary.
    static thread_local std::unique_ptr<gsl_rng, void(*)(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 std::mutex mutex_;
    static std::mutex init_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 : public 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:
    typedef unsigned long argument_type;

    explicit DiscreteUniform(unsigned long n=0);

    unsigned long operator()(void) const;

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

    unsigned long int min(void) const;

    unsigned long int max(void) const;

  private:
    unsigned long range_;
  };


  class Geometric : public Discrete
  {
  public:
    Geometric(double p);

    /*
      \return a number from Geomtric distribution
    */
    unsigned long operator()(void) const;

    unsigned long operator()(double p) const;

  private:
    double p_;
  };


  class HyperGeometric : public Discrete
  {
  public:
    HyperGeometric(void);

    HyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);

    unsigned long int operator()(void) const;

    unsigned long int operator()(unsigned int n1, unsigned int n2,
                                 unsigned int t) const;
  private:
    unsigned int n1_;
    unsigned int n2_;
    unsigned int t_;
  };

  class NegativeHyperGeometric : public Discrete
  {
  public:
    NegativeHyperGeometric(void);

    NegativeHyperGeometric(unsigned int n1, unsigned int n2, unsigned int t);

    unsigned long int operator()(void) const;

    unsigned long int operator()(unsigned int n1, unsigned int n2,
                                 unsigned int t) const;
  private:
    unsigned int n1_;
    unsigned int n2_;
    unsigned int t_;
  };


  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)
  {
    DiscreteUniform rnd;
    std::random_shuffle(first, last, rnd);
  }

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

#endif