Line data Source code
1 : // SPDX-License-Identifier: BSD-3-Clause
2 : // Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
3 : /// @file
4 : #pragma once
5 :
6 : #include <cmath>
7 :
8 : #include "scipp/common/numeric.h"
9 : #include "scipp/common/span.h"
10 : #include "scipp/core/dtype.h"
11 :
12 : namespace scipp::core {
13 :
14 : /// A value/variance pair with operators that propagate uncertainties.
15 : ///
16 : /// This is intended for small T such as double, float, and int. It is the
17 : /// central implementation of uncertainty propagation in scipp, for built-in
18 : /// operations as well as custom operations using one of the transform
19 : /// functions. Since T is assumed to be small it is copied into the class and
20 : /// extracted later. See also ValuesAndVariances.
21 : template <class T> struct ValueAndVariance {
22 : using value_type = std::remove_cv_t<T>;
23 :
24 : /// This constructor is essential to prevent warnings about narrowing the
25 : /// types in initializer lists used below
26 : template <class T1, class T2>
27 8029911 : constexpr ValueAndVariance(T1 t1, T2 t2) noexcept : value(t1), variance(t2) {}
28 2118 : constexpr ValueAndVariance(const double val) noexcept
29 2118 : : value(val), variance(0.0) {}
30 : T value;
31 : T variance;
32 :
33 : template <class T2>
34 0 : constexpr auto &operator=(const ValueAndVariance<T2> other) noexcept {
35 0 : value = other.value;
36 0 : variance = other.variance;
37 0 : return *this;
38 : }
39 :
40 0 : template <class T2> constexpr auto &operator=(const T2 other) noexcept {
41 0 : value = other;
42 0 : variance = 0.0;
43 0 : return *this;
44 : }
45 :
46 63691 : template <class T2> constexpr auto &operator+=(const T2 other) noexcept {
47 63691 : return *this = *this + other;
48 : }
49 0 : template <class T2> constexpr auto &operator-=(const T2 other) noexcept {
50 0 : return *this = *this - other;
51 : }
52 84 : template <class T2> constexpr auto &operator*=(const T2 other) noexcept {
53 84 : return *this = *this * other;
54 : }
55 0 : template <class T2> constexpr auto &operator/=(const T2 other) noexcept {
56 0 : return *this = *this / other;
57 : }
58 :
59 : template <class T2>
60 0 : explicit constexpr operator ValueAndVariance<T2>() const noexcept {
61 0 : return {static_cast<T2>(value), static_cast<T2>(variance)};
62 : }
63 : };
64 :
65 : template <class T>
66 84 : constexpr auto operator-(const ValueAndVariance<T> a) noexcept {
67 84 : return ValueAndVariance{-a.value, a.variance};
68 : }
69 :
70 : template <class B, class E>
71 9 : constexpr auto pow(const ValueAndVariance<B> base, const E exponent) noexcept {
72 9 : const auto pow_1 = numeric::pow(base.value, exponent - 1);
73 9 : const auto var_factor = std::abs(exponent) * pow_1;
74 9 : return ValueAndVariance{pow_1 * base.value,
75 9 : var_factor * var_factor * base.variance};
76 : }
77 :
78 0 : template <class T> constexpr auto sqrt(const ValueAndVariance<T> a) noexcept {
79 : using std::sqrt;
80 0 : return ValueAndVariance{sqrt(a.value),
81 0 : static_cast<T>(0.25 * (a.variance / a.value))};
82 : }
83 :
84 1695 : template <class T> constexpr auto abs(const ValueAndVariance<T> a) noexcept {
85 : using std::abs;
86 1695 : return ValueAndVariance{abs(a.value), a.variance};
87 : }
88 :
89 0 : template <typename T> constexpr auto exp(const ValueAndVariance<T> a) noexcept {
90 : using std::exp;
91 0 : const auto val = exp(a.value);
92 0 : return ValueAndVariance(val, val * val * a.variance);
93 : }
94 :
95 0 : template <typename T> constexpr auto log(const ValueAndVariance<T> a) noexcept {
96 : using std::log;
97 0 : return ValueAndVariance(log(a.value), a.variance / (a.value * a.value));
98 : }
99 :
100 : template <typename T>
101 0 : constexpr auto log10(const ValueAndVariance<T> a) noexcept {
102 : using std::log10;
103 0 : const auto x = a.value * std::log(static_cast<T>(10.0L));
104 0 : return ValueAndVariance(log10(a.value), a.variance / (x * x));
105 : }
106 :
107 1612 : template <class T> constexpr auto isnan(const ValueAndVariance<T> a) noexcept {
108 : using numeric::isnan;
109 1612 : return isnan(a.value);
110 : }
111 :
112 2 : template <class T> constexpr auto isinf(const ValueAndVariance<T> a) noexcept {
113 : using numeric::isinf;
114 2 : return isinf(a.value);
115 : }
116 :
117 : template <class T>
118 42 : constexpr auto isfinite(const ValueAndVariance<T> a) noexcept {
119 : using numeric::isfinite;
120 42 : return isfinite(a.value);
121 : }
122 :
123 : template <class T>
124 0 : constexpr auto signbit(const ValueAndVariance<T> a) noexcept {
125 0 : return numeric::signbit(a.value);
126 : }
127 :
128 : template <class T>
129 2 : constexpr auto isposinf(const ValueAndVariance<T> a) noexcept {
130 2 : return numeric::isinf(a.value) && !signbit(a);
131 : }
132 :
133 : template <class T>
134 2 : constexpr auto isneginf(const ValueAndVariance<T> a) noexcept {
135 2 : return numeric::isinf(a.value) && signbit(a);
136 : }
137 :
138 : template <class T1, class T2>
139 63706 : constexpr auto operator+(const ValueAndVariance<T1> a,
140 : const ValueAndVariance<T2> b) noexcept {
141 63706 : return ValueAndVariance{a.value + b.value, a.variance + b.variance};
142 : }
143 : template <class T1, class T2>
144 2 : constexpr auto operator-(const ValueAndVariance<T1> a,
145 : const ValueAndVariance<T2> b) noexcept {
146 2 : return ValueAndVariance{a.value - b.value, a.variance + b.variance};
147 : }
148 : template <class T1, class T2>
149 8 : constexpr auto operator*(const ValueAndVariance<T1> a,
150 : const ValueAndVariance<T2> b) noexcept {
151 8 : return ValueAndVariance{a.value * b.value,
152 8 : a.variance * b.value * b.value +
153 8 : b.variance * a.value * a.value};
154 : }
155 : template <class T1, class T2>
156 14 : constexpr auto operator/(const ValueAndVariance<T1> a,
157 : const ValueAndVariance<T2> b) noexcept {
158 : return ValueAndVariance{
159 14 : a.value / b.value,
160 14 : (a.variance + b.variance * (a.value * a.value) / (b.value * b.value)) /
161 14 : (b.value * b.value)};
162 : }
163 :
164 : template <class T1, class T2>
165 0 : constexpr auto operator+(const ValueAndVariance<T1> a, const T2 b) noexcept {
166 0 : return ValueAndVariance{a.value + b, a.variance};
167 : }
168 : template <class T1, class T2>
169 17 : constexpr auto operator+(const T2 a, const ValueAndVariance<T1> b) noexcept {
170 17 : return ValueAndVariance{a + b.value, b.variance};
171 : }
172 : template <class T1, class T2>
173 1655 : constexpr auto operator-(const ValueAndVariance<T1> a, const T2 b) noexcept {
174 1655 : return ValueAndVariance{a.value - b, a.variance};
175 : }
176 : template <class T1, class T2>
177 0 : constexpr auto operator-(const T2 a, const ValueAndVariance<T1> b) noexcept {
178 0 : return ValueAndVariance{a - b.value, b.variance};
179 : }
180 : template <class T1, class T2>
181 359 : constexpr auto operator*(const ValueAndVariance<T1> a, const T2 b) noexcept {
182 : // Order of operations to prevent integer overflow of `b*b`.
183 359 : return ValueAndVariance{a.value * b, a.variance * b * b};
184 : }
185 : template <class T1, class T2>
186 132 : constexpr auto operator*(const T1 a, const ValueAndVariance<T2> b) noexcept {
187 : // Order of operations to prevent integer overflow of `a*a`.
188 132 : return ValueAndVariance{a * b.value, b.variance * a * a};
189 : }
190 : template <class T1, class T2>
191 0 : constexpr auto operator/(const ValueAndVariance<T1> a, const T2 b) noexcept {
192 0 : const auto denominator = [b]() {
193 : if constexpr (std::is_integral_v<T2>)
194 : // Prevent int overflow.
195 0 : return static_cast<T1>(b) * static_cast<T1>(b);
196 : else
197 0 : return b * b;
198 0 : }();
199 0 : return ValueAndVariance{a.value / b, a.variance / denominator};
200 : }
201 : template <class T1, class T2>
202 0 : constexpr auto operator/(const T1 a, const ValueAndVariance<T2> b) noexcept {
203 0 : return ValueAndVariance{a / b.value, b.variance * a * a /
204 0 : (b.value * b.value) /
205 0 : (b.value * b.value)};
206 : }
207 :
208 : // Comparison operators. Note that all of these IGNORE VARIANCES
209 : template <class A, class B>
210 11 : constexpr auto operator==(const ValueAndVariance<A> a,
211 : const ValueAndVariance<B> b) noexcept {
212 11 : return a.value == b.value;
213 : }
214 : template <class A, class B>
215 0 : constexpr auto operator==(const ValueAndVariance<A> a, const B b) noexcept {
216 0 : return a.value == b;
217 : }
218 : template <class A, class B>
219 0 : constexpr auto operator==(const A a, const ValueAndVariance<B> b) noexcept {
220 0 : return a == b.value;
221 : }
222 :
223 : template <class A, class B>
224 2 : constexpr auto operator!=(const ValueAndVariance<A> a,
225 : const ValueAndVariance<B> b) noexcept {
226 2 : return !(a == b);
227 : }
228 :
229 : template <class A, class B>
230 0 : constexpr auto operator!=(const ValueAndVariance<A> a, const B b) noexcept {
231 0 : return !(a == b);
232 : }
233 :
234 : template <class A, class B>
235 0 : constexpr auto operator!=(const A a, const ValueAndVariance<B> b) noexcept {
236 0 : return !(a == b);
237 : }
238 :
239 : template <class A, class B>
240 36 : constexpr auto operator<(const ValueAndVariance<A> a,
241 : const ValueAndVariance<B> b) noexcept {
242 36 : return a.value < b.value;
243 : }
244 : template <class A, class B>
245 53 : constexpr auto operator<(const ValueAndVariance<A> a, const B b) noexcept {
246 53 : return a.value < b;
247 : }
248 : template <class A, class B>
249 0 : constexpr auto operator<(const A a, const ValueAndVariance<B> b) noexcept {
250 0 : return a < b.value;
251 : }
252 :
253 : template <class A, class B>
254 2 : constexpr auto operator<=(const ValueAndVariance<A> a,
255 : const ValueAndVariance<B> b) noexcept {
256 2 : return a.value <= b.value;
257 : }
258 : template <class A, class B>
259 1656 : constexpr auto operator<=(const ValueAndVariance<A> a, const B b) noexcept {
260 1656 : return a.value <= b;
261 : }
262 : template <class A, class B>
263 0 : constexpr auto operator<=(const A a, const ValueAndVariance<B> b) noexcept {
264 0 : return a <= b.value;
265 : }
266 :
267 : template <class A, class B>
268 9 : constexpr auto operator>(const ValueAndVariance<A> a,
269 : const ValueAndVariance<B> b) noexcept {
270 9 : return a.value > b.value;
271 : }
272 : template <class A, class B>
273 9 : constexpr auto operator>(const ValueAndVariance<A> a, const B b) noexcept {
274 9 : return a.value > b;
275 : }
276 : template <class A, class B>
277 0 : constexpr auto operator>(const A a, const ValueAndVariance<B> b) noexcept {
278 0 : return a > b.value;
279 : }
280 :
281 : template <class A, class B>
282 2 : constexpr auto operator>=(const ValueAndVariance<A> a,
283 : const ValueAndVariance<B> b) noexcept {
284 2 : return a.value >= b.value;
285 : }
286 : template <class A, class B>
287 0 : constexpr auto operator>=(const ValueAndVariance<A> a, const B b) noexcept {
288 0 : return a.value >= b;
289 : }
290 : template <class A, class B>
291 0 : constexpr auto operator>=(const A a, const ValueAndVariance<B> b) noexcept {
292 0 : return a >= b.value;
293 : }
294 : // end comparison operators
295 :
296 : template <class T>
297 48 : constexpr auto min(const ValueAndVariance<T> a,
298 : const ValueAndVariance<T> b) noexcept {
299 48 : return a.value < b.value ? a : b;
300 : }
301 : template <class T>
302 714 : constexpr auto max(const ValueAndVariance<T> a,
303 : const ValueAndVariance<T> b) noexcept {
304 714 : return a.value > b.value ? a : b;
305 : }
306 :
307 : /// Deduction guide for class ValueAndVariances. Using decltype to deal with
308 : /// potential mixed-type val and var arguments arising in binary operations
309 : /// between, e.g., double and float.
310 : template <class T1, class T2>
311 : ValueAndVariance(const T1 &val, const T2 &var)
312 : -> ValueAndVariance<decltype(T1() + T2())>;
313 : template <class T>
314 : ValueAndVariance(const scipp::span<T> &val, const scipp::span<T> &var)
315 : -> ValueAndVariance<scipp::span<T>>;
316 :
317 : template <class T> struct is_ValueAndVariance : std::false_type {};
318 : template <class T>
319 : struct is_ValueAndVariance<ValueAndVariance<T>> : std::true_type {};
320 : template <class T>
321 : inline constexpr bool is_ValueAndVariance_v = is_ValueAndVariance<T>::value;
322 :
323 : } // namespace scipp::core
324 :
325 : namespace scipp {
326 : using core::is_ValueAndVariance_v;
327 : using core::ValueAndVariance;
328 : } // namespace scipp
|
