Chebyshev
Unit testing for scientific software
Loading...
Searching...
No Matches
estimator.h
Go to the documentation of this file.
1
5
6#ifndef CHEBYSHEV_ESTIMATOR
7#define CHEBYSHEV_ESTIMATOR
8
9
10#include <functional>
11#include <cmath>
12
13#include "../core/common.h"
14#include "../core/random.h"
15#include "./prec_structures.h"
16
17
18namespace chebyshev {
19namespace prec {
20
21
23 namespace estimator {
24
25
29 template<typename FloatType = double>
30 inline auto quadrature1D() {
31
32 return [](
36
37 if(options.domain.size() != 1) {
38 throw std::runtime_error(
39 "estimator::quadrature1D only works on mono-dimensional domains"
40 );
41 }
42
43 interval domain = options.domain[0];
44
45 FloatType sum = 0;
46 FloatType sumSqr = 0;
47 FloatType sumAbs = 0;
48 FloatType max = 0;
49
50 const FloatType length = domain.length();
51 const FloatType dx = length / options.iterations;
54
55 FloatType diff = std::abs(funcApprox(domain.a) - funcExpected(domain.a));
56
57 sum += diff;
58 sumSqr += diff * diff;
59 sumAbs += std::abs(funcExpected(domain.a));
60 max = diff;
61
62 for (unsigned int i = 1; i < options.iterations; ++i) {
63
64 x = domain.a + i * dx;
65 diff = std::abs(funcApprox(x) - funcExpected(x));
66
67 if(diff > max)
68 max = diff;
69
70 if(i % 2 == 0)
71 coeff = 2;
72 else
73 coeff = 4;
74
75 sum += coeff * diff;
76 sumSqr += coeff * diff * diff;
78 }
79
80 diff = std::abs(funcApprox(domain.b) - funcExpected(domain.b));
81
82 sum += diff;
83 sumSqr += diff * diff;
84 sumAbs += std::abs(funcExpected(domain.b));
85
86 if(diff > max)
87 max = diff;
88
90 res.absErr = sum;
91 res.maxErr = max;
92 res.meanErr = (sum * dx / 3.0) / length;
93 res.rmsErr = std::sqrt((sumSqr * dx / 3.0) / length);
94 res.relErr = std::abs((sum * dx / 3.0) / (sumAbs * dx / 3.0));
95
96 return res;
97 };
98 }
99
100
109 template <
110 typename IntType = int,
111 typename ReturnType = IntType
112 >
113 inline auto discrete1D() {
114
115 // Return a one-dimensional discrete estimator as a lambda function
116 return [](
117 std::function<IntType(ReturnType)> funcApprox,
118 std::function<IntType(ReturnType)> funcExpected,
120
121 if(options.domain.size() != 1) {
122 throw std::runtime_error(
123 "estimator::discrete1D only works on mono-dimensional domains"
124 );
125 }
126
127 IntType extreme1 = IntType(std::ceil(options.domain[0].a));
128 IntType extreme2 = IntType(std::floor(options.domain[0].b));
129
132
133 long double maxErr = 0;
134 long double sumDiff = 0;
135 long double sumSqr = 0;
136 long double sumAbs = 0;
138
139 for (IntType n = lower; n <= upper; ++n) {
140
143
144 const long double diff = (long double) resExpected > resApprox ?
146
147 maxErr = std::max(maxErr, diff);
148 sumDiff += diff;
149 sumSqr += diff * diff;
150 sumAbs += std::abs((long double) funcExpected(n));
151 totalPoints++;
152 }
153
155 res.absErr = sumAbs;
156 res.maxErr = maxErr;
157 res.meanErr = totalPoints > 0 ? (sumDiff / totalPoints) : 0;
158 res.rmsErr = totalPoints > 0 ? (std::sqrt(sumSqr) / totalPoints) : 0;
159 res.relErr = sumDiff / sumAbs;
160 return res;
161 };
162 }
163
164
168 template<typename FloatType = double>
169 inline auto montecarlo1D(std::shared_ptr<random::random_context> rnd_ctx) {
170
171 // Return a one-dimensional Monte Carlo estimator as a lambda function
172 return [rnd_ctx](
176
177 random::random_source rnd = rnd_ctx->get_rnd();
178
179 if(options.domain.size() != 1) {
180 throw std::runtime_error(
181 "estimator::montecarlo1D only works on mono-dimensional domains"
182 );
183 }
184
185 FloatType sum = 0;
186 FloatType sumSqr = 0;
187 FloatType sumAbs = 0;
188 FloatType max = 0;
189 const FloatType length = options.domain[0].length();
190
191 for (int i = 0; i < options.iterations; ++i) {
192
193 FloatType x = rnd.uniform(options.domain[0].a, options.domain[0].b);
194 const FloatType diff = std::abs(funcApprox(x) - funcExpected(x));
195
196 max = std::max(max, diff);
197 sum += diff;
198 sumSqr += diff * diff;
199 sumAbs += std::abs(funcExpected(x));
200 }
201
203 res.maxErr = max;
204 res.meanErr = sum / options.iterations;
205 res.absErr = sum * (length / options.iterations);
206 res.rmsErr = sumSqr * (length / options.iterations);
207 res.relErr = sum / sumAbs;
208
209 return res;
210 };
211 }
212
213
226 template <
227 typename FloatType = double,
228 typename Vector = std::vector<FloatType>
229 >
230 inline auto montecarlo(
231 std::shared_ptr<random::random_context> rnd_ctx, unsigned int dimensions) {
232
233 // Return an n-dimensional Monte Carlo estimator as a lambda function
234 return [rnd_ctx, dimensions](
235 std::function<FloatType(Vector)> funcApprox,
236 std::function<FloatType(Vector)> fExpected,
238
239 random::random_source rnd = rnd_ctx->get_rnd();
240
241 if(options.domain.size() != dimensions) {
242 throw std::runtime_error(
243 "The estimation domain's dimension does not match "
244 "the instantiated number of dimensions "
245 "in estimator::montecarlo"
246 );
247 }
248
249 FloatType sum = 0;
250 FloatType sumSqr = 0;
251 FloatType sumAbs = 0;
252 FloatType max = 0;
253
254 // Compute the measure of a multi-interval
255 FloatType volume = 1;
256 for (interval k : options.domain)
257 volume *= k.length();
258
260
261 for (int i = 0; i < options.iterations; ++i) {
262
263 rnd.uniform(x, options.domain);
264
265 const FloatType diff = std::abs(funcApprox(x) - funcExpected(x));
266
267 if(max < diff)
268 max = diff;
269
270 sum += diff;
271 sumSqr += diff * diff;
272 sumAbs += std::abs(funcExpected(x));
273 }
274
276 res.maxErr = max;
277 res.meanErr = sum / options.iterations;
278 res.absErr = sum * (volume / options.iterations);
279 res.rmsErr = sumSqr * (volume / options.iterations);
280 res.relErr = sum / sumAbs;
281
282 return res;
283 };
284 }
285 }
286}}
287
288
289#endif
auto quadrature1D()
Use Simpson's quadrature scheme to approximate error integrals for univariate real functions (endofun...
Definition estimator.h:30
auto montecarlo(std::shared_ptr< random::random_context > rnd_ctx, unsigned int dimensions)
Use crude Monte Carlo integration to approximate error integrals for multivariate real functions.
Definition estimator.h:230
auto montecarlo1D(std::shared_ptr< random::random_context > rnd_ctx)
Use crude Monte Carlo integration to approximate error integrals for univariate real functions.
Definition estimator.h:169
auto discrete1D()
Use a discrete estimator over a lattice of points, here implemented in one dimension,...
Definition estimator.h:113
General namespace of the framework.
Definition benchmark.h:22
constexpr FloatType get_nan()
Get a quiet NaN of the specified floating point type.
Definition common.h:65
std::function< Type(Type)> EndoFunction
An endofunction is a function which has the same type of input and output, such as a real function of...
Definition common.h:60
Structures for precision testing.
A structure holding the options for precision estimation.
Definition prec_structures.h:77
A structure holding the result of precision estimation.
Definition prec_structures.h:25
long double maxErr
Estimated maximum absolute error on interval.
Definition prec_structures.h:37
long double absErr
Estimated absolute error on interval.
Definition prec_structures.h:49
An interval on the real numbers.
Definition interval.h:16
long double length() const
Returns the length of the interval.
Definition interval.h:33
long double a
Lower extreme of the interval.
Definition interval.h:19
long double b
Upper extreme of the interval.
Definition interval.h:22
A source of pseudorandom numbers.
Definition random.h:39
long double uniform(long double a, long double b)
Generate a uniformly distributed random number.
Definition random.h:76