Random.cpp

// Random.cpp is based on Random.h
// Mersenne Twister random number generator -- a C++ class Random
// 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.

// Parts of this file are modified beginning in 29.10.09 for adaption in PXL.
// Parts of this file are modified beginning in 10.02.12 for adaption in CRPropa.

#include "crpropa/Random.h"

#include "crpropa/base64.h"

#include <cstdio>

namespace crpropa {

Random::Random(const uint32_t& oneSeed) {
	seed(oneSeed);
}

Random::Random(uint32_t * const bigSeed, const uint32_t seedLength) {
	seed(bigSeed, seedLength);
}

Random::Random() {
	seed();
}

double Random::rand() {
	return double(randInt()) * (1.0 / 4294967295.0);
}

double Random::rand(const double& n) {
	return rand() * n;
}

double Random::randExc() {
	return double(randInt()) * (1.0 / 4294967296.0);
}

double Random::randExc(const double& n) {
	return randExc() * n;
}

double Random::randDblExc() {
	return (double(randInt()) + 0.5) * (1.0 / 4294967296.0);
}

double Random::randDblExc(const double& n) {
	return randDblExc() * n;
}

double Random::rand53() {
	uint32_t a = randInt() >> 5, b = randInt() >> 6;
	return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0); // by Isaku Wada
}

// Return a real number from a normal (Gaussian) distribution with given
// mean and variance by Box-Muller method
double Random::randNorm(const double& mean, const double& variance) {
	double r = sqrt(-2.0 * log(1.0 - randDblExc())) * variance;
	double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
	return mean + r * cos(phi);
}

double Random::randUniform(double min, double max) {
	return min + (max - min) * rand();
}

double Random::randRayleigh(double sigma) {
	return sigma * sqrt(-2.0 * log(1 - rand()));
}

double Random::randFisher(double kappa) {
	return acos(1. + 1. / kappa * log(1 - rand() * (1 - exp(-2 * kappa))));
}

size_t Random::randBin(const std::vector<float> &cdf) {
	std::vector<float>::const_iterator it = std::lower_bound(cdf.begin(),
			cdf.end(), rand() * cdf.back());
	return it - cdf.begin();
}

size_t Random::randBin(const std::vector<double> &cdf) {
	std::vector<double>::const_iterator it = std::lower_bound(cdf.begin(),
			cdf.end(), rand() * cdf.back());
	return it - cdf.begin();
}

Vector3d Random::randVector() {
	double z = randUniform(-1.0, 1.0);
	double t = randUniform(-1.0 * M_PI, M_PI);
	double r = sqrt(1 - z * z);
	return Vector3d(r * cos(t), r * sin(t), z);
}

Vector3d Random::randVectorAroundMean(const Vector3d &meanDirection,
		double angle) {
	Vector3d axis = meanDirection.cross(randVector());
	Vector3d v = meanDirection;
	return v.getRotated(axis, angle);
}

Vector3d Random::randFisherVector(const Vector3d &meanDirection, double kappa) {
	return randVectorAroundMean(meanDirection, randFisher(kappa));
}

Vector3d Random::randConeVector(const Vector3d &meanDirection, double angularRadius) {
        const double theta = acos(randUniform(1, cos(angularRadius)));
        return randVectorAroundMean(meanDirection, theta);
}

Vector3d Random::randVectorLamberts() {
	// random vector following Lamberts cosine law (https://en.wikipedia.org/wiki/Lambert%27s_cosine_law)
	// for a surface element with normal vector pointing in positive z-axis (0, 0, 1)
	double phi = randUniform(-1.0 * M_PI, M_PI);
	double theta = M_PI / 2.0 - acos(sqrt(randUniform(0, 1)));
	return Vector3d(cos(phi) * cos(theta), sin(phi) * cos(theta), sin(theta));
}

Vector3d Random::randVectorLamberts(const Vector3d &normalVector) {
	// random vector following Lamberts cosine law for a surface element described by normalVector
	Vector3d vLambertz = randVectorLamberts();
	// find rotation axis that rotates the z-axis to the normalVector of the surface element
	Vector3d axis = normalVector.cross(Vector3d(0, 0, 1));
	if (axis.getR() < std::numeric_limits<double>::epsilon()) {
		axis = Vector3d(0, 0, 1);
	}
	double angle = normalVector.getAngleTo(Vector3d(0, 0, 1));
	// rotate the random Lamberts vector from z-axis to respective surface element
	return vLambertz.getRotated(axis / axis.getR(), -angle);
}

Vector3d Random::randomInterpolatedPosition(const Vector3d &a, const Vector3d &b) {
	return a + rand() * (b - a);
}

