1 /*
2 ** Zabbix
3 ** Copyright (C) 2001-2021 Zabbix SIA
4 **
5 ** This program is free software; you can redistribute it and/or modify
6 ** it under the terms of the GNU General Public License as published by
7 ** the Free Software Foundation; either version 2 of the License, or
8 ** (at your option) any later version.
9 **
10 ** This program is distributed in the hope that it will be useful,
11 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
12 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 ** GNU General Public License for more details.
14 **
15 ** You should have received a copy of the GNU General Public License
16 ** along with this program; if not, write to the Free Software
17 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 **/
19
20 #include "common.h"
21 #include "log.h"
22
23 #include "zbxalgo.h"
24
25 #define DB_INFINITY (1e12 - 1e-4)
26
27 #define ZBX_MATH_EPSILON (1e-6)
28
29 #define ZBX_IS_NAN(x) ((x) != (x))
30
31 #define ZBX_VALID_MATRIX(m) (0 < (m)->rows && 0 < (m)->columns && NULL != (m)->elements)
32 #define ZBX_MATRIX_EL(m, row, col) ((m)->elements[(row) * (m)->columns + (col)])
33 #define ZBX_MATRIX_ROW(m, row) ((m)->elements + (row) * (m)->columns)
34
35 typedef struct
36 {
37 int rows;
38 int columns;
39 double *elements;
40 }
41 zbx_matrix_t;
42
zbx_matrix_struct_alloc(zbx_matrix_t ** pm)43 static void zbx_matrix_struct_alloc(zbx_matrix_t **pm)
44 {
45 *pm = (zbx_matrix_t *)zbx_malloc(*pm, sizeof(zbx_matrix_t));
46
47 (*pm)->rows = 0;
48 (*pm)->columns = 0;
49 (*pm)->elements = NULL;
50 }
51
zbx_matrix_alloc(zbx_matrix_t * m,int rows,int columns)52 static int zbx_matrix_alloc(zbx_matrix_t *m, int rows, int columns)
53 {
54 if (0 >= rows || 0 >= columns)
55 goto error;
56
57 m->rows = rows;
58 m->columns = columns;
59
60 m->elements = (double *)zbx_malloc(m->elements, sizeof(double) * rows * columns);
61
62 return SUCCEED;
63 error:
64 THIS_SHOULD_NEVER_HAPPEN;
65 return FAIL;
66 }
67
zbx_matrix_free(zbx_matrix_t * m)68 static void zbx_matrix_free(zbx_matrix_t *m)
69 {
70 if (NULL != m)
71 zbx_free(m->elements);
72
73 zbx_free(m);
74 }
75
zbx_matrix_copy(zbx_matrix_t * dest,zbx_matrix_t * src)76 static int zbx_matrix_copy(zbx_matrix_t *dest, zbx_matrix_t *src)
77 {
78 if (!ZBX_VALID_MATRIX(src))
79 goto error;
80
81 if (SUCCEED != zbx_matrix_alloc(dest, src->rows, src->columns))
82 return FAIL;
83
84 memcpy(dest->elements, src->elements, sizeof(double) * src->rows * src->columns);
85 return SUCCEED;
86 error:
87 THIS_SHOULD_NEVER_HAPPEN;
88 return FAIL;
89 }
90
zbx_identity_matrix(zbx_matrix_t * m,int n)91 static int zbx_identity_matrix(zbx_matrix_t *m, int n)
92 {
93 int i, j;
94
95 if (SUCCEED != zbx_matrix_alloc(m, n, n))
96 return FAIL;
97
98 for (i = 0; i < n; i++)
99 for (j = 0; j < n; j++)
100 ZBX_MATRIX_EL(m, i, j) = (i == j ? 1.0 : 0.0);
101
102 return SUCCEED;
103 }
104
zbx_transpose_matrix(zbx_matrix_t * m,zbx_matrix_t * r)105 static int zbx_transpose_matrix(zbx_matrix_t *m, zbx_matrix_t *r)
106 {
107 int i, j;
108
109 if (!ZBX_VALID_MATRIX(m))
110 goto error;
111
112 if (SUCCEED != zbx_matrix_alloc(r, m->columns, m->rows))
113 return FAIL;
114
115 for (i = 0; i < r->rows; i++)
116 for (j = 0; j < r->columns; j++)
117 ZBX_MATRIX_EL(r, i, j) = ZBX_MATRIX_EL(m, j, i);
118
119 return SUCCEED;
120 error:
121 THIS_SHOULD_NEVER_HAPPEN;
122 return FAIL;
123 }
124
zbx_matrix_swap_rows(zbx_matrix_t * m,int r1,int r2)125 static void zbx_matrix_swap_rows(zbx_matrix_t *m, int r1, int r2)
126 {
127 double tmp;
128 int i;
129
130 for (i = 0; i < m->columns; i++)
131 {
132 tmp = ZBX_MATRIX_EL(m, r1, i);
133 ZBX_MATRIX_EL(m, r1, i) = ZBX_MATRIX_EL(m, r2, i);
134 ZBX_MATRIX_EL(m, r2, i) = tmp;
135 }
136 }
137
zbx_matrix_divide_row_by(zbx_matrix_t * m,int row,double denominator)138 static void zbx_matrix_divide_row_by(zbx_matrix_t *m, int row, double denominator)
139 {
140 int i;
141
142 for (i = 0; i < m->columns; i++)
143 ZBX_MATRIX_EL(m, row, i) /= denominator;
144 }
145
zbx_matrix_add_rows_with_factor(zbx_matrix_t * m,int dest,int src,double factor)146 static void zbx_matrix_add_rows_with_factor(zbx_matrix_t *m, int dest, int src, double factor)
147 {
148 int i;
149
150 for (i = 0; i < m->columns; i++)
151 ZBX_MATRIX_EL(m, dest, i) += ZBX_MATRIX_EL(m, src, i) * factor;
152 }
153
zbx_inverse_matrix(zbx_matrix_t * m,zbx_matrix_t * r)154 static int zbx_inverse_matrix(zbx_matrix_t *m, zbx_matrix_t *r)
155 {
156 zbx_matrix_t *l = NULL;
157 double pivot, factor, det;
158 int i, j, k, n, res;
159
160 if (!ZBX_VALID_MATRIX(m) || m->rows != m->columns)
161 goto error;
162
163 n = m->rows;
164
165 if (1 == n)
166 {
167 if (SUCCEED != zbx_matrix_alloc(r, 1, 1))
168 return FAIL;
169
170 if (0.0 == ZBX_MATRIX_EL(m, 0, 0))
171 {
172 zabbix_log(LOG_LEVEL_DEBUG, "matrix is singular");
173 res = FAIL;
174 goto out;
175 }
176
177 ZBX_MATRIX_EL(r, 0, 0) = 1.0 / ZBX_MATRIX_EL(m, 0, 0);
178 return SUCCEED;
179 }
180
181 if (2 == n)
182 {
183 if (SUCCEED != zbx_matrix_alloc(r, 2, 2))
184 return FAIL;
185
186 if (0.0 == (det = ZBX_MATRIX_EL(m, 0, 0) * ZBX_MATRIX_EL(m, 1, 1) -
187 ZBX_MATRIX_EL(m, 0, 1) * ZBX_MATRIX_EL(m, 1, 0)))
188 {
189 zabbix_log(LOG_LEVEL_DEBUG, "matrix is singular");
190 res = FAIL;
191 goto out;
192 }
193
194 ZBX_MATRIX_EL(r, 0, 0) = ZBX_MATRIX_EL(m, 1, 1) / det;
195 ZBX_MATRIX_EL(r, 0, 1) = -ZBX_MATRIX_EL(m, 0, 1) / det;
196 ZBX_MATRIX_EL(r, 1, 0) = -ZBX_MATRIX_EL(m, 1, 0) / det;
197 ZBX_MATRIX_EL(r, 1, 1) = ZBX_MATRIX_EL(m, 0, 0) / det;
198 return SUCCEED;
199 }
200
201 if (SUCCEED != zbx_identity_matrix(r, n))
202 return FAIL;
203
204 zbx_matrix_struct_alloc(&l);
205
206 if (SUCCEED != (res = zbx_matrix_copy(l, m)))
207 goto out;
208
209 /* Gauss-Jordan elimination with partial (row) pivoting */
210 for (i = 0; i < n; i++)
211 {
212 k = i;
213 pivot = ZBX_MATRIX_EL(l, i, i);
214
215 for (j = i; j < n; j++)
216 {
217 if (fabs(ZBX_MATRIX_EL(l, j, i)) > fabs(pivot))
218 {
219 k = j;
220 pivot = ZBX_MATRIX_EL(l, j, i);
221 }
222 }
223
224 if (0.0 == pivot)
225 {
226 zabbix_log(LOG_LEVEL_DEBUG, "matrix is singular");
227 res = FAIL;
228 goto out;
229 }
230
231 if (k != i)
232 {
233 zbx_matrix_swap_rows(l, i, k);
234 zbx_matrix_swap_rows(r, i, k);
235 }
236
237 for (j = i + 1; j < n; j++)
238 {
239 if (0.0 != (factor = -ZBX_MATRIX_EL(l, j, i) / ZBX_MATRIX_EL(l, i, i)))
240 {
241 zbx_matrix_add_rows_with_factor(l, j, i, factor);
242 zbx_matrix_add_rows_with_factor(r, j, i, factor);
243 }
244 }
245 }
246
247 for (i = n - 1; i > 0; i--)
248 {
249 for (j = 0; j < i; j++)
250 {
251 if (0.0 != (factor = -ZBX_MATRIX_EL(l, j, i) / ZBX_MATRIX_EL(l, i, i)))
252 {
253 zbx_matrix_add_rows_with_factor(l, j, i, factor);
254 zbx_matrix_add_rows_with_factor(r, j, i, factor);
255 }
256 }
257 }
258
259 for (i = 0; i < n; i++)
260 zbx_matrix_divide_row_by(r, i, ZBX_MATRIX_EL(l, i, i));
261
262 res = SUCCEED;
263 out:
264 zbx_matrix_free(l);
265 return res;
266 error:
267 THIS_SHOULD_NEVER_HAPPEN;
268 return FAIL;
269 }
270
zbx_matrix_mult(zbx_matrix_t * left,zbx_matrix_t * right,zbx_matrix_t * result)271 static int zbx_matrix_mult(zbx_matrix_t *left, zbx_matrix_t *right, zbx_matrix_t *result)
272 {
273 double element;
274 int i, j, k;
275
276 if (!ZBX_VALID_MATRIX(left) || !ZBX_VALID_MATRIX(right) || left->columns != right->rows)
277 goto error;
278
279 if (SUCCEED != zbx_matrix_alloc(result, left->rows, right->columns))
280 return FAIL;
281
282 for (i = 0; i < result->rows; i++)
283 {
284 for (j = 0; j < result->columns; j++)
285 {
286 element = 0;
287
288 for (k = 0; k < left->columns; k++)
289 element += ZBX_MATRIX_EL(left, i, k) * ZBX_MATRIX_EL(right, k, j);
290
291 ZBX_MATRIX_EL(result, i, j) = element;
292 }
293 }
294
295 return SUCCEED;
296 error:
297 THIS_SHOULD_NEVER_HAPPEN;
298 return FAIL;
299 }
300
zbx_least_squares(zbx_matrix_t * independent,zbx_matrix_t * dependent,zbx_matrix_t * coefficients)301 static int zbx_least_squares(zbx_matrix_t *independent, zbx_matrix_t *dependent, zbx_matrix_t *coefficients)
302 {
303 /* |<----------to_be_inverted---------->| */
304 /* coefficients = inverse( transpose( independent ) * independent ) * transpose( independent ) * dependent */
305 /* |<------------------left_part------------------>| |<-----------right_part----------->| */
306 /* we change order of matrix multiplication to reduce operation count and memory usage */
307 zbx_matrix_t *independent_transposed = NULL, *to_be_inverted = NULL, *left_part = NULL, *right_part = NULL;
308 int res;
309
310 zbx_matrix_struct_alloc(&independent_transposed);
311 zbx_matrix_struct_alloc(&to_be_inverted);
312 zbx_matrix_struct_alloc(&left_part);
313 zbx_matrix_struct_alloc(&right_part);
314
315 if (SUCCEED != (res = zbx_transpose_matrix(independent, independent_transposed)))
316 goto out;
317
318 if (SUCCEED != (res = zbx_matrix_mult(independent_transposed, independent, to_be_inverted)))
319 goto out;
320
321 if (SUCCEED != (res = zbx_inverse_matrix(to_be_inverted, left_part)))
322 goto out;
323
324 if (SUCCEED != (res = zbx_matrix_mult(independent_transposed, dependent, right_part)))
325 goto out;
326
327 if (SUCCEED != (res = zbx_matrix_mult(left_part, right_part, coefficients)))
328 goto out;
329
330 out:
331 zbx_matrix_free(independent_transposed);
332 zbx_matrix_free(to_be_inverted);
333 zbx_matrix_free(left_part);
334 zbx_matrix_free(right_part);
335 return res;
336 }
337
zbx_fill_dependent(double * x,int n,zbx_fit_t fit,zbx_matrix_t * m)338 static int zbx_fill_dependent(double *x, int n, zbx_fit_t fit, zbx_matrix_t *m)
339 {
340 int i;
341
342 if (FIT_LINEAR == fit || FIT_POLYNOMIAL == fit || FIT_LOGARITHMIC == fit)
343 {
344 if (SUCCEED != zbx_matrix_alloc(m, n, 1))
345 return FAIL;
346
347 for (i = 0; i < n; i++)
348 ZBX_MATRIX_EL(m, i, 0) = x[i];
349 }
350 else if (FIT_EXPONENTIAL == fit || FIT_POWER == fit)
351 {
352 if (SUCCEED != zbx_matrix_alloc(m, n, 1))
353 return FAIL;
354
355 for (i = 0; i < n; i++)
356 {
357 if (0.0 >= x[i])
358 {
359 zabbix_log(LOG_LEVEL_DEBUG, "data contains negative or zero values");
360 return FAIL;
361 }
362
363 ZBX_MATRIX_EL(m, i, 0) = log(x[i]);
364 }
365 }
366
367 return SUCCEED;
368 }
369
zbx_fill_independent(double * t,int n,zbx_fit_t fit,int k,zbx_matrix_t * m)370 static int zbx_fill_independent(double *t, int n, zbx_fit_t fit, int k, zbx_matrix_t *m)
371 {
372 double element;
373 int i, j;
374
375 if (FIT_LINEAR == fit || FIT_EXPONENTIAL == fit)
376 {
377 if (SUCCEED != zbx_matrix_alloc(m, n, 2))
378 return FAIL;
379
380 for (i = 0; i < n; i++)
381 {
382 ZBX_MATRIX_EL(m, i, 0) = 1.0;
383 ZBX_MATRIX_EL(m, i, 1) = t[i];
384 }
385 }
386 else if (FIT_LOGARITHMIC == fit || FIT_POWER == fit)
387 {
388 if (SUCCEED != zbx_matrix_alloc(m, n, 2))
389 return FAIL;
390
391 for (i = 0; i < n; i++)
392 {
393 ZBX_MATRIX_EL(m, i, 0) = 1.0;
394 ZBX_MATRIX_EL(m, i, 1) = log(t[i]);
395 }
396 }
397 else if (FIT_POLYNOMIAL == fit)
398 {
399 if (k > n - 1)
400 k = n - 1;
401
402 if (SUCCEED != zbx_matrix_alloc(m, n, k+1))
403 return FAIL;
404
405 for (i = 0; i < n; i++)
406 {
407 element = 1.0;
408
409 for (j = 0; j < k; j++)
410 {
411 ZBX_MATRIX_EL(m, i, j) = element;
412 element *= t[i];
413 }
414
415 ZBX_MATRIX_EL(m, i, k) = element;
416 }
417 }
418
419 return SUCCEED;
420 }
421
zbx_regression(double * t,double * x,int n,zbx_fit_t fit,int k,zbx_matrix_t * coefficients)422 static int zbx_regression(double *t, double *x, int n, zbx_fit_t fit, int k, zbx_matrix_t *coefficients)
423 {
424 zbx_matrix_t *independent = NULL, *dependent = NULL;
425 int res;
426
427 zbx_matrix_struct_alloc(&independent);
428 zbx_matrix_struct_alloc(&dependent);
429
430 if (SUCCEED != (res = zbx_fill_independent(t, n, fit, k, independent)))
431 goto out;
432
433 if (SUCCEED != (res = zbx_fill_dependent(x, n, fit, dependent)))
434 goto out;
435
436 if (SUCCEED != (res = zbx_least_squares(independent, dependent, coefficients)))
437 goto out;
438
439 out:
440 zbx_matrix_free(independent);
441 zbx_matrix_free(dependent);
442 return res;
443 }
444
zbx_polynomial_value(double t,zbx_matrix_t * coefficients)445 static double zbx_polynomial_value(double t, zbx_matrix_t *coefficients)
446 {
447 double pow = 1.0, res = 0.0;
448 int i;
449
450 for (i = 0; i < coefficients->rows; i++, pow *= t)
451 res += ZBX_MATRIX_EL(coefficients, i, 0) * pow;
452
453 return res;
454 }
455
zbx_polynomial_antiderivative(double t,zbx_matrix_t * coefficients)456 static double zbx_polynomial_antiderivative(double t, zbx_matrix_t *coefficients)
457 {
458 double pow = t, res = 0.0;
459 int i;
460
461 for (i = 0; i < coefficients->rows; i++, pow *= t)
462 res += ZBX_MATRIX_EL(coefficients, i, 0) * pow / (i + 1);
463
464 return res;
465 }
466
zbx_derive_polynomial(zbx_matrix_t * polynomial,zbx_matrix_t * derivative)467 static int zbx_derive_polynomial(zbx_matrix_t *polynomial, zbx_matrix_t *derivative)
468 {
469 int i;
470
471 if (!ZBX_VALID_MATRIX(polynomial))
472 goto error;
473
474 if (SUCCEED != zbx_matrix_alloc(derivative, (polynomial->rows > 1 ? polynomial->rows - 1 : 1), 1))
475 return FAIL;
476
477 for (i = 1; i < polynomial->rows; i++)
478 ZBX_MATRIX_EL(derivative, i - 1, 0) = ZBX_MATRIX_EL(polynomial, i, 0) * i;
479
480 if (1 == i)
481 ZBX_MATRIX_EL(derivative, 0, 0) = 0.0;
482
483 return SUCCEED;
484 error:
485 THIS_SHOULD_NEVER_HAPPEN;
486 return FAIL;
487 }
488
zbx_polynomial_roots(zbx_matrix_t * coefficients,zbx_matrix_t * roots)489 static int zbx_polynomial_roots(zbx_matrix_t *coefficients, zbx_matrix_t *roots)
490 {
491 #define Re(z) (z)[0]
492 #define Im(z) (z)[1]
493
494 #define ZBX_COMPLEX_MULT(z1, z2, tmp) \
495 do \
496 { \
497 Re(tmp) = Re(z1) * Re(z2) - Im(z1) * Im(z2); \
498 Im(tmp) = Re(z1) * Im(z2) + Im(z1) * Re(z2); \
499 Re(z1) = Re(tmp); \
500 Im(z1) = Im(tmp); \
501 } \
502 while(0)
503
504 #define ZBX_MAX_ITERATIONS 200
505
506 zbx_matrix_t *denominator_multiplicands = NULL, *updates = NULL;
507 double z[2], mult[2], denominator[2], zpower[2], polynomial[2], highest_degree_coefficient,
508 lower_bound, upper_bound, radius, max_update, min_distance, residual, temp;
509 int i, j, degree, first_nonzero, res, iteration = 0, roots_ok = 0, root_init = 0;
510
511 if (!ZBX_VALID_MATRIX(coefficients))
512 goto error;
513
514 degree = coefficients->rows - 1;
515 highest_degree_coefficient = ZBX_MATRIX_EL(coefficients, degree, 0);
516
517 while (0.0 == highest_degree_coefficient && 0 < degree)
518 highest_degree_coefficient = ZBX_MATRIX_EL(coefficients, --degree, 0);
519
520 if (0 == degree)
521 {
522 /* please check explicitly for an attempt to solve equation 0 == 0 */
523 if (0.0 == highest_degree_coefficient)
524 goto error;
525
526 return SUCCEED;
527 }
528
529 if (1 == degree)
530 {
531 if (SUCCEED != zbx_matrix_alloc(roots, 1, 2))
532 return FAIL;
533
534 Re(ZBX_MATRIX_ROW(roots, 0)) = -ZBX_MATRIX_EL(coefficients, 0, 0) / ZBX_MATRIX_EL(coefficients, 1, 0);
535 Im(ZBX_MATRIX_ROW(roots, 0)) = 0.0;
536
537 return SUCCEED;
538 }
539
540 if (2 == degree)
541 {
542 if (SUCCEED != zbx_matrix_alloc(roots, 2, 2))
543 return FAIL;
544
545 if (0.0 < (temp = ZBX_MATRIX_EL(coefficients, 1, 0) * ZBX_MATRIX_EL(coefficients, 1, 0) -
546 4 * ZBX_MATRIX_EL(coefficients, 2, 0) * ZBX_MATRIX_EL(coefficients, 0, 0)))
547 {
548 temp = (0 < ZBX_MATRIX_EL(coefficients, 1, 0) ?
549 -ZBX_MATRIX_EL(coefficients, 1, 0) - sqrt(temp) :
550 -ZBX_MATRIX_EL(coefficients, 1, 0) + sqrt(temp));
551 Re(ZBX_MATRIX_ROW(roots, 0)) = 0.5 * temp / ZBX_MATRIX_EL(coefficients, 2, 0);
552 Re(ZBX_MATRIX_ROW(roots, 1)) = 2.0 * ZBX_MATRIX_EL(coefficients, 0, 0) / temp;
553 Im(ZBX_MATRIX_ROW(roots, 0)) = Im(ZBX_MATRIX_ROW(roots, 1)) = 0.0;
554 }
555 else
556 {
557 Re(ZBX_MATRIX_ROW(roots, 0)) = Re(ZBX_MATRIX_ROW(roots, 1)) =
558 -0.5 * ZBX_MATRIX_EL(coefficients, 1, 0) / ZBX_MATRIX_EL(coefficients, 2, 0);
559 Im(ZBX_MATRIX_ROW(roots, 0)) = -(Im(ZBX_MATRIX_ROW(roots, 1)) = 0.5 * sqrt(-temp)) /
560 ZBX_MATRIX_EL(coefficients, 2, 0);
561 }
562
563 return SUCCEED;
564 }
565
566 zbx_matrix_struct_alloc(&denominator_multiplicands);
567 zbx_matrix_struct_alloc(&updates);
568
569 if (SUCCEED != zbx_matrix_alloc(roots, degree, 2) ||
570 SUCCEED != zbx_matrix_alloc(denominator_multiplicands, degree, 2) ||
571 SUCCEED != zbx_matrix_alloc(updates, degree, 2))
572 {
573 res = FAIL;
574 goto out;
575 }
576
577 /* if n lower coefficients are zeros, zero is a root of multiplicity n */
578 for (first_nonzero = 0; 0.0 == ZBX_MATRIX_EL(coefficients, first_nonzero, 0); first_nonzero++)
579 Re(ZBX_MATRIX_ROW(roots, first_nonzero)) = Im(ZBX_MATRIX_ROW(roots, first_nonzero)) = 0.0;
580
581 /* compute bounds for the roots */
582 upper_bound = lower_bound = 1.0;
583
584 for (i = first_nonzero; i < degree; i++)
585 {
586 if (upper_bound < fabs(ZBX_MATRIX_EL(coefficients, i, 0) / highest_degree_coefficient))
587 upper_bound = fabs(ZBX_MATRIX_EL(coefficients, i, 0) / highest_degree_coefficient);
588
589 if (lower_bound < fabs(ZBX_MATRIX_EL(coefficients, i + 1, 0) /
590 ZBX_MATRIX_EL(coefficients, first_nonzero, 0)))
591 lower_bound = fabs(ZBX_MATRIX_EL(coefficients, i + 1, 0) /
592 ZBX_MATRIX_EL(coefficients, first_nonzero, 0));
593 }
594
595 radius = 1.0 / lower_bound;
596
597 /* Weierstrass (Durand-Kerner) method */
598 while (ZBX_MAX_ITERATIONS >= ++iteration && !roots_ok)
599 {
600 if (0 == root_init)
601 {
602 if (radius <= upper_bound)
603 {
604 for (i = 0; i < degree - first_nonzero; i++)
605 {
606 Re(ZBX_MATRIX_ROW(roots, i)) = radius * cos((2.0 * M_PI * (i + 0.25)) /
607 (degree - first_nonzero));
608 Im(ZBX_MATRIX_ROW(roots, i)) = radius * sin((2.0 * M_PI * (i + 0.25)) /
609 (degree - first_nonzero));
610 }
611
612 radius *= 2.0;
613 }
614 else
615 root_init = 1;
616 }
617
618 roots_ok = 1;
619 max_update = 0.0;
620 min_distance = HUGE_VAL;
621
622 for (i = first_nonzero; i < degree; i++)
623 {
624 Re(z) = Re(ZBX_MATRIX_ROW(roots, i));
625 Im(z) = Im(ZBX_MATRIX_ROW(roots, i));
626
627 /* subtract from z every one of denominator_multiplicands and multiplicate them */
628 Re(denominator) = highest_degree_coefficient;
629 Im(denominator) = 0.0;
630
631 for (j = first_nonzero; j < degree; j++)
632 {
633 if (j == i)
634 continue;
635
636 temp = (ZBX_MATRIX_EL(roots, i, 0) - ZBX_MATRIX_EL(roots, j, 0)) *
637 (ZBX_MATRIX_EL(roots, i, 0) - ZBX_MATRIX_EL(roots, j, 0)) +
638 (ZBX_MATRIX_EL(roots, i, 1) - ZBX_MATRIX_EL(roots, j, 1)) *
639 (ZBX_MATRIX_EL(roots, i, 1) - ZBX_MATRIX_EL(roots, j, 1));
640 if (temp < min_distance)
641 min_distance = temp;
642
643 Re(ZBX_MATRIX_ROW(denominator_multiplicands, j)) = Re(z) - Re(ZBX_MATRIX_ROW(roots, j));
644 Im(ZBX_MATRIX_ROW(denominator_multiplicands, j)) = Im(z) - Im(ZBX_MATRIX_ROW(roots, j));
645 ZBX_COMPLEX_MULT(denominator, ZBX_MATRIX_ROW(denominator_multiplicands, j), mult);
646 }
647
648 /* calculate complex value of polynomial for z */
649 Re(zpower) = 1.0;
650 Im(zpower) = 0.0;
651 Re(polynomial) = ZBX_MATRIX_EL(coefficients, first_nonzero, 0);
652 Im(polynomial) = 0.0;
653
654 for (j = first_nonzero + 1; j <= degree; j++)
655 {
656 ZBX_COMPLEX_MULT(zpower, z, mult);
657 Re(polynomial) += Re(zpower) * ZBX_MATRIX_EL(coefficients, j, 0);
658 Im(polynomial) += Im(zpower) * ZBX_MATRIX_EL(coefficients, j, 0);
659 }
660
661 /* check how good root approximation is */
662 residual = fabs(Re(polynomial)) + fabs(Im(polynomial));
663 roots_ok = roots_ok && (ZBX_MATH_EPSILON > residual);
664
665 /* divide polynomial value by denominator */
666 if (0.0 != (temp = Re(denominator) * Re(denominator) + Im(denominator) * Im(denominator)))
667 {
668 Re(ZBX_MATRIX_ROW(updates, i)) = (Re(polynomial) * Re(denominator) +
669 Im(polynomial) * Im(denominator)) / temp;
670 Im(ZBX_MATRIX_ROW(updates, i)) = (Im(polynomial) * Re(denominator) -
671 Re(polynomial) * Im(denominator)) / temp;
672 }
673 else /* Denominator is zero iff two or more root approximations are equal. */
674 /* Since root approximations are initially different their equality means that they */
675 /* converged to a multiple root (hopefully) and no updates are required in this case. */
676 {
677 Re(ZBX_MATRIX_ROW(updates, i)) = Im(ZBX_MATRIX_ROW(updates, i)) = 0.0;
678 }
679
680 temp = ZBX_MATRIX_EL(updates, i, 0) * ZBX_MATRIX_EL(updates, i, 0) +
681 ZBX_MATRIX_EL(updates, i, 1) * ZBX_MATRIX_EL(updates, i, 1);
682
683 if (temp > max_update)
684 max_update = temp;
685 }
686
687 if (max_update > radius * radius && 0 == root_init)
688 continue;
689 else
690 root_init = 1;
691
692 for (i = first_nonzero; i < degree; i++)
693 {
694 Re(ZBX_MATRIX_ROW(roots, i)) -= Re(ZBX_MATRIX_ROW(updates, i));
695 Im(ZBX_MATRIX_ROW(roots, i)) -= Im(ZBX_MATRIX_ROW(updates, i));
696 }
697 }
698
699 if (0 == roots_ok)
700 {
701 zabbix_log(LOG_LEVEL_DEBUG, "polynomial root finding problem is ill-defined");
702 res = FAIL;
703 }
704 else
705 res = SUCCEED;
706
707 out:
708 zbx_matrix_free(denominator_multiplicands);
709 zbx_matrix_free(updates);
710 return res;
711 error:
712 THIS_SHOULD_NEVER_HAPPEN;
713 return FAIL;
714
715 #undef ZBX_MAX_ITERATIONS
716
717 #undef Re
718 #undef Im
719 }
720
zbx_polynomial_minmax(double now,double time,zbx_mode_t mode,zbx_matrix_t * coefficients,double * result)721 static int zbx_polynomial_minmax(double now, double time, zbx_mode_t mode, zbx_matrix_t *coefficients,
722 double *result)
723 {
724 zbx_matrix_t *derivative = NULL, *derivative_roots = NULL;
725 double min, max, tmp;
726 int i, res;
727
728 if (!ZBX_VALID_MATRIX(coefficients))
729 goto error;
730
731 zbx_matrix_struct_alloc(&derivative);
732 zbx_matrix_struct_alloc(&derivative_roots);
733
734 if (SUCCEED != (res = zbx_derive_polynomial(coefficients, derivative)))
735 goto out;
736
737 if (SUCCEED != (res = zbx_polynomial_roots(derivative, derivative_roots)))
738 goto out;
739
740 /* choose min and max among now, now + time and derivative roots inbetween (these are potential local extrema) */
741 /* we ignore imaginary part of roots, this means that more calculations will be made, */
742 /* but result will not be affected and we wont need a boundary on minimal imaginary part that differs from zero */
743
744 min = zbx_polynomial_value(now, coefficients);
745 tmp = zbx_polynomial_value(now + time, coefficients);
746
747 if (tmp < min)
748 {
749 max = min;
750 min = tmp;
751 }
752 else
753 max = tmp;
754
755 for (i = 0; i < derivative_roots->rows; i++)
756 {
757 tmp = ZBX_MATRIX_EL(derivative_roots, i, 0);
758
759 if (tmp < now || tmp > now + time)
760 continue;
761
762 tmp = zbx_polynomial_value(tmp, coefficients);
763
764 if (tmp < min)
765 min = tmp;
766 else if (tmp > max)
767 max = tmp;
768 }
769
770 if (MODE_MAX == mode)
771 *result = max;
772 else if (MODE_MIN == mode)
773 *result = min;
774 else if (MODE_DELTA == mode)
775 *result = max - min;
776 else
777 THIS_SHOULD_NEVER_HAPPEN;
778
779 out:
780 zbx_matrix_free(derivative);
781 zbx_matrix_free(derivative_roots);
782 return res;
783 error:
784 THIS_SHOULD_NEVER_HAPPEN;
785 return FAIL;
786 }
787
zbx_polynomial_timeleft(double now,double threshold,zbx_matrix_t * coefficients,double * result)788 static int zbx_polynomial_timeleft(double now, double threshold, zbx_matrix_t *coefficients, double *result)
789 {
790 zbx_matrix_t *shifted_coefficients = NULL, *roots = NULL;
791 double tmp;
792 int i, res, no_root = 1;
793
794 if (!ZBX_VALID_MATRIX(coefficients))
795 goto error;
796
797 zbx_matrix_struct_alloc(&shifted_coefficients);
798 zbx_matrix_struct_alloc(&roots);
799
800 if (SUCCEED != (res = zbx_matrix_copy(shifted_coefficients, coefficients)))
801 goto out;
802
803 ZBX_MATRIX_EL(shifted_coefficients, 0, 0) -= threshold;
804
805 if (SUCCEED != (res = zbx_polynomial_roots(shifted_coefficients, roots)))
806 goto out;
807
808 /* choose the closest root right from now or set result to -1 otherwise */
809 /* if zbx_polynomial_value(tmp) is not close enough to zero it must be a complex root and must be skipped */
810
811 for (i = 0; i < roots->rows; i++)
812 {
813 tmp = ZBX_MATRIX_EL(roots, i, 0);
814
815 if (no_root)
816 {
817 if (tmp > now && ZBX_MATH_EPSILON > fabs(zbx_polynomial_value(tmp, shifted_coefficients)))
818 {
819 no_root = 0;
820 *result = tmp;
821 }
822 }
823 else if (now < tmp && tmp < *result &&
824 ZBX_MATH_EPSILON > fabs(zbx_polynomial_value(tmp, shifted_coefficients)))
825 {
826 *result = tmp;
827 }
828 }
829
830 if (no_root)
831 *result = DB_INFINITY;
832 else
833 *result -= now;
834
835 out:
836 zbx_matrix_free(shifted_coefficients);
837 zbx_matrix_free(roots);
838 return res;
839 error:
840 THIS_SHOULD_NEVER_HAPPEN;
841 return FAIL;
842 }
843
zbx_calculate_value(double t,zbx_matrix_t * coefficients,zbx_fit_t fit,double * value)844 static int zbx_calculate_value(double t, zbx_matrix_t *coefficients, zbx_fit_t fit, double *value)
845 {
846 if (!ZBX_VALID_MATRIX(coefficients))
847 goto error;
848
849 if (FIT_LINEAR == fit)
850 *value = ZBX_MATRIX_EL(coefficients, 0, 0) + ZBX_MATRIX_EL(coefficients, 1, 0) * t;
851 else if (FIT_POLYNOMIAL == fit)
852 *value = zbx_polynomial_value(t, coefficients);
853 else if (FIT_EXPONENTIAL == fit)
854 *value = exp(ZBX_MATRIX_EL(coefficients, 0, 0) + ZBX_MATRIX_EL(coefficients, 1, 0) * t);
855 else if (FIT_LOGARITHMIC == fit)
856 *value = ZBX_MATRIX_EL(coefficients, 0, 0) + ZBX_MATRIX_EL(coefficients, 1, 0) * log(t);
857 else if (FIT_POWER == fit)
858 *value = exp(ZBX_MATRIX_EL(coefficients, 0, 0) + ZBX_MATRIX_EL(coefficients, 1, 0) * log(t));
859 else
860 goto error;
861
862 return SUCCEED;
863 error:
864 THIS_SHOULD_NEVER_HAPPEN;
865 return FAIL;
866 }
867
zbx_fit_code(char * fit_str,zbx_fit_t * fit,unsigned * k,char ** error)868 int zbx_fit_code(char *fit_str, zbx_fit_t *fit, unsigned *k, char **error)
869 {
870 if ('\0' == *fit_str || 0 == strcmp(fit_str, "linear"))
871 {
872 *fit = FIT_LINEAR;
873 *k = 0;
874 }
875 else if (0 == strncmp(fit_str, "polynomial", strlen("polynomial")))
876 {
877 *fit = FIT_POLYNOMIAL;
878
879 if (SUCCEED != is_uint_range(fit_str + strlen("polynomial"), k, 1, 6))
880 {
881 *error = zbx_strdup(*error, "polynomial degree is invalid");
882 return FAIL;
883 }
884 }
885 else if (0 == strcmp(fit_str, "exponential"))
886 {
887 *fit = FIT_EXPONENTIAL;
888 *k = 0;
889 }
890 else if (0 == strcmp(fit_str, "logarithmic"))
891 {
892 *fit = FIT_LOGARITHMIC;
893 *k = 0;
894 }
895 else if (0 == strcmp(fit_str, "power"))
896 {
897 *fit = FIT_POWER;
898 *k = 0;
899 }
900 else
901 {
902 *error = zbx_strdup(*error, "invalid 'fit' parameter");
903 return FAIL;
904 }
905
906 return SUCCEED;
907 }
908
zbx_mode_code(char * mode_str,zbx_mode_t * mode,char ** error)909 int zbx_mode_code(char *mode_str, zbx_mode_t *mode, char **error)
910 {
911 if ('\0' == *mode_str || 0 == strcmp(mode_str, "value"))
912 {
913 *mode = MODE_VALUE;
914 }
915 else if (0 == strcmp(mode_str, "max"))
916 {
917 *mode = MODE_MAX;
918 }
919 else if (0 == strcmp(mode_str, "min"))
920 {
921 *mode = MODE_MIN;
922 }
923 else if (0 == strcmp(mode_str, "delta"))
924 {
925 *mode = MODE_DELTA;
926 }
927 else if (0 == strcmp(mode_str, "avg"))
928 {
929 *mode = MODE_AVG;
930 }
931 else
932 {
933 *error = zbx_strdup(*error, "invalid 'mode' parameter");
934 return FAIL;
935 }
936
937 return SUCCEED;
938 }
939
zbx_log_expression(double now,zbx_fit_t fit,int k,zbx_matrix_t * coeffs)940 static void zbx_log_expression(double now, zbx_fit_t fit, int k, zbx_matrix_t *coeffs)
941 {
942 /* x is item value, t is time in seconds counted from now */
943 if (FIT_LINEAR == fit)
944 {
945 zabbix_log(LOG_LEVEL_DEBUG, "fitted expression is: x = (" ZBX_FS_DBL ") + (" ZBX_FS_DBL ") * (" ZBX_FS_DBL " + t)",
946 ZBX_MATRIX_EL(coeffs, 0, 0), ZBX_MATRIX_EL(coeffs, 1, 0), now);
947 }
948 else if (FIT_POLYNOMIAL == fit)
949 {
950 char *polynomial = NULL;
951 size_t alloc, offset;
952
953 while (0 <= k)
954 {
955 zbx_snprintf_alloc(&polynomial, &alloc, &offset, "(" ZBX_FS_DBL ") * (" ZBX_FS_DBL " + t) ^ %d",
956 ZBX_MATRIX_EL(coeffs, k, 0), now, k);
957
958 if (0 < k--)
959 zbx_snprintf_alloc(&polynomial, &alloc, &offset, " + ");
960 }
961
962 zabbix_log(LOG_LEVEL_DEBUG, "fitted expression is: x = %s", polynomial);
963
964 zbx_free(polynomial);
965 }
966 else if (FIT_EXPONENTIAL == fit)
967 {
968 zabbix_log(LOG_LEVEL_DEBUG, "fitted expression is: x = (" ZBX_FS_DBL ") * exp( (" ZBX_FS_DBL ") * (" ZBX_FS_DBL " + t) )",
969 exp(ZBX_MATRIX_EL(coeffs, 0, 0)), ZBX_MATRIX_EL(coeffs, 1, 0), now);
970 }
971 else if (FIT_LOGARITHMIC == fit)
972 {
973 zabbix_log(LOG_LEVEL_DEBUG, "fitted expression is: x = (" ZBX_FS_DBL ") + (" ZBX_FS_DBL ") * log(" ZBX_FS_DBL " + t)",
974 ZBX_MATRIX_EL(coeffs, 0, 0), ZBX_MATRIX_EL(coeffs, 1, 0), now);
975 }
976 else if (FIT_POWER == fit)
977 {
978 zabbix_log(LOG_LEVEL_DEBUG, "fitted expression is: x = (" ZBX_FS_DBL ") * (" ZBX_FS_DBL " + t) ^ (" ZBX_FS_DBL ")",
979 exp(ZBX_MATRIX_EL(coeffs, 0, 0)), now, ZBX_MATRIX_EL(coeffs, 1, 0));
980 }
981 else
982 THIS_SHOULD_NEVER_HAPPEN;
983 }
984
zbx_forecast(double * t,double * x,int n,double now,double time,zbx_fit_t fit,unsigned k,zbx_mode_t mode)985 double zbx_forecast(double *t, double *x, int n, double now, double time, zbx_fit_t fit, unsigned k, zbx_mode_t mode)
986 {
987 zbx_matrix_t *coefficients = NULL;
988 double left, right, result;
989 int res;
990
991 if (1 == n)
992 {
993 if (MODE_VALUE == mode || MODE_MAX == mode || MODE_MIN == mode || MODE_AVG == mode)
994 return x[0];
995
996 if (MODE_DELTA == mode)
997 return 0.0;
998
999 THIS_SHOULD_NEVER_HAPPEN;
1000 return ZBX_MATH_ERROR;
1001 }
1002
1003 zbx_matrix_struct_alloc(&coefficients);
1004
1005 if (SUCCEED != (res = zbx_regression(t, x, n, fit, k, coefficients)))
1006 goto out;
1007
1008 zbx_log_expression(now, fit, (int)k, coefficients);
1009
1010 if (MODE_VALUE == mode)
1011 {
1012 res = zbx_calculate_value(now + time, coefficients, fit, &result);
1013 goto out;
1014 }
1015
1016 if (0.0 == time)
1017 {
1018 if (MODE_MAX == mode || MODE_MIN == mode || MODE_AVG == mode)
1019 {
1020 res = zbx_calculate_value(now + time, coefficients, fit, &result);
1021 }
1022 else if (MODE_DELTA == mode)
1023 {
1024 result = 0.0;
1025 res = SUCCEED;
1026 }
1027 else
1028 {
1029 THIS_SHOULD_NEVER_HAPPEN;
1030 res = FAIL;
1031 }
1032
1033 goto out;
1034 }
1035
1036 if (FIT_LINEAR == fit || FIT_EXPONENTIAL == fit || FIT_LOGARITHMIC == fit || FIT_POWER == fit)
1037 {
1038 /* fit is monotone, therefore maximum and minimum are either at now or at now + time */
1039 if (SUCCEED != zbx_calculate_value(now, coefficients, fit, &left) ||
1040 SUCCEED != zbx_calculate_value(now + time, coefficients, fit, &right))
1041 {
1042 res = FAIL;
1043 goto out;
1044 }
1045
1046 if (MODE_MAX == mode)
1047 {
1048 result = (left > right ? left : right);
1049 }
1050 else if (MODE_MIN == mode)
1051 {
1052 result = (left < right ? left : right);
1053 }
1054 else if (MODE_DELTA == mode)
1055 {
1056 result = (left > right ? left - right : right - left);
1057 }
1058 else if (MODE_AVG == mode)
1059 {
1060 if (FIT_LINEAR == fit)
1061 {
1062 result = 0.5 * (left + right);
1063 }
1064 else if (FIT_EXPONENTIAL == fit)
1065 {
1066 result = (right - left) / time / ZBX_MATRIX_EL(coefficients, 1, 0);
1067 }
1068 else if (FIT_LOGARITHMIC == fit)
1069 {
1070 result = right + ZBX_MATRIX_EL(coefficients, 1, 0) *
1071 (log(1.0 + time / now) * now / time - 1.0);
1072 }
1073 else if (FIT_POWER == fit)
1074 {
1075 if (-1.0 != ZBX_MATRIX_EL(coefficients, 1, 0))
1076 result = (right * (now + time) - left * now) / time /
1077 (ZBX_MATRIX_EL(coefficients, 1, 0) + 1.0);
1078 else
1079 result = exp(ZBX_MATRIX_EL(coefficients, 0, 0)) * log(1.0 + time / now) / time;
1080 }
1081 else
1082 {
1083 THIS_SHOULD_NEVER_HAPPEN;
1084 res = FAIL;
1085 goto out;
1086 }
1087 }
1088 else
1089 {
1090 THIS_SHOULD_NEVER_HAPPEN;
1091 res = FAIL;
1092 goto out;
1093 }
1094
1095 res = SUCCEED;
1096 }
1097 else if (FIT_POLYNOMIAL == fit)
1098 {
1099 if (MODE_MAX == mode || MODE_MIN == mode || MODE_DELTA == mode)
1100 {
1101 res = zbx_polynomial_minmax(now, time, mode, coefficients, &result);
1102 }
1103 else if (MODE_AVG == mode)
1104 {
1105 result = (zbx_polynomial_antiderivative(now + time, coefficients) -
1106 zbx_polynomial_antiderivative(now, coefficients)) / time;
1107 res = SUCCEED;
1108 }
1109 else
1110 {
1111 THIS_SHOULD_NEVER_HAPPEN;
1112 res = FAIL;
1113 }
1114 }
1115 else
1116 {
1117 THIS_SHOULD_NEVER_HAPPEN;
1118 res = FAIL;
1119 }
1120
1121 out:
1122 zbx_matrix_free(coefficients);
1123
1124 if (SUCCEED != res)
1125 {
1126 result = ZBX_MATH_ERROR;
1127 }
1128 else if (ZBX_IS_NAN(result))
1129 {
1130 zabbix_log(LOG_LEVEL_DEBUG, "numerical error");
1131 result = ZBX_MATH_ERROR;
1132 }
1133 else if (DB_INFINITY < result)
1134 {
1135 result = DB_INFINITY;
1136 }
1137 else if (-DB_INFINITY > result)
1138 {
1139 result = -DB_INFINITY;
1140 }
1141
1142 return result;
1143 }
1144
zbx_timeleft(double * t,double * x,int n,double now,double threshold,zbx_fit_t fit,unsigned k)1145 double zbx_timeleft(double *t, double *x, int n, double now, double threshold, zbx_fit_t fit, unsigned k)
1146 {
1147 zbx_matrix_t *coefficients = NULL;
1148 double current, result;
1149 int res;
1150
1151 if (1 == n)
1152 return (x[0] == threshold ? 0.0 : DB_INFINITY);
1153
1154 zbx_matrix_struct_alloc(&coefficients);
1155
1156 if (SUCCEED != (res = zbx_regression(t, x, n, fit, k, coefficients)))
1157 goto out;
1158
1159 zbx_log_expression(now, fit, (int)k, coefficients);
1160
1161 if (SUCCEED != (res = zbx_calculate_value(now, coefficients, fit, ¤t)))
1162 {
1163 goto out;
1164 }
1165 else if (current == threshold)
1166 {
1167 result = 0.0;
1168 goto out;
1169 }
1170
1171 if (FIT_LINEAR == fit)
1172 {
1173 result = (threshold - ZBX_MATRIX_EL(coefficients, 0, 0)) / ZBX_MATRIX_EL(coefficients, 1, 0) - now;
1174 }
1175 else if (FIT_POLYNOMIAL == fit)
1176 {
1177 res = zbx_polynomial_timeleft(now, threshold, coefficients, &result);
1178 }
1179 else if (FIT_EXPONENTIAL == fit)
1180 {
1181 result = (log(threshold) - ZBX_MATRIX_EL(coefficients, 0, 0)) / ZBX_MATRIX_EL(coefficients, 1, 0) - now;
1182 }
1183 else if (FIT_LOGARITHMIC == fit)
1184 {
1185 result = exp((threshold - ZBX_MATRIX_EL(coefficients, 0, 0)) / ZBX_MATRIX_EL(coefficients, 1, 0)) - now;
1186 }
1187 else if (FIT_POWER == fit)
1188 {
1189 result = exp((log(threshold) - ZBX_MATRIX_EL(coefficients, 0, 0)) / ZBX_MATRIX_EL(coefficients, 1, 0))
1190 - now;
1191 }
1192 else
1193 {
1194 THIS_SHOULD_NEVER_HAPPEN;
1195 res = FAIL;
1196 }
1197
1198 out:
1199 if (SUCCEED != res)
1200 {
1201 result = ZBX_MATH_ERROR;
1202 }
1203 else if (ZBX_IS_NAN(result))
1204 {
1205 zabbix_log(LOG_LEVEL_DEBUG, "numerical error");
1206 result = ZBX_MATH_ERROR;
1207 }
1208 else if (0.0 > result || DB_INFINITY < result)
1209 {
1210 result = DB_INFINITY;
1211 }
1212
1213 zbx_matrix_free(coefficients);
1214 return result;
1215 }
1216