1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7 #ifndef VARRAY_H
8 #define VARRAY_H
9
10 #include <iostream> // cerr
11 #include <vector>
12 #include <cstddef>
13 #include <cfloat>
14 #include <cassert>
15 #include "compare.h"
16
17 // compare
18 // clang-format off
19 const auto is_not_equal = [](auto a, auto b) noexcept { return (a < b || b < a); };
20 const auto is_equal = [](auto a, auto b) noexcept { return !(a < b || b < a); };
21 const auto dbl_is_equal = [](auto a, auto b) noexcept { return (std::isnan(a) || std::isnan(b) ? (std::isnan(a) && std::isnan(b)) : is_equal(a, b)); };
22 // clang-format on
23
24 // unary operators
25 // clang-format off
26 const auto unary_op_abs = [](const double x) noexcept { return std::fabs(x); };
27 const auto unary_op_int = [](const double x) noexcept { return (int)(x); };
28 const auto unary_op_nint = [](const double x) noexcept { return std::round(x); };
29 const auto unary_op_sqr = [](const double x) noexcept { return x * x; };
30 const auto unary_op_nop = [](const double x) noexcept { return x; };
31 const auto unary_op_reci = [](const double x) noexcept { return 1.0 / x; };
32 const auto unary_op_not = [](const double x) noexcept { return is_equal(x, 0.0); };
33 const auto unary_op_exp = [](const double x) noexcept { return std::exp(x); };
34 const auto unary_op_log = [](const double x) noexcept { return std::log(x); };
35 const auto unary_op_log10 = [](const double x) noexcept { return std::log10(x); };
36 const auto unary_op_sin = [](const double x) noexcept { return std::sin(x); };
37 const auto unary_op_cos = [](const double x) noexcept { return std::cos(x); };
38 const auto unary_op_tan = [](const double x) noexcept { return std::tan(x); };
39 const auto unary_op_asin = [](const double x) noexcept { return std::asin(x); };
40 const auto unary_op_acos = [](const double x) noexcept { return std::acos(x); };
41 const auto unary_op_atan = [](const double x) noexcept { return std::atan(x); };
42 // clang-format on
43
44 // binary operators
45 // clang-format off
46 const auto binary_op_LT = [](const double x, const double y) noexcept { return x < y; };
47 const auto binary_op_GT = [](const double x, const double y) noexcept { return x > y; };
48 const auto binary_op_LE = [](const double x, const double y) noexcept { return x <= y; };
49 const auto binary_op_GE = [](const double x, const double y) noexcept { return x >= y; };
50 const auto binary_op_NE = [](const double x, const double y) noexcept { return is_not_equal(x, y); };
51 const auto binary_op_EQ = [](const double x, const double y) noexcept { return is_equal(x, y); };
52 const auto binary_op_LEG = [](const double x, const double y) noexcept { return (x < y) ? -1.0 : (x > y); };
53 const auto binary_op_AND = [](const double x, const double y) noexcept { return is_not_equal(x, 0.0) && is_not_equal(y, 0.0); };
54 const auto binary_op_OR = [](const double x, const double y) noexcept { return is_not_equal(x, 0.0) || is_not_equal(y, 0.0); };
55 const auto binary_op_POW = [](const double x, const double y) noexcept { return std::pow(x, y); };
56 const auto binary_op_ADD = [](const double x, const double y) noexcept { return x + y; };
57 const auto binary_op_SUB = [](const double x, const double y) noexcept { return x - y; };
58 const auto binary_op_MUL = [](const double x, const double y) noexcept { return x * y; };
59 const auto binary_op_DIV = [](const double x, const double y) noexcept { return x / y; };
60 // clang-format on
61
62 //#define CHECK_UNUSED_VECTOR 1
63
64 #ifdef CHECK_UNUSED_VECTOR
65 template <typename T>
66 class
67 #ifdef WARN_UNUSED
68 [[gnu::warn_unused]]
69 #endif
70 CheckVector
71 {
72 public:
73 T *ptr;
74 size_t m_count;
75
CheckVector()76 CheckVector() : ptr(dummy), m_count(0) { }
CheckVector(size_t count)77 explicit CheckVector(size_t count) : ptr(dummy), m_count(count) { }
CheckVector(size_t count,const T & value)78 CheckVector(size_t count, const T& value) : ptr(dummy), m_count(count) { ptr[0] = value; }
79
begin()80 T * begin() noexcept { return &ptr[0]; }
end()81 T * end() noexcept { return &ptr[0] + 1; }
begin()82 const T * begin() const noexcept { return &ptr[0]; }
end()83 const T * end() const noexcept { return &ptr[0] + 1; }
84
empty()85 bool empty() const { return true; }
size()86 size_t size() const { return m_count; }
clear()87 void clear() { }
shrink_to_fit()88 void shrink_to_fit() { }
89
resize(size_t count)90 void resize(size_t count) { m_count = count; }
resize(size_t count,const T & value)91 void resize(size_t count, const T& value) { m_count = count; ptr[0] = value; }
92
reserve(size_t new_cap)93 void reserve(size_t new_cap) { (void)new_cap; }
push_back(const T & value)94 void push_back(const T& value) { ptr[0] = value; }
assign(size_t count,const T & value)95 void assign(size_t count, const T& value) { (void)count; ptr[0] = value; }
96
data()97 T * data() noexcept { return ptr; }
data()98 const T * data() const noexcept { return ptr; }
99
100 T& operator[](size_t pos) { (void)pos; return ptr[0]; }
101 const T & operator[](size_t pos) const { (void)pos; return ptr[0]; }
102 CheckVector& operator=(const CheckVector& other) { *this = other; ptr = dummy; m_count = 0; return *this; }
103 CheckVector& operator=(CheckVector&& other) { *this = other; ptr = dummy; m_count = 0; return *this; }
104 CheckVector& operator=(const std::vector<T>& other) { (void)other; ptr = dummy; m_count = 0; return *this; }
105 CheckVector& operator=(std::vector<T>&& other) { (void)other; ptr = dummy; m_count = 0; return *this; }
106 // Copy constructor
CheckVector(const CheckVector & v2)107 CheckVector(const CheckVector &v2) { m_count = 0; ptr = dummy; ptr[0] = v2.ptr[0]; }
CheckVector(const std::vector<T> & v2)108 CheckVector(const std::vector<T> &v2) { (void)v2; m_count = 0; ptr = dummy; }
109
110 bool operator==(const CheckVector& other) { (void)other; return true; }
111
112 private:
113 T dummy[1] {};
114 };
115
116 template <typename T>
117 using Varray = CheckVector<T>;
118
119 #else
120
121 template <typename T>
122 using Varray = std::vector<T>;
123
124 #endif
125
126 template <typename T>
127 using Varray2D = Varray<Varray<T>>;
128
129 template <typename T>
130 using Varray3D = Varray2D<Varray<T>>;
131
132 template <typename T>
133 using Varray4D = Varray3D<Varray<T>>;
134
135
136 struct MinMax
137 {
138 double min = DBL_MAX;
139 double max = -DBL_MAX;
140 size_t n = 0;
MinMaxMinMax141 MinMax() {};
MinMaxMinMax142 MinMax(double rmin, double rmax, size_t rn) : min(rmin), max(rmax), n(rn) {};
MinMaxMinMax143 MinMax(double rmin, double rmax) : min(rmin), max(rmax), n(0) {};
144 };
145
146 struct MinMaxSum : MinMax
147 {
148 double sum = 0.0;
MinMaxSumMinMaxSum149 MinMaxSum() {};
MinMaxSumMinMaxSum150 MinMaxSum(double rmin, double rmax, double rsum, size_t rn) : sum(rsum) { min = rmin; max = rmax; n = rn; };
MinMaxSumMinMaxSum151 MinMaxSum(double rmin, double rmax, double rsum) : sum(rsum) { min = rmin; max = rmax; n = 0; };
152 };
153
154 struct MinMaxMean : MinMax
155 {
156 double mean = 0.0;
MinMaxMeanMinMaxMean157 MinMaxMean() {};
MinMaxMeanMinMaxMean158 MinMaxMean(double rmin, double rmax, double rmean, size_t rn) : mean(rmean) { min = rmin; max = rmax; n = rn; };
MinMaxMeanMinMaxMean159 MinMaxMean(double rmin, double rmax, double rmean) : mean(rmean) { min = rmin; max = rmax; n = 0; };
160 };
161
162
163 template <typename T>
164 inline void
varray_free(Varray<T> & v)165 varray_free(Varray<T> &v)
166 {
167 v.clear();
168 v.shrink_to_fit();
169 }
170
171 #define varrayResize(p, s) varray_resize((p), (s), __FILE__, __LINE__)
172 #define varrayResizeInit(p, s, v) varray_resize((p), (s), (v), __FILE__, __LINE__)
173
174 template <typename T>
175 inline void
varray_resize(Varray<T> & v,size_t count,const char * file,int line)176 varray_resize(Varray<T> &v, size_t count, const char *file, int line)
177 {
178 try
179 {
180 v.resize(count);
181 }
182 catch (const std::exception &e)
183 {
184 std::cerr << "Exception caught when trying to allocate " << count << " vector elements: "
185 << e.what() << " in " << file << ":" << line << '\n';
186 throw;
187 }
188 }
189
190 template <typename T>
191 inline void
varray_resize(Varray<T> & v,size_t count,T value,const char * file,int line)192 varray_resize(Varray<T> &v, size_t count, T value, const char *file, int line)
193 {
194 try
195 {
196 v.resize(count, value);
197 }
198 catch (const std::exception &e)
199 {
200 std::cerr << "Exception caught when trying to allocate " << count << " vector elements: "
201 << e.what() << " in " << file << ":" << line << '\n';
202 throw;
203 }
204 }
205
206 /*
207 template <class T, size_t ROW, size_t COL>
208 using Matrix = std::array<std::array<T, COL>, ROW>;
209 */
210
211 // clang-format off
212 #define MADDMN(x, y) ((DBL_IS_EQUAL((x), missval1) || DBL_IS_EQUAL((y), missval2)) ? missval1 : (x) + (y))
213 #define MSUBMN(x, y) ((DBL_IS_EQUAL((x), missval1) || DBL_IS_EQUAL((y), missval2)) ? missval1 : (x) - (y))
214 #define MMULMN(x, y) ((DBL_IS_EQUAL((x), 0) || DBL_IS_EQUAL((y), 0)) ? 0 : (DBL_IS_EQUAL((x), missval1) || DBL_IS_EQUAL((y), missval2)) ? missval1 : (x) * (y))
215 #define MDIVMN(x, y) ((DBL_IS_EQUAL((x), missval1) || DBL_IS_EQUAL((y), missval2) || DBL_IS_EQUAL((y), 0)) ? missval1 : (x) / (y))
216 #define MPOWMN(x, y) ((DBL_IS_EQUAL((x), missval1) || DBL_IS_EQUAL((y), missval2)) ? missval1 : std::pow((x), (y)))
217 #define MSQRTMN(x) ((DBL_IS_EQUAL((x), missval1) || (x) < 0) ? missval1 : std::sqrt(x))
218
219 #define ADD(x, y) ((x) + (y))
220 #define SUB(x, y) ((x) - (y))
221 #define MUL(x, y) ((x) * (y))
222 #define DIV(x, y) (IS_EQUAL((y), 0.) ? missval1 : (x) / (y))
223 #define POW(x, y) std::pow((x), (y))
224 #define SQRT(x) std::sqrt(x)
225
226 #define ADDM(x, y) ((IS_EQUAL((x), missval1) || IS_EQUAL((y), missval2)) ? missval1 : (x) + (y))
227 #define SUBM(x, y) ((IS_EQUAL((x), missval1) || IS_EQUAL((y), missval2)) ? missval1 : (x) - (y))
228 #define MULM(x, y) ((IS_EQUAL((x), 0) || IS_EQUAL((y), 0)) ? 0 : (IS_EQUAL((x), missval1) || IS_EQUAL((y), missval2)) ? missval1 : (x) * (y))
229 #define DIVM(x, y) ((IS_EQUAL((x), missval1) || IS_EQUAL((y), missval2) || IS_EQUAL((y), 0)) ? missval1 : (x) / (y))
230 #define POWM(x, y) ((IS_EQUAL((x), missval1) || IS_EQUAL((y), missval2)) ? missval1 : std::pow((x), (y)))
231 #define SQRTM(x) ((IS_EQUAL((x), missval1) || (x) < 0) ? missval1 : std::sqrt(x))
232
233 #define ADDMN(x, y) FADDMN(x, y, missval1, missval2)
234 #define SUBMN(x, y) FSUBMN(x, y, missval1, missval2)
235 #define MULMN(x, y) FMULMN(x, y, missval1, missval2)
236 #define DIVMN(x, y) FDIVMN(x, y, missval1, missval2)
237 #define POWMN(x, y) FPOWMN(x, y, missval1, missval2)
238 #define SQRTMN(x) FSQRTMN(x, missval1)
239
FADDMN(const double x,const double y,const double missval1,const double missval2)240 static inline double FADDMN(const double x, const double y, const double missval1, const double missval2)
241 {
242 return MADDMN(x, y);
243 }
FSUBMN(const double x,const double y,const double missval1,const double missval2)244 static inline double FSUBMN(const double x, const double y, const double missval1, const double missval2)
245 {
246 return MSUBMN(x, y);
247 }
FMULMN(const double x,const double y,const double missval1,const double missval2)248 static inline double FMULMN(const double x, const double y, const double missval1, const double missval2)
249 {
250 return MMULMN(x, y);
251 }
FDIVMN(const double x,const double y,const double missval1,const double missval2)252 static inline double FDIVMN(const double x, const double y, const double missval1, const double missval2)
253 {
254 return MDIVMN(x, y);
255 }
FPOWMN(const double x,const double y,const double missval1,const double missval2)256 static inline double FPOWMN(const double x, const double y, const double missval1, const double missval2)
257 {
258 return MPOWMN(x, y);
259 }
FSQRTMN(const double x,const double missval1)260 static inline double FSQRTMN(const double x, const double missval1)
261 {
262 return MSQRTMN(x);
263 }
264 // clang-format on
265
266 const char *fpe_errstr(int fpeRaised);
267
268 template <typename T>
269 MinMaxSum varray_min_max_sum(size_t len, const Varray<T> &v, MinMaxSum mms);
270
271 template <typename T>
272 MinMaxSum varray_min_max_sum_mv(size_t len, const Varray<T> &v, T missval, MinMaxSum mms);
273
274 template <typename T>
275 MinMaxMean varray_min_max_mean(size_t len, const Varray<T> &v);
276
277 template <typename T>
278 MinMaxMean varray_min_max_mean_mv(size_t len, const Varray<T> &v, T missval);
279
280 template <typename T>
281 MinMax array_min_max_mask(size_t len, const T *array, const Varray<int> &mask);
282
283 void array_add_array(size_t len, double *array1, const double *array2);
284 void array_add_array_mv(size_t len, double *array1, const double *array2, double missval);
285
286 template <typename T>
287 void
varray_fill(const size_t len,T * v,const T value)288 varray_fill(const size_t len, T *v, const T value)
289 {
290 for (size_t i = 0; i < len; ++i) v[i] = value;
291 }
292
293 template <typename T>
294 void
varray_fill(Varray<T> & v,const T value)295 varray_fill(Varray<T> &v, const T value)
296 {
297 varray_fill(v.size(), v.data(), value);
298 }
299
300 template <typename T>
301 void
varray_fill(const size_t len,Varray<T> & v,const T value)302 varray_fill(const size_t len, Varray<T> &v, const T value)
303 {
304 varray_fill(len, v.data(), value);
305 }
306
307 template <typename T1, typename T2>
308 void
array_copy(const size_t len,const T1 * array1,T2 * array2)309 array_copy(const size_t len, const T1 *array1, T2 *array2)
310 {
311 for (size_t i = 0; i < len; ++i) array2[i] = array1[i];
312 }
313
314 template <typename T1, typename T2>
315 void
varray_copy(const size_t len,const Varray<T1> & v1,Varray<T2> & v2)316 varray_copy(const size_t len, const Varray<T1> &v1, Varray<T2> &v2)
317 {
318 for (size_t i = 0; i < len; ++i) v2[i] = v1[i];
319 }
320
321 template <typename T, class UnaryOperation>
varray_transform(const Varray<T> & vIn,Varray<T> & vOut,UnaryOperation unary_op)322 void varray_transform(const Varray<T> &vIn, Varray<T> &vOut, UnaryOperation unary_op)
323 {
324 assert(vIn.size() > 0);
325 assert(vOut.size() > 0);
326 assert(vOut.size() <= vIn.size());
327
328 const auto len = vIn.size();
329 for (size_t i = 0; i < len; ++i) vOut[i] = unary_op(vIn[i]);
330 }
331
332 template <typename T>
333 size_t array_num_mv(size_t len, const T *array, T missval);
334
335 template <typename T>
336 size_t varray_num_mv(size_t len, const Varray<T> &v, T missval);
337
338 template <typename T>
339 MinMax varray_min_max(size_t len, const Varray<T> &v);
340
341 template <typename T>
342 MinMax varray_min_max(size_t len, const T *array);
343
344 template <typename T>
345 MinMax varray_min_max(const Varray<T> &v);
346
347 template <typename T>
348 MinMax varray_min_max_mv(size_t len, const Varray<T> &v, T missval);
349
350 template <typename T>
351 MinMax varray_min_max_mv(size_t len, const T *array, T missval);
352
353 template <typename T>
354 T varray_min(size_t len, const Varray<T> &v);
355
356 template <typename T>
357 T varray_max(size_t len, const Varray<T> &v);
358
359 template <typename T>
360 T varray_range(size_t len, const Varray<T> &v);
361
362 template <typename T>
363 T varray_min_mv(size_t len, const Varray<T> &v, T missval);
364
365 template <typename T>
366 T varray_max_mv(size_t len, const Varray<T> &v, T missval);
367
368 template <typename T>
369 T varray_range_mv(size_t len, const Varray<T> &v, T missval);
370
371 template <typename T>
372 double varray_sum(size_t len, const Varray<T> &v);
373
374 template <typename T>
375 double varray_sum_mv(size_t len, const Varray<T> &v, T missval);
376
377 template <typename T>
378 double varray_mean(size_t len, const Varray<T> &v);
379
380 template <typename T>
381 double varray_mean_mv(size_t len, const Varray<T> &v, T missval);
382
383 template <typename T>
384 double varray_weighted_mean(size_t len, const Varray<T> &v, const Varray<double> &w, T missval);
385
386 template <typename T>
387 double varray_weighted_mean_mv(size_t len, const Varray<T> &v, const Varray<double> &w, T missval);
388
389 template <typename T>
390 double varray_avg_mv(size_t len, const Varray<T> &v, T missval);
391
392 template <typename T>
393 double varray_weighted_avg_mv(size_t len, const Varray<T> &v, const Varray<double> &w, T missval);
394
395 template <typename T>
396 double varray_var(size_t len, const Varray<T> &v, size_t nmiss, T missval);
397
398 template <typename T>
399 double varray_var_1(size_t len, const Varray<T> &v, size_t nmiss, T missval);
400
401 template <typename T>
402 double varray_weighted_var(size_t len, const Varray<T> &v, const Varray<double> &w, size_t nmiss, T missval);
403
404 template <typename T>
405 double varray_weighted_var_1(size_t len, const Varray<T> &v, const Varray<double> &w, size_t nmiss, T missval);
406
407 template <typename T>
408 double varray_skew(size_t len, const Varray<T> &v, size_t nmiss, T missval);
409
410 template <typename T>
411 double varray_kurt(size_t len, const Varray<T> &v, size_t nmiss, T missval);
412
413 template <typename T>
414 double varray_median(size_t len, const Varray<T> &v, size_t nmiss, T missval);
415
416 #endif // VARRAY_H
417