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>


namespace itpp
{


typedef DSFMT_19937_RNG Random_Generator;


void RNG_reset(unsigned int seed);
void RNG_reset();
void RNG_randomize();
void RNG_get_state(ivec &state);
void RNG_set_state(const ivec &state);

class 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 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 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 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 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 Gamma_RNG
{
public:
  Gamma_RNG(double a = 1.0, double b = 1.0): alpha(a), beta(b) {}
  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:
  double alpha;
  double beta;
  Random_Generator RNG;
  Normal_RNG NRNG;
};

class 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 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(0.0), sigma(1.0), norm_factor(1.0 / std::sqrt(2.0)) {}
  void setup(std::complex<double> mean, double variance) {
    m = mean;
    sigma = std::sqrt(variance);
  }
  void get_setup(std::complex<double> &mean, double &variance) {
    mean = m;
    variance = sigma * sigma;
  }
  std::complex<double> operator()() { return sigma*sample() + m; }
  cvec operator()(int n) {
    cvec temp(n);
    sample_vector(n, temp);
    temp *= sigma;
    temp += m;
    return temp;
  }
  cmat operator()(int h, int w) {
    cmat temp(h, w);
    sample_matrix(h, w, temp);
    temp *= sigma;
    temp += m;
    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:
  std::complex<double> m;
  double sigma;
  const double norm_factor;
  Normal_RNG nRNG;
};

class 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 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 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 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

---