master
c 520 lines 16.7 KB
Raw
1 // SPDX-License-Identifier: GPL-3.0-or-later
2
3 #include "../libnetdata.h"
4
5 NETDATA_DOUBLE default_single_exponential_smoothing_alpha = 0.1;
6
7 void log_series_to_stderr(NETDATA_DOUBLE *series, size_t entries, NETDATA_DOUBLE result, const char *msg) {
8 const NETDATA_DOUBLE *value, *end = &series[entries];
9
10 fprintf(stderr, "%s of %zu entries [ ", msg, entries);
11 for(value = series; value < end ;value++) {
12 if(value != series) fprintf(stderr, ", ");
13 fprintf(stderr, "%" NETDATA_DOUBLE_MODIFIER, *value);
14 }
15 fprintf(stderr, " ] results in " NETDATA_DOUBLE_FORMAT "\n", result);
16 }
17
18 // --------------------------------------------------------------------------------------------------------------------
19
20 inline NETDATA_DOUBLE sum_and_count(const NETDATA_DOUBLE *series, size_t entries, size_t *count) {
21 const NETDATA_DOUBLE *value, *end = &series[entries];
22 NETDATA_DOUBLE sum = 0;
23 size_t c = 0;
24
25 for(value = series; value < end ; value++) {
26 if(netdata_double_isnumber(*value)) {
27 sum += *value;
28 c++;
29 }
30 }
31
32 if(unlikely(!c)) sum = NAN;
33 if(likely(count)) *count = c;
34
35 return sum;
36 }
37
38 inline NETDATA_DOUBLE sum(const NETDATA_DOUBLE *series, size_t entries) {
39 return sum_and_count(series, entries, NULL);
40 }
41
42 inline NETDATA_DOUBLE average(const NETDATA_DOUBLE *series, size_t entries) {
43 size_t count = 0;
44 NETDATA_DOUBLE sum = sum_and_count(series, entries, &count);
45
46 if(unlikely(!count)) return NAN;
47 return sum / (NETDATA_DOUBLE)count;
48 }
49
50 // --------------------------------------------------------------------------------------------------------------------
51
52 // periods up to this size use a stack buffer to avoid heap overhead
53 #define MOVING_AVERAGE_STACK_PERIOD 32
54
55 NETDATA_DOUBLE moving_average(const NETDATA_DOUBLE *series, size_t entries, size_t period) {
56 // Keep the zero-period fast path while making the rolling window size
57 // unambiguously non-zero for the arithmetic below.
58 const size_t window = period ? period : 1;
59
60 if(unlikely(period == 0))
61 return 0.0;
62
63 size_t i, count;
64 NETDATA_DOUBLE sum = 0, avg = 0;
65
66 NETDATA_DOUBLE stack_buf[MOVING_AVERAGE_STACK_PERIOD];
67 NETDATA_DOUBLE *heap_p = NULL;
68 NETDATA_DOUBLE *p;
69
70 if(window <= MOVING_AVERAGE_STACK_PERIOD) {
71 memset(stack_buf, 0, window * sizeof(*stack_buf));
72 p = stack_buf;
73 } else {
74 heap_p = callocz(window, sizeof(*heap_p));
75 p = heap_p;
76 }
77
78 for(i = 0, count = 0; i < entries; i++) {
79 NETDATA_DOUBLE value = series[i];
80 if(unlikely(!netdata_double_isnumber(value))) continue;
81 size_t slot = count % window;
82
83 if(unlikely(count < window)) {
84 sum += value;
85 avg = (count == window - 1) ? sum / (NETDATA_DOUBLE)window : 0;
86 }
87 else {
88 sum = sum - p[slot] + value;
89 avg = sum / (NETDATA_DOUBLE)window;
90 }
91
92 p[slot] = value;
93 count++;
94 }
95
96 freez(heap_p);
97 return avg;
98 }
99
100 static int statistical_unittest_assert_close(const char *name, NETDATA_DOUBLE expected, NETDATA_DOUBLE actual) {
101 if(ABS(expected - actual) <= 0.000001)
102 return 0;
103
104 fprintf(stderr, "statistical_unittest: %s failed, expected " NETDATA_DOUBLE_FORMAT ", got " NETDATA_DOUBLE_FORMAT "\n",
105 name, expected, actual);
106 return 1;
107 }
108
109 int statistical_unittest(void) {
110 int errors = 0;
111
112 NETDATA_DOUBLE series[] = { 1, 2, 3, 4, 5 };
113 NETDATA_DOUBLE heap_series[MOVING_AVERAGE_STACK_PERIOD + 8];
114
115 for(size_t i = 0; i < sizeof(heap_series) / sizeof(heap_series[0]); i++)
116 heap_series[i] = (NETDATA_DOUBLE)(i + 1);
117
118 errors += statistical_unittest_assert_close("moving_average(period=0)", 0.0,
119 moving_average(series, sizeof(series) / sizeof(series[0]), 0));
120 errors += statistical_unittest_assert_close("moving_average(stack path)", 4.0,
121 moving_average(series, sizeof(series) / sizeof(series[0]), 3));
122 errors += statistical_unittest_assert_close("moving_average(heap path)",
123 ((NETDATA_DOUBLE)(sizeof(heap_series) / sizeof(heap_series[0])) + 1.0) / 2.0,
124 moving_average(heap_series, sizeof(heap_series) / sizeof(heap_series[0]),
125 sizeof(heap_series) / sizeof(heap_series[0])));
126
127 if(errors)
128 fprintf(stderr, "statistical_unittest: %d errors found\n", errors);
129 else
130 fprintf(stderr, "statistical_unittest: all tests passed\n");
131
132 return errors ? 1 : 0;
133 }
134
135 // --------------------------------------------------------------------------------------------------------------------
136
137 static int qsort_compare(const void *a, const void *b) {
138 NETDATA_DOUBLE *p1 = (NETDATA_DOUBLE *)a, *p2 = (NETDATA_DOUBLE *)b;
139 NETDATA_DOUBLE n1 = *p1, n2 = *p2;
140
141 if(unlikely(isnan(n1) || isnan(n2))) {
142 if(isnan(n1) && !isnan(n2)) return -1;
143 if(!isnan(n1) && isnan(n2)) return 1;
144 return 0;
145 }
146 if(unlikely(isinf(n1) || isinf(n2))) {
147 if(!isinf(n1) && isinf(n2)) return -1;
148 if(isinf(n1) && !isinf(n2)) return 1;
149 return 0;
150 }
151
152 if(unlikely(n1 < n2)) return -1;
153 if(unlikely(n1 > n2)) return 1;
154 return 0;
155 }
156
157 inline void sort_series(NETDATA_DOUBLE *series, size_t entries) {
158 qsort(series, entries, sizeof(NETDATA_DOUBLE), qsort_compare);
159 }
160
161 inline NETDATA_DOUBLE *copy_series(const NETDATA_DOUBLE *series, size_t entries) {
162 NETDATA_DOUBLE *copy = mallocz(sizeof(NETDATA_DOUBLE) * entries);
163 memcpy(copy, series, sizeof(NETDATA_DOUBLE) * entries);
164 return copy;
165 }
166
167 NETDATA_DOUBLE percentile_on_sorted_series(const NETDATA_DOUBLE *series, size_t entries, double percentile) {
168 if (unlikely(entries == 0)) return NAN;
169 if (unlikely(entries == 1)) return series[0];
170
171 // Clamp percentile between 0.0 and 1.0
172 percentile = fmax(0.0, fmin(1.0, percentile));
173
174 // Compute fractional index
175 NETDATA_DOUBLE index = percentile * (NETDATA_DOUBLE)(entries - 1);
176 size_t low_idx = (size_t)floor(index);
177 size_t high_idx = (size_t)ceil(index);;
178
179 // If index is an integer or at the last element, return directly
180 if (high_idx >= entries || low_idx == high_idx || considered_equal_ndd(index, (NETDATA_DOUBLE)low_idx))
181 return series[low_idx];
182
183 // Linear interpolation
184 NETDATA_DOUBLE weight = index - (NETDATA_DOUBLE)low_idx;
185 return series[low_idx] + weight * (series[high_idx] - series[low_idx]);
186 }
187
188 NETDATA_DOUBLE median_on_sorted_series(const NETDATA_DOUBLE *series, size_t entries) {
189 return percentile_on_sorted_series(series, entries, 0.5);
190 }
191
192 NETDATA_DOUBLE median(const NETDATA_DOUBLE *series, size_t entries) {
193 if(unlikely(entries == 0)) return NAN;
194 if(unlikely(entries == 1)) return series[0];
195
196 if(unlikely(entries == 2))
197 return (series[0] + series[1]) / 2;
198
199 NETDATA_DOUBLE *copy = copy_series(series, entries);
200 sort_series(copy, entries);
201
202 NETDATA_DOUBLE avg = median_on_sorted_series(copy, entries);
203
204 freez(copy);
205 return avg;
206 }
207
208 // --------------------------------------------------------------------------------------------------------------------
209
210 NETDATA_DOUBLE moving_median(const NETDATA_DOUBLE *series, size_t entries, size_t period) {
211 if(entries <= period)
212 return median(series, entries);
213
214 NETDATA_DOUBLE *data = copy_series(series, entries);
215
216 size_t i;
217 for(i = period; i < entries; i++) {
218 data[i - period] = median(&series[i - period], period);
219 }
220
221 NETDATA_DOUBLE avg = median(data, entries - period);
222 freez(data);
223 return avg;
224 }
225
226 // --------------------------------------------------------------------------------------------------------------------
227
228 // http://stackoverflow.com/a/15150143/4525767
229 NETDATA_DOUBLE running_median_estimate(const NETDATA_DOUBLE *series, size_t entries) {
230 NETDATA_DOUBLE median = 0.0f;
231 NETDATA_DOUBLE average = 0.0f;
232 size_t i;
233
234 for(i = 0; i < entries ; i++) {
235 NETDATA_DOUBLE value = series[i];
236 if(unlikely(!netdata_double_isnumber(value))) continue;
237
238 average += ( value - average ) * 0.1f; // rough running average.
239 median += copysignndd( average * 0.01, value - median );
240 }
241
242 return median;
243 }
244
245 // --------------------------------------------------------------------------------------------------------------------
246
247 NETDATA_DOUBLE standard_deviation(const NETDATA_DOUBLE *series, size_t entries) {
248 if(unlikely(entries == 0)) return NAN;
249 if(unlikely(entries == 1)) return series[0];
250
251 const NETDATA_DOUBLE *value, *end = &series[entries];
252 size_t count;
253 NETDATA_DOUBLE sum;
254
255 for(count = 0, sum = 0, value = series ; value < end ;value++) {
256 if(likely(netdata_double_isnumber(*value))) {
257 count++;
258 sum += *value;
259 }
260 }
261
262 if(unlikely(count == 0)) return NAN;
263 if(unlikely(count == 1)) return sum;
264
265 NETDATA_DOUBLE average = sum / (NETDATA_DOUBLE)count;
266
267 for(count = 0, sum = 0, value = series ; value < end ;value++) {
268 if(netdata_double_isnumber(*value)) {
269 count++;
270 sum += powndd(*value - average, 2);
271 }
272 }
273
274 if(unlikely(count == 0)) return NAN;
275 if(unlikely(count == 1)) return average;
276
277 NETDATA_DOUBLE variance = sum / (NETDATA_DOUBLE)(count); // remove -1 from count to have a population stddev
278 NETDATA_DOUBLE stddev = sqrtndd(variance);
279 return stddev;
280 }
281
282 // --------------------------------------------------------------------------------------------------------------------
283
284 NETDATA_DOUBLE single_exponential_smoothing(const NETDATA_DOUBLE *series, size_t entries, NETDATA_DOUBLE alpha) {
285 if(unlikely(entries == 0))
286 return NAN;
287
288 if(unlikely(isnan(alpha)))
289 alpha = default_single_exponential_smoothing_alpha;
290
291 const NETDATA_DOUBLE *value = series, *end = &series[entries];
292 NETDATA_DOUBLE level = (1.0 - alpha) * (*value);
293
294 for(value++ ; value < end; value++) {
295 if(likely(netdata_double_isnumber(*value)))
296 level = alpha * (*value) + (1.0 - alpha) * level;
297 }
298
299 return level;
300 }
301
302 NETDATA_DOUBLE single_exponential_smoothing_reverse(const NETDATA_DOUBLE *series, size_t entries, NETDATA_DOUBLE alpha) {
303 if(unlikely(entries == 0))
304 return NAN;
305
306 if(unlikely(isnan(alpha)))
307 alpha = default_single_exponential_smoothing_alpha;
308
309 const NETDATA_DOUBLE *value = &series[entries -1];
310 NETDATA_DOUBLE level = (1.0 - alpha) * (*value);
311
312 for(value++ ; value >= series; value--) {
313 if(likely(netdata_double_isnumber(*value)))
314 level = alpha * (*value) + (1.0 - alpha) * level;
315 }
316
317 return level;
318 }
319
320 // --------------------------------------------------------------------------------------------------------------------
321
322 // http://grisha.org/blog/2016/02/16/triple-exponential-smoothing-forecasting-part-ii/
323 NETDATA_DOUBLE double_exponential_smoothing(const NETDATA_DOUBLE *series, size_t entries,
324 NETDATA_DOUBLE alpha,
325 NETDATA_DOUBLE beta,
326 NETDATA_DOUBLE *forecast) {
327 if(unlikely(entries == 0))
328 return NAN;
329
330 NETDATA_DOUBLE level, trend;
331
332 if(unlikely(isnan(alpha)))
333 alpha = 0.3;
334
335 if(unlikely(isnan(beta)))
336 beta = 0.05;
337
338 level = series[0];
339
340 if(likely(entries > 1))
341 trend = series[1] - series[0];
342 else
343 trend = 0;
344
345 const NETDATA_DOUBLE *value = series;
346 for(value++ ; value >= series; value--) {
347 if(likely(netdata_double_isnumber(*value))) {
348 NETDATA_DOUBLE last_level = level;
349 level = alpha * *value + (1.0 - alpha) * (level + trend);
350 trend = beta * (level - last_level) + (1.0 - beta) * trend;
351
352 }
353 }
354
355 if(forecast)
356 *forecast = level + trend;
357
358 return level;
359 }
360
361 // --------------------------------------------------------------------------------------------------------------------
362
363 /*
364 * Based on th R implementation
365 *
366 * a: level component
367 * b: trend component
368 * s: seasonal component
369 *
370 * Additive:
371 *
372 * Yhat[t+h] = a[t] + h * b[t] + s[t + 1 + (h - 1) mod p],
373 * a[t] = α (Y[t] - s[t-p]) + (1-α) (a[t-1] + b[t-1])
374 * b[t] = β (a[t] - a[t-1]) + (1-β) b[t-1]
375 * s[t] = γ (Y[t] - a[t]) + (1-γ) s[t-p]
376 *
377 * Multiplicative:
378 *
379 * Yhat[t+h] = (a[t] + h * b[t]) * s[t + 1 + (h - 1) mod p],
380 * a[t] = α (Y[t] / s[t-p]) + (1-α) (a[t-1] + b[t-1])
381 * b[t] = β (a[t] - a[t-1]) + (1-β) b[t-1]
382 * s[t] = γ (Y[t] / a[t]) + (1-γ) s[t-p]
383 */
384 static int __HoltWinters(
385 const NETDATA_DOUBLE *series,
386 int entries, // start_time + h
387
388 NETDATA_DOUBLE alpha, // alpha parameter of Holt-Winters Filter.
389 NETDATA_DOUBLE
390 beta, // beta parameter of Holt-Winters Filter. If set to 0, the function will do exponential smoothing.
391 NETDATA_DOUBLE
392 gamma, // gamma parameter used for the seasonal component. If set to 0, an non-seasonal model is fitted.
393
394 const int *seasonal,
395 const int *period,
396 const NETDATA_DOUBLE *a, // Start value for level (a[0]).
397 const NETDATA_DOUBLE *b, // Start value for trend (b[0]).
398 NETDATA_DOUBLE *s, // Vector of start values for the seasonal component (s_1[0] ... s_p[0])
399
400 /* return values */
401 NETDATA_DOUBLE *SSE, // The final sum of squared errors achieved in optimizing
402 NETDATA_DOUBLE *level, // Estimated values for the level component (size entries - t + 2)
403 NETDATA_DOUBLE *trend, // Estimated values for the trend component (size entries - t + 2)
404 NETDATA_DOUBLE *season // Estimated values for the seasonal component (size entries - t + 2)
405 )
406 {
407 if(unlikely(entries < 4))
408 return 0;
409
410 int start_time = 2;
411
412 NETDATA_DOUBLE res = 0, xhat = 0, stmp = 0;
413 int i, i0, s0;
414
415 /* copy start values to the beginning of the vectors */
416 level[0] = *a;
417 if(beta > 0) trend[0] = *b;
418 if(gamma > 0) memcpy(season, s, *period * sizeof(NETDATA_DOUBLE));
419
420 for(i = start_time - 1; i < entries; i++) {
421 /* indices for period i */
422 i0 = i - start_time + 2;
423 s0 = i0 + *period - 1;
424
425 /* forecast *for* period i */
426 xhat = level[i0 - 1] + (beta > 0 ? trend[i0 - 1] : 0);
427 stmp = gamma > 0 ? season[s0 - *period] : (*seasonal != 1);
428 if (*seasonal == 1)
429 xhat += stmp;
430 else
431 xhat *= stmp;
432
433 /* Sum of Squared Errors */
434 res = series[i] - xhat;
435 *SSE += res * res;
436
437 /* estimate of level *in* period i */
438 if (*seasonal == 1)
439 level[i0] = alpha * (series[i] - stmp)
440 + (1 - alpha) * (level[i0 - 1] + trend[i0 - 1]);
441 else
442 level[i0] = alpha * (series[i] / stmp)
443 + (1 - alpha) * (level[i0 - 1] + trend[i0 - 1]);
444
445 /* estimate of trend *in* period i */
446 if (beta > 0)
447 trend[i0] = beta * (level[i0] - level[i0 - 1])
448 + (1 - beta) * trend[i0 - 1];
449
450 /* estimate of seasonal component *in* period i */
451 if (gamma > 0) {
452 if (*seasonal == 1)
453 season[s0] = gamma * (series[i] - level[i0])
454 + (1 - gamma) * stmp;
455 else
456 season[s0] = gamma * (series[i] / level[i0])
457 + (1 - gamma) * stmp;
458 }
459 }
460
461 return 1;
462 }
463
464 NETDATA_DOUBLE holtwinters(const NETDATA_DOUBLE *series, size_t entries,
465 NETDATA_DOUBLE alpha,
466 NETDATA_DOUBLE beta,
467 NETDATA_DOUBLE gamma,
468 NETDATA_DOUBLE *forecast) {
469 if(unlikely(isnan(alpha)))
470 alpha = 0.3;
471
472 if(unlikely(isnan(beta)))
473 beta = 0.05;
474
475 if(unlikely(isnan(gamma)))
476 gamma = 0;
477
478 int seasonal = 0;
479 int period = 0;
480 NETDATA_DOUBLE a0 = series[0];
481 NETDATA_DOUBLE b0 = 0;
482 NETDATA_DOUBLE s[] = {};
483
484 NETDATA_DOUBLE errors = 0.0;
485 size_t nb_computations = entries;
486 NETDATA_DOUBLE *estimated_level = callocz(nb_computations, sizeof(NETDATA_DOUBLE));
487 NETDATA_DOUBLE *estimated_trend = callocz(nb_computations, sizeof(NETDATA_DOUBLE));
488 NETDATA_DOUBLE *estimated_season = callocz(nb_computations, sizeof(NETDATA_DOUBLE));
489
490 int ret = __HoltWinters(
491 series,
492 (int)entries,
493 alpha,
494 beta,
495 gamma,
496 &seasonal,
497 &period,
498 &a0,
499 &b0,
500 s,
501 &errors,
502 estimated_level,
503 estimated_trend,
504 estimated_season
505 );
506
507 NETDATA_DOUBLE value = estimated_level[nb_computations - 1];
508
509 if(forecast)
510 *forecast = 0.0;
511
512 freez(estimated_level);
513 freez(estimated_trend);
514 freez(estimated_season);
515
516 if(!ret)
517 return 0.0;
518
519 return value;
520 }