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
8 #include <algorithm>
9 #include <cdi.h>
10
11 #include "cdo_options.h"
12 #include "process_int.h"
13 #include "field.h"
14
15 // clang-format off
__anoncf6460190102(auto a, auto b) 16 auto memtype_is_float_float = [](auto a, auto b) { return (a == MemType::Float && b == MemType::Float); };
__anoncf6460190202(auto a, auto b) 17 auto memtype_is_double_float = [](auto a, auto b) { return (a == MemType::Double && b == MemType::Float); };
18 // auto memtype_is_float_double = [](auto a, auto b) { return (a == MemType::Float && b == MemType::Double); };
19 // auto memtype_is_double_double = [](auto a, auto b) { return (a == MemType::Double && b == MemType::Double); };
20
21 // arithmetic
__anoncf6460190302(const double a, const double b) 22 auto arith_add = [](const double a, const double b) { return a + b; };
__anoncf6460190402(const double a, const double b) 23 auto arith_sub = [](const double a, const double b) { return a - b; };
__anoncf6460190502(const double a, const double b) 24 auto arith_mul = [](const double a, const double b) { return a * b; };
__anoncf6460190602(const double a, const double b) 25 auto arith_div = [](const double a, const double b) { return a / b; };
__anoncf6460190702(const auto a, const auto b) 26 auto arith_min = [](const auto a, const auto b) { return (b > a) ? a : b; };
__anoncf6460190802(const auto a, const auto b) 27 auto arith_max = [](const auto a, const auto b) { return (b < a) ? a : b; };
__anoncf6460190902(const double a, const double b) 28 auto arith_sumq = [](const double a, const double b) { return a + b * b; };
29 // auto arith_atan2 = [](const double a, const double b) { return std::atan2(a, b); };
30
31 // arithmetic with missing values
32 auto arith_add_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460190a02(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 33 { a = (isEQ(a, mv_a) || isEQ(b, mv_b)) ? mv_a : arith_add(a, b); };
34 auto arith_sub_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460190b02(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 35 { a = (isEQ(a, mv_a) || isEQ(b, mv_b)) ? mv_a : arith_sub(a, b); };
36 auto arith_mul_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460190c02(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 37 { a = (isEQ(a, 0) || isEQ(b, 0)) ? 0 : ((isEQ(a, mv_a) || isEQ(b, mv_b)) ? mv_a : arith_mul(a, b)); };
38 auto arith_div_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460190d02(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 39 { a = (isEQ(a, mv_a) || isEQ(b, mv_b) || isEQ(b, 0)) ? mv_a : arith_div(a, b); };
40 auto arith_min_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460190e02(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 41 { a = isEQ(b, mv_b) ? a : (isEQ(a, mv_a) ? b : arith_min(a, b)); };
42 auto arith_max_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460190f02(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 43 { a = isEQ(b, mv_b) ? a : (isEQ(a, mv_a) ? b : arith_max(a, b)); };
44 auto arith_sum_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460191002(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 45 { if (!isEQ(b, mv_b)) a = (isEQ(a, mv_a) ? b : arith_add(a, b)); };
46 auto arith_sumq_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460191102(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 47 { if (!isEQ(b, mv_b)) a = (isEQ(a, mv_a) ? arith_mul(b, b) : arith_sumq(a, b)); };
48 /*
49 auto arith_atan2_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
50 { a = (isEQ(a, mv_a) || isEQ(b, mv_b)) ? mv_a : arith_atan2(a, b); };
51 */
52 auto arith_vinit_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460191202(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 53 { a = !isEQ(b, mv_b); (void)mv_a; };
54 auto arith_vincr_mv = [](auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ)
__anoncf6460191302(auto &a, const auto mv_a, const auto b, const auto mv_b, auto isEQ) 55 { if (!isEQ(b, mv_b)) a = a + 1; (void)mv_a; };
56 // clang-format on
57
58 template <typename T1, typename T2>
59 static void
varray2_div(const size_t len,Varray<T1> & v1,const Varray<T2> & v2,const double missval1)60 varray2_div(const size_t len, Varray<T1> &v1, const Varray<T2> &v2, const double missval1)
61 {
62 assert(len > 0);
63 assert(v1.size() > 0);
64 assert(v2.size() > 0);
65 assert(len <= v1.size());
66 assert(len <= v2.size());
67
68 for (size_t i = 0; i < len; ++i) v1[i] = is_equal(v2[i], 0.0) ? missval1 : arith_div(v1[i], v2[i]);
69 }
70
71 template <typename T1, typename T2, typename FUNC>
72 static void
varray2_arith(const size_t len,Varray<T1> & v1,const Varray<T2> & v2,FUNC func)73 varray2_arith(const size_t len, Varray<T1> &v1, const Varray<T2> &v2, FUNC func)
74 {
75 assert(len > 0);
76 assert(v1.size() > 0);
77 assert(v2.size() > 0);
78 assert(len <= v1.size());
79 assert(len <= v2.size());
80
81 for (size_t i = 0; i < len; ++i) v1[i] = func(v1[i], v2[i]);
82 }
83
84 template <typename T1, typename T2, typename FUNC>
85 static void
varray2_arith_mv(const size_t len,Varray<T1> & v1,const Varray<T2> & v2,const double missval1,const double missval2,FUNC func)86 varray2_arith_mv(const size_t len, Varray<T1> &v1, const Varray<T2> &v2, const double missval1, const double missval2, FUNC func)
87 {
88 assert(len > 0);
89 assert(v1.size() > 0);
90 assert(v2.size() > 0);
91 assert(len <= v1.size());
92 assert(len <= v2.size());
93
94 const T1 mv1 = missval1;
95 const T2 mv2 = missval2;
96
97 if (std::isnan(missval1) || std::isnan(missval2))
98 {
99 #ifdef _OPENMP
100 #pragma omp parallel for default(shared)
101 #endif
102 for (size_t i = 0; i < len; ++i) func(v1[i], mv1, v2[i], mv2, dbl_is_equal);
103 }
104 else
105 {
106 #ifdef _OPENMP
107 #pragma omp parallel for default(shared)
108 #endif
109 for (size_t i = 0; i < len; ++i) func(v1[i], mv1, v2[i], mv2, is_equal);
110 }
111 }
112
113 // init valid values
114 void
field2_vinit(Field & field1,const Field & field2)115 field2_vinit(Field &field1, const Field &field2)
116 {
117 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
118
119 if (memtype_is_float_float(field1.memType, field2.memType))
120 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_vinit_mv);
121 else if (memtype_is_double_float(field1.memType, field2.memType))
122 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_vinit_mv);
123 else
124 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_vinit_mv);
125
126 field1.nmiss = field2.nmiss;
127 }
128
129 // increment valid values
130 void
field2_vincr(Field & field1,const Field & field2)131 field2_vincr(Field &field1, const Field &field2)
132 {
133 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
134
135 if (memtype_is_float_float(field1.memType, field2.memType))
136 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_vincr_mv);
137 else if (memtype_is_double_float(field1.memType, field2.memType))
138 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_vincr_mv);
139 else
140 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_vincr_mv);
141
142 field1.nmiss = field2.nmiss;
143 }
144
145 void
field2_add(Field & field1,const Field & field2)146 field2_add(Field &field1, const Field &field2)
147 {
148 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
149
150 if (field1.nmiss || field2.nmiss)
151 {
152 if (memtype_is_float_float(field1.memType, field2.memType))
153 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_add_mv);
154 else if (memtype_is_double_float(field1.memType, field2.memType))
155 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_add_mv);
156 else
157 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_add_mv);
158
159 field_num_mv(field1);
160 }
161 else
162 {
163 if (memtype_is_float_float(field1.memType, field2.memType))
164 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_add);
165 else if (memtype_is_double_float(field1.memType, field2.memType))
166 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_add);
167 else
168 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_add);
169 }
170 }
171
172 void
field2_sum(Field & field1,const Field & field2)173 field2_sum(Field &field1, const Field &field2)
174 {
175 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
176
177 if (field1.nmiss || field2.nmiss)
178 {
179 if (memtype_is_float_float(field1.memType, field2.memType))
180 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_sum_mv);
181 else if (memtype_is_double_float(field1.memType, field2.memType))
182 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_sum_mv);
183 else
184 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_sum_mv);
185
186 field_num_mv(field1);
187 }
188 else
189 {
190 if (memtype_is_float_float(field1.memType, field2.memType))
191 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_add);
192 else if (memtype_is_double_float(field1.memType, field2.memType))
193 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_add);
194 else
195 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_add);
196 }
197 }
198
199 void
field2_sumw(Field & field1,const Field & field2,double w)200 field2_sumw(Field &field1, const Field &field2, double w)
201 {
202 const auto missval1 = field1.missval;
203 const auto missval2 = field2.missval;
204 auto &array1 = field1.vec_d;
205 const auto &array2 = field2.vec_d;
206
207 auto len = field1.size;
208 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
209
210 if (field1.nmiss || field2.nmiss)
211 {
212 for (size_t i = 0; i < len; i++)
213 if (!DBL_IS_EQUAL(array2[i], missval2))
214 {
215 if (!DBL_IS_EQUAL(array1[i], missval1))
216 array1[i] += w * array2[i];
217 else
218 array1[i] = w * array2[i];
219 }
220
221 field1.nmiss = field_num_miss(field1);
222 }
223 else
224 {
225 #ifdef _OPENMP
226 #pragma omp parallel for default(none) shared(len, array1, array2, w)
227 #endif
228 for (size_t i = 0; i < len; i++) array1[i] += w * array2[i];
229 }
230 }
231
232 /*
233 * Compute the occurrence of values in field, if they do not equal refval.
234 * This can be used to compute the lengths of multiple periods in a timeseries.
235 * Missing field values are handled like refval, i.e. they stop a running period.
236 * If there is missing data in the occurence field, missing fields values do not
237 * change anything (they do not start a non-period by setting occurrence to zero).
238 */
239 void
field2_sumtr(Field & occur,const Field & field,const double refval)240 field2_sumtr(Field &occur, const Field &field, const double refval)
241 {
242 auto omissval = occur.missval;
243 auto fmissval = field.missval;
244 auto &oarray = occur.vec_d;
245 const auto &farray = field.vec_d;
246
247 auto len = occur.size;
248 if (len != field.size) cdo_abort("Fields have different size (%s)", __func__);
249
250 if (occur.nmiss || field.nmiss)
251 {
252 #ifdef _OPENMP
253 #pragma omp parallel for default(shared) schedule(static)
254 #endif
255 for (size_t i = 0; i < len; i++)
256 if (!DBL_IS_EQUAL(farray[i], fmissval))
257 {
258 if (!DBL_IS_EQUAL(oarray[i], omissval))
259 oarray[i] = (DBL_IS_EQUAL(farray[i], refval)) ? 0.0 : oarray[i] + 1.0;
260 else
261 oarray[i] = (DBL_IS_EQUAL(farray[i], refval)) ? 0.0 : 1.0;
262 }
263 else
264 {
265 if (!DBL_IS_EQUAL(oarray[i], omissval)) oarray[i] = 0.0;
266 }
267
268 occur.nmiss = field_num_miss(occur);
269 }
270 else
271 {
272 #ifdef _OPENMP
273 #pragma omp parallel for default(shared)
274 #endif
275 for (size_t i = 0; i < len; i++) oarray[i] = (DBL_IS_EQUAL(farray[i], refval)) ? 0.0 : oarray[i] + 1.0;
276 }
277 }
278
279 void
field2_sumq(Field & field1,const Field & field2)280 field2_sumq(Field &field1, const Field &field2)
281 {
282 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
283
284 if (field1.nmiss || field2.nmiss)
285 {
286 if (memtype_is_float_float(field1.memType, field2.memType))
287 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_sumq_mv);
288 else if (memtype_is_double_float(field1.memType, field2.memType))
289 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_sumq_mv);
290 else
291 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_sumq_mv);
292
293 field_num_mv(field1);
294 }
295 else
296 {
297 if (memtype_is_float_float(field1.memType, field2.memType))
298 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_sumq);
299 else if (memtype_is_double_float(field1.memType, field2.memType))
300 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_sumq);
301 else
302 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_sumq);
303 }
304 }
305
306 void
field2_sumsumq(Field & field1,Field & field2,const Field & field3)307 field2_sumsumq(Field &field1, Field &field2, const Field &field3)
308 {
309 field2_sumq(field2, field3);
310 field2_sum(field1, field3);
311 }
312
313 void
field2_sumqw(Field & field1,const Field & field2,const double w)314 field2_sumqw(Field &field1, const Field &field2, const double w)
315 {
316 const auto missval1 = field1.missval;
317 const auto missval2 = field2.missval;
318 auto &array1 = field1.vec_d;
319 const auto &array2 = field2.vec_d;
320
321 const auto len = field1.size;
322 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
323
324 if (field1.nmiss || field2.nmiss)
325 {
326 for (size_t i = 0; i < len; i++)
327 if (!DBL_IS_EQUAL(array2[i], missval2))
328 {
329 if (!DBL_IS_EQUAL(array1[i], missval1))
330 array1[i] += w * array2[i] * array2[i];
331 else
332 array1[i] = w * array2[i] * array2[i];
333 }
334
335 field1.nmiss = field_num_miss(field1);
336 }
337 else
338 {
339 for (size_t i = 0; i < len; i++) array1[i] += w * array2[i] * array2[i];
340 }
341 }
342
343 void
field2_sub(Field & field1,const Field & field2)344 field2_sub(Field &field1, const Field &field2)
345 {
346 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
347
348 if (field1.nmiss || field2.nmiss)
349 {
350 if (memtype_is_float_float(field1.memType, field2.memType))
351 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_sub_mv);
352 else if (memtype_is_double_float(field1.memType, field2.memType))
353 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_sub_mv);
354 else
355 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_sub_mv);
356
357 field_num_mv(field1);
358 }
359 else
360 {
361 if (memtype_is_float_float(field1.memType, field2.memType))
362 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_sub);
363 else if (memtype_is_double_float(field1.memType, field2.memType))
364 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_sub);
365 else
366 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_sub);
367 }
368 }
369
370 void
field2_mul(Field & field1,const Field & field2)371 field2_mul(Field &field1, const Field &field2)
372 {
373 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
374
375 if (field1.nmiss || field2.nmiss)
376 {
377 if (memtype_is_float_float(field1.memType, field2.memType))
378 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_mul_mv);
379 else if (memtype_is_double_float(field1.memType, field2.memType))
380 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_mul_mv);
381 else
382 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_mul_mv);
383
384 field_num_mv(field1);
385 }
386 else
387 {
388 if (memtype_is_float_float(field1.memType, field2.memType))
389 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_mul);
390 else if (memtype_is_double_float(field1.memType, field2.memType))
391 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_mul);
392 else
393 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_mul);
394 }
395 }
396
397 void
field2_div(Field & field1,const Field & field2)398 field2_div(Field &field1, const Field &field2)
399 {
400 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
401
402 if (field1.nmiss || field2.nmiss)
403 {
404 if (memtype_is_float_float(field1.memType, field2.memType))
405 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_div_mv);
406 else if (memtype_is_double_float(field1.memType, field2.memType))
407 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_div_mv);
408 else
409 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_div_mv);
410
411 field_num_mv(field1);
412 }
413 else
414 {
415 if (memtype_is_float_float(field1.memType, field2.memType))
416 varray2_div(field1.size, field1.vec_f, field2.vec_f, field1.missval);
417 else if (memtype_is_double_float(field1.memType, field2.memType))
418 varray2_div(field1.size, field1.vec_d, field2.vec_f, field1.missval);
419 else
420 varray2_div(field1.size, field1.vec_d, field2.vec_d, field1.missval);
421
422 field_num_mv(field1);
423 }
424 }
425
426 void
field2_atan2(Field & field1,const Field & field2)427 field2_atan2(Field &field1, const Field &field2)
428 {
429 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
430
431 if (memtype_is_float_float(field1.memType, field2.memType))
432 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_div_mv);
433 else if (memtype_is_double_float(field1.memType, field2.memType))
434 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_div_mv);
435 else
436 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_div_mv);
437
438 field_num_mv(field1);
439 }
440
441 void
field2_set_miss(Field & field1,const Field & field2)442 field2_set_miss(Field &field1, const Field &field2)
443 {
444 const auto missval1 = field1.missval;
445 auto &array1 = field1.vec_d;
446 const auto &array2 = field2.vec_d;
447
448 const auto len = field1.size;
449 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
450
451 if (field1.nmiss)
452 {
453 for (size_t i = 0; i < len; i++) array1[i] = DBL_IS_EQUAL(array1[i], missval1) ? array2[i] : array1[i];
454
455 field1.nmiss = field_num_miss(field1);
456 }
457 }
458
459 void
field2_min(Field & field1,const Field & field2)460 field2_min(Field &field1, const Field &field2)
461 {
462 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
463
464 if (field1.nmiss || field2.nmiss)
465 {
466 if (memtype_is_float_float(field1.memType, field2.memType))
467 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_min_mv);
468 else if (memtype_is_double_float(field1.memType, field2.memType))
469 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_min_mv);
470 else
471 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_min_mv);
472
473 field_num_mv(field1);
474 }
475 else
476 {
477 if (memtype_is_float_float(field1.memType, field2.memType))
478 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_min);
479 else if (memtype_is_double_float(field1.memType, field2.memType))
480 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_min);
481 else
482 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_min);
483 }
484 }
485
486 void
field2_max(Field & field1,const Field & field2)487 field2_max(Field &field1, const Field &field2)
488 {
489 if (field1.size != field2.size) cdo_abort("Fields have different size (%s)", __func__);
490
491 if (field1.nmiss || field2.nmiss)
492 {
493 if (memtype_is_float_float(field1.memType, field2.memType))
494 varray2_arith_mv(field1.size, field1.vec_f, field2.vec_f, field1.missval, field2.missval, arith_max_mv);
495 else if (memtype_is_double_float(field1.memType, field2.memType))
496 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_f, field1.missval, field2.missval, arith_max_mv);
497 else
498 varray2_arith_mv(field1.size, field1.vec_d, field2.vec_d, field1.missval, field2.missval, arith_max_mv);
499
500 field_num_mv(field1);
501 }
502 else
503 {
504 if (memtype_is_float_float(field1.memType, field2.memType))
505 varray2_arith(field1.size, field1.vec_f, field2.vec_f, arith_max);
506 else if (memtype_is_double_float(field1.memType, field2.memType))
507 varray2_arith(field1.size, field1.vec_d, field2.vec_f, arith_max);
508 else
509 varray2_arith(field1.size, field1.vec_d, field2.vec_d, arith_max);
510 }
511 }
512
513 void
field2_maxmin(Field & field1,Field & field2,const Field & field3)514 field2_maxmin(Field &field1, Field &field2, const Field &field3)
515 {
516 field2_min(field2, field3);
517 field2_max(field1, field3);
518 }
519
__anoncf6460191402(auto &v1, auto &v2, const auto v3, const auto idx) 520 auto field2_set_index = [](auto &v1, auto &v2, const auto v3, const auto idx) {
521 v2 = v3;
522 v1 = idx;
523 };
524
525 template <typename T>
526 void
field2_minidx(size_t nmiss3,size_t len,double missval3,Field & field1,Field & field2,const Varray<T> & v3,int idx)527 field2_minidx(size_t nmiss3, size_t len, double missval3, Field &field1, Field &field2, const Varray<T> &v3, int idx)
528 {
529 const auto missval1 = field1.missval;
530 const auto missval2 = field2.missval;
531 auto &v1 = field1.vec_d;
532 auto &v2 = field2.vec_d;
533
534 if (field2.nmiss || nmiss3)
535 {
536 for (size_t i = 0; i < len; i++)
537 {
538 if (DBL_IS_EQUAL(v3[i], missval3))
539 {
540 if (DBL_IS_EQUAL(v2[i], missval2)) v1[i] = missval1;
541 }
542 else if (v3[i] < v2[i] || DBL_IS_EQUAL(v2[i], missval2))
543 {
544 field2_set_index(v1[i], v2[i], v3[i], idx);
545 }
546 }
547
548 field_num_mv(field1);
549 field_num_mv(field2);
550 }
551 else
552 {
553 for (size_t i = 0; i < len; i++)
554 {
555 if (v3[i] < v2[i]) field2_set_index(v1[i], v2[i], v3[i], idx);
556 }
557 }
558 }
559
560 void
field2_minidx(Field & field1,Field & field2,const Field & field3,int idx)561 field2_minidx(Field &field1, Field &field2, const Field &field3, int idx)
562 {
563 if (field1.size != field3.size) cdo_abort("Fields have different size (%s)", __func__);
564 if (field2.size != field3.size) cdo_abort("Fields have different size (%s)", __func__);
565
566 if (field3.memType == MemType::Float)
567 field2_minidx(field3.nmiss, field3.size, field3.missval, field1, field2, field3.vec_f, idx);
568 else
569 field2_minidx(field3.nmiss, field3.size, field3.missval, field1, field2, field3.vec_d, idx);
570 }
571
572 template <typename T>
573 void
field2_maxidx(size_t nmiss3,size_t len,double missval3,Field & field1,Field & field2,const Varray<T> & v3,int idx)574 field2_maxidx(size_t nmiss3, size_t len, double missval3, Field &field1, Field &field2, const Varray<T> &v3, int idx)
575 {
576 const auto missval1 = field1.missval;
577 const auto missval2 = field2.missval;
578 auto &v1 = field1.vec_d;
579 auto &v2 = field2.vec_d;
580
581 if (field2.nmiss || nmiss3)
582 {
583 for (size_t i = 0; i < len; i++)
584 {
585 if (DBL_IS_EQUAL(v3[i], missval3))
586 {
587 if (DBL_IS_EQUAL(v2[i], missval2)) v1[i] = missval1;
588 }
589 else if (v3[i] > v2[i] || DBL_IS_EQUAL(v2[i], missval2))
590 {
591 field2_set_index(v1[i], v2[i], v3[i], idx);
592 }
593 }
594
595 field_num_mv(field1);
596 field_num_mv(field2);
597 }
598 else
599 {
600 for (size_t i = 0; i < len; i++)
601 {
602 if (v3[i] > v2[i]) field2_set_index(v1[i], v2[i], v3[i], idx);
603 }
604 }
605 }
606
607 void
field2_maxidx(Field & field1,Field & field2,const Field & field3,int idx)608 field2_maxidx(Field &field1, Field &field2, const Field &field3, int idx)
609 {
610 if (field1.size != field3.size) cdo_abort("Fields have different size (%s)", __func__);
611 if (field2.size != field3.size) cdo_abort("Fields have different size (%s)", __func__);
612
613 if (field3.memType == MemType::Float)
614 field2_maxidx(field3.nmiss, field3.size, field3.missval, field1, field2, field3.vec_f, idx);
615 else
616 field2_maxidx(field3.nmiss, field3.size, field3.missval, field1, field2, field3.vec_d, idx);
617 }
618
619 static size_t
field_set_nmiss(const size_t len,Varray<double> & v,double missval)620 field_set_nmiss(const size_t len, Varray<double> &v, double missval)
621 {
622 size_t nmiss = 0;
623
624 if (std::isnan(missval))
625 {
626 for (size_t i = 0; i < len; i++)
627 if (DBL_IS_EQUAL(v[i], missval) || v[i] < 0)
628 {
629 v[i] = missval;
630 nmiss++;
631 }
632 }
633 else
634 {
635 for (size_t i = 0; i < len; i++)
636 if (IS_EQUAL(v[i], missval) || v[i] < 0)
637 {
638 v[i] = missval;
639 nmiss++;
640 }
641 }
642
643 return nmiss;
644 }
645
646 void
field2_var(Field & field1,const Field & field2,const Field & field3,const int divisor)647 field2_var(Field &field1, const Field &field2, const Field &field3, const int divisor)
648 {
649 const auto missval1 = field1.missval;
650 const auto missval2 = field2.missval;
651 auto &array1 = field1.vec_d;
652 const auto &array2 = field2.vec_d;
653 const auto &array3 = field3.vec_d;
654
655 const auto len = field1.size;
656 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
657
658 if (field1.nmiss || field2.nmiss)
659 {
660 for (size_t i = 0; i < len; i++)
661 {
662 const auto temp = DIVMN(MULMN(array1[i], array1[i]), array3[i]);
663 array1[i] = DIVMN(SUBMN(array2[i], temp), array3[i] - divisor);
664 if (array1[i] < 0 && array1[i] > -1.e-5) array1[i] = 0;
665 }
666 }
667 else
668 {
669 for (size_t i = 0; i < len; i++)
670 {
671 const auto temp = DIV(MUL(array1[i], array1[i]), array3[i]);
672 array1[i] = DIV(SUB(array2[i], temp), array3[i] - divisor);
673 if (array1[i] < 0 && array1[i] > -1.e-5) array1[i] = 0;
674 }
675 }
676
677 field1.nmiss = field_set_nmiss(len, array1, missval1);
678 }
679
680 void
field2_std(Field & field1,const Field & field2,const Field & field3,const int divisor)681 field2_std(Field &field1, const Field &field2, const Field &field3, const int divisor)
682 {
683 const auto missval1 = field1.missval;
684 auto &array1 = field1.vec_d;
685
686 const auto len = field1.size;
687 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
688
689 field2_var(field1, field2, field3, divisor);
690
691 size_t nmiss = 0;
692 for (size_t i = 0; i < len; i++)
693 if (DBL_IS_EQUAL(array1[i], missval1) || array1[i] < 0)
694 {
695 array1[i] = missval1;
696 nmiss++;
697 }
698 else
699 {
700 array1[i] = IS_NOT_EQUAL(array1[i], 0) ? std::sqrt(array1[i]) : 0;
701 }
702 field1.nmiss = nmiss;
703 }
704
705 void
fieldc_var(Field & field1,const Field & field2,const int nsets,const int divisor)706 fieldc_var(Field &field1, const Field &field2, const int nsets, const int divisor)
707 {
708 const auto nsetx = nsets - divisor;
709 const auto missval1 = field1.missval;
710 const auto missval2 = field2.missval;
711 auto &array1 = field1.vec_d;
712 const auto &array2 = field2.vec_d;
713 double temp;
714
715 const auto len = field1.size;
716 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
717
718 if (nsetx == 0)
719 {
720 for (size_t i = 0; i < len; i++) array1[i] = missval1;
721 }
722 else if (field1.nmiss || field2.nmiss)
723 {
724 for (size_t i = 0; i < len; i++)
725 {
726 temp = MULMN(array1[i], array1[i]) / nsets;
727 array1[i] = SUBMN(array2[i], temp) / nsetx;
728 if (array1[i] < 0 && array1[i] > -1.e-5) array1[i] = 0;
729 }
730 }
731 else
732 {
733 for (size_t i = 0; i < len; i++)
734 {
735 temp = MUL(array1[i], array1[i]) / nsets;
736 array1[i] = SUB(array2[i], temp) / nsetx;
737 if (array1[i] < 0 && array1[i] > -1.e-5) array1[i] = 0;
738 }
739 }
740
741 field1.nmiss = field_set_nmiss(len, array1, missval1);
742 }
743
744 void
fieldc_std(Field & field1,const Field & field2,const int nsets,const int divisor)745 fieldc_std(Field &field1, const Field &field2, const int nsets, const int divisor)
746 {
747 const auto missval1 = field1.missval;
748 auto &array1 = field1.vec_d;
749
750 const auto len = field1.size;
751 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
752
753 fieldc_var(field1, field2, nsets, divisor);
754
755 size_t nmiss = 0;
756 for (size_t i = 0; i < len; i++)
757 if (DBL_IS_EQUAL(array1[i], missval1) || array1[i] < 0)
758 {
759 array1[i] = missval1;
760 nmiss++;
761 }
762 else
763 {
764 array1[i] = IS_NOT_EQUAL(array1[i], 0) ? std::sqrt(array1[i]) : 0;
765 }
766
767 field1.nmiss = nmiss;
768 }
769
770 void
field2_moq(Field & field1,const Field & field2)771 field2_moq(Field &field1, const Field &field2)
772 {
773 const auto missval1 = field1.missval;
774 const auto missval2 = field2.missval;
775 auto &array1 = field1.vec_d;
776 const auto &array2 = field2.vec_d;
777
778 const auto len = field1.size;
779 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
780
781 if (field2.nmiss)
782 {
783 for (size_t i = 0; i < len; i++)
784 if (!DBL_IS_EQUAL(array2[i], missval2))
785 array1[i] = array2[i] * array2[i];
786 else
787 array1[i] = missval1;
788
789 field1.nmiss = field_num_miss(field1);
790 }
791 else
792 {
793 for (size_t i = 0; i < len; i++) array1[i] = array2[i] * array2[i];
794 }
795 }
796
797 void
field2_moqw(Field & field1,const Field & field2,const double w)798 field2_moqw(Field &field1, const Field &field2, const double w)
799 {
800 const auto missval1 = field1.missval;
801 const auto missval2 = field2.missval;
802 auto &array1 = field1.vec_d;
803 const auto &array2 = field2.vec_d;
804
805 const auto len = field1.size;
806 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
807
808 if (field2.nmiss)
809 {
810 for (size_t i = 0; i < len; i++)
811 if (!DBL_IS_EQUAL(array2[i], missval2))
812 array1[i] = w * array2[i] * array2[i];
813 else
814 array1[i] = missval1;
815
816 field1.nmiss = field_num_miss(field1);
817 }
818 else
819 {
820 for (size_t i = 0; i < len; i++) array1[i] = w * array2[i] * array2[i];
821 }
822 }
823
824 /**
825 * Counts the number of nonmissing values. The result of the operation
826 * is computed according to the following rules:
827 *
828 * field1 field2 result
829 * a b a + 1
830 * a miss a
831 * miss b 1
832 * miss miss miss
833 *
834 * @param field1 the 1st input field, also holds the result
835 * @param field2 the 2nd input field
836 */
837 void
field2_count(Field & field1,const Field & field2)838 field2_count(Field &field1, const Field &field2)
839 {
840 const auto missval1 = field1.missval;
841 const auto missval2 = field2.missval;
842 auto &array1 = field1.vec_d;
843 const auto &array2 = field2.vec_d;
844
845 const auto len = field1.size;
846 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
847
848 if (field1.nmiss || field2.nmiss)
849 {
850 for (size_t i = 0; i < len; i++)
851 if (!DBL_IS_EQUAL(array2[i], missval2))
852 {
853 if (!DBL_IS_EQUAL(array1[i], missval1))
854 array1[i] += 1.0;
855 else
856 array1[i] = 1.0;
857 }
858
859 field1.nmiss = field_num_miss(field1);
860 }
861 else
862 {
863 for (size_t i = 0; i < len; i++) array1[i] += 1.0;
864 }
865 }
866
867 void
field2_function(Field & field1,const Field & field2,int function)868 field2_function(Field &field1, const Field &field2, int function)
869 {
870 // clang-format off
871 switch (function)
872 {
873 case FieldFunc_Add: field2_add(field1, field2); break;
874 case FieldFunc_Min: field2_min(field1, field2); break;
875 case FieldFunc_Max: field2_max(field1, field2); break;
876 case FieldFunc_Sum: field2_sum(field1, field2); break;
877 case FieldFunc_Mean: field2_sum(field1, field2); break;
878 case FieldFunc_Avg: field2_add(field1, field2); break;
879 case FieldFunc_Sub: field2_sub(field1, field2); break;
880 case FieldFunc_Mul: field2_mul(field1, field2); break;
881 case FieldFunc_Div: field2_div(field1, field2); break;
882 case FieldFunc_Atan2: field2_atan2(field1, field2); break;
883 case FieldFunc_Setmiss: field2_set_miss(field1, field2); break;
884 default: cdo_abort("%s: function %d not implemented!", __func__, function);
885 }
886 // clang-format on
887 }
888