random.h

Go to the documentation of this file.

#ifndef RANDOM_H
#define RANDOM_H

#include <itpp/base/random_dsfmt.h>
#include <itpp/base/operators.h>
#include <itpp/itexports.h>

namespace itpp
{


ITPP_EXPORT void GlobalRNG_reset(unsigned int seed);

ITPP_EXPORT void GlobalRNG_reset();

ITPP_EXPORT unsigned int GlobalRNG_get_local_seed();

ITPP_EXPORT void GlobalRNG_randomize();

ITPP_EXPORT void GlobalRNG_get_state(ivec &state);

ITPP_EXPORT void GlobalRNG_set_state(const ivec &state);




ITPP_EXPORT void RNG_reset(unsigned int seed);

ITPP_EXPORT void RNG_reset();

ITPP_EXPORT void RNG_randomize();

ITPP_EXPORT void RNG_get_state(ivec &state);

ITPP_EXPORT void RNG_set_state(const ivec &state);



class ITPP_EXPORT Random_Generator
{
  typedef random_details::ActiveDSFMT DSFMT;
public:
  Random_Generator(): _dsfmt(random_details::lc_get()) {
    if(!random_details::lc_is_initialized()) {
      _dsfmt.init_gen_rand(GlobalRNG_get_local_seed());
      random_details::lc_mark_initialized();
    }
  }

  //Constructor using a certain seed - do not want to provide it. Use free-standing functions to change per-thread local seed
  //Random_Generator(unsigned int seed) { init_gen_rand(seed); random_details::lc_mark_initialized();}

  //provide wrappers for DSFMT algorithm functions

  //Set the seed to a semi-random value (based on hashed time and clock)  - do not want to provide it. Use free-standing functions to change per-thread local seed
  //void randomize(){RNG_randomize();}

  //Reset the generator with the same seed as used last time  - do not want to provide it. Use free-standing functions to change per-thread local seed
  //void reset() {RNG_reset();}

  //Initialise the generator with a new seed (\sa init_gen_rand())  - do not want to provide it. Use free-standing functions to change per-thread local seed
  //void reset(unsigned int seed) { RNG_reset(seed); }

  //Resume the state of the generator from a previously saved ivec  - do not want to provide it. Use free-standing functions to change per-thread local seed
  //void set_state(const ivec &state) {RNG_set_state(state);}

  //Return current state of generator in the form of ivec  - do not want to provide it. Use free-standing functions to change per-thread local seed
  //ivec get_state() const {ivec ret; RNG_get_state(ret); return ret; }

  double random_01() { return genrand_open_open(); }
  double random_01_lclosed() { return genrand_close_open(); }
  double random_01_rclosed() { return genrand_open_close(); }
  uint32_t random_int() { return genrand_uint32(); }

  uint32_t genrand_uint32() { return _dsfmt.genrand_uint32(); }

  double genrand_close1_open2() { return _dsfmt.genrand_close1_open2(); }

  double genrand_close_open() { return genrand_close1_open2() - 1.0; }

  double genrand_open_close() { return 2.0 - genrand_close1_open2(); }

  double genrand_open_open() { return _dsfmt.genrand_open_open();}

private:
  DSFMT _dsfmt;
};

class ITPP_EXPORT Bernoulli_RNG
{
public:
  Bernoulli_RNG(double prob) { setup(prob); }
  Bernoulli_RNG() { p = 0.5; }
  void setup(double prob) {
    it_assert(prob >= 0.0 && prob <= 1.0, "The Bernoulli source probability "
              "must be between 0 and 1");
    p = prob;
  }
  double get_setup() const { return p; }
  bin operator()() { return sample(); }
  bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; }
  bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; }
  bin sample() { return RNG.genrand_close_open() < p ? bin(1) : bin(0); }
  void sample_vector(int size, bvec &out) {
    out.set_size(size, false);
    for(int i = 0; i < size; i++) out(i) = sample();
  }
  void sample_matrix(int rows, int cols, bmat &out) {
    out.set_size(rows, cols, false);
    for(int i = 0; i < rows * cols; i++) out(i) = sample();
  }
protected:
private:
  double p;
  Random_Generator RNG;
};

class ITPP_EXPORT I_Uniform_RNG
{
public:
  I_Uniform_RNG(int min = 0, int max = 1);
  void setup(int min, int max);
  void get_setup(int &min, int &max) const;
  int operator()() { return sample(); }
  ivec operator()(int n);
  imat operator()(int h, int w);
  int sample() {
    return floor_i(RNG.genrand_close_open() * (hi - lo + 1)) + lo;
  }
private:
  int lo;
  int hi;
  Random_Generator RNG;
};

