/*!
 * \file
 * \brief Definition of classes for random number generators
 * \author Tony Ottosson and Adam Piatyszek
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2010  (see AUTHORS file for a list of contributors)
 *
 * This program 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 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 *
 * -------------------------------------------------------------------------
 */

#ifndef RANDOM_H
#define RANDOM_H

#include <itpp/base/operators.h>
#include <ctime>


namespace itpp
{

//! \addtogroup randgen

/*!
 * \brief Base class for random (stochastic) sources.
 * \ingroup randgen
 *
 * The Random_Generator class is based on the MersenneTwister MTRand
 * class code in version 1.0 (15 May 2003) by Richard J. Wagner. See
 * http://www-personal.engin.umich.edu/~wagnerr/MersenneTwister.html
 * for details.
 *
 * Here are the original release notes copied from the
 * \c MersenneTwister.h file:
 *
 * \verbatim
 * Mersenne Twister random number generator -- a C++ class MTRand Based on
 * code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus Richard J.
 * Wagner v1.0 15 May 2003 rjwagner@writeme.com
 *
 * The Mersenne Twister is an algorithm for generating random numbers. It
 * was designed with consideration of the flaws in various other generators.
 * The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
 * are far greater. The generator is also fast; it avoids multiplication and
 * division, and it benefits from caches and pipelines. For more information
 * see the inventors' web page at
 * http://www.math.keio.ac.jp/~matumoto/emt.html

 * Reference:
 * M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
 * Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions
 * on Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp.
 * 3-30.
 *
 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
 * Copyright (C) 2000 - 2003, Richard J. Wagner
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 *   1. Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *
 *   2. Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *
 *   3. The names of its contributors may not be used to endorse or promote
 *      products derived from this software without specific prior written
 *      permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * The original code included the following notice:
 *
 *     When you use this, send an email to: matumoto@math.keio.ac.jp with an
 *     appropriate reference to your work.
 *
 * It would be nice to CC: rjwagner@writeme.com and
 * Cokus@math.washington.edu when you write.
 * \endverbatim
 */
class Random_Generator
{
public:
  //! Construct a new Random_Generator.
  Random_Generator() { if (!initialized) reset(4357U); }
  //! Construct Random_Generator object using \c seed
  Random_Generator(unsigned int seed) { reset(seed); }
  //! Set the seed to a semi-random value (based on hashed time and clock).
  void randomize() { reset(hash(time(0), clock())); }
  //! Reset the source. The same sequance will be generated as the last time.
  void reset() { initialize(last_seed); reload(); initialized = true; }
  //! Reset the source after setting the seed to seed.
  void reset(unsigned int seed) { last_seed = seed; reset(); }

  //! Return a uniformly distributed [0,2^32-1] integer.
  unsigned int random_int() {
    if (left == 0) reload();
    --left;

    register unsigned int s1;
    s1 = *pNext++;
    s1 ^= (s1 >> 11);
    s1 ^= (s1 <<  7) & 0x9d2c5680U;
    s1 ^= (s1 << 15) & 0xefc60000U;
    return (s1 ^(s1 >> 18));
  }

  //! Return a uniformly distributed (0,1) value.
  double random_01() { return (random_int() + 0.5) * (1.0 / 4294967296.0); }
  //! Return a uniformly distributed [0,1) value.
  double random_01_lclosed() { return random_int() * (1.0 / 4294967296.0); }
  //! Return a uniformly distributed [0,1] value.
  double random_01_closed() { return random_int() * (1.0 / 4294967295.0); }
  //! Return a uniformly distributed [0,1) value in 53-bit resolution.
  double random53_01_lclosed() {
    return ((random_int() >> 5) * 67108864.0 + (random_int() >> 6))
           * (1.0 / 9007199254740992.0); // by Isaku Wada
  }

