/*---------------------------------------------------------------------------*
* IT++ *
*---------------------------------------------------------------------------*
* Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger, *
* Tobias Ringström, and Jonas Samuelsson. *
* *
* Permission to use, copy, modify, and distribute this software and its *
* documentation under the terms of the GNU General Public License is hereby *
* granted. No representations are made about the suitability of this *
* software for any purpose. It is provided "as is" without expressed or *
* implied warranty. See the GNU General Public License for more details. *
*---------------------------------------------------------------------------*/
/*!
\file
\brief Implementation of classes for random number generators.
\author Tony Ottosson
1.14
2004/09/23 14:26:13
*/
#include <ctime>
#ifdef _MSC_VER
# include <process.h>
# define getpid _getpid
#else
# include <unistd.h>
#endif
#include "base/binary.h"
#include "base/matfunc.h"
#include "base/scalfunc.h"
#include "base/random.h"
namespace itpp {
int Random_Generator::left = 0;
const unsigned int Random_Generator::MAXINT = 0xffffffff; // largest value from randInt()
const unsigned int Random_Generator::MAGIC = 0x9908b0dfU; // magic constant
unsigned int Random_Generator::state[624];
unsigned int *Random_Generator::pNext;
//! Variable used to ensure proper seed initialization
static bool __Random_Generator_seed_is_initialized = false;
///////////////////////////////////////////////
// Random_Generator
///////////////////////////////////////////////
Random_Generator::Random_Generator()
{
if (!__Random_Generator_seed_is_initialized){
lastSeed = 4357U;
reset();
__Random_Generator_seed_is_initialized = true;
}
}
void Random_Generator::randomize()
{
lastSeed = time(0) * getpid(); // not good enough if randomize is used within a short time-interval
reset();
}
void Random_Generator::get_state(ivec &out_state)
{
out_state.set_size(625, false);
for (int i=0; i<624; i++)
out_state(i) = state[i];
out_state(624) = left; // the number of elements left in state before reload
//saved_pNext = pNext;
}
void Random_Generator::set_state(ivec &new_state)
{
it_assert(new_state.size()==625, "Random_Generator::set_state(): Not a valid state vector");
for (int i=0; i<624; i++)
state[i] = new_state(i);
left = new_state(624);
pNext = &state[624-left];
}
// Set the seed of the Global Random Number Generator
void RNG_reset(unsigned long seed)
{
Random_Generator RNG;
RNG.reset(seed);
}
// Set the seed of the Global Random Number Generator to the same as last time
void RNG_reset()
{
Random_Generator RNG;
RNG.reset();
}
// Set a random seed for the Global Random Number Generator
void RNG_randomize()
{
Random_Generator RNG;
RNG.randomize();
}
// Save current full state of generator in memory
void RNG_get_state(ivec &state)
{
Random_Generator RNG;
RNG.get_state(state);
}
// Resume the state saved in memory
void RNG_set_state(ivec &state)
{
Random_Generator RNG;
RNG.set_state(state);
}
///////////////////////////////////////////////
// I_Uniform_RNG
///////////////////////////////////////////////
I_Uniform_RNG::I_Uniform_RNG(int min, int max)
{
setup(min, max);
}
void I_Uniform_RNG::setup(int min, int max)
{
if (min <= max) {
lo = min;
hi = max;
}
else {
lo = max;
hi = min;
}
}
void I_Uniform_RNG::get_setup(int &min, int &max) const
{
min = lo;
max = hi;
}
ivec I_Uniform_RNG::operator()(int n)
{
ivec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
imat I_Uniform_RNG::operator()(int h, int w)
{
imat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
///////////////////////////////////////////////
// Uniform_RNG
///////////////////////////////////////////////
Uniform_RNG::Uniform_RNG(double min, double max)
{
setup(min, max);
}
void Uniform_RNG::setup(double min, double max)
{
if (min <= max) {
lo_bound = min;
hi_bound = max;
}
else {
lo_bound = max;
hi_bound = min;
}
}
void Uniform_RNG::get_setup(double &min, double &max) const
{
min = lo_bound;
max = hi_bound;
}
///////////////////////////////////////////////
// Exp_RNG
///////////////////////////////////////////////
Exponential_RNG::Exponential_RNG(double lambda)
{
setup(lambda);
}
vec Exponential_RNG::operator()(int n)
{
vec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
mat Exponential_RNG::operator()(int h, int w)
{
mat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
///////////////////////////////////////////////
// Normal_RNG
///////////////////////////////////////////////
void Normal_RNG::get_setup(double &meanval, double &variance) const
{
meanval = mean;
variance = sigma*sigma;
}
///////////////////////////////////////////////
// Laplace_RNG
///////////////////////////////////////////////
Laplace_RNG::Laplace_RNG(double meanval, double variance)
{
setup(meanval, variance);
}
void Laplace_RNG::setup(double meanval, double variance)
{
mean = meanval;
var = variance;
}
void Laplace_RNG::get_setup(double &meanval, double &variance) const
{
meanval = mean;
variance = var;
}
vec Laplace_RNG::operator()(int n)
{
vec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
mat Laplace_RNG::operator()(int h, int w)
{
mat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
///////////////////////////////////////////////
// AR1_Normal_RNG
///////////////////////////////////////////////
AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho)
{
mean = meanval;
var = variance;
r = rho;
mem = 0.0;
factr = -2.0 * var * (1.0 - rho*rho);
odd = true;
}
void AR1_Normal_RNG::setup(double meanval, double variance, double rho)
{
mean = meanval;
var = variance;
r = rho;
factr = -2.0 * var * (1.0 - rho*rho);
mem = 0.0;
odd = true;
}
void AR1_Normal_RNG::get_setup(double &meanval, double &variance, double &rho) const
{
meanval = mean;
variance = var;
rho = r;
}
vec AR1_Normal_RNG::operator()(int n)
{
vec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
mat AR1_Normal_RNG::operator()(int h, int w)
{
mat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
void AR1_Normal_RNG::reset()
{
mem = 0.0;
}
///////////////////////////////////////////////
// Weibull_RNG
///////////////////////////////////////////////
Weibull_RNG::Weibull_RNG(double lambda, double beta)
{
setup(lambda, beta);
}
void Weibull_RNG::setup(double lambda, double beta)
{
l=lambda;
b=beta;
mean = gamma(1.0 + 1.0 / b) / l;
var = gamma(1.0+2.0/b)/(l*l) - mean;
}
vec Weibull_RNG::operator()(int n)
{
vec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
mat Weibull_RNG::operator()(int h, int w)
{
mat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
///////////////////////////////////////////////
// Rayleigh_RNG
///////////////////////////////////////////////
Rayleigh_RNG::Rayleigh_RNG(double sigma)
{
setup(sigma);
}
vec Rayleigh_RNG::operator()(int n)
{
vec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
mat Rayleigh_RNG::operator()(int h, int w)
{
mat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
///////////////////////////////////////////////
// Rice_RNG
///////////////////////////////////////////////
Rice_RNG::Rice_RNG(double lambda, double beta)
{
setup(lambda, beta);
}
vec Rice_RNG::operator()(int n)
{
vec vv(n);
for (int i=0; i<n; i++)
vv(i) = sample();
return vv;
}
mat Rice_RNG::operator()(int h, int w)
{
mat mm(h,w);
int i,j;
for (i=0; i<h; i++)
for (j=0; j<w; j++)
mm(i,j) = sample();
return mm;
}
// -------------------------------------------------------------------------------
} //namespace itpp
syntax highlighted by Code2HTML, v. 0.9.1