There are two popular variants of the Holt-Winters forecasting method; RRDtool

supports the "additive" method, which means that seasonal variation is simply
added to the baseline. For our application, it would be more appropriate to use
the "multiplicative" Holt-Winters method, where seasonal variation is a
coefficient multiplied by the baseline.  Quick example to illustrate the
difference: if the average doubles season-over-season, the additive method
would predict the delta between min and max to be constant, whereas the
multiplicative method would predict the delta to double as well.

Attached is a patch against trunk to support the multiplicative method.  I've
done this with a new consolidation function, MHWPREDICT, which is essentially
interchangeable with HWPREDICT. There is a noticeable improvement in prediction
deviations for certain types of functions; the attachments show HWPREDICT and
MHWPREDICT predictions for a function with an x*sin(x) component.

Because HWPREDICT and MHWPREDICT differ only in their equations, I've factored
out their math into rrd_hw_math.c. The appropriate smoothing functions are
passed to the update functions in a container of function pointers, which are
called where appropriate. Thus the additive and multiplicative methods use the
same update functions, and the right equations are evaluated without having
flag checks everywhere. This approach, I think, makes the algorithms quite
clear, with minimal duplicate code.

I have moved update_hwpredict, update_seasonal, update_devpredict,
update_devseasonal, and update_failures into a separate file, rrd_hw_update.c,
with some slight refactoring related to rrd_hw_math.c. I ran some
regression tests against trunk to make sure I didn't break anything with
the existing HWPREDICT code.

MHWPREDICT uses the same deviation smoothing and failure detection algorithms
as HWPREDICT.

Some helpful references on the multiplicative Holt-Winters method:

http://www.it.iitb.ac.in/~praj/acads/seminar/04329008_ExponentialSmoothing.pdf
(a student's quick overview of additive vs. multiplicative HW)

http://ideas.repec.org/p/msh/ebswps/1999-1.html (paper on variations to the
multiplicative Holt-Winters, including variance calculations; FYI, my
implementation uses "Model 1")

My employer and the owner of this patch (IMVU, Inc.) is happy to license it
under the same terms as RRDtool, i.e. give it to the project.
-- Evan Miller emiller imvu.com

git-svn-id: svn://svn.oetiker.ch/rrdtool/trunk/program@1125 a5681a0c-68f1-0310-ab6d-d61299d08faa