double Random::randPowerLaw(double index, double min, double max) {
	if ((min < 0) || (max < min)) {
		throw std::runtime_error(
				"Power law distribution only possible for 0 <= min <= max");
	}
	//check for index -1!
	if ((std::abs(index + 1.0)) < std::numeric_limits<double>::epsilon()) {
		double part1 = log(max);
		double part2 = log(min);
		return exp((part1 - part2) * rand() + part2);
	} else {
		double part1 = pow(max, index + 1);
		double part2 = pow(min, index + 1);
		double ex = 1 / (index + 1);
		return pow((part1 - part2) * rand() + part2, ex);
	}
}

double Random::randBrokenPowerLaw(double index1, double index2,
		double breakpoint, double min, double max) {
	if ((min <= 0) || (max < min)) {
		throw std::runtime_error(
				"Power law distribution only possible for 0 < min <= max");
	}
	if (min >= breakpoint) {
		return this->randPowerLaw(index2, min, max);
	} else if (max <= breakpoint) {
		return this->randPowerLaw(index2, min, max);
	} else {
		double intPL1;
		// check if index1 = -1
		if ((std::abs(index1 + 1.0)) < std::numeric_limits<double>::epsilon()) {
			intPL1 = log(breakpoint / min);
		} else {
			intPL1 = (pow(breakpoint, index1 + 1) - pow(min, index1 + 1))
					/ (index1 + 1);
		}
		double intPL2;
		// check if index2 = -1
		if ((std::abs(index2 + 1.0)) < std::numeric_limits<double>::epsilon()) {
			intPL2 = log(max / breakpoint) * pow(breakpoint, index1 - index2);
		} else {
			intPL2 = (pow(max, index2 + 1) - pow(breakpoint, index2 + 1))
					* pow(breakpoint, index1 - index2) / (index2 + 1);
		}
		if (rand() > intPL1 / (intPL1 + intPL2))
			return randPowerLaw(index2, breakpoint, max);
		else
			return randPowerLaw(index1, min, breakpoint);
	}
}

double Random::randExponential() {
	double dum;
	do {
		dum = rand();
	} while (dum < std::numeric_limits<double>::epsilon());
	return -1.0 * log(dum);
}

uint32_t Random::randInt() {
	if (left == 0)
		reload();
	--left;

	uint32_t s1;
	s1 = *pNext++;
	s1 ^= (s1 >> 11);
	s1 ^= (s1 << 7) & 0x9d2c5680UL;
	s1 ^= (s1 << 15) & 0xefc60000UL;
	return (s1 ^ (s1 >> 18));
}

uint32_t Random::randInt(const uint32_t& n) {
// Find which bits are used in n
// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
	uint32_t used = n;
	used |= used >> 1;
	used |= used >> 2;
	used |= used >> 4;
	used |= used >> 8;
	used |= used >> 16;

// Draw numbers until one is found in [0,n]
	uint32_t i;
	do
		i = randInt() & used; // toss unused bits to shorten search
	while (i > n);
	return i;
}


uint64_t Random::randInt64()
{
	int64_t a = randInt();
	int64_t b = randInt();
	return (b + a << 32);
}


uint64_t Random::randInt64(const uint64_t &n)
{
	uint64_t used = n;
	used |= used >> 1;
	used |= used >> 2;
	used |= used >> 4;
	used |= used >> 8;
	used |= used >> 16;
	used |= used >> 32;

	// Draw numbers until one is found in [0,n]
	uint64_t i;
	do
		i = randInt64() & used; // toss unused bits to shorten search
	while (i > n);
	return i;
}



void Random::seed(const uint32_t oneSeed) {
	initial_seed.resize(1);
	initial_seed[0] = oneSeed;
	initialize(oneSeed);
	reload();
}

void Random::seed(uint32_t * const bigSeed, const uint32_t seedLength) {

	initial_seed.resize(seedLength);
	for (size_t i =0; i< seedLength; i++)
	{
		initial_seed[i] = bigSeed[i];
	}

	initialize(19650218UL);
	int i = 1;
	uint32_t j = 0;
	int k = (N > seedLength ? N : seedLength);
	for (; k; --k) {
		state[i] = state[i]
				^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL);
		state[i] += (bigSeed[j] & 0xffffffffUL) + j;
		state[i] &= 0xffffffffUL;
		++i;
		++j;
		if (i >= N) {
			state[0] = state[N - 1];
			i = 1;
		}
		if (j >= seedLength)
			j = 0;
	}
	for (k = N - 1; k; --k) {
		state[i] = state[i]
				^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL);
		state[i] -= i;
		state[i] &= 0xffffffffUL;
		++i;
		if (i >= N) {
			state[0] = state[N - 1];
			i = 1;
		}
	}
	state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
	reload();
}

void Random::seed() {
// First try getting an array from /dev/urandom
	FILE* urandom = std::fopen("/dev/urandom", "rb");
	if (urandom) {
		uint32_t bigSeed[N];
		uint32_t *s = bigSeed;
		int i = N;
		bool success = true;
		while (success && i--)
			success = std::fread(s++, sizeof(uint32_t), 1, urandom) != 0;
		std::fclose(urandom);
		if (success) {
			seed(bigSeed, N);
			return;
		}
	}

// Was not successful, so use time() and clock() instead
	seed(hash(time(NULL), clock()));
}


