1 #pragma once
2
3 /// \file Rng.h
4 /// \brief Random number generators
5 /// \author Pavel Sevecek (sevecek at sirrah.troja.mff.cuni.cz)
6 /// \date 2016-2021
7
8 #include "objects/containers/Array.h"
9 #include "objects/geometry/Box.h"
10 #include "objects/wrappers/AutoPtr.h"
11 #include <random>
12
13 NAMESPACE_SPH_BEGIN
14
15 /// \brief Random number generator with uniform distribution.
16 class UniformRng : public Noncopyable {
17 private:
18 std::mt19937_64 mt;
19
20 /// \note While std::mt19937_64 will always generate the same sequence of numbers (given the same seed, of
21 /// cource), the std::uniform_real_distribution might produce different result for different compilers, or
22 /// even for different versions of the compiler.
23 std::uniform_real_distribution<Float> dist;
24
25 public:
26 explicit UniformRng(const int seed = 1234) {
27 mt = std::mt19937_64(seed);
28 }
29
UniformRng(UniformRng && other)30 UniformRng(UniformRng&& other)
31 : mt(other.mt)
32 , dist(other.dist) {}
33
operator()34 Float operator()(const int UNUSED(s) = 0) {
35 return dist(mt);
36 }
37 };
38
39
40 /// \brief Random number generator used in code SPH5 of Benz & Asphaug (1994).
41 ///
42 /// Reimplemented for reproducibility of results.
43 class BenzAsphaugRng : public Noncopyable {
44 private:
45 const int seed;
46 const int im1 = 2147483563;
47 const int im2 = 2147483399;
48 const Float am = 1._f / im1;
49 const int imm1 = im1 - 1;
50 const int ia1 = 40014;
51 const int ia2 = 40692;
52 const int iq1 = 53668;
53 const int iq2 = 52774;
54 const int ir1 = 12211;
55 const int ir2 = 3791;
56 const Float eps = 1.2e-7_f;
57 const Float rnmx = 1._f - eps;
58
59 static constexpr int ntab = 32;
60 StaticArray<int, ntab> iv;
61 int iy = 0;
62 int idum2 = 123456789;
63 int idum;
64
65 public:
66 explicit BenzAsphaugRng(const int seed);
67
68 BenzAsphaugRng(BenzAsphaugRng&& other);
69
70 Float operator()(const int s = 0);
71 };
72
73 /// \brief Quasi-random number generator.
74 class HaltonQrng : public Noncopyable {
75 private:
76 static const int dimension = 6;
77 StaticArray<int, dimension> primes; /// \todo extend in case more is necessary
78 StaticArray<int, dimension> c;
79
80 public:
81 HaltonQrng();
82
83 HaltonQrng(HaltonQrng&& other);
84
85 Float operator()(const int s);
86 };
87
88
89 /// \brief Polymorphic holder allowing to store any RNG (type erasure).
90 class IRng : public Polymorphic {
91 public:
92 /// Generates a random number.
93 virtual Float operator()(const int s = 0) = 0;
94 };
95
96 template <typename TRng>
97 class RngWrapper : public IRng {
98 private:
99 TRng rng;
100
101 public:
102 template <typename... TArgs>
RngWrapper(TArgs &&...args)103 explicit RngWrapper(TArgs&&... args)
104 : rng(std::forward<TArgs>(args)...) {}
105
operator()106 virtual Float operator()(const int s) override {
107 return rng(s);
108 }
109 };
110
111 template <typename TRng, typename... TArgs>
makeRng(TArgs &&...args)112 AutoPtr<RngWrapper<TRng>> makeRng(TArgs&&... args) {
113 return makeAuto<RngWrapper<TRng>>(std::forward<TArgs>(args)...);
114 }
115
116 /// \brief Generates a random number from normal distribution, using Box-Muller algorithm.
117 ///
118 /// Could be optimized (it discards the second independent value), but it's the algorith easiest to implement.
119 template <typename TRng>
sampleNormalDistribution(TRng & rng,const Float mu,const Float sigma)120 INLINE Float sampleNormalDistribution(TRng& rng, const Float mu, const Float sigma) {
121 static const Float epsilon = std::numeric_limits<Float>::min();
122
123 Float u1, u2;
124 do {
125 u1 = rng();
126 u2 = rng();
127 } while (u1 <= epsilon);
128
129 const Float z1 = sqrt(-2._f * log(u1)) * cos(2._f * PI * u2);
130 SPH_ASSERT(isReal(z1));
131 return z1 * sigma + mu;
132 }
133
134 /// \brief Generates a random number from exponential distribution.
135 ///
136 /// Uses inverse transform sampling.
137 template <typename TRng>
sampleExponentialDistribution(TRng & rng,const Float lambda)138 INLINE Float sampleExponentialDistribution(TRng& rng, const Float lambda) {
139 static const Float epsilon = std::numeric_limits<Float>::min();
140
141 Float u;
142 do {
143 u = rng();
144 } while (u <= epsilon);
145 return -log(rng()) / lambda;
146 }
147
148 /// \brief Generates a random integer from Poisson distribution.
149 ///
150 /// Uses the Knuth's algorithm.
151 template <typename TRng>
samplePoissonDistribution(TRng & rng,const Float lambda)152 INLINE Size samplePoissonDistribution(TRng& rng, const Float lambda) {
153 const Float l = exp(-lambda);
154 Size k = 0;
155 Float p = 1;
156 do {
157 k = k + 1;
158 const Float u = rng();
159 p = p * u;
160 } while (p > l);
161 return k - 1;
162 }
163
164 /// \brief Generates a random position on a unit sphere
165 template <typename TRng>
sampleUnitSphere(TRng & rng)166 INLINE Vector sampleUnitSphere(TRng& rng) {
167 const Float phi = rng() * 2._f * PI;
168 const Float z = rng() * 2._f - 1._f;
169 const Float u = sqrt(1._f - sqr(z));
170
171 return Vector(u * cos(phi), u * sin(phi), z);
172 }
173
174
175 /// \brief Generates a random number from a generic distribution, using rejection sampling.
176 ///
177 /// Note that this function may be very inefficient and should be used only if the distribution cannot be
178 /// sampled with an explicit method.
179 /// \param rng Random number generator
180 /// \param range Minimal and maximal generated value
181 /// \param upperBound Upper bound for the values returned by the functor
182 /// \param func Probability distribution function. Does not have to be normalized.
183 template <typename TRng, typename TFunc>
sampleDistribution(TRng & rng,const Interval & range,const Float upperBound,const TFunc & func)184 INLINE Float sampleDistribution(TRng& rng, const Interval& range, const Float upperBound, const TFunc& func) {
185 while (true) {
186 const Float x = range.lower() + rng() * range.size();
187 const Float y = rng() * upperBound;
188 const Float pdf = func(x);
189 SPH_ASSERT(pdf >= 0._f && pdf < upperBound, pdf);
190 if (y < pdf) {
191 return x;
192 }
193 }
194 }
195
196 /// \brief Generates a random vector from a generic distribution, using rejection sampling.
197 ///
198 /// \param rng Random number generator
199 /// \param box Extents of generated vectors
200 /// \param upperBound Upper bound for the values returned by the functor
201 /// \param func Probability distribution function. Does not have to be normalized.
202 template <typename TRng, typename TFunc>
sampleDistribution(TRng & rng,const Box & box,const Float upperBound,const TFunc & func)203 INLINE Vector sampleDistribution(TRng& rng, const Box& box, const Float upperBound, const TFunc& func) {
204 while (true) {
205 const Vector r = box.lower() + Vector(rng(), rng(), rng()) * box.size();
206 const Float y = rng() * upperBound;
207 const Float pdf = func(r);
208 SPH_ASSERT(pdf >= 0._f && pdf < upperBound, pdf);
209 if (y < pdf) {
210 return r;
211 }
212 }
213 }
214
215 NAMESPACE_SPH_END
216