array.cc 19.4 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
17

Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
#ifdef HAVE_CONFIG_H
Uwe Schulzweida's avatar
Uwe Schulzweida committed
19
#include "config.h"  // restrict
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
21
#endif

22
23
#include <float.h>
#include <fenv.h>
24
#include <assert.h>
25

26
#include "compare.h"
27
#include "array.h"
28
#include "cimdOmp.h"
29
30

//#pragma STDC FENV_ACCESS ON
31

32
33
const char *
fpe_errstr(int fpeRaised)
34
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
  const char *errstr = nullptr;
36

37
38
39
40
41
42
43
  // clang-format off
  if      (fpeRaised & FE_DIVBYZERO) errstr = "division by zero";
  else if (fpeRaised & FE_INEXACT)   errstr = "inexact result";
  else if (fpeRaised & FE_INVALID)   errstr = "invalid result";
  else if (fpeRaised & FE_OVERFLOW)  errstr = "overflow";
  else if (fpeRaised & FE_UNDERFLOW) errstr = "underflow";
  // clang-format on
44
45
46
47

  return errstr;
}

48
MinMaxVal
49
varrayMinMaxMV(const size_t len, const double *array, const double missval, size_t &nvals)
50
{
51
52
  double vmin = DBL_MAX;
  double vmax = -DBL_MAX;
53

54
  nvals = 0;
55
  for (size_t i = 0; i < len; ++i)
56
    {
57
      if (!DBL_IS_EQUAL(array[i], missval))
58
        {
59
60
          if (array[i] < vmin) vmin = array[i];
          if (array[i] > vmax) vmax = array[i];
61
62
63
          nvals++;
        }
    }
64

65
  return MinMaxVal(vmin, vmax);
66
}
67

68
MinMaxVal
69
70
71
72
varrayMinMaxMV(const size_t len, const double *array, const double missval)
{
  size_t nvals;
  return varrayMinMaxMV(len, array, missval, nvals);
73
74
}

75
MinMaxVal
76
77
78
79
80
81
varrayMinMaxMV(const size_t len, const Varray<double> &array, const double missval)
{
  size_t nvals;
  return varrayMinMaxMV(len, array.data(), missval, nvals);
}

