include/algorithms/random.h

0001 #pragma once
0002
0003 #include <algorithm>
0004 #include <concepts>
0005 #include <cstdint>
0006 #include <functional>
0007 #include <mutex>
0008 #include <random>
0009
0010 #include <algorithms/detail/random.h>
0011 #include <algorithms/logger.h>
0012 #include <algorithms/service.h>
0013
0014 namespace algorithms {
0015
0016 // Random Engine callback function:
0017 //   - Signature: std::function<std::vector<value_type>(size_t N)> --> generates a vector
0018 //     of N numbers
0019 //   - RandomEngineCB is required to return N random numbers between 0 and
0020 //     std::numeric_limits<uint_fast64_t>::max()
0021 //   - RandomEngineCB is responsible to deal with possible simultaneous access by multiple
0022 //     Generator instances (required to be thread-safe).
0023 using RandomEngineCB = detail::CachedBitGenerator::GenFunc;
0024
0025 // thread-safe generator front-end. Requires that the underlying random engine used by
0026 // the RandomSvc is thread-safe.
0027 class Generator {
0028 public:
0029   Generator(const RandomEngineCB& gen, const size_t cache_size) : m_gen{gen, cache_size} {}
0030
0031   template <std::integral Int = int> Int uniform_int(const Int min = Int{0}, const Int max = Int{1}) const {
0032     std::uniform_int_distribution<Int> d{min, max};
0033     std::lock_guard<std::mutex> lock{m_mutex};
0034     return d(m_gen);
0035   }
0036   template <std::floating_point Float = double> Float uniform_double(const Float min = Float{0}, const Float max = Float{1}) const {
0037     std::uniform_real_distribution<Float> d{min, max};
0038     std::lock_guard<std::mutex> lock{m_mutex};
0039     return d(m_gen);
0040   }
0041   template <std::integral Int = int> Int poisson(const Int mean = Int{1}) const {
0042     std::poisson_distribution<Int> d{mean};
0043     std::lock_guard<std::mutex> lock{m_mutex};
0044     return d(m_gen);
0045   }
0046   template <std::floating_point Float = double> Float exponential(const Float lambda = Float{1}) const {
0047     std::exponential_distribution<Float> d{lambda};
0048     std::lock_guard<std::mutex> lock{m_mutex};
0049     return d(m_gen);
0050   }
0051   template <std::floating_point Float = double> Float gaussian(const Float mu = Float{0}, const Float sigma = Float{1}) const {
0052     std::normal_distribution<Float> d{mu, sigma};
0053     std::lock_guard<std::mutex> lock{m_mutex};
0054     return d(m_gen);
0055   }
0056
0057 private:
0058   mutable detail::CachedBitGenerator m_gen;
0059   mutable std::mutex m_mutex;
0060 };
0061
0062 // Random service that creates multiple Generators that are linked to a single random
0063 // engine. The Generators are safe to be used in parallel as long as the Engine itself is
0064 // thread-safe (this is a hard requirement for MT). The Generators avoid unnecesary locking
0065 // by running off a auto-refreshing cached random sequence.
0066 class RandomSvc : public LoggedService<RandomSvc> {
0067 public:
0068   using value_type = detail::CachedBitGenerator::result_type;
0069
0070   Generator generator() { return {m_gen, m_cache_size}; }
0071 // FIXME fix the CMake setup so these are properly found in Gaudi
0072 #if 0
0073   void init();
0074   void init(const RandomEngineCB& gen);
0075 #endif
0076   void init() {
0077     if (m_seed.hasValue()) {
0078       info() << "Custom random seed requested: " << m_seed << endmsg;
0079       m_gen = createEngine(m_seed);
0080     }
0081   }
0082   void init(const RandomEngineCB& gen) {
0083     info() << "Loading external generator function." << endmsg;
0084     m_gen = gen;
0085     if (m_seed.hasValue()) {
0086       warning() << "Custom random seed request ignored when using external generator function"
0087                 << endmsg;
0088     }
0089   }
0090
0091 #if 0
0092 private:
0093   RandomEngineCB createEngine(const size_t seed = 1);
0094 #endif
0095   RandomEngineCB createEngine(const size_t seed = 1) {
0096     return [=](const size_t size) {
0097       static std::mutex m;
0098       static std::mt19937_64 gen{seed};
0099       std::lock_guard<std::mutex> lock{m};
0100       std::vector<value_type> ret(size);
0101       std::generate(ret.begin(), ret.end(), gen);
0102       return ret;
0103     };
0104   }
0105   // end of FIXME
0106
0107 private:
0108   RandomEngineCB m_gen{createEngine()};
0109   Property<size_t> m_seed{this, "seed", "Random seed for the internal random engine"};
0110   Property<size_t> m_cache_size{this, "cacheSize", 1024, "Cache size for each generator instance"};
0111   std::mutex m_mutex;
0112
0113   ALGORITHMS_DEFINE_LOGGED_SERVICE(RandomSvc)
0114 };
0115
0116 } // namespace algorithms