Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
dataset.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_DATASET_H
7#define THEORETICA_DATASET_H
8
9#ifndef THEORETICA_NO_PRINT
10#include <sstream>
11#include <ostream>
12#endif
13
14#include "./core_traits.h"
15#include "./constants.h"
16#include "./error.h"
17
18
19namespace theoretica {
20
21
22 // Operations on datasets and generic ordered sets of numbers
23 // The Vector type must have size() and operator[]() methods
24 // (e.g. std::vector<real> and vec<real>)
25
26
28 template<typename Vector, enable_vector<Vector> = true>
29 inline auto product(const Vector& X) {
30
31 if(!X.size()) {
32 TH_MATH_ERROR("product", X.size(), INVALID_ARGUMENT);
34 }
35
36 auto res = X[0];
37 for(unsigned int i = 1; i < X.size(); i++)
38 res *= X[i];
39
40 return res;
41 }
42
43
45 template<typename Vector>
46 inline auto product_sum(const Vector& X, const Vector& Y) {
47
48 if(X.size() != Y.size() || !X.size()) {
49 TH_MATH_ERROR("product_sum", X.size(), INVALID_ARGUMENT);
50 return nan();
51 }
52
53 auto res = X[0] * Y[0];
54 for(unsigned int i = 1; i < X.size(); i++)
55 res += X[i] * Y[i];
56
57 return res;
58 }
59
60
62 template<typename Vector>
63 inline auto product_sum_squares(const Vector& X, const Vector& Y) {
64
65 if(X.size() != Y.size() || !X.size()) {
66 TH_MATH_ERROR("product_sum_squares", X.size(), INVALID_ARGUMENT);
67 return nan();
68 }
69
70 auto res = (X[0] * X[0]) * (Y[0] * Y[0]);
71
72 for(unsigned int i = 1; i < X.size(); i++)
73 res += (X[i] * X[i]) * (Y[i] * Y[i]);
74
75 return res;
76 }
77
78
80 template<typename Vector>
81 inline auto product_sum(const Vector& X, const Vector& Y, const Vector& Z) {
82
83 if(X.size() != Y.size() || X.size() != Z.size() || !X.size()) {
84 TH_MATH_ERROR("product_sum", X.size(), INVALID_ARGUMENT);
85 return nan();
86 }
87
88 auto res = X[0] * Y[0] * Z[0];
89 for(unsigned int i = 1; i < X.size(); i++)
90 res += X[i] * Y[i] * Z[i];
91
92 return res;
93 }
94
95
97 template<typename Vector>
98 inline auto quotient_sum(const Vector& X, const Vector& Y) {
99
100 if(X.size() != Y.size()) {
101 TH_MATH_ERROR("quotient_sum", X.size(), INVALID_ARGUMENT);
102 return nan();
103 }
104
105 if(abs(Y[0]) < MACH_EPSILON) {
106 TH_MATH_ERROR("quotient_sum", Y[0], DIV_BY_ZERO);
107 return nan();
108 }
109
110 auto res = X[0] / Y[0];
111 for(unsigned int i = 1; i < X.size(); i++) {
112
113 if(abs(Y[i]) < MACH_EPSILON) {
114 TH_MATH_ERROR("quotient_sum", Y[i], DIV_BY_ZERO);
115 return nan();
116 }
117
118 res += X[i] / Y[i];
119 }
120
121 return res;
122 }
123
124
126 template<typename Vector>
127 inline auto sum_squares(const Vector& X) {
128
129 if(!X.size()) {
130 TH_MATH_ERROR("sum_squares", X.size(), INVALID_ARGUMENT);
131 return nan();
132 }
133
134 auto res = X[0] * X[0];
135 for(unsigned int i = 1; i < X.size(); i++)
136 res += X[i] * X[i];
137
138 return res;
139 }
140
141
147 template<typename Vector>
148 inline real sum_compensated(const Vector& X) {
149
150 // Total sum
151 real sum = 0;
152
153 // Correction term
154 real corr = 0;
155
156 for (unsigned int i = 0; i < X.size(); i++) {
157
158 const real temp = sum + X[i];
159
160 // Sort the two addends to preserve bits
161 corr += (abs(sum) >= abs(X[i]))
162 ? ((sum - temp) + X[i])
163 : ((X[i] - temp) + sum);
164
165 sum = temp;
166 }
167
168 // Apply correction
169 return sum + corr;
170 }
171
172
181 template<typename Vector>
183 const Vector& X, size_t begin = 0,
184 size_t end = 0, size_t base_size = 128) {
185
186 if(end == 0)
187 end = X.size();
188
189 real sum = 0;
190
191 // Base case with given size (defaults to 128)
192 if((end - begin) <= base_size) {
193
194 for (size_t i = begin; i < end; ++i)
195 sum += X[i];
196
197 } else {
198
199 // Recursive sum of two halves
200 const size_t m = (end - begin) / 2;
201 const size_t cutoff = begin + m;
202
203 sum = sum_pairwise(X, begin, cutoff, base_size)
204 + sum_pairwise(X, cutoff, end, base_size);
205 }
206
207 return sum;
208 }
209
210
215 template <
216 typename Vector,
217 std::enable_if_t<has_real_elements<Vector>::value> = true
218 >
219 inline auto sum(const Vector& X) {
220 return sum_pairwise(X);
221 }
222
223
227 template<typename Vector>
228 inline auto sum(const Vector& X) {
229
230 // Use pairwise sum to reduce floating point errors.
232 return sum_pairwise(X);
233
234 auto res = X[0];
235 for(unsigned int i = 1; i < X.size(); i++)
236 res += X[i];
237
238 return res;
239 }
240
241
249 template<typename Vector, typename Function>
250 inline Vector& apply(Function f, Vector& X) {
251
252 for (unsigned int i = 0; i < X.size(); i++)
253 X[i] = f(X[i]);
254
255 return X;
256 }
257
258
265 template<typename Vector1, typename Vector2 = Vector1, typename Function>
266 inline Vector2& map(Function f, const Vector1& src, Vector2& dest) {
267
268 if(src.size() != dest.size()) {
269 TH_MATH_ERROR("th::map", dest.size(), INVALID_ARGUMENT);
270 dest = Vector2(nan());
271 return dest;
272 }
273
274 for (unsigned int i = 0; i < src.size(); i++)
275 dest[i] = f(src[i]);
276
277 return dest;
278 }
279
280
286 template<typename Vector2, typename Vector1, typename Function>
287 inline Vector2 map(Function f, const Vector1& X) {
288
289 Vector2 res;
290 res.resize(X.size());
291
292 for (unsigned int i = 0; i < X.size(); i++)
293 res[i] = f(X[i]);
294
295 return res;
296 }
297
298
304 template<typename Vector, typename Function>
305 inline Vector map(Function f, const Vector& X) {
307 }
308
309
311 template<typename Vector1, typename Vector2, typename Vector3 = Vector1>
312 inline Vector3 concatenate(const Vector1& v1, const Vector2& v2) {
313
314 Vector3 res;
315 res.resize(v1.size() + v2.size());
316 const unsigned int offset = v1.size();
317
318 for (unsigned int i = 0; i < offset; ++i)
319 res[i] = v1[i];
320
321 for (unsigned int i = 0; i < v2.size(); ++i)
322 res[i + offset] = v2[i];
323
324 return res;
325 }
326
327
329 template<typename Vector>
330 inline auto max(const Vector& X) {
331
333
334 if(!X.size()) {
335 TH_MATH_ERROR("max", X.size(), INVALID_ARGUMENT);
336 return Type(nan());
337 }
338
339 auto curr = X[0];
340
341 for (unsigned int i = 1; i < X.size(); ++i)
342 if(X[i] > curr)
343 curr = X[i];
344
345 return curr;
346 }
347
348
350 template<typename Vector>
351 inline auto min(const Vector& X) {
352
354
355 if(!X.size()) {
356 TH_MATH_ERROR("min", X.size(), INVALID_ARGUMENT);
357 return Type(nan());
358 }
359
360 auto curr = X[0];
361
362 for (unsigned int i = 1; i < X.size(); ++i)
363 if(X[i] < curr)
364 curr = X[i];
365
366 return curr;
367 }
368
369
370 // Different types of means
371
372
374 template<typename Dataset>
375 inline real arithmetic_mean(const Dataset& data) {
376
377 if(!data.size()) {
378 TH_MATH_ERROR("arithmetic_mean", data.size(), DIV_BY_ZERO);
379 return nan();
380 }
381
382 // Sum of x_i / N
383 return sum(data) / (real) data.size();
384 }
385
386
388 template<typename Dataset>
389 inline real harmonic_mean(const Dataset& data) {
390
391 if(!data.size()) {
392 TH_MATH_ERROR("harmonic_mean", data.size(), DIV_BY_ZERO);
393 return nan();
394 }
395
396 real res = 0;
397
398 for (unsigned int i = 0; i < data.size(); ++i) {
399
400 if(data[i] == 0) {
401 TH_MATH_ERROR("harmonic_mean", data[i], DIV_BY_ZERO);
402 return nan();
403 }
404
405 res += 1.0 / data[i];
406 }
407
408 return static_cast<real>(data.size()) / res;
409 }
410
411
414 template<typename Dataset>
415 inline real geometric_mean(const Dataset& data) {
416 return root(product(data), data.size());
417 }
418
419
422 template<typename Dataset1, typename Dataset2>
423 inline real weighted_mean(const Dataset1& data, const Dataset2& weights) {
424
425 // Sum of x_i * w_i / Sum of w_i
426 return product_sum(data, weights) / sum(weights);
427 }
428
429
432 template<typename Dataset>
433 inline real quadratic_mean(const Dataset& data) {
434
435 if(!data.size()) {
436 TH_MATH_ERROR("quadratic_mean", data.size(), INVALID_ARGUMENT);
437 return nan();
438 }
439
440 return sqrt(sum_squares(data) / data.size());
441 }
442}
443
444#endif
Mathematical constants and default algorithm parameters.
#define TH_CONSTIF
Enable constexpr in if statements if C++17 is supported.
Definition constants.h:169
Fundamental type traits.
Error handling.
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:219
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
double real
A real number, defined as a floating point type.
Definition constants.h:198
real weighted_mean(const Dataset1 &data, const Dataset2 &weights)
Compute the weighted mean of a set of values <data> and <weights> must have the same size.
Definition dataset.h:423
real arithmetic_mean(const Dataset &data)
Compute the arithmetic mean of a set of values.
Definition dataset.h:375
auto min(const Vector &X)
Finds the minimum value inside a dataset.
Definition dataset.h:351
real harmonic_mean(const Dataset &data)
Compute the harmonic mean of a set of values.
Definition dataset.h:389
real sum_compensated(const Vector &X)
Compute the sum of a set of values using the compensated Neumaier-Kahan-Babushka summation algorithm ...
Definition dataset.h:148
dual2 sqrt(dual2 x)
Compute the square root of a second order dual number.
Definition dual2_functions.h:54
Vector3 concatenate(const Vector1 &v1, const Vector2 &v2)
Concatenate two datasets to form a single one.
Definition dataset.h:312
Vector2 & map(Function f, const Vector1 &src, Vector2 &dest)
Get a new vector obtained by applying the function element-wise.
Definition dataset.h:266
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
real sum_pairwise(const Vector &X, size_t begin=0, size_t end=0, size_t base_size=128)
Compute the sum of a set of values using pairwise summation to reduce round-off error.
Definition dataset.h:182
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
auto sum(const Vector &X)
Compute the sum of a vector of real values using pairwise summation to reduce round-off error.
Definition dataset.h:219
real geometric_mean(const Dataset &data)
Compute the geometric mean of a set of values as .
Definition dataset.h:415
auto max(const Vector &X)
Finds the maximum value inside a dataset.
Definition dataset.h:330
auto product_sum(const Vector &X, const Vector &Y)
Sum the products of two sets of values.
Definition dataset.h:46
auto product_sum_squares(const Vector &X, const Vector &Y)
Sum the products of the squares of two sets of data.
Definition dataset.h:63
real quadratic_mean(const Dataset &data)
Compute the quadratic mean (Root Mean Square) of a set of values .
Definition dataset.h:433
auto quotient_sum(const Vector &X, const Vector &Y)
Sum the quotients of two sets of values.
Definition dataset.h:98
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207
Vector & apply(Function f, Vector &X)
Apply a function to a set of values element-wise.
Definition dataset.h:250
real root(real x, int n)
Compute the n-th root of x.
Definition real_analysis.h:812
auto product(const Vector &X)
Compute the product of a set of values.
Definition dataset.h:29
real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
auto sum_squares(const Vector &X)
Sum the squares of a set of values.
Definition dataset.h:127