version 3.10-dev
random.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3//
4// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
12#ifndef DUMUX_COMMON_RANDOM_HH
13#define DUMUX_COMMON_RANDOM_HH
14
15#include <random>
16#include <type_traits>
17#include <cstdint>
18
19namespace Dumux {
20
28template<class Scalar = double>
30{
31 struct Parameters
32 {
33 Parameters(Scalar a, Scalar b)
34 : a_(a), b_(b) {}
35
36 Scalar a() const { return a_; }
37 Scalar b() const { return b_; }
38 private:
39 Scalar a_, b_;
40 };
41public:
42 using param_type = Parameters;
43 using result_type = Scalar;
44
45 explicit SimpleUniformDistribution(Scalar min, Scalar max = 1.0)
46 : offset_(min)
47 , size_(max-min)
48 {}
49
50 explicit SimpleUniformDistribution(const Parameters& p)
51 : SimpleUniformDistribution(p.a(), p.b())
52 {}
53
56 {}
57
58 void param(const Parameters& p)
59 {
60 offset_ = p.a();
61 size_ = p.b()-p.a();
62 }
63
64 Parameters param() const
65 { return { offset_, offset_+size_ }; }
66
67 Scalar a() const { return offset_; }
68 Scalar b() const { return offset_ + size_; }
69
70 template<class Generator,
71 typename std::enable_if_t<std::is_same_v<typename Generator::result_type, std::uint_fast32_t>, int> = 0>
72 Scalar operator()(Generator& gen)
73 { return offset_ + size_*(0x1.0p-32 * gen()); }
74
75private:
76 Scalar offset_;
77 Scalar size_;
78};
79
87template<class Scalar = double>
89{
90 struct Parameters
91 {
92 Parameters(Scalar m, Scalar s)
93 : m_(m), s_(s) {}
94
95 Scalar m() const { return m_; }
96 Scalar s() const { return s_; }
97 private:
98 Scalar m_, s_;
99 };
100public:
101 using param_type = Parameters;
102 using result_type = Scalar;
103
104 explicit SimpleNormalDistribution(Scalar mean, Scalar stddev = 1.0)
105 : mean_(mean)
106 , stddev_(stddev)
107 {}
108
109 explicit SimpleNormalDistribution(const Parameters& p)
110 : SimpleNormalDistribution(p.m(), p.s())
111 {}
112
115 {}
116
117 void param(const Parameters& p)
118 {
119 mean_ = p.m();
120 stddev_ = p.s();
121 }
122
123 Parameters param() const
124 { return { mean_, stddev_ }; }
125
126 Scalar m() const { return mean_; }
127 Scalar s() const { return stddev_; }
128
129 template<class Generator>
130 Scalar operator()(Generator& gen)
131 {
132 if (isCached_)
133 {
134 isCached_ = false;
135 return cachedValue_;
136 }
137
138 // Box-Mueller transform (https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform)
139 const auto [u1, u2] = generateUniformPair_(gen);
140
141 using std::sqrt; using std::log;
142 using std::cos; using std::sin;
143 constexpr Scalar twoPi = 2.0 * M_PI;
144
145 const Scalar magnitude = stddev_ * sqrt(-2.0 * log(u1));
146 const Scalar z0 = magnitude * cos(twoPi * u2) + mean_;
147 const Scalar z1 = magnitude * sin(twoPi * u2) + mean_;
148 cachedValue_ = z0;
149 isCached_ = true;
150 return z1;
151 }
152
153private:
154 template<class Generator,
155 typename std::enable_if_t<std::is_same_v<typename Generator::result_type, std::uint_fast32_t>, int> = 0>
156 auto generateUniformPair_(Generator& gen)
157 {
158 // biased uniform number generator (0,1)
159 // https://www.pcg-random.org/posts/bounded-rands.html
160 constexpr Scalar eps = std::numeric_limits<Scalar>::epsilon();
161 Scalar u1 = 0.0, u2 = 0.0;
162 do {
163 u1 = 0x1.0p-32 * gen();
164 u2 = 0x1.0p-32 * gen();
165 } while (u1 <= eps);
166 return std::make_pair(u1, u2);
167 }
168
169 Scalar mean_;
170 Scalar stddev_;
171 bool isCached_ = false;
172 Scalar cachedValue_ = {};
173};
174
181template<class Scalar = double>
183{
184 using Parameters = typename SimpleNormalDistribution<Scalar>::param_type;
185public:
186 using param_type = Parameters;
187 using result_type = Scalar;
188
189 explicit SimpleLogNormalDistribution(Scalar mean, Scalar stddev = 1.0)
190 : normal_(mean, stddev)
191 {}
192
193 explicit SimpleLogNormalDistribution(const Parameters& p)
194 : SimpleLogNormalDistribution(p.mean, p.stddev)
195 {}
196
199 {}
200
201 void param(const Parameters& p)
202 { normal_.param(p); }
203
204 Parameters param() const
205 { return normal_.param(); }
206
207 Scalar m() const { return normal_.m(); }
208 Scalar s() const { return normal_.s(); }
209
210 template<class Generator>
211 Scalar operator()(Generator& gen)
212 {
213 using std::exp;
214 return exp(normal_(gen));
215 }
216
217private:
219};
220
221} // end namespace Dumux
222
223#endif
A simple log-normal distribution.
Definition: random.hh:183
Scalar operator()(Generator &gen)
Definition: random.hh:211
Scalar m() const
Definition: random.hh:207
SimpleLogNormalDistribution()
Definition: random.hh:197
Parameters param_type
Definition: random.hh:186
Parameters param() const
Definition: random.hh:204
SimpleLogNormalDistribution(const Parameters &p)
Definition: random.hh:193
SimpleLogNormalDistribution(Scalar mean, Scalar stddev=1.0)
Definition: random.hh:189
void param(const Parameters &p)
Definition: random.hh:201
Scalar s() const
Definition: random.hh:208
Scalar result_type
Definition: random.hh:187
A simple normal distribution based on a biased uniform number generator and the Box-Mueller transform...
Definition: random.hh:89
Scalar operator()(Generator &gen)
Definition: random.hh:130
Parameters param_type
Definition: random.hh:101
void param(const Parameters &p)
Definition: random.hh:117
SimpleNormalDistribution(const Parameters &p)
Definition: random.hh:109
SimpleNormalDistribution(Scalar mean, Scalar stddev=1.0)
Definition: random.hh:104
Scalar result_type
Definition: random.hh:102
Parameters param() const
Definition: random.hh:123
Scalar m() const
Definition: random.hh:126
SimpleNormalDistribution()
Definition: random.hh:113
Scalar s() const
Definition: random.hh:127
A simple uniform distribution based on a biased uniform number generator.
Definition: random.hh:30
Scalar b() const
Definition: random.hh:68
Scalar a() const
Definition: random.hh:67
SimpleUniformDistribution()
Definition: random.hh:54
void param(const Parameters &p)
Definition: random.hh:58
SimpleUniformDistribution(Scalar min, Scalar max=1.0)
Definition: random.hh:45
SimpleUniformDistribution(const Parameters &p)
Definition: random.hh:50
Parameters param() const
Definition: random.hh:64
Scalar operator()(Generator &gen)
Definition: random.hh:72
Parameters param_type
Definition: random.hh:42
Scalar result_type
Definition: random.hh:43
Definition: adapt.hh:17