RandHelper.h

Go to the documentation of this file.

/****************************************************************************/
// Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
// Copyright (C) 2005-2024 German Aerospace Center (DLR) and others.
// This program and the accompanying materials are made available under the
// terms of the Eclipse Public License 2.0 which is available at
// https://www.eclipse.org/legal/epl-2.0/
// This Source Code may also be made available under the following Secondary
// Licenses when the conditions for such availability set forth in the Eclipse
// Public License 2.0 are satisfied: GNU General Public License, version 2
// or later which is available at
// https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
/****************************************************************************/
//
/****************************************************************************/
#pragma once
#include <config.h>
 
#include <cassert>
#include <vector>
#include <map>
#include <random>
#include <sstream>
#include <iostream>
#include <algorithm>
 
 
// ===========================================================================
// class declaration
// ===========================================================================
 
class OptionsCont;
 
// ===========================================================================
// class definitions
// ===========================================================================
class XoShiRo256PlusPlus {
public:
    inline void seed(const uint64_t seed)     {
        this->state[0] = splitmix64(splitmix64(seed));
        this->state[1] = splitmix64(this->state[0]);
        this->state[2] = splitmix64(this->state[1]);
        this->state[3] = splitmix64(this->state[2]);
    }
 
    inline uint64_t operator()()     {
        const uint64_t result = rotl64(this->state[0] + this->state[3], 23) + this->state[0];
        const uint64_t t = this->state[1] << 17;
        this->state[2] ^= this->state[0];
        this->state[3] ^= this->state[1];
        this->state[1] ^= this->state[2];
        this->state[0] ^= this->state[3];
        this->state[2] ^= t;
        this->state[3] = rotl64(this->state[3], 45);
        return result;
    }
 
    void discard(unsigned long long skip) {
        while (skip-- > 0) {
            (*this)();
        }
    }
 
    friend std::ostream& operator<<(std::ostream& os, const XoShiRo256PlusPlus& r) {
        os << r.state[0] << r.state[1] << r.state[2] << r.state[3];
        return os;
    }
 
    friend std::istream& operator>>(std::istream& is, XoShiRo256PlusPlus& r) {
        is >> r.state[0] >> r.state[1] >> r.state[2] >> r.state[3];
        return is;
    }
 
 
private:
    static inline uint64_t rotl64(const uint64_t x, const int k) {
        return (x << k) | (x >> (64 - k));
    }
 
    static inline uint64_t splitmix64(const uint64_t seed) {
        uint64_t z = (seed + 0x9e3779b97f4a7c15);
        z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
        z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
        return z ^ (z >> 31);
    }
 
    uint64_t state[4];
 
};
 
 
//class SumoRNG : public XoShiRo256PlusPlus {
class SumoRNG : public std::mt19937 {
public:
    SumoRNG(const std::string& _id) : id(_id) {}
 
    unsigned long long int count = 0;
    std::string id;
};
 
 
class RandHelper {
 
public:
    static void insertRandOptions(OptionsCont& oc);
 
    static void initRand(SumoRNG* which = nullptr, const bool random = false, const int seed = 23423);
 
    static void initRandGlobal(SumoRNG* which = nullptr);
 
    static double rand(SumoRNG* rng = nullptr);
 
    static inline double rand(double maxV, SumoRNG* rng = nullptr) {
        return maxV * rand(rng);
    }
 
    static inline double rand(double minV, double maxV, SumoRNG* rng = nullptr) {
        return minV + (maxV - minV) * rand(rng);
    }
 
    static inline int rand(int maxV, SumoRNG* rng = nullptr) {
        if (rng == nullptr) {
            rng = &myRandomNumberGenerator;
        }
        unsigned int usedBits = maxV - 1;
        usedBits |= usedBits >> 1;
        usedBits |= usedBits >> 2;
        usedBits |= usedBits >> 4;
        usedBits |= usedBits >> 8;
        usedBits |= usedBits >> 16;
 
        // Draw numbers until one is found in [0, maxV-1]
        int result;
        do {
            result = (*rng)() & usedBits;
            rng->count++;
        } while (result >= maxV);
        return result;
    }
 
    static inline int rand(int minV, int maxV, SumoRNG* rng = nullptr) {
        return minV + rand(maxV - minV, rng);
    }
 
    static inline long long int rand(long long int maxV, SumoRNG* rng = nullptr) {
        if (maxV <= std::numeric_limits<int>::max()) {
            return rand((int)maxV, rng);
        }
        if (rng == nullptr) {
            rng = &myRandomNumberGenerator;
        }
        unsigned long long int usedBits = maxV - 1;
        usedBits |= usedBits >> 1;
        usedBits |= usedBits >> 2;
        usedBits |= usedBits >> 4;
        usedBits |= usedBits >> 8;
        usedBits |= usedBits >> 16;
        usedBits |= usedBits >> 32;
 
        // Draw numbers until one is found in [0, maxV-1]
        long long int result;
        do {
            result = (((unsigned long long int)(*rng)() << 32) | (*rng)()) & usedBits;    // toss unused bits to shorten search
            rng->count += 2;
        } while (result >= maxV);
        return result;
    }
 
    static inline long long int rand(long long int minV, long long int maxV, SumoRNG* rng = nullptr) {
        return minV + rand(maxV - minV, rng);
    }
 
    static double randNorm(double mean, double variance, SumoRNG* rng = nullptr);
 
    static double randExp(double rate, SumoRNG* rng = nullptr);
 
    template<class T>
    static inline const T&
    getRandomFrom(const std::vector<T>& v, SumoRNG* rng = nullptr) {
        assert(v.size() > 0);
        return v[rand((int)v.size(), rng)];
    }
 
    static std::string saveState(SumoRNG* rng = nullptr) {
        if (rng == nullptr) {
            rng = &myRandomNumberGenerator;
        }
        std::ostringstream oss;
        if (rng->count < 1000000) { // TODO make this configurable
            oss << rng->count;
        } else {
            oss << (*rng);
        }
        return oss.str();
    }
 
    static void loadState(const std::string& state, SumoRNG* rng = nullptr) {
        if (rng == nullptr) {
            rng = &myRandomNumberGenerator;
        }
        std::istringstream iss(state);
        if (state.size() < 10) {
            iss >> rng->count;
            rng->discard(rng->count);
        } else {
            iss >> (*rng);
        }
    }
 
    template<class T>
    static void shuffle(std::vector<T>& v, SumoRNG* rng = nullptr) {
        for (int i = (int)(v.size() - 1); i > 0; --i) {
            std::swap(*(v.begin() + i), *(v.begin() + rand(i, rng)));
        }
    }
 
protected:
    static SumoRNG myRandomNumberGenerator;
 
};