  //! Save current full state of generator in memory
  void get_state(ivec &out_state);
  //! Resume the state saved in memory. Clears memory.
  void set_state(ivec &new_state);

private:
  //! initialization flag
  static bool initialized;
  //! seed used for initialisation
  unsigned int last_seed;
  //! internal state
  static unsigned int state[624];
  //! next value to get from state
  static unsigned int *pNext;
  //! number of values left before reload needed
  static int left;

  /*!
   * \brief Initialize generator state with seed.
   * See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
   * \note In previous versions, most significant bits (MSBs) of the seed
   * affect only MSBs of the state array. Modified 9 Jan 2002 by Makoto
   * Matsumoto.
   */
  void initialize(unsigned int seed) {
    register unsigned int *s = state;
    register unsigned int *r = state;
    register int i = 1;
    *s++ = seed & 0xffffffffU;
    for (; i < 624; ++i) {
      *s++ = (1812433253U * (*r ^(*r >> 30)) + i) & 0xffffffffU;
      r++;
    }
  }

  /*!
   * \brief Generate N new values in state.
   * Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
   */
  void reload() {
    register unsigned int *p = state;
    register int i;
    for (i = 624 - 397; i--; ++p)
      *p = twist(p[397], p[0], p[1]);
    for (i = 397; --i; ++p)
      *p = twist(p[397-624], p[0], p[1]);
    *p = twist(p[397-624], p[0], state[0]);

    left = 624, pNext = state;
  }
  //!
  unsigned int hiBit(const unsigned int& u) const { return u & 0x80000000U; }
  //!
  unsigned int loBit(const unsigned int& u) const { return u & 0x00000001U; }
  //!
  unsigned int loBits(const unsigned int& u) const { return u & 0x7fffffffU; }
  //!
  unsigned int mixBits(const unsigned int& u, const unsigned int& v) const
  { return hiBit(u) | loBits(v); }