void Random::initialize(const uint32_t seed) {
	uint32_t *s = state;
	uint32_t *r = state;
	int i = 1;
	*s++ = seed & 0xffffffffUL;
	for (; i < N; ++i) {
		*s++ = (1812433253UL * (*r ^ (*r >> 30)) + i) & 0xffffffffUL;
		r++;
	}
}

void Random::reload() {
	uint32_t *p = state;
	int i;
	for (i = N - M; i--; ++p)
		*p = twist(p[M], p[0], p[1]);
	for (i = M; --i; ++p)
		*p = twist(p[M - N], p[0], p[1]);
	*p = twist(p[M - N], p[0], state[0]);

	left = N, pNext = state;
}

uint32_t Random::hash(time_t t, clock_t c) {
	static uint32_t differ = 0; // guarantee time-based seeds will change

	uint32_t h1 = 0;
	unsigned char *p = (unsigned char *) &t;
	for (size_t i = 0; i < sizeof(t); ++i) {
		h1 *= std::numeric_limits<unsigned char>::max() + 2U;
		h1 += p[i];
	}
	uint32_t h2 = 0;
	p = (unsigned char *) &c;
	for (size_t j = 0; j < sizeof(c); ++j) {
		h2 *= std::numeric_limits<unsigned char>::max() + 2U;
		h2 += p[j];
	}
	return (h1 + differ++) ^ h2;
}

void Random::save(uint32_t* saveArray) const {
	uint32_t *sa = saveArray;
	const uint32_t *s = state;
	int i = N;
	for (; i--; *sa++ = *s++) {
	}
	*sa = left;
}

const std::vector<uint32_t> &Random::getSeed() const
{
	return initial_seed;
}

void Random::load(uint32_t * const loadArray) {
	uint32_t *s = state;
	uint32_t *la = loadArray;
	int i = N;
	for (; i--; *s++ = *la++) {
	}
	left = *la;
	pNext = &state[N - left];
}

std::ostream& operator<<(std::ostream& os, const Random& mtrand) {
	const uint32_t *s = mtrand.state;
	int i = mtrand.N;
	for (; i--; os << *s++ << "\t") {
	}
	return os << mtrand.left;
}

std::istream& operator>>(std::istream& is, Random& mtrand) {
	uint32_t *s = mtrand.state;
	int i = mtrand.N;
	for (; i--; is >> *s++) {
	}
	is >> mtrand.left;
	mtrand.pNext = &mtrand.state[mtrand.N - mtrand.left];
	return is;
}

#ifdef _OPENMP
#include <omp.h>
#include <stdexcept>

// see http://stackoverflow.com/questions/8051108/using-the-openmp-threadprivate-directive-on-static-instances-of-c-stl-types
const static int MAX_THREAD = 256;

struct RANDOM_TLS_ITEM {
	Random r;
	char padding[(sizeof(Random) / 64 + 1) * 64 - sizeof(Random)];
};

#ifdef _MSC_VER
__declspec(align(64)) static RANDOM_TLS_ITEM _tls[MAX_THREAD];
#else
__attribute__ ((aligned(64))) static RANDOM_TLS_ITEM _tls[MAX_THREAD];
#endif

Random &Random::instance() {
	int i = omp_get_thread_num();
	if (i >= MAX_THREAD)
	throw std::runtime_error("crpropa::Random: more than MAX_THREAD threads!");
	return _tls[i].r;
}

void Random::seedThreads(const uint32_t oneSeed) {
	for(size_t i = 0; i < MAX_THREAD; ++i)
	_tls[i].r.seed(oneSeed + i);
}

std::vector< std::vector<uint32_t> > Random::getSeedThreads()
{
	std::vector< std::vector<uint32_t> > seeds;
	for(size_t i = 0; i < omp_get_num_threads(); ++i)
		seeds.push_back(_tls[i].r.getSeed() ); 
	return seeds;
}

#else
static Random _random;
Random &Random::instance() {
	return _random;
}
void Random::seedThreads(const uint32_t oneSeed) {
	_random.seed(oneSeed);
}
std::vector< std::vector<uint32_t> > Random::getSeedThreads()
{
	std::vector< std::vector<uint32_t> > seeds;
		seeds.push_back(_random.getSeed() ); 
	return seeds;
}
#endif

const std::string Random::getSeed_base64() const
{
	return Base64::encode((unsigned char*) &initial_seed[0], sizeof(initial_seed[0]) * initial_seed.size() / sizeof(unsigned char));
}

void Random::seed(const std::string &b64Seed)
{
	std::string decoded_data = Base64::decode(b64Seed);
	size_t seedSize = decoded_data.size() * sizeof(decoded_data[0]) / sizeof(uint32_t);
	seed((uint32_t*)decoded_data.c_str(), seedSize );
}

} // namespace crpropa