Line data    Source code

       1             : /****************************************************************************/
       2             : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
       3             : // Copyright (C) 2005-2024 German Aerospace Center (DLR) and others.
       4             : // This program and the accompanying materials are made available under the
       5             : // terms of the Eclipse Public License 2.0 which is available at
       6             : // https://www.eclipse.org/legal/epl-2.0/
       7             : // This Source Code may also be made available under the following Secondary
       8             : // Licenses when the conditions for such availability set forth in the Eclipse
       9             : // Public License 2.0 are satisfied: GNU General Public License, version 2
      10             : // or later which is available at
      11             : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
      12             : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
      13             : /****************************************************************************/
      14             : /// @file    RandHelper.h
      15             : /// @author  Daniel Krajzewicz
      16             : /// @author  Michael Behrisch
      17             : /// @author  Jakob Erdmann
      18             : /// @date    Fri, 29.04.2005
      19             : ///
      20             : //
      21             : /****************************************************************************/
      22             : #pragma once
      23             : #include <config.h>
      24             : 
      25             : #include <cassert>
      26             : #include <vector>
      27             : #include <map>
      28             : #include <random>
      29             : #include <sstream>
      30             : #include <iostream>
      31             : #include <algorithm>
      32             : 
      33             : 
      34             : // ===========================================================================
      35             : // class declaration
      36             : // ===========================================================================
      37             : 
      38             : class OptionsCont;
      39             : 
      40             : // ===========================================================================
      41             : // class definitions
      42             : // ===========================================================================
      43             : /**
      44             :  * @class XoShiRo256PlusPlus
      45             :  * @brief A random number generator as proposed here https://prng.di.unimi.it/xoshiro256plusplus.c
      46             :  */
      47             : class XoShiRo256PlusPlus {
      48             : public:
      49             :     inline void seed(const uint64_t seed)     {
      50             :         this->state[0] = splitmix64(splitmix64(seed));
      51             :         this->state[1] = splitmix64(this->state[0]);
      52             :         this->state[2] = splitmix64(this->state[1]);
      53             :         this->state[3] = splitmix64(this->state[2]);
      54             :     }
      55             : 
      56             :     inline uint64_t operator()()     {
      57             :         const uint64_t result = rotl64(this->state[0] + this->state[3], 23) + this->state[0];
      58             :         const uint64_t t = this->state[1] << 17;
      59             :         this->state[2] ^= this->state[0];
      60             :         this->state[3] ^= this->state[1];
      61             :         this->state[1] ^= this->state[2];
      62             :         this->state[0] ^= this->state[3];
      63             :         this->state[2] ^= t;
      64             :         this->state[3] = rotl64(this->state[3], 45);
      65             :         return result;
      66             :     }
      67             : 
      68             :     void discard(unsigned long long skip) {
      69             :         while (skip-- > 0) {
      70             :             (*this)();
      71             :         }
      72             :     }
      73             : 
      74             :     friend std::ostream& operator<<(std::ostream& os, const XoShiRo256PlusPlus& r) {
      75             :         os << r.state[0] << r.state[1] << r.state[2] << r.state[3];
      76             :         return os;
      77             :     }
      78             : 
      79             :     friend std::istream& operator>>(std::istream& is, XoShiRo256PlusPlus& r) {
      80             :         is >> r.state[0] >> r.state[1] >> r.state[2] >> r.state[3];
      81             :         return is;
      82             :     }
      83             : 
      84             : 
      85             : private:
      86             :     static inline uint64_t rotl64(const uint64_t x, const int k) {
      87             :         return (x << k) | (x >> (64 - k));
      88             :     }
      89             : 
      90             :     static inline uint64_t splitmix64(const uint64_t seed) {
      91             :         uint64_t z = (seed + 0x9e3779b97f4a7c15);
      92             :         z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
      93             :         z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
      94             :         return z ^ (z >> 31);
      95             :     }
      96             : 
      97             :     uint64_t state[4];
      98             : 
      99             : };
     100             : 
     101             : 
     102             : //class SumoRNG : public XoShiRo256PlusPlus {
     103     6844225 : class SumoRNG : public std::mt19937 {
     104             : public:
     105     5167148 :     SumoRNG(const std::string& _id) : id(_id) {}
     106             : 
     107             :     unsigned long long int count = 0;
     108             :     std::string id;
     109             : };
     110             : 
     111             : 
     112             : /**
     113             :  * @class RandHelper
     114             :  * @brief Utility functions for using a global, resetable random number generator
     115             :  */
     116             : class RandHelper {
     117             : 
     118             : public:
     119             :     /// @brief Initialises the given options container with random number options
     120             :     static void insertRandOptions(OptionsCont& oc);
     121             : 
     122             :     /// @brief Initialises the random number generator with hardware randomness or seed
     123             :     static void initRand(SumoRNG* which = nullptr, const bool random = false, const int seed = 23423);
     124             : 
     125             :     /// @brief Reads the given random number options and initialises the random number generator in accordance
     126             :     static void initRandGlobal(SumoRNG* which = nullptr);
     127             : 
     128             :     /// @brief Returns a random real number in [0, 1)
     129             :     static double rand(SumoRNG* rng = nullptr);
     130             : 
     131             :     /// @brief Returns a random real number in [0, maxV)
     132             :     static inline double rand(double maxV, SumoRNG* rng = nullptr) {
     133    50957878 :         return maxV * rand(rng);
     134             :     }
     135             : 
     136             :     /// @brief Returns a random real number in [minV, maxV)
     137             :     static inline double rand(double minV, double maxV, SumoRNG* rng = nullptr) {
     138      235514 :         return minV + (maxV - minV) * rand(rng);
     139             :     }
     140             : 
     141             :     /// @brief Returns a random integer in [0, maxV-1]
     142   151669906 :     static inline int rand(int maxV, SumoRNG* rng = nullptr) {
     143   151730351 :         if (rng == nullptr) {
     144             :             rng = &myRandomNumberGenerator;
     145             :         }
     146   151669906 :         unsigned int usedBits = maxV - 1;
     147   151669906 :         usedBits |= usedBits >> 1;
     148   151669906 :         usedBits |= usedBits >> 2;
     149   151669906 :         usedBits |= usedBits >> 4;
     150   151669906 :         usedBits |= usedBits >> 8;
     151   151669906 :         usedBits |= usedBits >> 16;
     152             : 
     153             :         // Draw numbers until one is found in [0, maxV-1]
     154             :         int result;
     155             :         do {
     156   151955017 :             result = (*rng)() & usedBits;
     157   151955017 :             rng->count++;
     158   151887133 :         } while (result >= maxV);
     159   151669906 :         return result;
     160             :     }
     161             : 
     162             :     /// @brief Returns a random integer in [minV, maxV-1]
     163             :     static inline int rand(int minV, int maxV, SumoRNG* rng = nullptr) {
     164       89175 :         return minV + rand(maxV - minV, rng);
     165             :     }
     166             : 
     167             :     /// @brief Returns a random 64 bit integer in [0, maxV-1]
     168       14272 :     static inline long long int rand(long long int maxV, SumoRNG* rng = nullptr) {
     169       14272 :         if (maxV <= std::numeric_limits<int>::max()) {
     170       14272 :             return rand((int)maxV, rng);
     171             :         }
     172           0 :         if (rng == nullptr) {
     173             :             rng = &myRandomNumberGenerator;
     174             :         }
     175           0 :         unsigned long long int usedBits = maxV - 1;
     176           0 :         usedBits |= usedBits >> 1;
     177           0 :         usedBits |= usedBits >> 2;
     178           0 :         usedBits |= usedBits >> 4;
     179           0 :         usedBits |= usedBits >> 8;
     180           0 :         usedBits |= usedBits >> 16;
     181           0 :         usedBits |= usedBits >> 32;
     182             : 
     183             :         // Draw numbers until one is found in [0, maxV-1]
     184             :         long long int result;
     185             :         do {
     186           0 :             result = (((unsigned long long int)(*rng)() << 32) | (*rng)()) & usedBits;    // toss unused bits to shorten search
     187           0 :             rng->count += 2;
     188           0 :         } while (result >= maxV);
     189             :         return result;
     190             :     }
     191             : 
     192             :     /// @brief Returns a random 64 bit integer in [minV, maxV-1]
     193             :     static inline long long int rand(long long int minV, long long int maxV, SumoRNG* rng = nullptr) {
     194       14032 :         return minV + rand(maxV - minV, rng);
     195             :     }
     196             : 
     197             :     /// @brief Access to a random number from a normal distribution
     198             :     static double randNorm(double mean, double variance, SumoRNG* rng = nullptr);
     199             : 
     200             :     /// @brief Access to a random number from an exponential distribution
     201             :     static double randExp(double rate, SumoRNG* rng = nullptr);
     202             : 
     203             :     /// @brief Returns a random element from the given vector
     204             :     template<class T>
     205             :     static inline const T&
     206             :     getRandomFrom(const std::vector<T>& v, SumoRNG* rng = nullptr) {
     207             :         assert(v.size() > 0);
     208   151495383 :         return v[rand((int)v.size(), rng)];
     209             :     }
     210             : 
     211             :     /// @brief save rng state to string
     212         639 :     static std::string saveState(SumoRNG* rng = nullptr) {
     213         639 :         if (rng == nullptr) {
     214             :             rng = &myRandomNumberGenerator;
     215             :         }
     216         639 :         std::ostringstream oss;
     217         639 :         if (rng->count < 1000000) { // TODO make this configurable
     218             :             oss << rng->count;
     219             :         } else {
     220           0 :             oss << (*rng);
     221             :         }
     222         639 :         return oss.str();
     223         639 :     }
     224             : 
     225             :     /// @brief load rng state from string
     226         639 :     static void loadState(const std::string& state, SumoRNG* rng = nullptr) {
     227         639 :         if (rng == nullptr) {
     228             :             rng = &myRandomNumberGenerator;
     229             :         }
     230         639 :         std::istringstream iss(state);
     231         639 :         if (state.size() < 10) {
     232         639 :             iss >> rng->count;
     233         639 :             rng->discard(rng->count);
     234             :         } else {
     235           0 :             iss >> (*rng);
     236             :         }
     237         639 :     }
     238             : 
     239             :     template<class T>
     240           0 :     static void shuffle(std::vector<T>& v, SumoRNG* rng = nullptr) {
     241           0 :         for (int i = (int)(v.size() - 1); i > 0; --i) {
     242           0 :             std::swap(*(v.begin() + i), *(v.begin() + rand(i, rng)));
     243             :         }
     244           0 :     }
     245             : 
     246             : protected:
     247             :     /// @brief the default random number generator to use
     248             :     static SumoRNG myRandomNumberGenerator;
     249             : 
     250             : };