  /*
   * ----------------------------------------------------------------------
   * --- ediap - 2007/01/17 ---
   * ----------------------------------------------------------------------
   * Wagners's implementation of the twist() function was as follows:
   *  { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfU); }
   * However, this code caused a warning/error under MSVC++, because
   * unsigned value loBit(s1) is being negated with `-' (c.f. bug report
   * [1635988]). I changed this to the same implementation as is used in
   * original C sources of Mersenne Twister RNG:
   *  #define MATRIX_A 0x9908b0dfUL
   *  #define UMASK 0x80000000UL
   *  #define LMASK 0x7fffffffUL
   *  #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
   *  #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
   * ----------------------------------------------------------------------
   */
  //!
  unsigned int twist(const unsigned int& m, const unsigned int& s0,
                     const unsigned int& s1) const
  { return m ^(mixBits(s0, s1) >> 1) ^(loBit(s1) ? 0x9908b0dfU : 0U); }
  //!
  unsigned int hash(time_t t, clock_t c);
};


//! \addtogroup randgen
//!@{

//! Set the seed of the Global Random Number Generator
void RNG_reset(unsigned int seed);
//! Set the seed of the Global Random Number Generator to the same as last reset/init
void RNG_reset();
//! Set a random seed for the Global Random Number Generator.
void RNG_randomize();
//! Save current full state of generator in memory
void RNG_get_state(ivec &state);
//! Resume the state saved in memory
void RNG_set_state(ivec &state);
//!@}

/*!
  \brief Bernoulli distribution
  \ingroup randgen
*/
class Bernoulli_RNG
{
public:
  //! Binary source with probability prob for a 1
  Bernoulli_RNG(double prob) { setup(prob); }
  //! Binary source with probability prob for a 1
  Bernoulli_RNG() { p = 0.5; }
  //! set the probability
  void setup(double prob) {
    it_assert(prob >= 0.0 && prob <= 1.0, "The bernoulli source probability must be between 0 and 1");
    p = prob;
  }
  //! return the probability
  double get_setup() const { return p; }
  //! Get one sample.
  bin operator()() { return sample(); }
  //! Get a sample vector.
  bvec operator()(int n) { bvec temp(n); sample_vector(n, temp); return temp; }
  //! Get a sample matrix.
  bmat operator()(int h, int w) { bmat temp(h, w); sample_matrix(h, w, temp); return temp; }
  //! Get a sample
  bin sample() { return bin(RNG.random_01() < p ? 1 : 0); }
  //! Get a sample vector.
  void sample_vector(int size, bvec &out) {
    out.set_size(size, false);
    for (int i = 0; i < size; i++) out(i) = sample();
  }
  //! Get a sample matrix.
  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;
};

/*!
  \brief Integer uniform distribution
  \ingroup randgen

  Example: Generation of random uniformly distributed integers in the interval [0,10].
  \code
  #include "itpp/sigproc.h"

  int main() {

  I_Uniform_RNG gen(0, 10);

  cout << gen() << endl; // prints a random integer
  cout << gen(10) << endl; // prints 10 random integers
  }
  \endcode
*/
class I_Uniform_RNG
{
public:
  //! constructor. Sets min and max values.
  I_Uniform_RNG(int min = 0, int max = 1);
  //! set min and max values
  void setup(int min, int max);
  //! get the parameters
  void get_setup(int &min, int &max) const;
  //! Get one sample.
  int operator()() { return sample(); }
  //! Get a sample vector.
  ivec operator()(int n);
  //! Get a sample matrix.
  imat operator()(int h, int w);
  //! Return a single value from this random generator
  int sample() { return (floor_i(RNG.random_01() * (hi - lo + 1)) + lo); }
protected:
private:
  //!
  int lo;
  //!
  int hi;
  //!
  Random_Generator RNG;
};

/*!
  \brief Uniform distribution
  \ingroup randgen
*/
class Uniform_RNG
{
public:
  //! Constructor. Set min, max and seed.
  Uniform_RNG(double min = 0, double max = 1.0);
  //! set min and max
  void setup(double min, double max);
  //! get parameters
  void get_setup(double &min, double &max) const;
  //! Get one sample.
  double operator()() { return (sample() * (hi_bound - lo_bound) + lo_bound); }
  //! Get a sample vector.
  vec operator()(int n) {
    vec temp(n);
    sample_vector(n, temp);
    temp *= hi_bound - lo_bound;
    temp += lo_bound;
    return temp;
  }
  //! Get a sample matrix.
  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;
  }
  //! Get a Uniformly distributed (0,1) sample
  double sample() {  return RNG.random_01(); }
  //! Get a Uniformly distributed (0,1) vector
  void sample_vector(int size, vec &out) {
    out.set_size(size, false);
    for (int i = 0; i < size; i++) out(i) = sample();
  }
  //! Get a Uniformly distributed (0,1) matrix
  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;
};

/*!
  \brief Exponential distribution
  \ingroup randgen
*/
class Exponential_RNG
{
public:
  //! constructor. Set lambda.
  Exponential_RNG(double lambda = 1.0);
  //! Set lambda
  void setup(double lambda) { l = lambda; }
  //! get lambda
  double get_setup() const;
  //! Get one sample.
  double operator()() { return sample(); }
  //! Get a sample vector.
  vec operator()(int n);
  //! Get a sample matrix.
  mat operator()(int h, int w);
protected:
private:
  //!
  double sample() {  return (-std::log(RNG.random_01()) / l); }
  //!
  double l;
  //!
  Random_Generator RNG;
};

/*!
 * \brief Normal distribution
 * \ingroup randgen
 *
 * Normal (Gaussian) random variables, using a simplified Ziggurat method.
 *
 * For details see the following arcticle: George Marsaglia, Wai Wan
 * Tsang, "The Ziggurat Method for Generating Random Variables", Journal
 * of Statistical Software, vol. 5 (2000), no. 8
 *
 * This implementation is based on the generator written by Jochen Voss
 * found at http://seehuhn.de/comp/ziggurat/, which is also included in
 * the GSL library (randlist/gauss.c).
 */
