From 08aece013cb577d803ee058400b79e64a0dcc9bb Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 14 Nov 2024 18:03:20 +0100 Subject: [PATCH] implement `RandomSeedSequence` --- .../ActsExamples/Framework/RandomNumbers.hpp | 188 +++++------------- .../Framework/src/Framework/RandomNumbers.cpp | 183 +++++++++++++++++ 2 files changed, 229 insertions(+), 142 deletions(-) diff --git a/Examples/Framework/include/ActsExamples/Framework/RandomNumbers.hpp b/Examples/Framework/include/ActsExamples/Framework/RandomNumbers.hpp index 2534117f2ec..44c6b53dde8 100644 --- a/Examples/Framework/include/ActsExamples/Framework/RandomNumbers.hpp +++ b/Examples/Framework/include/ActsExamples/Framework/RandomNumbers.hpp @@ -10,7 +10,8 @@ #include #include -#include + +#include namespace ActsExamples { struct AlgorithmContext; @@ -21,6 +22,44 @@ using RandomEngine = std::mt19937; ///< Mersenne Twister /// The seed type used in the framework. using RandomSeed = std::uint32_t; +class RandomSeedSequence { + public: + using result_type = std::uint_least32_t; + static constexpr std::size_t N = 4; + using container_type = boost::container::small_vector; + + RandomSeedSequence(); + explicit RandomSeedSequence(RandomSeed seed); + template + RandomSeedSequence(InputIterator first, InputIterator last) + : m_seed(first, last) {} + template + RandomSeedSequence(const InputRange& range) + : m_seed(range.begin(), range.end()) {} + RandomSeedSequence(std::initializer_list il) + : m_seed(il.begin(), il.end()) {} + + std::size_t size() const; + + container_type generate(std::size_t size) const; + + template + void generate(OutputIterator first, OutputIterator last) const { + container_type result = generate(static_cast(last - first)); + std::copy(result.begin(), result.end(), first); + } + + template + void param(OutputIterator first) const { + std::copy(m_seed.begin(), m_seed.end(), first); + } + + RandomSeedSequence append(RandomSeed seed) const; + + private: + container_type m_seed; +}; + /// Provide event and algorithm specific random number generator.s /// /// This provides local random number generators, allowing for @@ -40,6 +79,12 @@ class RandomNumbers { explicit RandomNumbers(const Config& cfg); + RandomSeedSequence createEventAlgorithmSeedSequence( + const AlgorithmContext& context) const; + + RandomEngine createEventAlgorithmEngine( + const AlgorithmContext& context) const; + /// Spawn an algorithm-local random number generator. To avoid inefficiencies /// and multiple uses of a given RNG seed, this should only be done once per /// Algorithm invocation, after what the generator object should be reused. @@ -59,145 +104,4 @@ class RandomNumbers { Config m_cfg; }; -/// Distribute entropy from a sequence of input values to a sequence of output -/// values. -/// -/// This function is copied from std::seed_seq::generate() and modified to avoid -/// the intermediate std::vector. -/// -/// @tparam InputIteratorBegin is the type of the input iterator begin -/// @tparam InputIteratorEnd is the type of the input iterator end -/// @tparam OutputIteratorBegin is the type of the output iterator begin -/// @tparam OutputIteratorEnd is the type of the output iterator end -/// -/// @param inputIteratorBegin is the begin of the input sequence -/// @param inputIteratorEnd is the end of the input sequence -/// @param outputIteratorBegin is the begin of the output sequence -/// @param outputIteratorEnd is the end of the output sequence -template -void distributeEntropy(InputIteratorBegin inputIteratorBegin, - InputIteratorEnd inputIteratorEnd, - OutputIteratorBegin outputIteratorBegin, - OutputIteratorEnd outputIteratorEnd) - requires((std::input_iterator && - std::random_access_iterator && - std::is_assignable_v) && - (std::input_iterator && - std::random_access_iterator && - std::is_assignable_v) && - (std::output_iterator && - std::random_access_iterator) && - (std::output_iterator && - std::random_access_iterator)) -{ - using result_type = std::uint_least32_t; - - if (inputIteratorBegin == inputIteratorEnd) { - return; - } - if (outputIteratorBegin == outputIteratorEnd) { - throw std::invalid_argument("outputIteratorBegin == outputIteratorEnd"); - } - - const auto __v_ = inputIteratorBegin; - const auto __v_size = - static_cast(inputIteratorEnd - inputIteratorBegin); - const auto __first = outputIteratorBegin; - const auto __last = outputIteratorEnd; - - std::fill(__first, __last, 0x8b8b8b8b); - const std::size_t __n = static_cast(__last - __first); - const std::size_t __s = __v_size; - const std::size_t __t = (__n >= 623) ? 11 - : (__n >= 68) ? 7 - : (__n >= 39) ? 5 - : (__n >= 7) ? 3 - : (__n - 1) / 2; - const std::size_t __p = (__n - __t) / 2; - const std::size_t __q = __p + __t; - const std::size_t __m = std::max(__s + 1, __n); - // __k = 0; - { - result_type __r = - 1664525 * _Tp(__first[0] ^ __first[__p] ^ __first[__n - 1]); - __first[__p] += __r; - __r += __s; - __first[__q] += __r; - __first[0] = __r; - } - // Initialize indexing terms used with if statements as an optimization to - // avoid calculating modulo n on every loop iteration for each term. - std::size_t __kmodn = 0; // __k % __n - std::size_t __k1modn = __n - 1; // (__k - 1) % __n - std::size_t __kpmodn = __p % __n; // (__k + __p) % __n - std::size_t __kqmodn = __q % __n; // (__k + __q) % __n - - for (std::size_t __k = 1; __k <= __s; ++__k) { - if (++__kmodn == __n) { - __kmodn = 0; - } - if (++__k1modn == __n) { - __k1modn = 0; - } - if (++__kpmodn == __n) { - __kpmodn = 0; - } - if (++__kqmodn == __n) { - __kqmodn = 0; - } - - result_type __r = - 1664525 * _Tp(__first[__kmodn] ^ __first[__kpmodn] ^ __first[__k1modn]); - __first[__kpmodn] += __r; - __r += __kmodn + __v_[__k - 1]; - __first[__kqmodn] += __r; - __first[__kmodn] = __r; - } - for (std::size_t __k = __s + 1; __k < __m; ++__k) { - if (++__kmodn == __n) { - __kmodn = 0; - } - if (++__k1modn == __n) { - __k1modn = 0; - } - if (++__kpmodn == __n) { - __kpmodn = 0; - } - if (++__kqmodn == __n) { - __kqmodn = 0; - } - - result_type __r = - 1664525 * _Tp(__first[__kmodn] ^ __first[__kpmodn] ^ __first[__k1modn]); - __first[__kpmodn] += __r; - __r += __kmodn; - __first[__kqmodn] += __r; - __first[__kmodn] = __r; - } - for (std::size_t __k = __m; __k < __m + __n; ++__k) { - if (++__kmodn == __n) { - __kmodn = 0; - } - if (++__k1modn == __n) { - __k1modn = 0; - } - if (++__kpmodn == __n) { - __kpmodn = 0; - } - if (++__kqmodn == __n) { - __kqmodn = 0; - } - - result_type __r = 1566083941 * _Tp(__first[__kmodn] + __first[__kpmodn] + - __first[__k1modn]); - __first[__kpmodn] ^= __r; - __r -= __kmodn; - __first[__kqmodn] ^= __r; - __first[__kmodn] = __r; - } -} - } // namespace ActsExamples diff --git a/Examples/Framework/src/Framework/RandomNumbers.cpp b/Examples/Framework/src/Framework/RandomNumbers.cpp index 4313a5cce77..4183f9b18c3 100644 --- a/Examples/Framework/src/Framework/RandomNumbers.cpp +++ b/Examples/Framework/src/Framework/RandomNumbers.cpp @@ -10,10 +10,193 @@ #include "ActsExamples/Framework/AlgorithmContext.hpp" +#include + +namespace { + +/// Distribute entropy from a sequence of input values to a sequence of output +/// values. +/// +/// This function is copied from std::seed_seq::generate() and modified to avoid +/// the intermediate std::vector. +/// +/// @tparam InputIteratorBegin is the type of the input iterator begin +/// @tparam InputIteratorEnd is the type of the input iterator end +/// @tparam OutputIteratorBegin is the type of the output iterator begin +/// @tparam OutputIteratorEnd is the type of the output iterator end +/// +/// @param inputIteratorBegin is the begin of the input sequence +/// @param inputIteratorEnd is the end of the input sequence +/// @param outputIteratorBegin is the begin of the output sequence +/// @param outputIteratorEnd is the end of the output sequence +template +void distributeEntropy(InputIteratorBegin inputIteratorBegin, + InputIteratorEnd inputIteratorEnd, + OutputIteratorBegin outputIteratorBegin, + OutputIteratorEnd outputIteratorEnd) + requires((std::input_iterator && + std::random_access_iterator) && + (std::input_iterator && + std::random_access_iterator) && + (std::output_iterator && + std::random_access_iterator) && + (std::output_iterator && + std::random_access_iterator)) +{ + using result_type = std::uint_least32_t; + + std::seed_seq __seq(inputIteratorBegin, inputIteratorEnd); + + if (inputIteratorBegin == inputIteratorEnd) { + return; + } + if (outputIteratorBegin == outputIteratorEnd) { + return; + } + + const auto __v_ = inputIteratorBegin; + const auto __v_size = + static_cast(inputIteratorEnd - inputIteratorBegin); + const auto __first = outputIteratorBegin; + const auto __last = outputIteratorEnd; + + auto _Tp = [](result_type __x) -> result_type { return __x ^ (__x >> 27); }; + + std::fill(__first, __last, 0x8b8b8b8b); + const std::size_t __n = static_cast(__last - __first); + const std::size_t __s = __v_size; + const std::size_t __t = (__n >= 623) ? 11 + : (__n >= 68) ? 7 + : (__n >= 39) ? 5 + : (__n >= 7) ? 3 + : (__n - 1) / 2; + const std::size_t __p = (__n - __t) / 2; + const std::size_t __q = __p + __t; + const std::size_t __m = std::max(__s + 1, __n); + // __k = 0; + { + result_type __r = + 1664525 * _Tp(__first[0] ^ __first[__p] ^ __first[__n - 1]); + __first[__p] += __r; + __r += __s; + __first[__q] += __r; + __first[0] = __r; + } + // Initialize indexing terms used with if statements as an optimization to + // avoid calculating modulo n on every loop iteration for each term. + std::size_t __kmodn = 0; // __k % __n + std::size_t __k1modn = __n - 1; // (__k - 1) % __n + std::size_t __kpmodn = __p % __n; // (__k + __p) % __n + std::size_t __kqmodn = __q % __n; // (__k + __q) % __n + + for (std::size_t __k = 1; __k <= __s; ++__k) { + if (++__kmodn == __n) { + __kmodn = 0; + } + if (++__k1modn == __n) { + __k1modn = 0; + } + if (++__kpmodn == __n) { + __kpmodn = 0; + } + if (++__kqmodn == __n) { + __kqmodn = 0; + } + + result_type __r = + 1664525 * _Tp(__first[__kmodn] ^ __first[__kpmodn] ^ __first[__k1modn]); + __first[__kpmodn] += __r; + __r += __kmodn + __v_[__k - 1]; + __first[__kqmodn] += __r; + __first[__kmodn] = __r; + } + for (std::size_t __k = __s + 1; __k < __m; ++__k) { + if (++__kmodn == __n) { + __kmodn = 0; + } + if (++__k1modn == __n) { + __k1modn = 0; + } + if (++__kpmodn == __n) { + __kpmodn = 0; + } + if (++__kqmodn == __n) { + __kqmodn = 0; + } + + result_type __r = + 1664525 * _Tp(__first[__kmodn] ^ __first[__kpmodn] ^ __first[__k1modn]); + __first[__kpmodn] += __r; + __r += __kmodn; + __first[__kqmodn] += __r; + __first[__kmodn] = __r; + } + for (std::size_t __k = __m; __k < __m + __n; ++__k) { + if (++__kmodn == __n) { + __kmodn = 0; + } + if (++__k1modn == __n) { + __k1modn = 0; + } + if (++__kpmodn == __n) { + __kpmodn = 0; + } + if (++__kqmodn == __n) { + __kqmodn = 0; + } + + result_type __r = 1566083941 * _Tp(__first[__kmodn] + __first[__kpmodn] + + __first[__k1modn]); + __first[__kpmodn] ^= __r; + __r -= __kmodn; + __first[__kqmodn] ^= __r; + __first[__kmodn] = __r; + } +} + +} // namespace + namespace ActsExamples { +RandomSeedSequence::RandomSeedSequence() = default; + +RandomSeedSequence::RandomSeedSequence(RandomSeed seed) { + m_seed = {seed}; +} + +std::size_t RandomSeedSequence::size() const { + return m_seed.size(); +} + +RandomSeedSequence::container_type RandomSeedSequence::generate( + std::size_t size) const { + container_type result(size, 0); + distributeEntropy(m_seed.begin(), m_seed.end(), result.begin(), result.end()); + return result; +} + +RandomSeedSequence RandomSeedSequence::append(RandomSeed seed) const { + RandomSeedSequence result; + result.m_seed = m_seed; + result.m_seed.push_back(seed); + return result; +} + RandomNumbers::RandomNumbers(const Config& cfg) : m_cfg(cfg) {} +RandomSeedSequence RandomNumbers::createEventAlgorithmSeedSequence( + const AlgorithmContext& context) const { + return RandomSeedSequence( + {m_cfg.seed, static_cast(context.eventNumber)}); +} + +RandomEngine RandomNumbers::createEventAlgorithmEngine( + const AlgorithmContext& context) const { + RandomSeedSequence seedSequence = createEventAlgorithmSeedSequence(context); + return RandomEngine(seedSequence); +} + RandomEngine RandomNumbers::spawnGenerator( const AlgorithmContext& context) const { return RandomEngine(generateSeed(context));