class ITPP_EXPORT Uniform_RNG
{
public:
  Uniform_RNG(double min = 0, double max = 1.0);
  void setup(double min, double max);
  void get_setup(double &min, double &max) const;
  double operator()() { return (sample() * (hi_bound - lo_bound) + lo_bound); }
  vec operator()(int n) {
    vec temp(n);
    sample_vector(n, temp);
    temp *= hi_bound - lo_bound;
    temp += lo_bound;
    return temp;
  }
  mat operator()(int h, int w) {
    mat temp(h, w);
    sample_matrix(h, w, temp);
    temp *= hi_bound - lo_bound;
    temp += lo_bound;
    return temp;
  }
  double sample() {  return RNG.genrand_close_open(); }
  void sample_vector(int size, vec &out) {
    out.set_size(size, false);
    for(int i = 0; i < size; i++) out(i) = sample();
  }
  void sample_matrix(int rows, int cols, mat &out) {
    out.set_size(rows, cols, false);
    for(int i = 0; i < rows * cols; i++) out(i) = sample();
  }
protected:
private:
  double lo_bound, hi_bound;
  Random_Generator RNG;
};

class ITPP_EXPORT Exponential_RNG
{
public:
  Exponential_RNG(double lambda = 1.0);
  void setup(double lambda) { l = lambda; }
  double get_setup() const;
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int h, int w);
private:
  double sample() { return (-std::log(RNG.genrand_open_close()) / l); }
  double l;
  Random_Generator RNG;
};

class ITPP_EXPORT Normal_RNG
{
public:
  Normal_RNG(double meanval, double variance):
    mean(meanval), sigma(std::sqrt(variance)) {}
  Normal_RNG(): mean(0.0), sigma(1.0) {}
  void setup(double meanval, double variance)
  { mean = meanval; sigma = std::sqrt(variance); }
  void get_setup(double &meanval, double &variance) const;
  double operator()() { return (sigma * sample() + mean); }
  vec operator()(int n) {
    vec temp(n);
    sample_vector(n, temp);
    temp *= sigma;
    temp += mean;
    return temp;
  }
  mat operator()(int h, int w) {
    mat temp(h, w);
    sample_matrix(h, w, temp);
    temp *= sigma;
    temp += mean;
    return temp;
  }
  double sample();

  void sample_vector(int size, vec &out) {
    out.set_size(size, false);
    for(int i = 0; i < size; i++) out(i) = sample();
  }

  void sample_matrix(int rows, int cols, mat &out) {
    out.set_size(rows, cols, false);
    for(int i = 0; i < rows * cols; i++) out(i) = sample();
  }
private:
  double mean, sigma;
  static const double ytab[128];
  static const unsigned int ktab[128];
  static const double wtab[128];
  static const double PARAM_R;
  Random_Generator RNG;
};

class ITPP_EXPORT Gamma_RNG
{
public:
  Gamma_RNG(double a = 1.0, double b = 1.0): alpha(a), beta(b) {init_state();}
  void setup(double a, double b) { alpha = a; beta = b; }
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int r, int c);
  double sample();
private:
  void init_state();
  double alpha;
  double beta;

  Random_Generator RNG;
  Normal_RNG NRNG;

  /* State variables - used in Gamma_Rng::sample()*/
  double _s, _s2, _d, _scale;
  double _q0, _b, _si, _c;
};

class ITPP_EXPORT Laplace_RNG
{
public:
  Laplace_RNG(double meanval = 0.0, double variance = 1.0);
  void setup(double meanval, double variance);
  void get_setup(double &meanval, double &variance) const;
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int h, int w);
  double sample() {
    double u = RNG.genrand_open_open();
    double l = sqrt_12var;
    if(u < 0.5)
      l *= std::log(2.0 * u);
    else
      l *= -std::log(2.0 * (1 - u));
    return (l + mean);
  }
private:
  double mean, var, sqrt_12var;
  Random_Generator RNG;
};

class ITPP_EXPORT Complex_Normal_RNG
{
public:
  Complex_Normal_RNG(std::complex<double> mean, double variance):
    norm_factor(1.0 / std::sqrt(2.0)) {
    setup(mean, variance);
  }
  Complex_Normal_RNG(): m_re(0.0), m_im(0.0), sigma(1.0), norm_factor(1.0 / std::sqrt(2.0)) {}
  void setup(std::complex<double> mean, double variance) {
    m_re = mean.real();
    m_im = mean.imag();
    sigma = std::sqrt(variance);
  }
  void get_setup(std::complex<double> &mean, double &variance) {
    mean = std::complex<double>(m_re, m_im);
    variance = sigma * sigma;
  }
  std::complex<double> operator()() { return sigma * sample() + std::complex<double>(m_re, m_im); }
  cvec operator()(int n) {
    cvec temp(n);
    sample_vector(n, temp);
    temp *= sigma;
    temp += std::complex<double>(m_re, m_im);
    return temp;
  }
  cmat operator()(int h, int w) {
    cmat temp(h, w);
    sample_matrix(h, w, temp);
    temp *= sigma;
    temp += std::complex<double>(m_re, m_im);
    return temp;
  }
  std::complex<double> sample() {
    double a = nRNG.sample() * norm_factor;
    double b = nRNG.sample() * norm_factor;
    return std::complex<double>(a, b);
  }