class Normal_RNG
{
public:
  //! Constructor. Set mean and variance.
  Normal_RNG(double meanval, double variance):
      mean(meanval), sigma(std::sqrt(variance)) {}
  //! Constructor. Set mean and variance.
  Normal_RNG(): mean(0.0), sigma(1.0) {}
  //! Set mean, and variance
  void setup(double meanval, double variance)
  { mean = meanval; sigma = std::sqrt(variance); }
  //! Get mean and variance
  void get_setup(double &meanval, double &variance) const;
  //! Get one sample.
  double operator()() { return (sigma*sample() + mean); }
  //! Get a sample vector.
  vec operator()(int n) {
    vec temp(n);
    sample_vector(n, temp);
    temp *= sigma;
    temp += mean;
    return temp;
  }
  //! Get a sample matrix.
  mat operator()(int h, int w) {
    mat temp(h, w);
    sample_matrix(h, w, temp);
    temp *= sigma;
    temp += mean;
    return temp;
  }
  //! Get a Normal distributed (0,1) sample
  double sample();

  //! Get a Normal distributed (0,1) vector
  void sample_vector(int size, vec &out) {
    out.set_size(size, false);
    for (int i = 0; i < size; i++) out(i) = sample();
  }

  //! Get a Normal distributed (0,1) matrix
  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;
};

/*!
  \brief Laplacian distribution
  \ingroup randgen
*/
class Laplace_RNG
{
public:
  //! Constructor. Set mean and variance.
  Laplace_RNG(double meanval = 0.0, double variance = 1.0);
  //! Set mean and variance
  void setup(double meanval, double variance);
  //! Get mean and variance
  void get_setup(double &meanval, double &variance) const;
  //! Get one sample.
  double operator()() { return sample(); }
  //! Get a sample vector.
  vec operator()(int n);
  //! Get a sample matrix.
  mat operator()(int h, int w);
  //! Returns a single sample
  double sample() {
    double u = RNG.random_01();
    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;
};

/*!
  \brief A Complex Normal Source
  \ingroup randgen
*/
class Complex_Normal_RNG
{
public:
  //! Constructor. Set mean and variance.
  Complex_Normal_RNG(std::complex<double> mean, double variance):
      norm_factor(1.0 / std::sqrt(2.0)) {
    setup(mean, variance);
  }
  //! Default constructor
  Complex_Normal_RNG(): m(0.0), sigma(1.0), norm_factor(1.0 / std::sqrt(2.0)) {}
  //! Set mean and variance
  void setup(std::complex<double> mean, double variance) {
    m = mean;
    sigma = std::sqrt(variance);
  }
  //! Get mean and variance
  void get_setup(std::complex<double> &mean, double &variance) {
    mean = m;
    variance = sigma * sigma;
  }
  //! Get one sample.
  std::complex<double> operator()() { return sigma*sample() + m; }
  //! Get a sample vector.
  cvec operator()(int n) {
    cvec temp(n);
    sample_vector(n, temp);
    temp *= sigma;
    temp += m;
    return temp;
  }
  //! Get a sample matrix.
  cmat operator()(int h, int w) {
    cmat temp(h, w);
    sample_matrix(h, w, temp);
    temp *= sigma;
    temp += m;
    return temp;
  }
  //! Get a Complex Normal (0,1) distributed sample
  std::complex<double> sample() {
    double a = nRNG.sample() * norm_factor;
    double b = nRNG.sample() * norm_factor;
    return std::complex<double>(a, b);
  }

  //! Get a Complex Normal (0,1) distributed vector
  void sample_vector(int size, cvec &out) {
    out.set_size(size, false);
    for (int i = 0; i < size; i++) out(i) = sample();
  }

  //! Get a Complex Normal (0,1) distributed matrix
  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();
  }

