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
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Sample/Correlation/IDistribution1DSampler.cpp
//! @brief Defines class interface IProfile1D, and children thereof.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Sample/Correlation/IDistribution1DSampler.h"
#include <numbers>
#include <random>
using std::numbers::pi;
IDistribution1DSampler::~IDistribution1DSampler() = default;
double Distribution1DCauchySampler::randomSample(int seed) const
{
// BornAgain Cauchy Distribution = std library Exponential distribution
std::random_device rd; // random device class instance
std::mt19937 gen(seed < 0 ? rd() : seed); // Standard mersenne_twister_engine seeded with rd()
std::exponential_distribution<double> expDist(m_lambda);
double value = expDist(gen);
std::bernoulli_distribution bernoulliDist(0.5);
bool sign = bernoulliDist(gen);
return sign ? value : -value;
}
double Distribution1DGaussSampler::randomSample(int seed) const
{
// BornAgain Gauss Distribution = std library Normal distribution
std::random_device rd;
std::mt19937 gen(seed < 0 ? rd() : seed);
std::normal_distribution<double> normalDist(m_mean, m_stddev);
return normalDist(gen);
}
double Distribution1DGateSampler::randomSample(int seed) const
{
// BornAgain Gate Distribution = std library Uniform distribution
std::random_device rd;
std::mt19937 gen(seed < 0 ? rd() : seed);
std::uniform_real_distribution<double> uniformDist(m_a, m_b);
return uniformDist(gen);
}
double Distribution1DTriangleSampler::randomSample(int seed) const
{
std::random_device rd;
std::mt19937 gen(seed < 0 ? rd() : seed);
// generate a cdf value between 0 and 1
std::uniform_real_distribution<> uniformDist(0.0, 1.0);
double cdf_value = uniformDist(gen);
// solve for x by inverting the cdf of Triangle Distribution
if (cdf_value <= 0.5)
return (-m_omega + m_omega * std::sqrt(2 * cdf_value));
return (m_omega - m_omega * std::sqrt(2 * (1 - cdf_value)));
}
double Distribution1DCosineSampler::randomSample(int seed) const
{
std::random_device rd;
std::mt19937 gen(seed < 0 ? rd() : seed);
// generate a cdf value between 0 and 1
std::uniform_real_distribution<> uniformDist(0.0, 1.0);
double cdf_value = uniformDist(gen);
// solve for x from the cdf of Cosine Distribution using Newton-Raphson method
double func = 0.0, funcDeriv = 0.0, x = 0.0;
// initial guess for x
if (cdf_value <= 0.5)
x = -m_omega / 2;
else
x = m_omega / 2;
bool convergedSoln = false;
while (!convergedSoln) {
func = x + m_omega / pi * std::sin(pi * x / m_omega) + m_omega * (1 - 2 * cdf_value);
funcDeriv = 1 + std::cos(pi * x / m_omega);
x = x - func / funcDeriv;
if (std::abs(func / funcDeriv) < 0.001)
convergedSoln = true;
}
return x;
}
|