diff --git a/doc/rrdcreate.pod b/doc/rrdcreate.pod
index 27ef702afd5097ca94e08e13c2b49cb01c3d0559..3689e2937aa09a158bebe6c20de60c4ca5cbfc4e 100644 (file)
--- a/doc/rrdcreate.pod
+++ b/doc/rrdcreate.pod
@@ -206,6 +206,10 @@ B<RRA:>I<HWPREDICT>B<:>I<rows>B<:>I<alpha>B<:>I<beta>B<:>I<seasonal period>[B<:>
 
 =item *
 
+B<RRA:>I<MHWPREDICT>B<:>I<rows>B<:>I<alpha>B<:>I<beta>B<:>I<seasonal period>[B<:>I<rra-num>]
+
+=item *
+
 B<RRA:>I<SEASONAL>B<:>I<seasonal period>B<:>I<gamma>B<:>I<rra-num>
 
 =item *
@@ -225,19 +229,32 @@ B<RRA:>I<FAILURES>B<:>I<rows>B<:>I<threshold>B<:>I<window length>B<:>I<rra-num>
 These B<RRAs> differ from the true consolidation functions in several ways.
 First, each of the B<RRA>s is updated once for every primary data point.
 Second, these B<RRAs> are interdependent. To generate real-time confidence
-bounds, a matched set of HWPREDICT, SEASONAL, DEVSEASONAL, and
-DEVPREDICT must exist. Generating smoothed values of the primary data points
-requires both a HWPREDICT B<RRA> and SEASONAL B<RRA>. Aberrant behavior
-detection requires FAILURES, HWPREDICT, DEVSEASONAL, and SEASONAL.
-
-The actual predicted, or smoothed, values are stored in the HWPREDICT
-B<RRA>. The predicted deviations are stored in DEVPREDICT (think a standard
-deviation which can be scaled to yield a confidence band). The FAILURES
-B<RRA> stores binary indicators. A 1 marks the indexed observation as
-failure; that is, the number of confidence bounds violations in the
-preceding window of observations met or exceeded a specified threshold. An
-example of using these B<RRAs> to graph confidence bounds and failures
-appears in L<rrdgraph>.
+bounds, a matched set of SEASONAL, DEVSEASONAL, DEVPREDICT, and either
+HWPREDICT or MHWPREDICT must exist. Generating smoothed values of the primary
+data points requires a SEASONAL B<RRA> and either an HWPREDICT or MHWPREDICT 
+B<RRA>. Aberrant behavior detection requires FAILURES, DEVSEASONAL, SEASONAL,
+and either HWPREDICT or MHWPREDICT.
+
+The predicted, or smoothed, values are stored in the HWPREDICT or MHWPREDICT
+B<RRA>. HWPREDICT and MHWPREDICT are actually two variations on the
+Holt-Winters method. They are interchangeable. Both attempt to decompose data
+into three components: a baseline, a trend, and a seasonal coefficient.
+HWPREDICT adds its seasonal coefficient to the baseline to form a prediction, whereas
+MHWPREDICT multiplies its seasonal coefficient by the baseline to form a
+prediction. The difference is noticeable when the baseline changes
+significantly in the course of a season; HWPREDICT will predict the seasonality
+to stay constant as the baseline changes, but MHWPREDICT will predict the
+seasonality to grow or shrink in proportion to the baseline. The proper choice
+of method depends on the thing being modeled. For simplicity, the rest of this
+discussion will refer to HWPREDICT, but MHWPREDICT may be substituted in its
+place.
+
+The predicted deviations are stored in DEVPREDICT (think a standard deviation
+which can be scaled to yield a confidence band). The FAILURES B<RRA> stores 
+binary indicators. A 1 marks the indexed observation as failure; that is, the 
+number of confidence bounds violations in the preceding window of observations 
+met or exceeded a specified threshold. An example of using these B<RRAs> to graph 
+confidence bounds and failures appears in L<rrdgraph>.
 
 The SEASONAL and DEVSEASONAL B<RRAs> store the seasonal coefficients for the
 Holt-Winters forecasting algorithm and the seasonal deviations, respectively.

diff --git a/doc/rrdtune.pod b/doc/rrdtune.pod
index 5e9029bf12a2772a4267c249262645e59fee2e7b..5a08394f55ecdb9bd89ab09ef6cf547ed1b55c69 100644 (file)
--- a/doc/rrdtune.pod
+++ b/doc/rrdtune.pod
@@ -121,14 +121,15 @@ B<RRA>. This parameter must be between 0 and 1.
 
 This option causes the aberrant behavior detection algorithm to reset
 for the specified data source; that is, forget all it is has learnt so far.
-Specifically, for the HWPREDICT B<RRA>, it sets the intercept and slope
-coefficients to unknown. For the SEASONAL B<RRA>, it sets all seasonal
+Specifically, for the HWPREDICT or MHWPREDICT B<RRA>, it sets the intercept and
+slope coefficients to unknown. For the SEASONAL B<RRA>, it sets all seasonal
 coefficients to unknown. For the DEVSEASONAL B<RRA>, it sets all seasonal
-deviation coefficients to unknown. For the FAILURES B<RRA>, it erases
-the violation history. Note that reset does not erase past predictions
-(the values of the HWPREDICT B<RRA>), predicted deviations (the values of the
-DEVPREDICT B<RRA>), or failure history (the values of the FAILURES B<RRA>).
-This option will function even if not all the listed B<RRAs> are present.
+deviation coefficients to unknown. For the FAILURES B<RRA>, it erases the
+violation history. Note that reset does not erase past predictions
+(the values of the HWPREDICT or MHWPREDICT B<RRA>), predicted deviations (the
+values of the DEVPREDICT B<RRA>), or failure history (the values of the 
+FAILURES B<RRA>).  This option will function even if not all the listed 
+B<RRAs> are present.
 
 Due to the implementation of this option, there is an indirect impact on
 other data sources in the RRD. A smoothing algorithm is applied to

diff --git a/src/Makefile.am b/src/Makefile.am
index 68780689db459c4b7546bd26fc58a29ec48b131f..6430f5a22ef83aa73e229323d4a7547f191aa302 100644 (file)
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -19,6 +19,8 @@ UPD_C_FILES =         \
        rrd_getopt1.c   \
        parsetime.c     \
        rrd_hw.c        \
+       rrd_hw_math.c   \
+       rrd_hw_update.c \
        rrd_diff.c      \
        rrd_format.c    \
        rrd_info.c      \
@@ -50,7 +52,8 @@ RRD_C_FILES =         \
 noinst_HEADERS = \
        unused.h \
        rrd_getopt.h parsetime.h \
-       rrd_format.h rrd_tool.h rrd_xport.h rrd.h rrd_hw.h rrd_rpncalc.h \
+       rrd_format.h rrd_tool.h rrd_xport.h rrd.h rrd_rpncalc.h \
+       rrd_hw.h rrd_hw_math.h rrd_hw_update.h \
        rrd_nan_inf.h fnv.h rrd_graph.h \
        rrd_is_thread_safe.h
 

diff --git a/src/rrd_create.c b/src/rrd_create.c
index 35a06b1d10165899ff7c505b5cc7248785b378bb..766b783460e205b7f7132b820b77a6298926209c 100644 (file)
--- a/src/rrd_create.c
+++ b/src/rrd_create.c
@@ -247,6 +247,7 @@ int rrd_create_r(
                     switch (cf_conv
                             (rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam)) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                         /* initialize some parameters */
                         rrd.rra_def[rrd.stat_head->rra_cnt].par[RRA_hw_alpha].
                             u_val = 0.1;
@@ -293,6 +294,7 @@ int rrd_create_r(
                     switch (cf_conv
                             (rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam)) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                     case CF_DEVSEASONAL:
                     case CF_SEASONAL:
                     case CF_DEVPREDICT:
@@ -316,6 +318,7 @@ int rrd_create_r(
                     switch (cf_conv
                             (rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam)) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                         rrd.rra_def[rrd.stat_head->rra_cnt].par[RRA_hw_alpha].
                             u_val = atof(token);
                         if (atof(token) <= 0.0 || atof(token) >= 1.0)
@@ -361,6 +364,7 @@ int rrd_create_r(
                     switch (cf_conv
                             (rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam)) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                         rrd.rra_def[rrd.stat_head->rra_cnt].par[RRA_hw_beta].
                             u_val = atof(token);
                         if (atof(token) < 0.0 || atof(token) > 1.0)
@@ -415,6 +419,7 @@ int rrd_create_r(
                             atoi(token) - 1;
                         break;
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                         /* length of the associated CF_SEASONAL and CF_DEVSEASONAL arrays. */
                         period = atoi(token);
                         if (period >
@@ -463,8 +468,10 @@ int rrd_create_r(
                     par[RRA_dependent_rra_idx].u_cnt, rrd.stat_head->rra_cnt);
 #endif
             /* should we create CF_SEASONAL, CF_DEVSEASONAL, and CF_DEVPREDICT? */
-            if (cf_conv(rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam) ==
-                CF_HWPREDICT
+            if ((cf_conv(rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam) ==
+                 CF_HWPREDICT
+                 || cf_conv(rrd.rra_def[rrd.stat_head->rra_cnt].cf_nam) ==
+                 CF_MHWPREDICT)
                 && rrd.rra_def[rrd.stat_head->rra_cnt].
                 par[RRA_dependent_rra_idx].u_cnt == rrd.stat_head->rra_cnt) {
 #ifdef DEBUG
@@ -671,6 +678,7 @@ int rrd_create_fn(
     for (i = 0; i < rrd->stat_head->rra_cnt; i++) {
         switch (cf_conv(rrd->rra_def[i].cf_nam)) {
         case CF_HWPREDICT:
+        case CF_MHWPREDICT:
             init_hwpredict_cdp(rrd->cdp_prep);
             break;
         case CF_SEASONAL:

diff --git a/src/rrd_dump.c b/src/rrd_dump.c
index 316814d349c7f62f4832594815375e7e6e2ebecd..eb9b4690e0fbb22f1788b95711695a80a8fa3499 100644 (file)
--- a/src/rrd_dump.c
+++ b/src/rrd_dump.c
@@ -176,6 +176,7 @@ int rrd_dump_r(
         fprintf(out_file, "\t\t<params>\n");
         switch (cf_conv(rrd.rra_def[i].cf_nam)) {
         case CF_HWPREDICT:
+        case CF_MHWPREDICT:
             fprintf(out_file, "\t\t<hw_alpha> %0.10e </hw_alpha>\n",
                     rrd.rra_def[i].par[RRA_hw_alpha].u_val);
             fprintf(out_file, "\t\t<hw_beta> %0.10e </hw_beta>\n",
@@ -255,6 +256,7 @@ int rrd_dump_r(
             }
             switch (cf_conv(rrd.rra_def[i].cf_nam)) {
             case CF_HWPREDICT:
+            case CF_MHWPREDICT:
                 value =
                     rrd.cdp_prep[i * rrd.stat_head->ds_cnt +
                                  ii].scratch[CDP_hw_intercept].u_val;

diff --git a/src/rrd_format.c b/src/rrd_format.c
index 6002495dcaa040c0828e807827085fba16a283f0..41dee7e0fdd26b7b0df6973d1acd540e0488aac9 100644 (file)
--- a/src/rrd_format.c
+++ b/src/rrd_format.c
@@ -75,6 +75,7 @@ enum cf_en cf_conv(
         converter(MAX, CF_MAXIMUM)
         converter(LAST, CF_LAST)
         converter(HWPREDICT, CF_HWPREDICT)
+        converter(MHWPREDICT, CF_MHWPREDICT)
         converter(DEVPREDICT, CF_DEVPREDICT)
         converter(SEASONAL, CF_SEASONAL)
         converter(DEVSEASONAL, CF_DEVSEASONAL)

diff --git a/src/rrd_format.h b/src/rrd_format.h
index c841a6d0dc3134559589540ac3287a757c2b4b59..c4c8690f171e8c7d477dfa068fba190482b24c2d 100644 (file)
--- a/src/rrd_format.h
+++ b/src/rrd_format.h
@@ -165,6 +165,7 @@ enum cf_en { CF_AVERAGE = 0,    /* data consolidation functions */
     CF_MAXIMUM,
     CF_LAST,
     CF_HWPREDICT,
+    CF_MHWPREDICT,
     /* An array of predictions using the seasonal 
      * Holt-Winters algorithm. Requires an RRA of type
      * CF_SEASONAL for this data source. */

diff --git a/src/rrd_graph.c b/src/rrd_graph.c
index 2d1785f40d8102f829aa8204dd7a82bb985812c7..3022227815392680370987f85351f28c023a9e07 100644 (file)
--- a/src/rrd_graph.c
+++ b/src/rrd_graph.c
@@ -717,6 +717,7 @@ void reduce_data(
                 else {
                     switch (cf) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                     case CF_DEVSEASONAL:
                     case CF_DEVPREDICT:
                     case CF_SEASONAL:
@@ -742,6 +743,7 @@ void reduce_data(
             } else {
                 switch (cf) {
                 case CF_HWPREDICT:
+                case CF_MHWPREDICT:
                 case CF_DEVSEASONAL:
                 case CF_DEVPREDICT:
                 case CF_SEASONAL:
@@ -1438,6 +1440,7 @@ int print_calc(
 
                     switch (im->gdes[i].cf) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                     case CF_DEVPREDICT:
                     case CF_DEVSEASONAL:
                     case CF_SEASONAL:

diff --git a/src/rrd_hw.c b/src/rrd_hw.c
index d9247c9ecafb7c07271c712978217906be19ce55..24cfb1709a6869a72913537228cf447759570fc2 100644 (file)
--- a/src/rrd_hw.c
+++ b/src/rrd_hw.c
@@ -8,139 +8,17 @@
 
 #include "rrd_tool.h"
 #include "rrd_hw.h"
+#include "rrd_hw_math.h"
+#include "rrd_hw_update.h"
+
+#define hw_dep_idx(rrd, rra_idx) rrd->rra_def[rra_idx].par[RRA_dependent_rra_idx].u_cnt
 
 /* #define DEBUG */
 
 /* private functions */
-unsigned long MyMod(
+static unsigned long MyMod(
     signed long val,
     unsigned long mod);
-int       update_hwpredict(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx);
-int       update_seasonal(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx,
-    rrd_value_t *seasonal_coef);
-int       update_devpredict(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx);
-int       update_devseasonal(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx,
-    rrd_value_t *seasonal_dev);
-int       update_failures(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx);
-
-int update_hwpredict(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx)
-{
-    rrd_value_t prediction, seasonal_coef;
-    unsigned long dependent_rra_idx, seasonal_cdp_idx;
-    unival   *coefs = rrd->cdp_prep[cdp_idx].scratch;
-    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
-
-    /* save coefficients from current prediction */
-    coefs[CDP_hw_last_intercept].u_val = coefs[CDP_hw_intercept].u_val;
-    coefs[CDP_hw_last_slope].u_val = coefs[CDP_hw_slope].u_val;
-    coefs[CDP_last_null_count].u_cnt = coefs[CDP_null_count].u_cnt;
-
-    /* retrieve the current seasonal coef */
-    dependent_rra_idx = current_rra->par[RRA_dependent_rra_idx].u_cnt;
-    seasonal_cdp_idx = dependent_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
-    if (dependent_rra_idx < rra_idx)
-        seasonal_coef =
-            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_last_seasonal].
-            u_val;
-    else
-        seasonal_coef =
-            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_seasonal].u_val;
-
-    /* compute the prediction */
-    if (isnan(coefs[CDP_hw_intercept].u_val)
-        || isnan(coefs[CDP_hw_slope].u_val)
-        || isnan(seasonal_coef)) {
-        prediction = DNAN;
-
-        /* bootstrap initialization of slope and intercept */
-        if (isnan(coefs[CDP_hw_intercept].u_val) &&
-            !isnan(coefs[CDP_scratch_idx].u_val)) {
-#ifdef DEBUG
-            fprintf(stderr, "Initialization of slope/intercept\n");
-#endif
-            coefs[CDP_hw_intercept].u_val = coefs[CDP_scratch_idx].u_val;
-            coefs[CDP_hw_last_intercept].u_val = coefs[CDP_scratch_idx].u_val;
-            /* initialize the slope to 0 */
-            coefs[CDP_hw_slope].u_val = 0.0;
-            coefs[CDP_hw_last_slope].u_val = 0.0;
-            /* initialize null count to 1 */
-            coefs[CDP_null_count].u_cnt = 1;
-            coefs[CDP_last_null_count].u_cnt = 1;
-        }
-        /* if seasonal coefficient is NA, then don't update intercept, slope */
-    } else {
-        prediction = coefs[CDP_hw_intercept].u_val +
-            (coefs[CDP_hw_slope].u_val) * (coefs[CDP_null_count].u_cnt)
-            + seasonal_coef;
-#ifdef DEBUG
-        fprintf(stderr, "computed prediction: %f\n", prediction);
-#endif
-        if (isnan(coefs[CDP_scratch_idx].u_val)) {
-            /* NA value, no updates of intercept, slope;
-             * increment the null count */
-            (coefs[CDP_null_count].u_cnt)++;
-        } else {
-#ifdef DEBUG
-            fprintf(stderr, "Updating intercept, slope\n");
-#endif
-            /* update the intercept */
-            coefs[CDP_hw_intercept].u_val =
-                (current_rra->par[RRA_hw_alpha].u_val) *
-                (coefs[CDP_scratch_idx].u_val - seasonal_coef) + (1 -
-                                                                  current_rra->
-                                                                  par
-                                                                  [RRA_hw_alpha].
-                                                                  u_val) *
-                (coefs[CDP_hw_intercept].u_val +
-                 (coefs[CDP_hw_slope].u_val) * (coefs[CDP_null_count].u_cnt));
-            /* update the slope */
-            coefs[CDP_hw_slope].u_val =
-                (current_rra->par[RRA_hw_beta].u_val) *
-                (coefs[CDP_hw_intercept].u_val -
-                 coefs[CDP_hw_last_intercept].u_val) + (1 -
-                                                        current_rra->
-                                                        par[RRA_hw_beta].
-                                                        u_val) *
-                (coefs[CDP_hw_slope].u_val);
-            /* reset the null count */
-            coefs[CDP_null_count].u_cnt = 1;
-        }
-    }
-
-    /* store the prediction for writing */
-    coefs[CDP_scratch_idx].u_val = prediction;
-    return 0;
-}
 
 int lookup_seasonal(
     rrd_t *rrd,
@@ -199,348 +77,6 @@ int lookup_seasonal(
     return -1;
 }
 
-int update_seasonal(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx,
-    rrd_value_t *seasonal_coef)
-{
-/* TODO: extract common if subblocks in the wake of I/O optimization */
-    rrd_value_t intercept, seasonal;
-    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
-    rra_def_t *hw_rra =
-        &(rrd->rra_def[current_rra->par[RRA_dependent_rra_idx].u_cnt]);
-    /* obtain cdp_prep index for HWPREDICT */
-    unsigned long hw_cdp_idx = (current_rra->par[RRA_dependent_rra_idx].u_cnt)
-        * (rrd->stat_head->ds_cnt) + ds_idx;
-    unival   *coefs = rrd->cdp_prep[hw_cdp_idx].scratch;
-
-    /* update seasonal coefficient in cdp prep areas */
-    seasonal = rrd->cdp_prep[cdp_idx].scratch[CDP_hw_seasonal].u_val;
-    rrd->cdp_prep[cdp_idx].scratch[CDP_hw_last_seasonal].u_val = seasonal;
-    rrd->cdp_prep[cdp_idx].scratch[CDP_hw_seasonal].u_val =
-        seasonal_coef[ds_idx];
-
-    /* update seasonal value for disk */
-    if (current_rra->par[RRA_dependent_rra_idx].u_cnt < rra_idx)
-        /* associated HWPREDICT has already been updated */
-        /* check for possible NA values */
-        if (isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
-            /* no update, store the old value unchanged,
-             * doesn't matter if it is NA */
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = seasonal;
-        } else if (isnan(coefs[CDP_hw_last_intercept].u_val)
-                   || isnan(coefs[CDP_hw_last_slope].u_val)) {
-            /* this should never happen, as HWPREDICT was already updated */
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = DNAN;
-        } else if (isnan(seasonal)) {
-            /* initialization: intercept is not currently being updated */
-#ifdef DEBUG
-            fprintf(stderr, "Initialization of seasonal coef %lu\n",
-                    rrd->rra_ptr[rra_idx].cur_row);
-#endif
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val
-                -= coefs[CDP_hw_last_intercept].u_val;
-        } else {
-            intercept = coefs[CDP_hw_intercept].u_val;
-#ifdef DEBUG
-            fprintf(stderr,
-                    "Updating seasonal, params: gamma %f, new intercept %f, old seasonal %f\n",
-                    current_rra->par[RRA_seasonal_gamma].u_val,
-                    intercept, seasonal);
-#endif
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
-                (current_rra->par[RRA_seasonal_gamma].u_val) *
-                (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val -
-                 intercept) + (1 -
-                               current_rra->par[RRA_seasonal_gamma].u_val) *
-                seasonal;
-    } else {
-        /* SEASONAL array is updated first, which means the new intercept
-         * hasn't be computed; so we compute it here. */
-
-        /* check for possible NA values */
-        if (isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
-            /* no update, simple store the old value unchanged,
-             * doesn't matter if it is NA */
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = seasonal;
-        } else if (isnan(coefs[CDP_hw_intercept].u_val)
-                   || isnan(coefs[CDP_hw_slope].u_val)) {
-            /* Initialization of slope and intercept will occur.
-             * force seasonal coefficient to 0. */
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = 0.0;
-        } else if (isnan(seasonal)) {
-            /* initialization: intercept will not be updated
-             * CDP_hw_intercept = CDP_hw_last_intercept; just need to 
-             * subtract this baseline value. */
-#ifdef DEBUG
-            fprintf(stderr, "Initialization of seasonal coef %lu\n",
-                    rrd->rra_ptr[rra_idx].cur_row);
-#endif
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val -=
-                coefs[CDP_hw_intercept].u_val;
-        } else {
-            /* Note that we must get CDP_scratch_idx from SEASONAL array, as CDP_scratch_idx
-             * for HWPREDICT array will be DNAN. */
-            intercept = (hw_rra->par[RRA_hw_alpha].u_val) *
-                (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val -
-                 seasonal)
-                + (1 -
-                   hw_rra->par[RRA_hw_alpha].u_val) *
-                (coefs[CDP_hw_intercept].u_val +
-                 (coefs[CDP_hw_slope].u_val) * (coefs[CDP_null_count].u_cnt));
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
-                (current_rra->par[RRA_seasonal_gamma].u_val) *
-                (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val -
-                 intercept) + (1 -
-                               current_rra->par[RRA_seasonal_gamma].u_val) *
-                seasonal;
-        }
-    }
-#ifdef DEBUG
-    fprintf(stderr, "seasonal coefficient set= %f\n",
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val);
-#endif
-    return 0;
-}
-
-int update_devpredict(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx)
-{
-    /* there really isn't any "update" here; the only reason this information
-     * is stored separately from DEVSEASONAL is to preserve deviation predictions
-     * for a longer duration than one seasonal cycle. */
-    unsigned long seasonal_cdp_idx =
-        (rrd->rra_def[rra_idx].par[RRA_dependent_rra_idx].u_cnt)
-        * (rrd->stat_head->ds_cnt) + ds_idx;
-
-    if (rrd->rra_def[rra_idx].par[RRA_dependent_rra_idx].u_cnt < rra_idx) {
-        /* associated DEVSEASONAL array already updated */
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val
-            =
-            rrd->cdp_prep[seasonal_cdp_idx].
-            scratch[CDP_last_seasonal_deviation].u_val;
-    } else {
-        /* associated DEVSEASONAL not yet updated */
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val
-            =
-            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_seasonal_deviation].
-            u_val;
-    }
-    return 0;
-}
-
-int update_devseasonal(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx,
-    rrd_value_t *seasonal_dev)
-{
-    rrd_value_t prediction = 0, seasonal_coef = DNAN;
-    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
-
-    /* obtain cdp_prep index for HWPREDICT */
-    unsigned long hw_rra_idx = current_rra->par[RRA_dependent_rra_idx].u_cnt;
-    unsigned long hw_cdp_idx = hw_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
-    unsigned long seasonal_cdp_idx;
-    unival   *coefs = rrd->cdp_prep[hw_cdp_idx].scratch;
-
-    rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].u_val =
-        rrd->cdp_prep[cdp_idx].scratch[CDP_seasonal_deviation].u_val;
-    /* retrieve the next seasonal deviation value, could be NA */
-    rrd->cdp_prep[cdp_idx].scratch[CDP_seasonal_deviation].u_val =
-        seasonal_dev[ds_idx];
-
-    /* retrieve the current seasonal_coef (not to be confused with the
-     * current seasonal deviation). Could make this more readable by introducing
-     * some wrapper functions. */
-    seasonal_cdp_idx =
-        (rrd->rra_def[hw_rra_idx].par[RRA_dependent_rra_idx].u_cnt)
-        * (rrd->stat_head->ds_cnt) + ds_idx;
-    if (rrd->rra_def[hw_rra_idx].par[RRA_dependent_rra_idx].u_cnt < rra_idx)
-        /* SEASONAL array already updated */
-        seasonal_coef =
-            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_last_seasonal].
-            u_val;
-    else
-        /* SEASONAL array not yet updated */
-        seasonal_coef =
-            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_seasonal].u_val;
-
-    /* compute the abs value of the difference between the prediction and
-     * observed value */
-    if (hw_rra_idx < rra_idx) {
-        /* associated HWPREDICT has already been updated */
-        if (isnan(coefs[CDP_hw_last_intercept].u_val) ||
-            isnan(coefs[CDP_hw_last_slope].u_val) || isnan(seasonal_coef)) {
-            /* one of the prediction values is uinitialized */
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = DNAN;
-            return 0;
-        } else {
-            prediction = coefs[CDP_hw_last_intercept].u_val +
-                (coefs[CDP_hw_last_slope].u_val) *
-                (coefs[CDP_last_null_count].u_cnt)
-                + seasonal_coef;
-        }
-    } else {
-        /* associated HWPREDICT has NOT been updated */
-        if (isnan(coefs[CDP_hw_intercept].u_val) ||
-            isnan(coefs[CDP_hw_slope].u_val) || isnan(seasonal_coef)) {
-            /* one of the prediction values is uinitialized */
-            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = DNAN;
-            return 0;
-        } else {
-            prediction = coefs[CDP_hw_intercept].u_val +
-                (coefs[CDP_hw_slope].u_val) * (coefs[CDP_null_count].u_cnt)
-                + seasonal_coef;
-        }
-    }
-
-    if (isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
-        /* no update, store existing value unchanged, doesn't
-         * matter if it is NA */
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
-            rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].u_val;
-    } else
-        if (isnan
-            (rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].
-             u_val)) {
-        /* initialization */
-#ifdef DEBUG
-        fprintf(stderr, "Initialization of seasonal deviation\n");
-#endif
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
-            fabs(prediction -
-                 rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val);
-    } else {
-        /* exponential smoothing update */
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
-            (rrd->rra_def[rra_idx].par[RRA_seasonal_gamma].u_val) *
-            fabs(prediction -
-                 rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)
-            + (1 -
-               rrd->rra_def[rra_idx].par[RRA_seasonal_gamma].u_val) *
-            (rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].
-             u_val);
-    }
-    return 0;
-}
-
-/* Check for a failure based on a threshold # of violations within the specified
- * window. */
-int update_failures(
-    rrd_t *rrd,
-    unsigned long cdp_idx,
-    unsigned long rra_idx,
-    unsigned long ds_idx,
-    unsigned short CDP_scratch_idx)
-{
-    /* detection of a violation depends on 3 RRAs:
-     * HWPREDICT, SEASONAL, and DEVSEASONAL */
-    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
-    unsigned long dev_rra_idx = current_rra->par[RRA_dependent_rra_idx].u_cnt;
-    rra_def_t *dev_rra = &(rrd->rra_def[dev_rra_idx]);
-    unsigned long hw_rra_idx = dev_rra->par[RRA_dependent_rra_idx].u_cnt;
-    rra_def_t *hw_rra = &(rrd->rra_def[hw_rra_idx]);
-    unsigned long seasonal_rra_idx = hw_rra->par[RRA_dependent_rra_idx].u_cnt;
-    unsigned long temp_cdp_idx;
-    rrd_value_t deviation = DNAN;
-    rrd_value_t seasonal_coef = DNAN;
-    rrd_value_t prediction = DNAN;
-    char      violation = 0;
-    unsigned short violation_cnt = 0, i;
-    char     *violations_array;
-
-    /* usual checks to determine the order of the RRAs */
-    temp_cdp_idx = dev_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
-    if (rra_idx < seasonal_rra_idx) {
-        /* DEVSEASONAL not yet updated */
-        deviation =
-            rrd->cdp_prep[temp_cdp_idx].scratch[CDP_seasonal_deviation].u_val;
-    } else {
-        /* DEVSEASONAL already updated */
-        deviation =
-            rrd->cdp_prep[temp_cdp_idx].scratch[CDP_last_seasonal_deviation].
-            u_val;
-    }
-    if (!isnan(deviation)) {
-
-        temp_cdp_idx = seasonal_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
-        if (rra_idx < seasonal_rra_idx) {
-            /* SEASONAL not yet updated */
-            seasonal_coef =
-                rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_seasonal].u_val;
-        } else {
-            /* SEASONAL already updated */
-            seasonal_coef =
-                rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_last_seasonal].
-                u_val;
-        }
-        /* in this code block, we know seasonal coef is not DNAN, because deviation is not
-         * null */
-
-        temp_cdp_idx = hw_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
-        if (rra_idx < hw_rra_idx) {
-            /* HWPREDICT not yet updated */
-            prediction =
-                rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_intercept].u_val +
-                (rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_slope].u_val)
-                * (rrd->cdp_prep[temp_cdp_idx].scratch[CDP_null_count].u_cnt)
-                + seasonal_coef;
-        } else {
-            /* HWPREDICT already updated */
-            prediction =
-                rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_last_intercept].
-                u_val +
-                (rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_last_slope].u_val)
-                *
-                (rrd->cdp_prep[temp_cdp_idx].scratch[CDP_last_null_count].
-                 u_cnt)
-                + seasonal_coef;
-        }
-
-        /* determine if the observed value is a violation */
-        if (!isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
-            if (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val >
-                prediction +
-                (current_rra->par[RRA_delta_pos].u_val) * deviation
-                || rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val <
-                prediction -
-                (current_rra->par[RRA_delta_neg].u_val) * deviation)
-                violation = 1;
-        } else {
-            violation = 1;  /* count DNAN values as violations */
-        }
-
-    }
-
-    /* determine if a failure has occurred and update the failure array */
-    violation_cnt = violation;
-    violations_array = (char *) ((void *) rrd->cdp_prep[cdp_idx].scratch);
-    for (i = current_rra->par[RRA_window_len].u_cnt; i > 1; i--) {
-        /* shift */
-        violations_array[i - 1] = violations_array[i - 2];
-        violation_cnt += violations_array[i - 1];
-    }
-    violations_array[0] = violation;
-
-    if (violation_cnt < current_rra->par[RRA_failure_threshold].u_cnt)
-        /* not a failure */
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = 0.0;
-    else
-        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = 1.0;
-
-    return (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val);
-}
-
 /* For the specified CDP prep area and the FAILURES RRA,
  * erase all history of past violations.
  */
@@ -693,14 +229,34 @@ int apply_smoother(
     free(working_average);
 
     if (cf_conv(rrd->rra_def[rra_idx].cf_nam) == CF_SEASONAL) {
+        rrd_value_t (
+    *init_seasonality) (
+    rrd_value_t seasonal_coef,
+    rrd_value_t intercept);
+
+        switch (cf_conv(rrd->rra_def[hw_dep_idx(rrd, rra_idx)].cf_nam)) {
+        case CF_HWPREDICT:
+            init_seasonality = hw_additive_init_seasonality;
+            break;
+        case CF_MHWPREDICT:
+            init_seasonality = hw_multiplicative_init_seasonality;
+            break;
+        default:
+            rrd_set_error("apply smoother: SEASONAL rra doesn't have "
+                          "valid dependency: %s",
+                          rrd->rra_def[hw_dep_idx(rrd, rra_idx)].cf_nam);
+            return -1;
+        }
+
         for (j = 0; j < row_length; ++j) {
             for (i = 0; i < row_count; ++i) {
-                rrd_values[i * row_length + j] -= baseline[j];
+                rrd_values[i * row_length + j] =
+                    init_seasonality(rrd_values[i * row_length + j],
+                                     baseline[j]);
             }
             /* update the baseline coefficient,
              * first, compute the cdp_index. */
-            offset = (rrd->rra_def[rra_idx].par[RRA_dependent_rra_idx].u_cnt)
-                * row_length + j;
+            offset = hw_dep_idx(rrd, rra_idx) * row_length + j;
             (rrd->cdp_prep[offset]).scratch[CDP_hw_intercept].u_val +=
                 baseline[j];
         }
@@ -775,6 +331,7 @@ void reset_aberrant_coefficients(
         cdp_idx = rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
         switch (cf_conv(rrd->rra_def[rra_idx].cf_nam)) {
         case CF_HWPREDICT:
+        case CF_MHWPREDICT:
             init_hwpredict_cdp(&(rrd->cdp_prep[cdp_idx]));
             break;
         case CF_SEASONAL:
@@ -850,31 +407,88 @@ int update_aberrant_CF(
     unsigned short CDP_scratch_idx,
     rrd_value_t *seasonal_coef)
 {
+    static hw_functions_t hw_multiplicative_functions = {
+        hw_multiplicative_calculate_prediction,
+        hw_multiplicative_calculate_intercept,
+        hw_calculate_slope,
+        hw_multiplicative_calculate_seasonality,
+        hw_multiplicative_init_seasonality,
+        hw_calculate_seasonal_deviation,
+        hw_init_seasonal_deviation,
+        1.0             // identity value
+    };
+
+    static hw_functions_t hw_additive_functions = {
+        hw_additive_calculate_prediction,
+        hw_additive_calculate_intercept,
+        hw_calculate_slope,
+        hw_additive_calculate_seasonality,
+        hw_additive_init_seasonality,
+        hw_calculate_seasonal_deviation,
+        hw_init_seasonal_deviation,
+        0.0             // identity value 
+    };
+
     rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = pdp_val;
     switch (current_cf) {
-    case CF_AVERAGE:
-    default:
-        return 0;
     case CF_HWPREDICT:
         return update_hwpredict(rrd, cdp_idx, rra_idx, ds_idx,
-                                CDP_scratch_idx);
+                                CDP_scratch_idx, &hw_additive_functions);
+    case CF_MHWPREDICT:
+        return update_hwpredict(rrd, cdp_idx, rra_idx, ds_idx,
+                                CDP_scratch_idx,
+                                &hw_multiplicative_functions);
     case CF_DEVPREDICT:
         return update_devpredict(rrd, cdp_idx, rra_idx, ds_idx,
                                  CDP_scratch_idx);
     case CF_SEASONAL:
-        return update_seasonal(rrd, cdp_idx, rra_idx, ds_idx, CDP_scratch_idx,
-                               seasonal_coef);
+        switch (cf_conv(rrd->rra_def[hw_dep_idx(rrd, rra_idx)].cf_nam)) {
+        case CF_HWPREDICT:
+            return update_seasonal(rrd, cdp_idx, rra_idx, ds_idx,
+                                   CDP_scratch_idx, seasonal_coef,
+                                   &hw_additive_functions);
+        case CF_MHWPREDICT:
+            return update_seasonal(rrd, cdp_idx, rra_idx, ds_idx,
+                                   CDP_scratch_idx, seasonal_coef,
+                                   &hw_multiplicative_functions);
+        default:
+            return -1;
+        }
     case CF_DEVSEASONAL:
-        return update_devseasonal(rrd, cdp_idx, rra_idx, ds_idx,
-                                  CDP_scratch_idx, seasonal_coef);
+        switch (cf_conv(rrd->rra_def[hw_dep_idx(rrd, rra_idx)].cf_nam)) {
+        case CF_HWPREDICT:
+            return update_devseasonal(rrd, cdp_idx, rra_idx, ds_idx,
+                                      CDP_scratch_idx, seasonal_coef,
+                                      &hw_additive_functions);
+        case CF_MHWPREDICT:
+            return update_devseasonal(rrd, cdp_idx, rra_idx, ds_idx,
+                                      CDP_scratch_idx, seasonal_coef,
+                                      &hw_multiplicative_functions);
+        default:
+            return -1;
+        }
     case CF_FAILURES:
-        return update_failures(rrd, cdp_idx, rra_idx, ds_idx,
-                               CDP_scratch_idx);
+        switch (cf_conv
+                (rrd->rra_def[hw_dep_idx(rrd, hw_dep_idx(rrd, rra_idx))].
+                 cf_nam)) {
+        case CF_HWPREDICT:
+            return update_failures(rrd, cdp_idx, rra_idx, ds_idx,
+                                   CDP_scratch_idx, &hw_additive_functions);
+        case CF_MHWPREDICT:
+            return update_failures(rrd, cdp_idx, rra_idx, ds_idx,
+                                   CDP_scratch_idx,
+                                   &hw_multiplicative_functions);
+        default:
+            return -1;
+        }
+    case CF_AVERAGE:
+    default:
+        return 0;
     }
     return -1;
 }
 
-unsigned long MyMod(
+static unsigned long MyMod(
     signed long val,
     unsigned long mod)
 {

diff --git a/src/rrd_hw_math.c b/src/rrd_hw_math.c
new file mode 100644 (file)
index 0000000..2b2c00f
--- /dev/null
+++ b/src/rrd_hw_math.c
@@ -0,0 +1,143 @@
+/*****************************************************************************
+ * rrd_hw_math.c  Math functions for Holt-Winters computations
+ *****************************************************************************/
+
+#include "rrd_hw_math.h"
+#include "rrd_config.h"
+
+/*****************************************************************************
+ * RRDtool supports both the additive and multiplicative Holt-Winters methods. 
+ * The additive method makes predictions by adding seasonality to the baseline, 
+ * whereas the multiplicative method multiplies the seasonality coefficient by 
+ * the baseline to make a prediction. This file contains all the differences
+ * between the additive and multiplicative methods, as well as a few math 
+ * functions common to them both.
+ ****************************************************************************/
+
+/*****************************************************************************
+ * Functions for additive Holt-Winters
+ *****************************************************************************/
+
+rrd_value_t hw_additive_calculate_prediction(
+    rrd_value_t intercept,
+    rrd_value_t slope,
+    int null_count,
+    rrd_value_t seasonal_coef)
+{
+    return intercept + slope * null_count + seasonal_coef;
+}
+
+rrd_value_t hw_additive_calculate_intercept(
+    rrd_value_t hw_alpha,
+    rrd_value_t observed,
+    rrd_value_t seasonal_coef,
+    unival *coefs)
+{
+    return hw_alpha * (observed - seasonal_coef)
+        + (1 - hw_alpha) * (coefs[CDP_hw_intercept].u_val
+                            +
+                            (coefs[CDP_hw_slope].u_val) *
+                            (coefs[CDP_null_count].u_cnt));
+}
+
+rrd_value_t hw_additive_calculate_seasonality(
+    rrd_value_t hw_gamma,
+    rrd_value_t observed,
+    rrd_value_t intercept,
+    rrd_value_t seasonal_coef)
+{
+    return hw_gamma * (observed - intercept)
+        + (1 - hw_gamma) * seasonal_coef;
+}
+
+rrd_value_t hw_additive_init_seasonality(
+    rrd_value_t seasonal_coef,
+    rrd_value_t intercept)
+{
+    return seasonal_coef - intercept;
+}
+
+/*****************************************************************************
+ * Functions for multiplicative Holt-Winters
+ *****************************************************************************/
+
+rrd_value_t hw_multiplicative_calculate_prediction(
+    rrd_value_t intercept,
+    rrd_value_t slope,
+    int null_count,
+    rrd_value_t seasonal_coef)
+{
+    return (intercept + slope * null_count) * seasonal_coef;
+}
+
+rrd_value_t hw_multiplicative_calculate_intercept(
+    rrd_value_t hw_alpha,
+    rrd_value_t observed,
+    rrd_value_t seasonal_coef,
+    unival *coefs)
+{
+    if (seasonal_coef <= 0) {
+        return DNAN;
+    }
+
+    return hw_alpha * (observed / seasonal_coef)
+        + (1 - hw_alpha) * (coefs[CDP_hw_intercept].u_val
+                            +
+                            (coefs[CDP_hw_slope].u_val) *
+                            (coefs[CDP_null_count].u_cnt));
+}
+
+rrd_value_t hw_multiplicative_calculate_seasonality(
+    rrd_value_t hw_gamma,
+    rrd_value_t observed,
+    rrd_value_t intercept,
+    rrd_value_t seasonal_coef)
+{
+    if (intercept <= 0) {
+        return DNAN;
+    }
+
+    return hw_gamma * (observed / intercept)
+        + (1 - hw_gamma) * seasonal_coef;
+}
+
+rrd_value_t hw_multiplicative_init_seasonality(
+    rrd_value_t seasonal_coef,
+    rrd_value_t intercept)
+{
+    if (intercept <= 0) {
+        return DNAN;
+    }
+
+    return seasonal_coef / intercept;
+}
+
+/*****************************************************************************
+ * Math functions common to additive and multiplicative Holt-Winters
+ *****************************************************************************/
+
+rrd_value_t hw_calculate_slope(
+    rrd_value_t hw_beta,
+    unival *coefs)
+{
+    return hw_beta * (coefs[CDP_hw_intercept].u_val -
+                      coefs[CDP_hw_last_intercept].u_val)
+        + (1 - hw_beta) * coefs[CDP_hw_slope].u_val;
+}
+
+rrd_value_t hw_calculate_seasonal_deviation(
+    rrd_value_t hw_gamma,
+    rrd_value_t prediction,
+    rrd_value_t observed,
+    rrd_value_t last)
+{
+    return hw_gamma * fabs(prediction - observed)
+        + (1 - hw_gamma) * last;
+}
+
+rrd_value_t hw_init_seasonal_deviation(
+    rrd_value_t prediction,
+    rrd_value_t observed)
+{
+    return fabs(prediction - observed);
+}

diff --git a/src/rrd_hw_math.h b/src/rrd_hw_math.h
new file mode 100644 (file)
index 0000000..3677b31
--- /dev/null
+++ b/src/rrd_hw_math.h
@@ -0,0 +1,132 @@
+/*****************************************************************************
+ * rrd_hw_math.h  Math functions for Holt-Winters computations
+ *****************************************************************************/
+
+#include "rrd.h"
+#include "rrd_format.h"
+
+/* since /usr/include/bits/mathcalls.h:265 defines gamma already */
+#define gamma hw_gamma
+
+/*****************************************************************************
+ * Functions for additive Holt-Winters
+ *****************************************************************************/
+
+rrd_value_t hw_additive_calculate_prediction(
+    rrd_value_t intercept,
+    rrd_value_t slope,
+    int null_count,
+    rrd_value_t seasonal_coef);
+
+rrd_value_t hw_additive_calculate_intercept(
+    rrd_value_t alpha,
+    rrd_value_t scratch,
+    rrd_value_t seasonal_coef,
+    unival *coefs);
+
+rrd_value_t hw_additive_calculate_seasonality(
+    rrd_value_t gamma,
+    rrd_value_t scratch,
+    rrd_value_t intercept,
+    rrd_value_t seasonal_coef);
+
+rrd_value_t hw_additive_init_seasonality(
+    rrd_value_t seasonal_coef,
+    rrd_value_t intercept);
+
+/*****************************************************************************
+ * Functions for multiplicative Holt-Winters
+ *****************************************************************************/
+
+rrd_value_t hw_multiplicative_calculate_prediction(
+    rrd_value_t intercept,
+    rrd_value_t slope,
+    int null_count,
+    rrd_value_t seasonal_coef);
+
+rrd_value_t hw_multiplicative_calculate_intercept(
+    rrd_value_t alpha,
+    rrd_value_t scratch,
+    rrd_value_t seasonal_coef,
+    unival *coefs);
+
+rrd_value_t hw_multiplicative_calculate_seasonality(
+    rrd_value_t gamma,
+    rrd_value_t scratch,
+    rrd_value_t intercept,
+    rrd_value_t seasonal_coef);
+
+rrd_value_t hw_multiplicative_init_seasonality(
+    rrd_value_t seasonal_coef,
+    rrd_value_t intercept);
+
+/*****************************************************************************
+ * Math functions common to additive and multiplicative Holt-Winters
+ *****************************************************************************/
+
+rrd_value_t hw_calculate_slope(
+    rrd_value_t beta,
+    unival *coefs);
+
+rrd_value_t hw_calculate_seasonal_deviation(
+    rrd_value_t gamma,
+    rrd_value_t prediction,
+    rrd_value_t observed,
+    rrd_value_t last);
+
+rrd_value_t hw_init_seasonal_deviation(
+    rrd_value_t prediction,
+    rrd_value_t observed);
+
+
+/* Function container */
+
+typedef struct hw_functions_t {
+    rrd_value_t (
+    *predict) (
+    rrd_value_t intercept,
+    rrd_value_t slope,
+    int null_count,
+    rrd_value_t seasonal_coef);
+
+    rrd_value_t (
+    *intercept) (
+    rrd_value_t alpha,
+    rrd_value_t observed,
+    rrd_value_t seasonal_coef,
+    unival *coefs);
+
+    rrd_value_t (
+    *slope)   (
+    rrd_value_t beta,
+    unival *coefs);
+
+    rrd_value_t (
+    *seasonality) (
+    rrd_value_t gamma,
+    rrd_value_t observed,
+    rrd_value_t intercept,
+    rrd_value_t seasonal_coef);
+
+    rrd_value_t (
+    *init_seasonality) (
+    rrd_value_t seasonal_coef,
+    rrd_value_t intercept);
+
+    rrd_value_t (
+    *seasonal_deviation) (
+    rrd_value_t gamma,
+    rrd_value_t prediction,
+    rrd_value_t observed,
+    rrd_value_t last);
+
+    rrd_value_t (
+    *init_seasonal_deviation) (
+    rrd_value_t prediction,
+    rrd_value_t observed);
+
+    rrd_value_t identity;
+} hw_functions_t;
+
+
+#undef gamma

diff --git a/src/rrd_hw_update.c b/src/rrd_hw_update.c
new file mode 100644 (file)
index 0000000..89d0c8e
--- /dev/null
+++ b/src/rrd_hw_update.c
@@ -0,0 +1,474 @@
+/*****************************************************************************
+ * rrd_hw_update.c  Functions for updating a Holt-Winters RRA
+ ****************************************************************************/
+
+#include "rrd_format.h"
+#include "rrd_config.h"
+#include "rrd_hw_math.h"
+#include "rrd_hw_update.h"
+
+static void init_slope_intercept(
+    unival *coefs,
+    unsigned short CDP_scratch_idx)
+{
+#ifdef DEBUG
+    fprintf(stderr, "Initialization of slope/intercept\n");
+#endif
+    coefs[CDP_hw_intercept].u_val = coefs[CDP_scratch_idx].u_val;
+    coefs[CDP_hw_last_intercept].u_val = coefs[CDP_scratch_idx].u_val;
+    /* initialize the slope to 0 */
+    coefs[CDP_hw_slope].u_val = 0.0;
+    coefs[CDP_hw_last_slope].u_val = 0.0;
+    /* initialize null count to 1 */
+    coefs[CDP_null_count].u_cnt = 1;
+    coefs[CDP_last_null_count].u_cnt = 1;
+}
+
+static int hw_is_violation(
+    rrd_value_t observed,
+    rrd_value_t prediction,
+    rrd_value_t deviation,
+    rrd_value_t delta_pos,
+    rrd_value_t delta_neg)
+{
+    return (observed > prediction + delta_pos * deviation
+            || observed < prediction - delta_neg * deviation);
+}
+
+int update_hwpredict(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    hw_functions_t * functions)
+{
+    rrd_value_t prediction;
+    unsigned long dependent_rra_idx, seasonal_cdp_idx;
+    unival   *coefs = rrd->cdp_prep[cdp_idx].scratch;
+    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
+
+    /* save coefficients from current prediction */
+    coefs[CDP_hw_last_intercept].u_val = coefs[CDP_hw_intercept].u_val;
+    coefs[CDP_hw_last_slope].u_val = coefs[CDP_hw_slope].u_val;
+    coefs[CDP_last_null_count].u_cnt = coefs[CDP_null_count].u_cnt;
+
+    /* retrieve the current seasonal coef */
+    dependent_rra_idx = current_rra->par[RRA_dependent_rra_idx].u_cnt;
+    seasonal_cdp_idx = dependent_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
+
+    rrd_value_t seasonal_coef = (dependent_rra_idx < rra_idx)
+        ? rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_last_seasonal].u_val
+        : rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_seasonal].u_val;
+
+    /* compute the prediction */
+    if (isnan(coefs[CDP_hw_intercept].u_val)
+        || isnan(coefs[CDP_hw_slope].u_val)
+        || isnan(seasonal_coef)) {
+        prediction = DNAN;
+
+        /* bootstrap initialization of slope and intercept */
+        if (isnan(coefs[CDP_hw_intercept].u_val) &&
+            !isnan(coefs[CDP_scratch_idx].u_val)) {
+            init_slope_intercept(coefs, CDP_scratch_idx);
+        }
+        /* if seasonal coefficient is NA, then don't update intercept, slope */
+    } else {
+        prediction = functions->predict(coefs[CDP_hw_intercept].u_val,
+                                        coefs[CDP_hw_slope].u_val,
+                                        coefs[CDP_null_count].u_cnt,
+                                        seasonal_coef);
+#ifdef DEBUG
+        fprintf(stderr,
+                "computed prediction: %f (intercept %f, slope %f, season %f)\n",
+                prediction, coefs[CDP_hw_intercept].u_val,
+                coefs[CDP_hw_slope].u_val, seasonal_coef);
+#endif
+        if (isnan(coefs[CDP_scratch_idx].u_val)) {
+            /* NA value, no updates of intercept, slope;
+             * increment the null count */
+            (coefs[CDP_null_count].u_cnt)++;
+        } else {
+            /* update the intercept */
+            coefs[CDP_hw_intercept].u_val =
+                functions->intercept(current_rra->par[RRA_hw_alpha].u_val,
+                                     coefs[CDP_scratch_idx].u_val,
+                                     seasonal_coef, coefs);
+
+            /* update the slope */
+            coefs[CDP_hw_slope].u_val =
+                functions->slope(current_rra->par[RRA_hw_beta].u_val, coefs);
+
+            /* reset the null count */
+            coefs[CDP_null_count].u_cnt = 1;
+#ifdef DEBUG
+            fprintf(stderr, "Updating intercept = %f, slope = %f\n",
+                    coefs[CDP_hw_intercept].u_val, coefs[CDP_hw_slope].u_val);
+#endif
+        }
+    }
+
+    /* store the prediction for writing */
+    coefs[CDP_scratch_idx].u_val = prediction;
+    return 0;
+}
+
+int update_seasonal(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    rrd_value_t *seasonal_coef,
+    hw_functions_t * functions)
+{
+/* TODO: extract common if subblocks in the wake of I/O optimization */
+    rrd_value_t intercept, seasonal;
+    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
+    rra_def_t *hw_rra =
+        &(rrd->rra_def[current_rra->par[RRA_dependent_rra_idx].u_cnt]);
+
+    /* obtain cdp_prep index for HWPREDICT */
+    unsigned long hw_cdp_idx = (current_rra->par[RRA_dependent_rra_idx].u_cnt)
+        * (rrd->stat_head->ds_cnt) + ds_idx;
+    unival   *coefs = rrd->cdp_prep[hw_cdp_idx].scratch;
+
+    /* update seasonal coefficient in cdp prep areas */
+    seasonal = rrd->cdp_prep[cdp_idx].scratch[CDP_hw_seasonal].u_val;
+    rrd->cdp_prep[cdp_idx].scratch[CDP_hw_last_seasonal].u_val = seasonal;
+    rrd->cdp_prep[cdp_idx].scratch[CDP_hw_seasonal].u_val =
+        seasonal_coef[ds_idx];
+
+    if (isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
+        /* no update, store the old value unchanged,
+         * doesn't matter if it is NA */
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = seasonal;
+        return 0;
+    }
+
+    /* update seasonal value for disk */
+    if (current_rra->par[RRA_dependent_rra_idx].u_cnt < rra_idx) {
+        /* associated HWPREDICT has already been updated */
+        /* check for possible NA values */
+        if (isnan(coefs[CDP_hw_last_intercept].u_val)
+            || isnan(coefs[CDP_hw_last_slope].u_val)) {
+            /* this should never happen, as HWPREDICT was already updated */
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = DNAN;
+        } else if (isnan(seasonal)) {
+            /* initialization: intercept is not currently being updated */
+#ifdef DEBUG
+            fprintf(stderr, "Initialization of seasonal coef %lu\n",
+                    rrd->rra_ptr[rra_idx].cur_row);
+#endif
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+                functions->init_seasonality(rrd->cdp_prep[cdp_idx].
+                                            scratch[CDP_scratch_idx].u_val,
+                                            coefs[CDP_hw_last_intercept].
+                                            u_val);
+        } else {
+            intercept = coefs[CDP_hw_intercept].u_val;
+
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+                functions->seasonality(current_rra->par[RRA_seasonal_gamma].
+                                       u_val,
+                                       rrd->cdp_prep[cdp_idx].
+                                       scratch[CDP_scratch_idx].u_val,
+                                       intercept, seasonal);
+#ifdef DEBUG
+            fprintf(stderr,
+                    "Updating seasonal = %f (params: gamma %f, new intercept %f, old seasonal %f)\n",
+                    rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val,
+                    current_rra->par[RRA_seasonal_gamma].u_val,
+                    intercept, seasonal);
+#endif
+        }
+    } else {
+        /* SEASONAL array is updated first, which means the new intercept
+         * hasn't be computed; so we compute it here. */
+
+        /* check for possible NA values */
+        if (isnan(coefs[CDP_hw_intercept].u_val)
+            || isnan(coefs[CDP_hw_slope].u_val)) {
+            /* Initialization of slope and intercept will occur.
+             * force seasonal coefficient to 0 or 1. */
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+                functions->identity;
+        } else if (isnan(seasonal)) {
+            /* initialization: intercept will not be updated
+             * CDP_hw_intercept = CDP_hw_last_intercept; just need to 
+             * subtract/divide by this baseline value. */
+#ifdef DEBUG
+            fprintf(stderr, "Initialization of seasonal coef %lu\n",
+                    rrd->rra_ptr[rra_idx].cur_row);
+#endif
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+                functions->init_seasonality(rrd->cdp_prep[cdp_idx].
+                                            scratch[CDP_scratch_idx].u_val,
+                                            coefs[CDP_hw_intercept].u_val);
+        } else {
+            /* Note that we must get CDP_scratch_idx from SEASONAL array, as CDP_scratch_idx
+             * for HWPREDICT array will be DNAN. */
+            intercept = functions->intercept(hw_rra->par[RRA_hw_alpha].u_val,
+                                             rrd->cdp_prep[cdp_idx].
+                                             scratch[CDP_scratch_idx].u_val,
+                                             seasonal, coefs);
+
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+                functions->seasonality(current_rra->par[RRA_seasonal_gamma].
+                                       u_val,
+                                       rrd->cdp_prep[cdp_idx].
+                                       scratch[CDP_scratch_idx].u_val,
+                                       intercept, seasonal);
+        }
+    }
+#ifdef DEBUG
+    fprintf(stderr, "seasonal coefficient set= %f\n",
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val);
+#endif
+    return 0;
+}
+
+int update_devpredict(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx)
+{
+    /* there really isn't any "update" here; the only reason this information
+     * is stored separately from DEVSEASONAL is to preserve deviation predictions
+     * for a longer duration than one seasonal cycle. */
+    unsigned long seasonal_cdp_idx =
+        (rrd->rra_def[rra_idx].par[RRA_dependent_rra_idx].u_cnt)
+        * (rrd->stat_head->ds_cnt) + ds_idx;
+
+    if (rrd->rra_def[rra_idx].par[RRA_dependent_rra_idx].u_cnt < rra_idx) {
+        /* associated DEVSEASONAL array already updated */
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val
+            =
+            rrd->cdp_prep[seasonal_cdp_idx].
+            scratch[CDP_last_seasonal_deviation].u_val;
+    } else {
+        /* associated DEVSEASONAL not yet updated */
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val
+            =
+            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_seasonal_deviation].
+            u_val;
+    }
+    return 0;
+}
+
+int update_devseasonal(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    rrd_value_t *seasonal_dev,
+    hw_functions_t * functions)
+{
+    rrd_value_t prediction = 0, seasonal_coef = DNAN;
+    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
+
+    /* obtain cdp_prep index for HWPREDICT */
+    unsigned long hw_rra_idx = current_rra->par[RRA_dependent_rra_idx].u_cnt;
+    unsigned long hw_cdp_idx = hw_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
+    unsigned long seasonal_cdp_idx;
+    unival   *coefs = rrd->cdp_prep[hw_cdp_idx].scratch;
+
+    rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].u_val =
+        rrd->cdp_prep[cdp_idx].scratch[CDP_seasonal_deviation].u_val;
+    /* retrieve the next seasonal deviation value, could be NA */
+    rrd->cdp_prep[cdp_idx].scratch[CDP_seasonal_deviation].u_val =
+        seasonal_dev[ds_idx];
+
+    /* retrieve the current seasonal_coef (not to be confused with the
+     * current seasonal deviation). Could make this more readable by introducing
+     * some wrapper functions. */
+    seasonal_cdp_idx =
+        (rrd->rra_def[hw_rra_idx].par[RRA_dependent_rra_idx].u_cnt)
+        * (rrd->stat_head->ds_cnt) + ds_idx;
+    if (rrd->rra_def[hw_rra_idx].par[RRA_dependent_rra_idx].u_cnt < rra_idx)
+        /* SEASONAL array already updated */
+        seasonal_coef =
+            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_last_seasonal].
+            u_val;
+    else
+        /* SEASONAL array not yet updated */
+        seasonal_coef =
+            rrd->cdp_prep[seasonal_cdp_idx].scratch[CDP_hw_seasonal].u_val;
+
+    /* compute the abs value of the difference between the prediction and
+     * observed value */
+    if (hw_rra_idx < rra_idx) {
+        /* associated HWPREDICT has already been updated */
+        if (isnan(coefs[CDP_hw_last_intercept].u_val) ||
+            isnan(coefs[CDP_hw_last_slope].u_val) || isnan(seasonal_coef)) {
+            /* one of the prediction values is uinitialized */
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = DNAN;
+            return 0;
+        } else {
+            prediction =
+                functions->predict(coefs[CDP_hw_last_intercept].u_val,
+                                   coefs[CDP_hw_last_slope].u_val,
+                                   coefs[CDP_last_null_count].u_cnt,
+                                   seasonal_coef);
+        }
+    } else {
+        /* associated HWPREDICT has NOT been updated */
+        if (isnan(coefs[CDP_hw_intercept].u_val) ||
+            isnan(coefs[CDP_hw_slope].u_val) || isnan(seasonal_coef)) {
+            /* one of the prediction values is uinitialized */
+            rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = DNAN;
+            return 0;
+        } else {
+            prediction = functions->predict(coefs[CDP_hw_intercept].u_val,
+                                            coefs[CDP_hw_slope].u_val,
+                                            coefs[CDP_null_count].u_cnt,
+                                            seasonal_coef);
+        }
+    }
+
+    if (isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
+        /* no update, store existing value unchanged, doesn't
+         * matter if it is NA */
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+            rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].u_val;
+    } else
+        if (isnan
+            (rrd->cdp_prep[cdp_idx].scratch[CDP_last_seasonal_deviation].
+             u_val)) {
+        /* initialization */
+#ifdef DEBUG
+        fprintf(stderr, "Initialization of seasonal deviation\n");
+#endif
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+            functions->init_seasonal_deviation(prediction,
+                                               rrd->cdp_prep[cdp_idx].
+                                               scratch[CDP_scratch_idx].
+                                               u_val);
+    } else {
+        /* exponential smoothing update */
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val =
+            functions->seasonal_deviation(rrd->rra_def[rra_idx].
+                                          par[RRA_seasonal_gamma].u_val,
+                                          prediction,
+                                          rrd->cdp_prep[cdp_idx].
+                                          scratch[CDP_scratch_idx].u_val,
+                                          rrd->cdp_prep[cdp_idx].
+                                          scratch
+                                          [CDP_last_seasonal_deviation].
+                                          u_val);
+    }
+    return 0;
+}
+
+/* Check for a failure based on a threshold # of violations within the specified
+ * window. */
+int update_failures(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    hw_functions_t * functions)
+{
+    /* detection of a violation depends on 3 RRAs:
+     * HWPREDICT, SEASONAL, and DEVSEASONAL */
+    rra_def_t *current_rra = &(rrd->rra_def[rra_idx]);
+    unsigned long dev_rra_idx = current_rra->par[RRA_dependent_rra_idx].u_cnt;
+    rra_def_t *dev_rra = &(rrd->rra_def[dev_rra_idx]);
+    unsigned long hw_rra_idx = dev_rra->par[RRA_dependent_rra_idx].u_cnt;
+    rra_def_t *hw_rra = &(rrd->rra_def[hw_rra_idx]);
+    unsigned long seasonal_rra_idx = hw_rra->par[RRA_dependent_rra_idx].u_cnt;
+    unsigned long temp_cdp_idx;
+    rrd_value_t deviation = DNAN;
+    rrd_value_t seasonal_coef = DNAN;
+    rrd_value_t prediction = DNAN;
+    char      violation = 0;
+    unsigned short violation_cnt = 0, i;
+    char     *violations_array;
+
+    /* usual checks to determine the order of the RRAs */
+    temp_cdp_idx = dev_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
+    if (rra_idx < seasonal_rra_idx) {
+        /* DEVSEASONAL not yet updated */
+        deviation =
+            rrd->cdp_prep[temp_cdp_idx].scratch[CDP_seasonal_deviation].u_val;
+    } else {
+        /* DEVSEASONAL already updated */
+        deviation =
+            rrd->cdp_prep[temp_cdp_idx].scratch[CDP_last_seasonal_deviation].
+            u_val;
+    }
+    if (!isnan(deviation)) {
+
+        temp_cdp_idx = seasonal_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
+        if (rra_idx < seasonal_rra_idx) {
+            /* SEASONAL not yet updated */
+            seasonal_coef =
+                rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_seasonal].u_val;
+        } else {
+            /* SEASONAL already updated */
+            seasonal_coef =
+                rrd->cdp_prep[temp_cdp_idx].scratch[CDP_hw_last_seasonal].
+                u_val;
+        }
+        /* in this code block, we know seasonal coef is not DNAN, because deviation is not
+         * null */
+
+        temp_cdp_idx = hw_rra_idx * (rrd->stat_head->ds_cnt) + ds_idx;
+        if (rra_idx < hw_rra_idx) {
+            /* HWPREDICT not yet updated */
+            prediction =
+                functions->predict(rrd->cdp_prep[temp_cdp_idx].
+                                   scratch[CDP_hw_intercept].u_val,
+                                   rrd->cdp_prep[temp_cdp_idx].
+                                   scratch[CDP_hw_slope].u_val,
+                                   rrd->cdp_prep[temp_cdp_idx].
+                                   scratch[CDP_null_count].u_cnt,
+                                   seasonal_coef);
+        } else {
+            /* HWPREDICT already updated */
+            prediction =
+                functions->predict(rrd->cdp_prep[temp_cdp_idx].
+                                   scratch[CDP_hw_last_intercept].u_val,
+                                   rrd->cdp_prep[temp_cdp_idx].
+                                   scratch[CDP_hw_last_slope].u_val,
+                                   rrd->cdp_prep[temp_cdp_idx].
+                                   scratch[CDP_last_null_count].u_cnt,
+                                   seasonal_coef);
+        }
+
+        /* determine if the observed value is a violation */
+        if (!isnan(rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val)) {
+            if (hw_is_violation
+                (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val,
+                 prediction, deviation, current_rra->par[RRA_delta_pos].u_val,
+                 current_rra->par[RRA_delta_neg].u_val)) {
+                violation = 1;
+            }
+        } else {
+            violation = 1;  /* count DNAN values as violations */
+        }
+
+    }
+
+    /* determine if a failure has occurred and update the failure array */
+    violation_cnt = violation;
+    violations_array = (char *) ((void *) rrd->cdp_prep[cdp_idx].scratch);
+    for (i = current_rra->par[RRA_window_len].u_cnt; i > 1; i--) {
+        /* shift */
+        violations_array[i - 1] = violations_array[i - 2];
+        violation_cnt += violations_array[i - 1];
+    }
+    violations_array[0] = violation;
+
+    if (violation_cnt < current_rra->par[RRA_failure_threshold].u_cnt)
+        /* not a failure */
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = 0.0;
+    else
+        rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val = 1.0;
+
+    return (rrd->cdp_prep[cdp_idx].scratch[CDP_scratch_idx].u_val);
+}

diff --git a/src/rrd_hw_update.h b/src/rrd_hw_update.h
new file mode 100644 (file)
index 0000000..e59e2db
--- /dev/null
+++ b/src/rrd_hw_update.h
@@ -0,0 +1,44 @@
+/*****************************************************************************
+ * rrd_hw_update.h  Functions for updating a Holt-Winters RRA
+ ****************************************************************************/
+
+int       update_hwpredict(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    hw_functions_t * functions);
+
+int       update_seasonal(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    rrd_value_t *seasonal_coef,
+    hw_functions_t * functions);
+
+int       update_devpredict(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx);
+
+int       update_devseasonal(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    rrd_value_t *seasonal_dev,
+    hw_functions_t * functions);
+
+int       update_failures(
+    rrd_t *rrd,
+    unsigned long cdp_idx,
+    unsigned long rra_idx,
+    unsigned long ds_idx,
+    unsigned short CDP_scratch_idx,
+    hw_functions_t * functions);

diff --git a/src/rrd_info.c b/src/rrd_info.c
index fe65acac9062a0625edf680db6509a2ffe1b042d..f1866bef03d9b007b77485a73290434ff8f9d089 100644 (file)
--- a/src/rrd_info.c
+++ b/src/rrd_info.c
@@ -192,6 +192,7 @@ info_t   *rrd_info_r(
 
         switch (current_cf) {
         case CF_HWPREDICT:
+        case CF_MHWPREDICT:
             info.u_val = rrd.rra_def[i].par[RRA_hw_alpha].u_val;
             cd = info_push(cd, sprintf_alloc("rra[%d].alpha", i), RD_I_VAL,
                            info);
@@ -231,6 +232,7 @@ info_t   *rrd_info_r(
         for (ii = 0; ii < rrd.stat_head->ds_cnt; ii++) {
             switch (current_cf) {
             case CF_HWPREDICT:
+            case CF_MHWPREDICT:
                 info.u_val =
                     rrd.cdp_prep[i * rrd.stat_head->ds_cnt +
                                  ii].scratch[CDP_hw_intercept].u_val;

diff --git a/src/rrd_restore.c b/src/rrd_restore.c
index 368f7bcfbb041ee6de6d4ead1ce5a2d414351114..68a037fe6139f1daf949c1489a83dc3e8c08340d 100644 (file)
--- a/src/rrd_restore.c
+++ b/src/rrd_restore.c
@@ -366,6 +366,7 @@ int xml2rrd(
             } else {
                 switch (cf_conv(rrd->rra_def[rra_index].cf_nam)) {
                 case CF_HWPREDICT:
+                case CF_MHWPREDICT:
                     read_tag(&ptr2, "hw_alpha", "%lf",
                              &(rrd->rra_def[rra_index].par[RRA_hw_alpha].
                                u_val));
@@ -456,6 +457,7 @@ int xml2rrd(
                                         i].scratch[CDP_secondary_val].u_val));
                     switch (cf_conv(rrd->rra_def[rra_index].cf_nam)) {
                     case CF_HWPREDICT:
+                    case CF_MHWPREDICT:
                         read_tag(&ptr2, "intercept", "%lf",
                                  &(rrd->
                                    cdp_prep[rrd->stat_head->ds_cnt *

diff --git a/src/rrd_tune.c b/src/rrd_tune.c
index 91f6409ada384e79093e9ea28d20b90ff92284ce..160b32a5e66a961845e9655e98e77018cca59a68 100644 (file)
--- a/src/rrd_tune.c
+++ b/src/rrd_tune.c
@@ -237,14 +237,20 @@ int rrd_tune(
             break;
         case 'x':
             if (set_hwarg(&rrd, CF_HWPREDICT, RRA_hw_alpha, optarg)) {
-                rrd_free(&rrd);
-                return -1;
+                if (set_hwarg(&rrd, CF_MHWPREDICT, RRA_hw_alpha, optarg)) {
+                    rrd_free(&rrd);
+                    return -1;
+                }
+                rrd_clear_error();
             }
             break;
         case 'y':
             if (set_hwarg(&rrd, CF_HWPREDICT, RRA_hw_beta, optarg)) {
-                rrd_free(&rrd);
-                return -1;
+                if (set_hwarg(&rrd, CF_MHWPREDICT, RRA_hw_beta, optarg)) {
+                    rrd_free(&rrd);
+                    return -1;
+                }
+                rrd_clear_error();
             }
             break;
         case 'z':

diff --git a/src/rrd_update.c b/src/rrd_update.c
index abe9704a2bf1080ba75ee112c47e6d9c04fe571a..1feecac58f8ba2862d68ba30badbf4f7ec866de6 100644 (file)
--- a/src/rrd_update.c
+++ b/src/rrd_update.c
@@ -1198,6 +1198,7 @@ int _rrd_update(
                                     u_val = seasonal_coef[ii];
                                 break;
                             case CF_HWPREDICT:
+                            case CF_MHWPREDICT:
                                 /* need to update the null_count and last_null_count.
                                  * even do this for non-DNAN pdp_temp because the
                                  * algorithm is not learning from batch updates. */