Program Listing for File Random.cpp

Return to documentation for file (src/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