-
Notifications
You must be signed in to change notification settings - Fork 41
/
RandomGen.cpp
182 lines (159 loc) · 5.42 KB
/
RandomGen.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
#include "RandomGen.hh"
#include <stdexcept>
using namespace std;
// Global static pointer used to ensure a single instance of the class.
RandomGen* RandomGen::m_pInstance = nullptr;
// Only allow one instance of class to be generated.
RandomGen* RandomGen::rndm() {
if (!m_pInstance) m_pInstance = new RandomGen;
return m_pInstance;
}
std::uint64_t splitmix64(std::uint64_t z) {
z += 0x9e3779b97f4a7c15;
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
return z ^ (z >> 31);
}
void RandomGen::SetSeed(uint64_t s) {
if (rng_locked) {
throw std::runtime_error("You can not change the seed because it is locked.");
}
uint64_t s1 = splitmix64(s);
rng = xoroshiro128plus64(s1, splitmix64(s1));
}
void RandomGen::LockSeed() {
rng_locked = true;
}
void RandomGen::UnlockSeed() {
rng_locked = false;
}
double RandomGen::rand_uniform() {
return (static_cast<double>(rng()) - xoroshiro128plus64_min) /
xoroshiro128plus64_minmax;
}
double RandomGen::rand_gauss(double mean, double sigma, bool zero_min) {
// std::normal_distribution<double> norm(mean, sigma);
// return norm(rng);
double u = rand_uniform(), v = rand_uniform();
double draw = mean + sigma * sqrt2 * sqrt(-log(u)) * cos(two_PI * v);
if (zero_min) {
return max(draw, 0.);
}
return draw;
}
double RandomGen::rand_zero_trunc_gauss(double mean, double sigma) {
double r = rand_gauss(mean, sigma, false);
while (r <= 0) {
r = rand_gauss(mean, sigma, false);
}
return r;
}
double RandomGen::FindNewMean(double sigma) {
// Follow https://en.wikipedia.org/wiki/Truncated_normal_distribution
double TruncGaussAlpha = -1. / sigma;
double LittlePhi_Alpha =
exp(-0.5 * TruncGaussAlpha * TruncGaussAlpha) / gcem::sqrt(2. * M_PI);
double BigPhi_Alpha = 0.5 * (1. + erf(TruncGaussAlpha / sqrt2));
return 1. + (LittlePhi_Alpha / (1. - BigPhi_Alpha)) * sigma;
}
double RandomGen::rand_exponential(double half_life, double t_min,
double t_max) {
double r = rand_uniform();
double r_min = 0;
double r_max = 1;
if (t_min >= 0) {
r_min = 1 - exp(-t_min * log2 / half_life);
}
if (t_max >= 0) {
r_max = 1 - exp(-t_max * log2 / half_life);
}
r = (r_max - r_min) * r + r_min;
return -log(1 - r) * half_life / log2;
}
double RandomGen::rand_skewGauss(double xi, double omega, double alpha) {
double delta = alpha / sqrt(1 + alpha * alpha);
double gamma1 = four_minus_PI_div_2 *
(pow(delta * sqrt2_div_PI, 3.) /
pow(1 - 2. * delta * delta / M_PI, 1.5)); // skewness
double muz = delta * sqrt2_div_PI;
double sigz = sqrt(1. - muz * muz);
double m_o;
if (alpha > 0.) {
m_o = muz - 0.5 * gamma1 * sigz - 0.5 * exp(-two_PI / alpha);
}
if (alpha < 0.) {
m_o = muz - 0.5 * gamma1 * sigz + 0.5 * exp(two_PI / alpha);
}
double mode = xi + omega * m_o;
// the height should be the value of the PDF at the mode
double height = exp(-0.5 * (pow((mode - xi) / omega, 2.))) /
(sqrt2_PI * omega) *
erfc(-1. * alpha * (mode - xi) / omega / sqrt2);
bool gotValue = false;
double minX = xi - 6. * omega;
double maxX =
xi + 6. * omega; // +/- 6sigma should be essentially +/- infinity
// can increase these for even better accuracy, at the
// cost of speed
double testX, testY, testProb;
while (gotValue == false) {
testX = minX + (maxX - minX) * RandomGen::rndm()->rand_uniform();
testY = height *
RandomGen::rndm()->rand_uniform(); // between 0 and peak height
// calculate the value of the skewGauss PDF at the test x-value
testProb = exp(-0.5 * (pow((testX - xi) / omega, 2.))) /
(sqrt2_PI * omega) *
erfc(-1. * alpha * (testX - xi) / omega / sqrt2);
if (testProb >= testY) {
gotValue = true;
}
}
return testX;
}
int RandomGen::poisson_draw(double mean) {
return std::poisson_distribution<int>(mean)(rng);
}
int64_t RandomGen::binom_draw(int64_t N0, double prob) {
return std::binomial_distribution<int64_t>(N0, prob)(rng);
}
int RandomGen::integer_range(int min, int max) {
return std::uniform_int_distribution<int>(min, max)(rng);
}
vector<double> RandomGen::VonNeumann(double xMin, double xMax, double yMin,
double yMax, double xTest, double yTest,
double fValue) {
vector<double> xyTry(3);
xyTry[0] = xTest;
xyTry[1] = yTest;
if (xyTry[1] > fValue) {
xyTry[0] = xMin + (xMax - xMin) * RandomGen::rndm()->rand_uniform();
xyTry[1] = yMin + (yMax - yMin) * RandomGen::rndm()->rand_uniform();
xyTry[2] = 1.;
} else
xyTry[2] = 0.;
return xyTry; // doing a vector means you can return 2 values at the same
// time
}
int RandomGen::SelectRanXeAtom() { // to determine the isotope of Xe
int A;
double isotope = rand_uniform() * 100.;
if (isotope > 0.000 && isotope <= 0.090)
A = 124;
else if (isotope > 0.090 && isotope <= 0.180)
A = 126;
else if (isotope > 0.180 && isotope <= 2.100)
A = 128;
else if (isotope > 2.100 && isotope <= 28.54)
A = 129;
else if (isotope > 28.54 && isotope <= 32.62)
A = 130;
else if (isotope > 32.62 && isotope <= 53.80)
A = 131;
else if (isotope > 53.80 && isotope <= 80.69)
A = 132;
else if (isotope > 80.69 && isotope <= 91.13)
A = 134;
else
A = 136;
return A;
}