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
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Base/Math/Functions.cpp
//! @brief Implements functions in namespace Math.
//!
//! @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 "Base/Math/Functions.h"
#include "Base/Util/Assert.h"
#include <chrono>
#include <cstring>
#include <fftw3.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_expint.h>
#include <limits>
#include <numbers>
#include <random>
using std::numbers::pi;
// ************************************************************************************************
// Various functions
// ************************************************************************************************
double Math::StandardNormal(double x)
{
return std::exp(-x * x / 2.0) / std::sqrt((2 * pi));
}
double Math::Gaussian(double x, double average, double std_dev)
{
return StandardNormal((x - average) / std_dev) / std_dev;
}
double Math::IntegratedGaussian(double x, double average, double std_dev)
{
if (std_dev <= std::numeric_limits<double>::min())
return Math::Heaviside(x - average);
double normalized_x = (x - average) / std_dev;
static double root2 = std::sqrt(2.0);
return (gsl_sf_erf(normalized_x / root2) + 1.0) / 2.0;
}
double Math::Heaviside(double x)
{
return x > 0.0 ? 1.0 : 0.0;
}
double Math::cot(double x)
{
return tan((pi / 2) - x);
}
double Math::sinc(double x) // Sin(x)/x
{
if (x == 0)
return 1;
return std::sin(x) / x;
}
complex_t Math::sinc(const complex_t z) // Sin(x)/x
{
// This is an exception from the rule that we must not test floating-point numbers for equality.
// For small non-zero arguments, sin(z) returns quite accurately z or z-z^3/6.
// There is no loss of precision in computing sin(z)/z.
// Therefore there is no need for an expensive test like abs(z)<eps.
if (z == complex_t(0., 0.))
return 1.0;
return std::sin(z) / z;
}
complex_t Math::tanhc(const complex_t z) // tanh(x)/x
{
if (std::abs(z) < std::numeric_limits<double>::epsilon())
return 1.0;
return std::tanh(z) / z;
}
double Math::Laue(const double x, size_t N)
{
static const double SQRT6DOUBLE_EPS = std::sqrt(6.0 * std::numeric_limits<double>::epsilon());
auto nd = static_cast<double>(N);
if (std::abs(nd * x) < SQRT6DOUBLE_EPS)
return nd;
double num = std::sin(nd * x);
double den = std::sin(x);
return num / den;
}
double Math::erf(double arg)
{
ASSERT(arg >= 0.0);
if (std::isinf(arg))
return 1.0;
return gsl_sf_erf(arg);
}
// ************************************************************************************************
// Random number generators
// ************************************************************************************************
/* currently unused
double Math::GenerateStandardNormalRandom() // using c++11 standard library
{
unsigned seed =
static_cast<unsigned>(std::chrono::system_clock::now().time_since_epoch().count());
std::mt19937 generator(seed);
std::normal_distribution<double> distribution(0.0, 1.0);
return distribution(generator);
}
double Math::GenerateNormalRandom(double average, double std_dev)
{
return GenerateStandardNormalRandom() * std_dev + average;
}
*/
double Math::GeneratePoissonRandom(double average) // using c++11 standard library
{
unsigned seed =
static_cast<unsigned>(std::chrono::system_clock::now().time_since_epoch().count());
std::mt19937 generator(seed);
if (average <= 0.0)
return 0.0;
if (average < 1000.0) { // Use std::poisson_distribution
std::poisson_distribution<int> distribution(average);
int sample = distribution(generator);
return (double)sample;
}
// Use normal approximation
std::normal_distribution<double> distribution(average, std::sqrt(average));
double sample = distribution(generator);
return std::max(0.0, sample);
}
int Math::GenerateNextSeed(unsigned seed)
{
std::mt19937 gen(seed);
std::uniform_int_distribution<int> distr(0, std::numeric_limits<int>::max());
return distr(gen);
}
|