  void sample_vector(int size, cvec &out) {
    out.set_size(size, false);
    for(int i = 0; i < size; i++) out(i) = sample();
  }

  void sample_matrix(int rows, int cols, cmat &out) {
    out.set_size(rows, cols, false);
    for(int i = 0; i < rows * cols; i++) out(i) = sample();
  }

  Complex_Normal_RNG & operator=(const Complex_Normal_RNG&) { return *this; }

private:
  double m_re;
  double m_im;
  double sigma;
  const double norm_factor;
  Normal_RNG nRNG;
};

class ITPP_EXPORT AR1_Normal_RNG
{
public:
  AR1_Normal_RNG(double meanval = 0.0, double variance = 1.0,
                 double rho = 0.0);
  void setup(double meanval, double variance, double rho);
  void get_setup(double &meanval, double &variance, double &rho) const;
  void reset();
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int h, int w);
private:
  double sample() {
    mem *= r;
    if(odd) {
      r1 = m_2pi * RNG.genrand_open_close();
      r2 = std::sqrt(factr * std::log(RNG.genrand_open_close()));
      mem += r2 * std::cos(r1);
    }
    else {
      mem += r2 * std::sin(r1);
    }
    odd = !odd;
    return (mem + mean);
  }
  double mem, r, factr, mean, var, r1, r2;
  bool odd;
  Random_Generator RNG;
};

typedef Normal_RNG Gauss_RNG;

typedef AR1_Normal_RNG AR1_Gauss_RNG;

class ITPP_EXPORT Weibull_RNG
{
public:
  Weibull_RNG(double lambda = 1.0, double beta = 1.0);
  void setup(double lambda, double beta);
  void get_setup(double &lambda, double &beta) { lambda = l; beta = b; }
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int h, int w);
private:
  double sample() {
    return (std::pow(-std::log(RNG.genrand_open_close()), 1.0 / b) / l);
  }
  double l, b;
  double mean, var;
  Random_Generator RNG;
};

class ITPP_EXPORT Rayleigh_RNG
{
public:
  Rayleigh_RNG(double sigma = 1.0);
  void setup(double sigma) { sig = sigma; }
  double get_setup() { return sig; }
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int h, int w);
private:
  double sample() {
    double s1 = nRNG.sample();
    double s2 = nRNG.sample();
    // s1 and s2 are N(0,1) and independent
    return (sig * std::sqrt(s1 * s1 + s2 * s2));
  }
  double sig;
  Normal_RNG nRNG;
};

class ITPP_EXPORT Rice_RNG
{
public:
  Rice_RNG(double sigma = 1.0, double v = 1.0);
  void setup(double sigma, double v) { sig = sigma; s = v; }
  void get_setup(double &sigma, double &v) { sigma = sig; v = s; }
  double operator()() { return sample(); }
  vec operator()(int n);
  mat operator()(int h, int w);
private:
  double sample() {
    double s1 = nRNG.sample() + s;
    double s2 = nRNG.sample();
    // s1 and s2 are N(0,1) and independent
    return (sig * std::sqrt(s1 * s1 + s2 * s2));
  }
  double sig, s;
  Normal_RNG nRNG;
};


inline bin randb(void) { Bernoulli_RNG src; return src.sample(); }
inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); }
inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; }
inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); }
inline bmat randb(int rows, int cols) { bmat temp; randb(rows, cols, temp); return temp; }

inline double randu(void) { Uniform_RNG src; return src.sample(); }
inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); }
inline vec randu(int size) { vec temp; randu(size, temp); return temp; }
inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); }
inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; }

inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); }
inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); }
inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); }

inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); }

inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); }

inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); }

inline double randn(void) { Normal_RNG src; return src.sample(); }
inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); }
inline vec randn(int size) { vec temp; randn(size, temp); return temp; }
inline void randn(int rows, int cols, mat &out)  { Normal_RNG src; src.sample_matrix(rows, cols, out); }
inline mat randn(int rows, int cols) { mat temp; randn(rows, cols, temp); return temp; }

inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); }
inline void randn_c(int size, cvec &out)  { Complex_Normal_RNG src; src.sample_vector(size, out); }
inline cvec randn_c(int size) { cvec temp; randn_c(size, temp); return temp; }
inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); }
inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; }


} // namespace itpp

#endif // #ifndef RANDOM_H