Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : // MOOSE includes
13 : #include "MooseError.h"
14 : #include "DataIO.h"
15 :
16 : // External library includes
17 : #include "randistrs.h"
18 :
19 : // STL includes
20 : #include <unordered_map>
21 :
22 : /**
23 : * A deterministic, indexable random number generator built on top of the
24 : * `randistrs` library.
25 : *
26 : * This class provides seed-able random generators that can produce
27 : * the n-th sample in a reproducible sequence without maintaining large
28 : * generator states between calls. It achieves this by keeping a small cached
29 : * generator state for efficient sequential access, while still allowing
30 : * deterministic random access to any index `n` via recomputation.
31 : *
32 : * Usage Example:
33 : *
34 : * ```cpp
35 : * MooseRandomStateless rng(12345);
36 : *
37 : * Real x = rng.rand(5); // 6th random number
38 : * unsigned int y = rng.randl(10); // 11th integer in full 32-bit range
39 : * unsigned int z = rng.randl(7, 0, 100); // 8th integer in [0, 100)
40 : *
41 : * rng.advance(100); // Advance cache by 100 for speedier access
42 : * Real x = rng.rand(0, 17); // The 117th number in the sequence
43 : * ```
44 : *
45 : * Thread Safety:
46 : * - Concurrent access is not thread-safe, since internal generator state is
47 : * cached and mutated during calls to `evaluate()`. If you need concurrency,
48 : * create separate objects per thread.
49 : *
50 : * Performance:
51 : * - Sequential access (e.g., `n = k, k+1, k+2, ...`) is O(1) per call.
52 : * - Small forward jumps reuse cached state efficiently.
53 : * - Backward jumps or large random jumps trigger full reseeding and
54 : * re-evaluation up to `n`, which is *O(n)* in the distance skipped.
55 : *
56 : * Design Overview:
57 : * - `MooseRandomStatless` owns several `Generator<T>` instances, each
58 : * parameterized by a distribution:
59 : * - `mts_ldrand`: Uniform real [0, 1)
60 : * - `mts_lrand`: Uniform unsigned int
61 : * - `rds_iuniform`: Uniform integer [lower, upper)
62 : * - `Generator<T>` encapsulates:
63 : * - An initial `mt_state` seeded by `mts_seed32new()`
64 : * - A cached current state and current index for sequential efficiency
65 : * - A stateless, deterministic evaluation method `evaluate(n)`
66 : * - The random distributions use the same algorithm and seeding semantics
67 : * as MOOSE's existing `MooseRandom` class.
68 : *
69 : * Implementation Notes:
70 : * - `evaluate()` is declared `const` but mutates internal state marked as
71 : * `mutable` to allow transparent caching; this is a deliberate trade-off
72 : * for performance.
73 : * - Each generator instance maintains its own `mt_state`; there is no global
74 : * RNG state shared across indices.
75 : * - There are separate generators for every instance of `rds_iuniform` with
76 : * unique ranges (upper - lower). This is because `rds_iuniform` increments
77 : * its state differently for each of these ranges. Creating separate
78 : * generators for each input range is unfortunate, but necessary for
79 : * stateless-ness.
80 : */
81 : class MooseRandomStateless
82 : {
83 : public:
84 : /// Deleted copy constructor and assignment to prevent deep-copy of mutable state
85 : MooseRandomStateless(const MooseRandomStateless &) = delete;
86 : MooseRandomStateless & operator=(const MooseRandomStateless &) = delete;
87 :
88 : /**
89 : * @brief Construct a new MooseRandomStateless with a given seed.
90 : *
91 : * Initializes internal generators for real and integer distributions.
92 : *
93 : * @param seed Base seed value for all RNGs in this handler
94 : */
95 94 : MooseRandomStateless(unsigned int seed)
96 94 : : _seed(seed), _rand_generator(mts_ldrand, seed), _randl_generator(mts_lrand, seed)
97 : {
98 94 : }
99 :
100 : /**
101 : * @brief Get the seed value used to initialize this handler's RNGs
102 : *
103 : * @return unsigned int
104 : */
105 6 : unsigned int getSeed() const { return _seed; }
106 :
107 : /**
108 : * @brief Get the amount that the RNGs have advanced by advance(count) calls.
109 : *
110 : * @return std::size_t
111 : */
112 6 : std::size_t getAdvanceCount() const { return _advance_count; }
113 :
114 : /**
115 : * @brief Return the n-th uniform Real random number in [0, 1).
116 : *
117 : * Uses the cached state for sequential access; otherwise, recomputes from seed.
118 : *
119 : * @param n 0-based index of the random number to generate
120 : * @return Random Real number in [0, 1)
121 : */
122 4088 : Real rand(std::size_t n) const { return _rand_generator.evaluate(n); }
123 :
124 : /**
125 : * @brief Return the n-th unsigned integer from the full 32-bit uniform range.
126 : * This is useful for generating seeds for new generators.
127 : *
128 : * @param n 0-based index of the random number to generate
129 : * @return Unsigned integer sampled uniformly in [0, 2^32)
130 : */
131 4030 : unsigned int randl(std::size_t n) const { return _randl_generator.evaluate(n); }
132 :
133 : /**
134 : * @brief Return the n-th bounded integer random number in [lower, upper).
135 : *
136 : * Uses a cached generator per unique range.
137 : *
138 : * @param n 0-based index of the random number to generate
139 : * @param lower Lower bound (inclusive)
140 : * @param upper Upper bound (exclusive)
141 : * @return Unsigned integer sampled uniformly in [lower, upper)
142 : */
143 4282 : unsigned int randl(std::size_t n, unsigned int lower, unsigned int upper) const
144 : {
145 : mooseAssert(upper >= lower, "randl: upper < lower");
146 4282 : const auto range = upper - lower;
147 :
148 4282 : auto [it, is_new] = _randlb_generators.try_emplace(
149 200872 : range, [&range](mt_state * state) { return rds_iuniform(state, 0, range); }, _seed);
150 :
151 4282 : if (is_new)
152 1424 : it->second.advance(_advance_count);
153 :
154 8564 : return lower + it->second.evaluate(n);
155 : }
156 :
157 : /**
158 : * @brief Advance all generators by the specified number.
159 : *
160 : * @param count The number of calls to advance.
161 : */
162 12 : void advance(std::size_t count)
163 : {
164 12 : _rand_generator.advance(count);
165 12 : _randl_generator.advance(count);
166 12 : for (auto & [range, gen] : _randlb_generators)
167 0 : gen.advance(count);
168 12 : _advance_count += count;
169 12 : }
170 :
171 : /**
172 : * @brief Template class wrapping a single random number generator function.
173 : *
174 : * @tparam T Output type of the RNG (e.g., `Real` or `unsigned int`)
175 : *
176 : * Each `Generator` stores:
177 : * - The RNG function pointer (from randistrs)
178 : * - An initial state (seeded once)
179 : * - A cached current state and index for fast sequential sampling
180 : */
181 : template <typename T>
182 : class Generator
183 : {
184 : public:
185 : /**
186 : * @brief Construct a new Generator with a function and seed.
187 : *
188 : * @param rng_func RNG function to apply to `mt_state` (e.g., `mts_ldrand`)
189 : * @param seed Seed used to initialize the base state
190 : */
191 1620 : Generator(std::function<T(mt_state *)> rng_func, unsigned int seed)
192 1620 : : _rng_func(std::move(rng_func))
193 : {
194 : // Create a new state with the given seed
195 1620 : mts_seed32new(&_initial_state, seed);
196 1620 : _current_state = _initial_state;
197 1620 : _current_index = 0;
198 1620 : }
199 :
200 : /**
201 : * @brief Evaluate the n-th random number in the sequence.
202 : *
203 : * Uses cached state when `n >= _current_index` for sequential access.
204 : * If `n < _current_index`, it resets to the initial seed and recomputes.
205 : *
206 : * @param n 0-based index in the random sequence
207 : * @return T The n-th random value
208 : */
209 13400 : T evaluate(std::size_t n) const
210 : {
211 13400 : if (n < _current_index)
212 : {
213 5774 : _current_state = _initial_state;
214 5774 : _current_index = 0;
215 : }
216 :
217 13400 : T val = T();
218 1089262 : for (; _current_index <= n; ++_current_index)
219 1075862 : val = _rng_func(&_current_state);
220 :
221 13400 : return val;
222 : }
223 :
224 : /**
225 : * @brief Advance the internal RNG state by a fixed count.
226 : *
227 : * This resets `_current_state` to `_initial_state` after advancing.
228 : *
229 : * @param count Number of RNG calls to skip
230 : */
231 1450 : void advance(std::size_t count)
232 : {
233 1450 : advance(_initial_state, count);
234 1450 : _current_state = _initial_state;
235 1450 : _current_index = 0;
236 1450 : }
237 :
238 : private:
239 : /// Advance a given RNG state `count` times.
240 1450 : void advance(mt_state & state, std::size_t count)
241 : {
242 34938 : for (std::size_t i = 0; i < count; ++i)
243 33488 : _rng_func(&state);
244 1450 : }
245 :
246 : /// RNG function pointer
247 : const std::function<T(mt_state *)> _rng_func;
248 : /// Base RNG state (seeded)
249 : mt_state _initial_state;
250 : /// Cached working RNG state
251 : mutable mt_state _current_state;
252 : /// Index of next number to be generated
253 : mutable std::size_t _current_index;
254 : };
255 :
256 : private:
257 : /// Seed shared by all internal generators
258 : const unsigned int _seed;
259 : /// Uniform Real [0, 1)
260 : Generator<Real> _rand_generator;
261 : /// Uniform uint32
262 : Generator<unsigned int> _randl_generator;
263 : /// Bounded uniform uint32 (indexed based on bounding range)
264 : mutable std::unordered_map<unsigned int, Generator<unsigned int>> _randlb_generators;
265 :
266 : /// The number of counts the generators have advanced
267 : /// This needs to be kept around for generation of new randlb generators
268 : std::size_t _advance_count = 0;
269 : };
270 :
271 : template <>
272 : inline void
273 6 : dataStore(std::ostream & stream, MooseRandomStateless & v, void * context)
274 : {
275 6 : std::size_t count = v.getAdvanceCount();
276 6 : storeHelper(stream, count, context);
277 6 : }
278 :
279 : template <>
280 : inline void
281 6 : dataLoad(std::istream & stream, MooseRandomStateless & v, void * context)
282 : {
283 : std::size_t count;
284 6 : loadHelper(stream, count, context);
285 6 : v.advance(count);
286 6 : }
287 :
288 : template <>
289 : inline void
290 6 : dataStore(std::ostream & stream, std::unique_ptr<MooseRandomStateless> & v, void * context)
291 : {
292 6 : unsigned int seed = v->getSeed();
293 6 : storeHelper(stream, seed, context);
294 6 : storeHelper(stream, *v, context);
295 6 : }
296 :
297 : template <>
298 : inline void
299 6 : dataLoad(std::istream & stream, std::unique_ptr<MooseRandomStateless> & v, void * context)
300 : {
301 : unsigned int seed;
302 6 : loadHelper(stream, seed, context);
303 6 : v = std::make_unique<MooseRandomStateless>(seed);
304 6 : loadHelper(stream, *v, context);
305 6 : }
|