82
template <typename T>
83
void
84
varrayMinMaxSum(const size_t len, const Varray<T> &array, double &rmin, double &rmax, double &rsum)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
{
86
  // rmin, rmax and rsum will be initialized in Info
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87

88
89
90
  /* 
  gnu 9.2.0: correct result
  intel 18.0.5: wrong result
91
#ifdef HAVE_OPENMP4
92
93
#pragma omp simd reduction(+:rsum) reduction(min:rmin) reduction(max:rmax)
#endif
94
  */
95
  for (size_t i = 0; i < len; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
    {
97
98
99
100
      const auto v = array[i];
      if (v < rmin) rmin = v;
      if (v > rmax) rmax = v;
      rsum += v;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
    }
102
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103

104
105
106
// Explicit instantiation
template void varrayMinMaxSum(const size_t len, const Varray<float> &array, double &rmin, double &rmax, double &rsum);
template void varrayMinMaxSum(const size_t len, const Varray<double> &array, double &rmin, double &rmax, double &rsum);
107

108
template <typename T>
109
size_t
110
varrayMinMaxSumMV(const size_t len, const Varray<T> &array, const T missval, double &rmin, double &rmax, double &rsum)
111
{
112
  // rmin, rmax and rsum will be initialized in Info
113
114

  size_t nvals = 0;
115
  for (size_t i = 0; i < len; ++i)
116
    {
117
118
      const auto v = array[i];
      if (!DBL_IS_EQUAL(v, missval))
119
        {
120
121
122
          if (v < rmin) rmin = v;
          if (v > rmax) rmax = v;
          rsum += v;
123
124
125
          nvals++;
        }
    }
126

Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
  if (nvals == 0 && IS_EQUAL(rmin, DBL_MAX)) rmin = missval;
128
  if (nvals == 0 && IS_EQUAL(rmax, -DBL_MAX)) rmax = missval;
129
130
131
132

  return nvals;
}

133
// Explicit instantiation
134
template size_t varrayMinMaxSumMV(const size_t len, const Varray<float> &array, const float missval, double &rmin, double &rmax, double &rsum);
135
template size_t varrayMinMaxSumMV(const size_t len, const Varray<double> &array, const double missval, double &rmin, double &rmax, double &rsum);
136

137
void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
varrayMinMaxMean(const size_t len, const Varray<double> &array, double &rmin, double &rmax, double &rmean)
139
{
140
141
142
  rmin = DBL_MAX;
  rmax = -DBL_MAX;
  double sum = 0.0;
143
  varrayMinMaxSum(len, array, rmin, rmax, sum);
144

145
  rmean = len ? sum / (double) len : 0;
146
147
}

148
size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
varrayMinMaxMeanMV(const size_t len, const Varray<double> &array, const double missval, double &rmin, double &rmax, double &rmean)
150
{
151
152
153
  rmin = DBL_MAX;
  rmax = -DBL_MAX;
  double sum = 0.0;
154
  size_t nvals = varrayMinMaxSumMV(len, array, missval, rmin, rmax, sum);
155

156
  rmean = nvals ? sum / (double) nvals : missval;
157
158

  return nvals;
159
160
}

161
void
162
arrayMinMaxMask(const size_t len, const double *array, const Varray<int> &mask, double &rmin, double &rmax)
163
{
164
165
  double zmin = DBL_MAX;
  double zmax = -DBL_MAX;
166

167
  if (!mask.empty())
168
    {
169
170
171
172
      for (size_t i = 0; i < len; ++i)
        {
          if (!mask[i])
            {
173
174
              if (array[i] < zmin) zmin = array[i];
              if (array[i] > zmax) zmax = array[i];
175
176
            }
        }
177
178
179
    }
  else
    {
180
181
      for (size_t i = 0; i < len; ++i)
        {
182
183
          if (array[i] < zmin) zmin = array[i];
          if (array[i] > zmax) zmax = array[i];
184
        }
185
186
    }

187
188
  rmin = zmin;
  rmax = zmax;
189
190
}

191
void
192
arrayAddArray(const size_t len, double *restrict array1, const double *restrict array2)
193
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
  //#ifdef  _OPENMP
195
196
  //#pragma omp parallel for default(none) shared(array1,array2)
  //#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
  for (size_t i = 0; i < len; ++i) array1[i] += array2[i];
198
}
199

200
void
201
arrayAddArrayMV(const size_t len, double *restrict array1, const double *restrict array2, const double missval)
202
{
203
  if (std::isnan(missval))
204
    {
205
206
207
208
209
210
211
212
      for (size_t i = 0; i < len; i++)
        if (!DBL_IS_EQUAL(array2[i], missval))
          {
            if (!DBL_IS_EQUAL(array1[i], missval))
              array1[i] += array2[i];
            else
              array1[i] = array2[i];
          }
213
214
215
    }
  else
    {
216
217
218
219
220
221
222
223
      for (size_t i = 0; i < len; i++)
        if (IS_NOT_EQUAL(array2[i], missval))
          {
            if (IS_NOT_EQUAL(array1[i], missval))
              array1[i] += array2[i];
            else
              array1[i] = array2[i];
          }
224
225
226
    }
}

227
size_t
228
arrayNumMV(const size_t len, const double *restrict array, const double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
230
231
{
  size_t nmiss = 0;

232
  if (std::isnan(missval))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
    {
234
235
      for (size_t i = 0; i < len; ++i)
        if (DBL_IS_EQUAL(array[i], missval)) nmiss++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
237
238
    }
  else
    {
239
240
      for (size_t i = 0; i < len; ++i)
        if (IS_EQUAL(array[i], missval)) nmiss++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
243
244
245
    }

  return nmiss;
}

246
template <typename T>
247
size_t
248
varrayNumMV(const size_t len, const Varray<T> &array, const T missval)
249
250
251
{
  size_t nmiss = 0;

252
  if (std::isnan(missval))
253
254
255
256
257
258
259
260
261
262
263
264
265
    {
      for (size_t i = 0; i < len; ++i)
        if (DBL_IS_EQUAL(array[i], missval)) nmiss++;
    }
  else
    {
      for (size_t i = 0; i < len; ++i)
        if (IS_EQUAL(array[i], missval)) nmiss++;
    }

  return nmiss;
}

266
// Explicit instantiation
267
template size_t varrayNumMV(const size_t len, const Varray<float> &array, const float missval);
268
269
template size_t varrayNumMV(const size_t len, const Varray<double> &array, const double missval);

270
MinMaxVal
Uwe Schulzweida's avatar
Uwe Schulzweida committed
271
272
273
274
275
varrayMinMax(const size_t len, const double *restrict array)
{
  double vmin = DBL_MAX;
  double vmax = -DBL_MAX;

276
277
278
#ifdef HAVE_OPENMP4
#pragma omp simd reduction(min:vmin) reduction(max:vmax)
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
280
281
282
283
284
  for (size_t i = 0; i < len; ++i)
    {
      if (array[i] < vmin) vmin = array[i];
      if (array[i] > vmax) vmax = array[i];
    }

285
  return MinMaxVal(vmin, vmax);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
287
}

288
MinMaxVal
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
291
292
293
varrayMinMax(const size_t len, const Varray<double> &v)
{
  return varrayMinMax(len, v.data());
}

294
template <typename T>
295
MinMaxVal
296
varrayMinMax(const Varray<T> &v)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
297
{
298
299
300
301
302
303
304
305
306
307
308
309
310
  double vmin = DBL_MAX;
  double vmax = -DBL_MAX;

  size_t len = v.size();
#ifdef HAVE_OPENMP4
#pragma omp simd reduction(min:vmin) reduction(max:vmax)
#endif
  for (size_t i = 0; i < len; ++i)
    {
      if (v[i] < vmin) vmin = v[i];
      if (v[i] > vmax) vmax = v[i];
    }

311
  return MinMaxVal(vmin, vmax);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
312
313
}

314
// Explicit instantiation
315
316
template MinMaxVal varrayMinMax(const Varray<float> &v);
template MinMaxVal varrayMinMax(const Varray<double> &v);
317

318
double
319
varrayMin(const size_t len, const Varray<double> &v)
320
{
321
  assert(v.size() >= len);
322

323
  double min = v[0];
324

325
  for (size_t i = 0; i < len; ++i)
326
    if (v[i] < min) min = v[i];
327
328
329
330

  return min;
}

331
double
332
varrayMax(const size_t len, const Varray<double> &v)
333
{
334
  assert(v.size() >= len);
335

336
  double max = v[0];
337

338
  for (size_t i = 0; i < len; ++i)
339
    if (v[i] > max) max = v[i];
340
341
342
343

  return max;
}

344
double
345
varrayRange(const size_t len, const Varray<double> &v)
346
{
347
  assert(v.size() >= len);
348

349
350
  double min = v[0];
  double max = v[0];
351

352
  for (size_t i = 0; i < len; ++i)
353
    {
354
355
      if (v[i] < min) min = v[i];
      if (v[i] > max) max = v[i];
356
357
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
358
  return (max - min);
359
360
}

361
double
362
varrayMinMV(const size_t len, const Varray<double> &v, const double missval)
363
{
364
  assert(v.size() >= len);
365
366
367

  double min = DBL_MAX;

368
  for (size_t i = 0; i < len; ++i)
369
370
    if (!DBL_IS_EQUAL(v[i], missval))
      if (v[i] < min) min = v[i];
371

372
  if (IS_EQUAL(min, DBL_MAX)) min = missval;
373
374
375
376

  return min;
}

377
double
378
varrayMaxMV(const size_t len, const Varray<double> &v, const double missval)
379
{
380
  assert(v.size() >= len);
381
382
383

  double max = -DBL_MAX;

384
  for (size_t i = 0; i < len; ++i)
385
386
    if (!DBL_IS_EQUAL(v[i], missval))
      if (v[i] > max) max = v[i];
387

388
  if (IS_EQUAL(max, -DBL_MAX)) max = missval;
389
390
391
392

  return max;
}

393
double
394
varrayRangeMV(const size_t len, const Varray<double> &v, const double missval)
395
{
396
  assert(v.size() >= len);
397

398
  double min = DBL_MAX;
399
400
  double max = -DBL_MAX;

401
  for (size_t i = 0; i < len; ++i)
402
    {
403
      if (!DBL_IS_EQUAL(v[i], missval))
404
        {
405
406
          if (v[i] < min) min = v[i];
          if (v[i] > max) max = v[i];
407
408
409
410
        }
    }

  double range;
411
  if (IS_EQUAL(min, DBL_MAX) && IS_EQUAL(max, -DBL_MAX))
412
413
    range = missval;
  else
414
415
    range = max - min;

416
417
418
  return range;
}

419
double
420
varraySum(const size_t len, const Varray<double> &v)
421
{
422
  assert(v.size() >= len);
423
424
425

  double sum = 0;

426
  for (size_t i = 0; i < len; ++i) sum += v[i];
427
428
429
430

  return sum;
}

431
double
432
varraySumMV(const size_t len, const Varray<double> &v, const double missval)
433
{
434
  assert(v.size() >= len);
435
436
437
438

  double sum = 0;
  size_t nvals = 0;

439
  if (std::isnan(missval))
440
    {
441
      for (size_t i = 0; i < len; ++i)
442
        if (!DBL_IS_EQUAL(v[i], missval))
443
          {
444
            sum += v[i];
445
446
447
448
449
            nvals++;
          }
    }
  else
    {
450
      for (size_t i = 0; i < len; ++i)
451
        if (IS_NOT_EQUAL(v[i], missval))
452
          {
453
            sum += v[i];
454
455
456
457
            nvals++;
          }
    }

458
  if (!nvals) sum = missval;
459
460
461

  return sum;
}
462

463
double
464
varrayMean(const size_t len, const Varray<double> &v)
465
{
466
  assert(v.size() >= len);
467

468
  const double sum = varraySum(len, v);
469

470
  return sum / len;
471
472
}

473
double
474
varrayMeanMV(const size_t len, const Varray<double> &v, const double missval)
475
{
476
  assert(v.size() >= len);
477
478
479

  double sum = 0, sumw = 0;

480
  for (size_t i = 0; i < len; ++i)
481
    if (!DBL_IS_EQUAL(v[i], missval))
482
      {
483
        sum += v[i];
484
485
486
487
488
489
490
        sumw += 1;
      }

  double missval1 = missval, missval2 = missval;
  return DIVMN(sum, sumw);
}

491
double
492
varrayWeightedMean(const size_t len, const Varray<double> &v, const Varray<double> &w, const double missval)
493
{
494
  assert(v.size() >= len);
495
  assert(w.size() >= len);
496

497
  double sum = 0, sumw = 0;
498

499
  for (size_t i = 0; i < len; ++i)
500
    {
501
      sum += w[i] * v[i];
502
      sumw += w[i];
503
504
    }

505
  return IS_EQUAL(sumw, 0.) ? missval : sum / sumw;
506
507
}

508
double
509
varrayWeightedMeanMV(const size_t len, const Varray<double> &v, const Varray<double> &w, const double missval)
510
{
511
  assert(v.size() >= len);
512
  assert(w.size() >= len);
513

514
  const double missval1 = missval, missval2 = missval;
515
516
  double sum = 0, sumw = 0;

517
  for (size_t i = 0; i < len; ++i)
518
    if (!DBL_IS_EQUAL(v[i], missval1) && !DBL_IS_EQUAL(w[i], missval1))
519
      {
520
        sum += w[i] * v[i];
521
522
523
524
        sumw += w[i];
      }

  return DIVMN(sum, sumw);
525
}
526

527
double
528
varrayAvgMV(const size_t len, const Varray<double> &v, const double missval)
529
{
530
  assert(v.size() >= len);
531

532
  const double missval1 = missval, missval2 = missval;
533
534
  double sum = 0, sumw = 0;

535
  for (size_t i = 0; i < len; ++i)
536
    {
537
      sum = ADDMN(sum, v[i]);
538
539
540
541
542
543
      sumw += 1;
    }

  return DIVMN(sum, sumw);
}

544
double
545
varrayWeightedAvgMV(const size_t len, const Varray<double> &v, const Varray<double> &w, const double missval)
546
{
547
  assert(v.size() >= len);
548
  assert(w.size() >= len);
549

550
  const double missval1 = missval, missval2 = missval;
551
552
  double sum = 0, sumw = 0;

553
554
  for (size_t i = 0; i < len; ++i)
    if (!DBL_IS_EQUAL(w[i], missval))
555
      {
556
        sum = ADDMN(sum, MULMN(w[i], v[i]));
557
558
559
560
561
        sumw = ADDMN(sumw, w[i]);
      }

  return DIVMN(sum, sumw);
}
562

Uwe Schulzweida's avatar
Uwe Schulzweida committed
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
static void
varrayPrevarsum0(size_t len, const Varray<double> &v, double &rsum, double &rsumw)
{
  rsum = 0;
  for (size_t i = 0; i < len; i++)
    {
      rsum += v[i];
    }
  rsumw = len;
}

static void
varrayPrevarsum0MV(size_t len, const Varray<double> &v, double missval, double &rsum, double &rsumw)
{
  rsum = rsumw = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        rsum += v[i];
        rsumw += 1;
      }
}

static void
587
588
589
590
591
592
593
594
595
596
597
598
599
varrayPrevarsum(size_t len, const Varray<double> &v,
                double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = 0;
  for (size_t i = 0; i < len; i++)
    {
      rsum += v[i];
      rsumq += v[i] * v[i];
    }
  rsumw = len;
  rsumwq = len;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
600
static void
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
varrayPrevarsumMV(size_t len, const Varray<double> &v, double missval,
                  double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = rsumw = rsumwq = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        rsum += v[i];
        rsumq += v[i] * v[i];
        rsumw += 1;
        rsumwq += 1;
      }
}

double
varrayVar(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayPrevarsumMV(len, v, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayPrevarsum(len, v, rsum, rsumw, rsumq, rsumwq);

  double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

double
varrayVar1(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayPrevarsumMV(len, v, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayPrevarsum(len, v, rsum, rsumw, rsumq, rsumwq);

  double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
static void
varrayWeightedPrevarsum(size_t len, const Varray<double> &v, const Varray<double> &w,
                        double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = rsumw = rsumwq = 0;
  for (size_t i = 0; i < len; i++)
    {
      rsum += w[i] * v[i];
      rsumq += w[i] * v[i] * v[i];
      rsumw += w[i];
      rsumwq += w[i] * w[i];
    }
}

static void
varrayWeightedPrevarsumMV(size_t len, const Varray<double> &v, const Varray<double> &w, double missval,
                          double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
  rsum = rsumq = rsumw = rsumwq = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval) && !DBL_IS_EQUAL(w[i], missval))
      {
        rsum += w[i] * v[i];
        rsumq += w[i] * v[i] * v[i];
        rsumw += w[i];
        rsumwq += w[i] * w[i];
      }
}

double
varrayWeightedVar(size_t len, const Varray<double> &v, const Varray<double> &w, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayWeightedPrevarsumMV(len, v, w, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayWeightedPrevarsum(len, v, w, rsum, rsumw, rsumq, rsumwq);

  double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}

double
varrayWeightedVar1(size_t len, const Varray<double> &v, const Varray<double> &w, size_t nmiss, double missval)
{
  double rsum, rsumw, rsumq, rsumwq;
  if (nmiss)
    varrayWeightedPrevarsumMV(len, v, w, missval, rsum, rsumw, rsumq, rsumwq);
  else
    varrayWeightedPrevarsum(len, v, w, rsum, rsumw, rsumq, rsumwq);

  double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
  if (rvar < 0 && rvar > -1.e-5) rvar = 0;

  return rvar;
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740

static void
varrayPrekurtsum(size_t len, const Varray<double> &v, const double mean,
                 double &rsum3w, double &rsum4w, double &rsum2diff, double &rsum4diff)
{
  rsum2diff = rsum4diff = 0;
  for (size_t i = 0; i < len; i++)
    {
      const auto vdiff = v[i] - mean;
      rsum2diff += vdiff * vdiff;
      rsum4diff += vdiff * vdiff * vdiff * vdiff;
    }
  rsum3w = len;
  rsum4w = len;
}

static void
varrayPrekurtsumMV(size_t len, const Varray<double> &v, double missval, const double mean,
                   double &rsum3w, double &rsum4w, double &rsum2diff, double &rsum4diff)
{
  rsum3w = rsum4w = rsum2diff = rsum4diff = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        const auto vdiff = v[i] - mean;
        rsum2diff += vdiff * vdiff;
        rsum4diff += vdiff * vdiff * vdiff * vdiff;
        rsum3w += 1;
        rsum4w += 1;
      }
}

double
varrayKurt(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum3w;  // 3rd moment variables
  double rsum4w;  // 4th moment variables
  double rsum2diff, rsum4diff;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
741
  double rsum, rsumw;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
742
743
  if (nmiss)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
744
      varrayPrevarsum0MV(len, v, missval, rsum, rsumw);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
745
746
747
748
      varrayPrekurtsumMV(len, v, missval, (rsum / rsumw), rsum3w, rsum4w, rsum2diff, rsum4diff);
    }
  else
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
749
      varrayPrevarsum0(len, v, rsum, rsumw);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
750
751
752
753
754
      varrayPrekurtsum(len, v, (rsum / rsumw), rsum3w, rsum4w, rsum2diff, rsum4diff);
    }

  if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum2diff, 0.0)) return missval;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
756
  double rkurt = ((rsum4diff / rsum3w) / std::pow(rsum2diff / rsum3w, 2)) - 3.0;
  if (rkurt < 0 && rkurt > -1.e-5) rkurt = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
757

Uwe Schulzweida's avatar
Uwe Schulzweida committed
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
  return rkurt;
}

static void
varrayPreskewsum(size_t len, const Varray<double> &v, const double mean,
                 double &rsum3w, double &rsum4w, double &rsum3diff, double &rsum2diff)
{
  rsum3diff = rsum2diff = 0;
  for (size_t i = 0; i < len; i++)
    {
      const auto vdiff = v[i] - mean;
      rsum3diff += vdiff * vdiff * vdiff;
      rsum2diff += vdiff * vdiff;
    }
  rsum3w = len;
  rsum4w = len;
}

static void
varrayPreskewsumMV(size_t len, const Varray<double> &v, double missval, const double mean,
                   double &rsum3w, double &rsum4w, double &rsum3diff, double &rsum2diff)
{
  rsum3w = rsum4w = rsum3diff = rsum2diff = 0;
  for (size_t i = 0; i < len; i++)
    if (!DBL_IS_EQUAL(v[i], missval))
      {
        const auto vdiff = v[i] - mean;
        rsum3diff += vdiff * vdiff * vdiff;
        rsum2diff += vdiff * vdiff;
        rsum3w += 1;
        rsum4w += 1;
      }
}

double
varraySkew(size_t len, const Varray<double> &v, size_t nmiss, double missval)
{
  double rsum3w;  // 3rd moment variables
  double rsum4w;  // 4th moment variables
  double rsum3diff, rsum2diff;
  double rsum, rsumw;
  if (nmiss)
    {
      varrayPrevarsum0MV(len, v, missval, rsum, rsumw);
      varrayPreskewsumMV(len, v, missval, (rsum / rsumw), rsum3w, rsum4w, rsum3diff, rsum2diff);
    }
  else
    {
      varrayPrevarsum0(len, v, rsum, rsumw);
      varrayPreskewsum(len, v, (rsum / rsumw), rsum3w, rsum4w, rsum3diff, rsum2diff);
    }

  if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum3w, 1.0) || IS_EQUAL(rsum2diff, 0.0)) return missval;

  double rskew = (rsum3diff / rsum3w) / std::pow((rsum2diff) / (rsum3w - 1.0), 1.5);
  if (rskew < 0 && rskew > -1.e-5) rskew = 0;

  return rskew;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816
}