VecGeom/base/RNG.h

0001 /// \file RNG.h
0002 /// \author Johannes de Fine Licht (johannes.definelicht@cern.ch)
0003
0004 #ifndef VECGEOM_BASE_RNG_H_
0005 #define VECGEOM_BASE_RNG_H_
0006
0007 #include "VecGeom/base/Global.h"
0008
0009 #include <cstdint>
0010 #include <random>
0011
0012 namespace vecgeom {
0013
0014 /**
0015  * @brief Thread-local singleton random number generator on CPU.
0016  */
0017 class RNG {
0018
0019 private:
0020   std::mt19937 rng;
0021   std::uniform_real_distribution<> uniform_dist;
0022
0023   VECGEOM_FORCE_INLINE
0024   Precision GetUniform() { return uniform_dist(rng); }
0025
0026 public:
0027   RNG() : rng(0), uniform_dist(0, 1) {}
0028
0029 public:
0030   using Seed_t = std::mt19937::result_type;
0031
0032   void seed(unsigned long seed_val) { rng.seed(static_cast<Seed_t>(seed_val)); }
0033
0034   /**
0035    * Build a deterministic seed for one logical stream starting from a base
0036    * seed. The stream id is part of the contract so callers can reproduce a
0037    * failing randomized test independent of which worker thread executes it.
0038    */
0039   static Seed_t MakeStreamSeed(unsigned long base_seed, unsigned long stream_id = 0)
0040   {
0041     std::uint64_t z = static_cast<std::uint64_t>(base_seed);
0042     z += 0x9e3779b97f4a7c15ULL + static_cast<std::uint64_t>(stream_id);
0043     z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
0044     z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
0045     z = z ^ (z >> 31);
0046     return static_cast<Seed_t>(z);
0047   }
0048
0049   /// Seed an arbitrary RNG instance from a deterministic base-seed/stream-id pair.
0050   static void SeedEngine(RNG &engine, unsigned long base_seed, unsigned long stream_id = 0)
0051   {
0052     engine.seed(MakeStreamSeed(base_seed, stream_id));
0053   }
0054
0055   /**
0056    * Access the calling thread's singleton RNG instance.
0057    */
0058   static RNG &Instance()
0059   {
0060     // The singleton is thread-local on CPU so randomized helpers can be used
0061     // safely in parallel without racing on one shared generator.
0062     static thread_local RNG instance;
0063     return instance;
0064   }
0065
0066   /// Seed the calling thread's RNG stream deterministically.
0067   static void SeedStream(unsigned long base_seed, unsigned long stream_id = 0)
0068   {
0069     SeedEngine(Instance(), base_seed, stream_id);
0070   }
0071
0072   /**
0073    * @return Uniformly distributed floating point number between 0 and 1 unless
0074    *         range arguments are passed.
0075    */
0076   Precision uniform(const Precision min = 0., const Precision max = 1.) { return min + (max - min) * GetUniform(); }
0077
0078   int Poisson(const Precision lambda)
0079   {
0080     int k                  = 0;
0081     const Precision target = exp(-lambda);
0082     Precision p            = GetUniform();
0083     while (p < target) {
0084       p *= GetUniform();
0085       ++k;
0086     }
0087     return k;
0088   }
0089
0090   // interface for ROOT compatibility
0091   Precision Gaus(Precision ave = 0.0, Precision sig = 1.0) { return Gauss(ave, sig); }
0092
0093   Precision Gauss(Precision ave = 0.0, Precision sig = 1.0)
0094   {
0095     Precision x1, x2, w;
0096
0097     do {
0098       x1 = 2.0 * GetUniform() - 1.0;
0099       x2 = 2.0 * GetUniform() - 1.0;
0100       w  = x1 * x1 + x2 * x2;
0101     } while (w >= 1.0);
0102
0103     w = std::sqrt((-2.0 * std::log(w)) / w);
0104     return ave + (x1 * w * sig);
0105   }
0106
0107   /**
0108    * Uniformly distributed array of floating point number between 0 and 1 unless
0109    *         range arguments are passed.
0110    */
0111   void uniform_array(size_t n, Precision *array, const Precision min = 0., const Precision max = 1.)
0112   {
0113     for (size_t i = 0; i < n; ++i) {
0114       array[i] = min + (max - min) * GetUniform();
0115     }
0116   }
0117
0118 private:
0119   RNG(RNG const &);
0120   RNG &operator=(RNG const &);
0121 };
0122
0123 } // namespace vecgeom
0124
0125 #endif // VECGEOM_BASE_RNG_H_