1 /*
2 * Copyright (c) 2013-2021, The OSKAR Developers.
3 * See the LICENSE file at the top-level directory of this distribution.
4 */
5
6 #include "mem/oskar_mem.h"
7
8 #include <math.h>
9 #include <float.h>
10
11 #ifdef __cplusplus
12 extern "C" {
13 #endif
14
15 /* Macro to update running statistics for mean and standard deviation using
16 * the method of Donald Knuth in "The Art of Computer Programming"
17 * vol 2, 3rd edition, page 232 */
18 #define RUNNING_STATS_KNUTH \
19 if (max && val > *max) \
20 { \
21 *max = val; \
22 } \
23 if (min && val < *min) \
24 { \
25 *min = val; \
26 } \
27 if (i == 0) \
28 { \
29 old_m = new_m = val; \
30 old_s = 0.0; \
31 } \
32 else \
33 { \
34 new_m = old_m + (val - old_m) / (i + 1); \
35 new_s = old_s + (val - old_m) * (val - new_m); \
36 old_m = new_m; \
37 old_s = new_s; \
38 }
39
oskar_mem_stats(const oskar_Mem * mem,size_t n,double * min,double * max,double * mean,double * std_dev,int * status)40 void oskar_mem_stats(const oskar_Mem* mem, size_t n, double* min, double* max,
41 double* mean, double* std_dev, int* status)
42 {
43 int type = 0;
44 size_t i = 0;
45 double old_m = 0.0, new_m = 0.0, old_s = 0.0, new_s = 0.0;
46
47 /* Check if safe to proceed. */
48 if (*status) return;
49
50 /* Check that data is in CPU accessible memory. */
51 if (oskar_mem_location(mem) != OSKAR_CPU)
52 {
53 *status = OSKAR_ERR_BAD_LOCATION;
54 return;
55 }
56
57 /* Check that the data type is single or double precision scalar. */
58 type = oskar_mem_type(mem);
59 if (oskar_type_is_complex(type) || oskar_type_is_matrix(type))
60 {
61 *status = OSKAR_ERR_BAD_DATA_TYPE;
62 return;
63 }
64
65 /* Initialise outputs. */
66 if (max) *max = -DBL_MAX;
67 if (min) *min = DBL_MAX;
68 if (mean) *mean = 0.0;
69 if (std_dev) *std_dev = 0.0;
70
71 /* Gather statistics. */
72 if (type == OSKAR_SINGLE)
73 {
74 double val = 0.0;
75 const float *data = 0;
76 data = oskar_mem_float_const(mem, status);
77 for (i = 0; i < n; ++i)
78 {
79 val = (double) data[i];
80 RUNNING_STATS_KNUTH
81 }
82 }
83 else if (type == OSKAR_DOUBLE)
84 {
85 double val = 0.0;
86 const double *data = 0;
87 data = oskar_mem_double_const(mem, status);
88 for (i = 0; i < n; ++i)
89 {
90 val = data[i];
91 RUNNING_STATS_KNUTH
92 }
93 }
94 else
95 {
96 *status = OSKAR_ERR_BAD_DATA_TYPE;
97 return;
98 }
99
100 /* Set mean and population standard deviation. */
101 if (mean) *mean = new_m;
102 if (std_dev) *std_dev = (n > 0) ? sqrt(new_s / n) : 0.0;
103 }
104
105 #ifdef __cplusplus
106 }
107 #endif
108