  //! Dummy assignment operator - MSVC++ warning C4512
  Complex_Normal_RNG & operator=(const Complex_Normal_RNG&) { return *this; }

private:
  std::complex<double> m;
  double sigma;
  const double norm_factor;
  Normal_RNG nRNG;
};

/*!
  \brief Filtered normal distribution
  \ingroup randgen
*/
class AR1_Normal_RNG
{
public:
  //! Constructor. Set mean, variance, and correlation.
  AR1_Normal_RNG(double meanval = 0.0, double variance = 1.0,
                 double rho = 0.0);
  //! Set mean, variance, and correlation
  void setup(double meanval, double variance, double rho);
  //! Get mean, variance and correlation
  void get_setup(double &meanval, double &variance, double &rho) const;
  //! Set memory contents to zero
  void reset();
  //! Get a single random sample
  double operator()() { return sample(); }
  //! Get a sample vector.
  vec operator()(int n);
  //! Get a sample matrix.
  mat operator()(int h, int w);
private:
  double sample() {
    mem *= r;
    if (odd) {
      r1 = m_2pi * RNG.random_01();
      r2 = std::sqrt(factr * std::log(RNG.random_01()));
      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;
};

/*!
  \brief Gauss_RNG is the same as Normal Source
  \ingroup randgen
*/
typedef Normal_RNG Gauss_RNG;

/*!
  \brief AR1_Gauss_RNG is the same as AR1_Normal_RNG
  \ingroup randgen
*/
typedef AR1_Normal_RNG AR1_Gauss_RNG;

/*!
  \brief Weibull distribution
  \ingroup randgen
*/
class Weibull_RNG
{
public:
  //! Constructor. Set lambda and beta.
  Weibull_RNG(double lambda = 1.0, double beta = 1.0);
  //! Set lambda, and beta
  void setup(double lambda, double beta);
  //! Get lambda and beta
  void get_setup(double &lambda, double &beta) { lambda = l; beta = b; }
  //! Get one sample.
  double operator()() { return sample(); }
  //! Get a sample vector.
  vec operator()(int n);
  //! Get a sample matrix.
  mat operator()(int h, int w);
private:
  double sample() {
    return (std::pow(-std::log(RNG.random_01()), 1.0 / b) / l);
  }
  double l, b;
  double mean, var;
  Random_Generator RNG;
};

/*!
  \brief Rayleigh distribution
  \ingroup randgen
*/
class Rayleigh_RNG
{
public:
  //! Constructor. Set sigma.
  Rayleigh_RNG(double sigma = 1.0);
  //! Set sigma
  void setup(double sigma) { sig = sigma; }
  //! Get sigma
  double get_setup() { return sig; }
  //! Get one sample.
  double operator()() { return sample(); }
  //! Get a sample vector.
  vec operator()(int n);
  //! Get a sample matrix.
  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;
};

/*!
  \brief Rice distribution
  \ingroup randgen
*/
class Rice_RNG
{
public:
  //! Constructor. Set sigma, and v (if v = 0, Rice -> Rayleigh).
  Rice_RNG(double sigma = 1.0, double v = 1.0);
  //! Set sigma, and v (if v = 0, Rice -> Rayleigh).
  void setup(double sigma, double v) { sig = sigma; s = v; }
  //! Get parameters
  void get_setup(double &sigma, double &v) { sigma = sig; v = s; }
  //! Get one sample
  double operator()() { return sample(); }
  //! Get a sample vector
  vec operator()(int n);
  //! Get a sample matrix
  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;
};

//! \addtogroup randgen
//!@{

//! Generates a random bit (equally likely 0s and 1s)
inline bin randb(void) { Bernoulli_RNG src; return src.sample(); }
//! Generates a random bit vector (equally likely 0s and 1s)
inline void randb(int size, bvec &out) { Bernoulli_RNG src; src.sample_vector(size, out); }
//! Generates a random bit vector (equally likely 0s and 1s)
inline bvec randb(int size) { bvec temp; randb(size, temp); return temp; }
//! Generates a random bit matrix (equally likely 0s and 1s)
inline void randb(int rows, int cols, bmat &out) { Bernoulli_RNG src; src.sample_matrix(rows, cols, out); }
//! Generates a random bit matrix (equally likely 0s and 1s)
inline bmat randb(int rows, int cols) { bmat temp; randb(rows, cols, temp); return temp; }

//! Generates a random uniform (0,1) number
inline double randu(void) { Uniform_RNG src; return src.sample(); }
//! Generates a random uniform (0,1) vector
inline void randu(int size, vec &out) { Uniform_RNG src; src.sample_vector(size, out); }
//! Generates a random uniform (0,1) vector
inline vec randu(int size) { vec temp; randu(size, temp); return temp; }
//! Generates a random uniform (0,1) matrix
inline void randu(int rows, int cols, mat &out) { Uniform_RNG src; src.sample_matrix(rows, cols, out); }
//! Generates a random uniform (0,1) matrix
inline mat randu(int rows, int cols) { mat temp; randu(rows, cols, temp); return temp; }

//! Generates a random integer in the interval [low,high]
inline int randi(int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(); }
//! Generates a random ivec with elements in the interval [low,high]
inline ivec randi(int size, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(size); }
//! Generates a random imat with elements in the interval [low,high]
inline imat randi(int rows, int cols, int low, int high) { I_Uniform_RNG src; src.setup(low, high); return src(rows, cols); }

//! Generates a random Rayleigh vector
inline vec randray(int size, double sigma = 1.0) { Rayleigh_RNG src; src.setup(sigma); return src(size); }

//! Generates a random Rice vector (See J.G. Poakis, "Digital Communications, 3rd ed." p.47)
inline vec randrice(int size, double sigma = 1.0, double s = 1.0) { Rice_RNG src; src.setup(sigma, s); return src(size); }

//! Generates a random complex Gaussian vector
inline vec randexp(int size, double lambda = 1.0) { Exponential_RNG src; src.setup(lambda); return src(size); }

//! Generates a random Gaussian (0,1) variable
inline double randn(void) { Normal_RNG src; return src.sample(); }
//! Generates a random Gaussian (0,1) vector
inline void randn(int size, vec &out) { Normal_RNG src; src.sample_vector(size, out); }
//! Generates a random Gaussian (0,1) vector
inline vec randn(int size) { vec temp; randn(size, temp); return temp; }
//! Generates a random Gaussian (0,1) matrix
inline void randn(int rows, int cols, mat &out)  { Normal_RNG src; src.sample_matrix(rows, cols, out); }
//! Generates a random Gaussian (0,1) matrix
inline mat randn(int rows, int cols) { mat temp; randn(rows, cols, temp); return temp; }

/*! \brief Generates a random complex Gaussian (0,1) variable

The real and imaginary parts are independent and have variances equal to 0.5
*/
inline std::complex<double> randn_c(void) { Complex_Normal_RNG src; return src.sample(); }
//! Generates a random complex Gaussian (0,1) vector
inline void randn_c(int size, cvec &out)  { Complex_Normal_RNG src; src.sample_vector(size, out); }
//! Generates a random complex Gaussian (0,1) vector
inline cvec randn_c(int size) { cvec temp; randn_c(size, temp); return temp; }
//! Generates a random complex Gaussian (0,1) matrix
inline void randn_c(int rows, int cols, cmat &out) { Complex_Normal_RNG src; src.sample_matrix(rows, cols, out); }
//! Generates a random complex Gaussian (0,1) matrix
inline cmat randn_c(int rows, int cols) { cmat temp; randn_c(rows, cols, temp); return temp; }

//!@}

} // namespace itpp

#endif // #ifndef RANDOM_H
