Added single exponential smoothing API

int ndpi_ses_init(struct ndpi_ses_struct *ses, double alpha, float significance);
int ndpi_ses_add_value(struct ndpi_ses_struct *ses, const u_int32_t _value, double *forecast, double *confidence_band);

4 files changed, 235 insertions, 29 deletions

diff --git a/example/ndpiReader.c b/example/ndpiReader.c
index aec78d778..3e9c54f8d 100644
--- a/example/ndpiReader.c
+++ b/example/ndpiReader.c
@@ -3722,20 +3722,53 @@ void rulesUnitTest() {
 void rsiUnitTest() {
   struct ndpi_rsi_struct s;
   unsigned int v[] = {
+    31,
+    87,
+    173,
+    213,
+    223,
+    230,
+    238,
+    245,
+    251,
+    151,
+    259,
+    261,
+    264,
+    264,
+    270,
+    273,
+    288,
+    288,
+    304,
+    304,
+    350,
+    384,
+    423,
+    439,
+    445,
+    445,
+    445,
+    445
+  };
+
+  unsigned int v1[] = {
+
     2227, 2219, 2208, 2217, 2218, 2213, 2223, 2243, 2224, 2229,
     2215, 2239, 2238, 2261, 2336, 2405, 2375, 2383, 2395, 2363,
     2382, 2387, 2365, 2319, 2310, 2333, 2268, 2310, 2240, 2217,
   };
   u_int i, n = sizeof(v) / sizeof(unsigned int);
-
+  u_int debug = 0;
+  
   assert(ndpi_alloc_rsi(&s, 8) == 0);
   
   for(i=0; i<n; i++) {
     float rsi = ndpi_rsi_add_value(&s, v[i]);
 
-#if 0
-    printf("%2d) RSI = %f\n", i, rsi);
-#endif
+
+    if(debug)
+      printf("%2d) RSI = %f\n", i, rsi);
   }
   
   ndpi_free_rsi(&s);
@@ -3796,7 +3829,6 @@ void hwUnitTest() {
 
 void hwUnitTest2() {
   struct ndpi_hw_struct hw;
-  u_int num_learning_points = 12;
   u_int8_t trace = 1;
   double v[] = {
     31.908466339111,
@@ -3814,7 +3846,7 @@ void hwUnitTest2() {
     264.78540039062,
     264.78540039062,
     270.47451782227,
-    273.3671875,
+    173.3671875,
     288.34222412109,
     288.34222412109,
     304.24795532227,
@@ -3828,13 +3860,20 @@ void hwUnitTest2() {
     445.05981445312,
     445.05981445312
   };
+  u_int num_learning_points = 1;
   u_int i, num = sizeof(v) / sizeof(double);
-  float alpha = 0.5, beta = 0.5, gamma = 0.1;
-
-  assert(ndpi_hw_init(&hw, num_learning_points, 0 /* 0=multiplicative, 1=additive */, alpha, beta, gamma, 0.05) == 0);
+  float alpha = 0.9, beta = 0.5, gamma = 1;
+  FILE *fd = fopen("/tmp/result.csv", "w");
   
-  if(trace)
+  assert(ndpi_hw_init(&hw, num_learning_points, 0 /* 0=multiplicative, 1=additive */,
+		      alpha, beta, gamma, 0.05) == 0);
+  
+  if(trace) {
     printf("\nHolt-Winters [alpha: %.1f][beta: %.1f][gamma: %.1f]\n", alpha, beta, gamma);
+
+    if(fd)
+      fprintf(fd, "index;value;prediction;lower;upper;anomaly\n");
+  }
   
   for(i=0; i<num; i++) {
     double prediction, confidence_band;
@@ -3843,13 +3882,94 @@ void hwUnitTest2() {
     
     lower = prediction - confidence_band, upper = prediction + confidence_band;
     
-    if(trace)
+    if(trace) {
       printf("%2u)\t%12.3f\t%.3f\t%12.3f\t%12.3f\t %s [%.3f]\n", i, v[i], prediction, lower, upper,
 	     ((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY",
 	     confidence_band);
+
+      if(fd)
+	fprintf(fd, "%u;%.0f;%.0f;%.0f;%.0f;%s\n",
+		i, v[i], prediction, lower, upper,
+		((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY");      
+    }
   }
-  
+
+  if(fd) fclose(fd);
+	   
   ndpi_hw_free(&hw);
+
+  //exit(0);
+}
+
+/* *********************************************** */
+
+void sesUnitTest() {
+  struct ndpi_ses_struct ses;
+  u_int8_t trace = 0;
+  double v[] = {
+    31.908466339111,
+    87.339714050293,
+    173.47660827637,
+    213.92568969727,
+    223.32124328613,
+    230.60134887695,
+    238.09457397461,
+    245.8137512207,
+    251.09228515625,
+    251.09228515625,
+    259.21997070312,
+    261.98754882812,
+    264.78540039062,
+    264.78540039062,
+    270.47451782227,
+    173.3671875,
+    288.34222412109,
+    288.34222412109,
+    304.24795532227,
+    304.24795532227,
+    350.92227172852,
+    384.54431152344,
+    423.25942993164,
+    439.43322753906,
+    445.05981445312,
+    445.05981445312,
+    445.05981445312,
+    445.05981445312
+  };
+  u_int num_learning_points = 1;
+  u_int i, num = sizeof(v) / sizeof(double);
+  float alpha = 0.9;
+  FILE *fd = fopen("/tmp/result.csv", "w");
+  
+  assert(ndpi_ses_init(&ses, alpha, 0.05) == 0);
+  
+  if(trace) {
+    printf("\nSingle Exponential Smoothing [alpha: %.1f]\n", alpha);
+
+    if(fd)
+      fprintf(fd, "index;value;prediction;lower;upper;anomaly\n");
+  }
+  
+  for(i=0; i<num; i++) {
+    double prediction, confidence_band;
+    double lower, upper;
+    int rc = ndpi_ses_add_value(&ses, v[i], &prediction, &confidence_band);
+    
+    lower = prediction - confidence_band, upper = prediction + confidence_band;
+    
+    if(trace) {
+      printf("%2u)\t%12.3f\t%.3f\t%12.3f\t%12.3f\t %s [%.3f]\n", i, v[i], prediction, lower, upper,
+	     ((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY",
+	     confidence_band);
+
+      if(fd)
+	fprintf(fd, "%u;%.0f;%.0f;%.0f;%.0f;%s\n",
+		i, v[i], prediction, lower, upper,
+		((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY");      
+    }
+  }
+
+  if(fd) fclose(fd);
 }
 
 /* *********************************************** */
@@ -3938,8 +4058,12 @@ int original_main(int argc, char **argv) {
       return(-1);
     }
 
-    // hwUnitTest2();
-    
+#ifdef HW_TEST
+    hwUnitTest2();
+#endif
+
+    sesUnitTest();
+
     /* Internal checks */    
     // binUnitTest();
     //hwUnitTest();
diff --git a/src/include/ndpi_api.h.in b/src/include/ndpi_api.h.in
index c12910609..cb5e45931 100644
--- a/src/include/ndpi_api.h.in
+++ b/src/include/ndpi_api.h.in
@@ -1474,6 +1474,11 @@ extern "C" {
 
   /* ******************************* */
 
+  int ndpi_ses_init(struct ndpi_ses_struct *ses, double alpha, float significance);
+  int ndpi_ses_add_value(struct ndpi_ses_struct *ses, const u_int32_t _value, double *forecast, double *confidence_band);
+
+  /* ******************************* */
+
   int   ndpi_jitter_init(struct ndpi_jitter_struct *hw, u_int16_t num_periods);
   void  ndpi_jitter_free(struct ndpi_jitter_struct *hw);
   float ndpi_jitter_add_value(struct ndpi_jitter_struct *s, const float value);
diff --git a/src/include/ndpi_typedefs.h b/src/include/ndpi_typedefs.h
index 35b3b95dd..de9a404c8 100644
--- a/src/include/ndpi_typedefs.h
+++ b/src/include/ndpi_typedefs.h
@@ -1663,6 +1663,17 @@ struct ndpi_hw_struct {
   double    *s;
 };
 
+struct ndpi_ses_struct {
+  struct {
+    double alpha, ro;
+  } params;
+
+  u_int32_t num_values;
+  double sum_square_error, last_forecast, last_value;
+};
+
+/* **************************************** */
+
 /* Prototype used to define custom DGA detection function */
 typedef int (*ndpi_custom_dga_predict_fctn)(const char* domain, int domain_length);
 
diff --git a/src/lib/ndpi_analyze.c b/src/lib/ndpi_analyze.c
index 3856e0f5f..19938b14f 100644
--- a/src/lib/ndpi_analyze.c
+++ b/src/lib/ndpi_analyze.c
@@ -852,6 +852,8 @@ void ndpi_free_rsi(struct ndpi_rsi_struct *s) {
 
 /* ************************************* */
 
+// #define DEBUG_RSI
+
 /*
   This function adds a new value and returns the computed RSI, or -1
   if there are too few points (< num_learning_values)
@@ -941,6 +943,7 @@ double ndpi_avg_inline(u_int32_t *v, u_int num) {
 }
 
 /* *********************************************************** */
+/* *********************************************************** */
 
 /*
   Initializes Holt-Winters with Confidence Interval
@@ -951,7 +954,7 @@ double ndpi_avg_inline(u_int32_t *v, u_int num) {
    additive     If set to 1 will use the Holt-Winters additive seasonal (should be the default), otherwise the multiplicative seasonal.
    alpha        Level: specifies the coefficient for the level smoothing. Range 0..1. The higher α, the faster the method forgets old values
    beta         Trend: specifies the coefficient for the trend smoothing. Range 0..1.
-   gamma        Seasonal: specifies the coefficient for the seasonal smoothing. Range 0..1. With gamma = 0, seaasonal correction is not used.
+   gamma        Seasonal: specifies the coefficient for the seasonal smoothing. Range 0..1. With gamma = 0, seasonal correction is not used.
 
    significance Significance level for the forecats sed for computing lower and upper bands. Range 0..1. Typical values 0.05 or less.
                 See https://en.wikipedia.org/wiki/Statistical_significance
@@ -978,8 +981,6 @@ int ndpi_hw_init(struct ndpi_hw_struct *hw,
   if((significance < 0) || (significance > 1)) significance = 0.05;
   hw->params.ro         = ndpi_normal_cdf_inverse(1 - (significance / 2.));
 
-  hw->params.ro = 1.95996398454005;
-  
   if((hw->y = (u_int32_t*)ndpi_calloc(hw->params.num_season_periods, sizeof(u_int32_t))) == NULL)
     return(-1);
 
@@ -1028,13 +1029,13 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
     u_int idx     = hw->num_values % hw->params.num_season_periods;
     double prev_u, prev_v, prev_s, value  = (double)_value;
     double sq, error;
-    
+
     if(hw->num_values == hw->params.num_season_periods) {
       double avg = ndpi_avg_inline(hw->y, hw->params.num_season_periods);
       u_int i;
 
       if(avg == 0) avg = 1; /* Avoid divisions by zero */
-      
+
       for(i=0; i<hw->params.num_season_periods; i++)
 	hw->s[i] = hw->y[i] / avg;
 
@@ -1043,7 +1044,7 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
 	hw->u = 0;
       else
 	hw->u = _value / hw->s[i];
-      
+
       hw->v = 0;
       ndpi_free(hw->y);
       hw->y = NULL;
@@ -1051,12 +1052,12 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
 
     idx     = hw->num_values % hw->params.num_season_periods;
     prev_u = hw->u, prev_v = hw->v, prev_s = hw->s[idx];
-     
+
     if(prev_s != 0)
       hw->u = ((hw->params.alpha * value) / prev_s)  + ( 1 - hw->params.alpha) * (hw->u + hw->v);
     else
       hw->u = 0; /* Avoid divisions by zero */
-    
+
     hw->v = (hw->params.beta   * (hw->u - prev_u)) + ((1 - hw->params.beta ) * hw->v);
 
     if(hw->u != 0)
@@ -1080,7 +1081,7 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
 	   *forecast, hw->sum_square_error,
 	   sq, *confidence_band);
 #endif
-    
+
     hw->num_values++, idx = (idx + 1) % hw->params.num_season_periods;
 
     return(1); /* We're in business: forecast is meaningful now */
@@ -1092,7 +1093,7 @@ int ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double
 
 /*
   Jitter calculator
-  
+
   Used to determine how noisy is a signal
 */
 
@@ -1100,7 +1101,7 @@ int ndpi_jitter_init(struct ndpi_jitter_struct *s, u_int16_t num_learning_values
   memset(s, 0, sizeof(struct ndpi_jitter_struct));
 
   if(num_learning_values < 2) num_learning_values = 2;
-  
+
   s->empty = 1, s->num_values = num_learning_values;
   s->observations = (float*)ndpi_calloc(num_learning_values, sizeof(float));
 
@@ -1108,7 +1109,7 @@ int ndpi_jitter_init(struct ndpi_jitter_struct *s, u_int16_t num_learning_values
     s->last_value = 0;
     return(0);
   } else
-    return(-1);  
+    return(-1);
 }
 
 /* ************************************* */
@@ -1132,7 +1133,7 @@ float ndpi_jitter_add_value(struct ndpi_jitter_struct *s, const float value) {
     s->observations[s->next_index] = val;
     s->jitter_total += val;
   }
-  
+
   s->last_value = value, s->next_index = (s->next_index + 1) % s->num_values;
   if(s->next_index == 0) s->jitter_ready = 1; /* We have completed one round */
 
@@ -1141,10 +1142,75 @@ float ndpi_jitter_add_value(struct ndpi_jitter_struct *s, const float value) {
 	 value, val, s->jitter_total,
 	 s->jitter_ready ? (s->jitter_total / s->num_values) : -1);
 #endif
-  
+
   if(!s->jitter_ready)
     return(-1); /* Too early */
-  else 
+  else
     return(s->jitter_total / s->num_values);
 }
 
+
+/* *********************************************************** */
+/* *********************************************************** */
+
+/*
+  Single Exponential Smoothing
+*/
+
+int ndpi_ses_init(struct ndpi_ses_struct *ses, double alpha, float significance) {
+  memset(ses, 0, sizeof(struct ndpi_ses_struct));
+
+  ses->params.alpha = alpha;
+
+  if((significance < 0) || (significance > 1)) significance = 0.05;
+  ses->params.ro         = ndpi_normal_cdf_inverse(1 - (significance / 2.));
+
+  return(0);
+}
+
+/* *********************************************************** */
+
+/*
+   Returns the forecast and the band (forecast +/- band are the upper and lower values)
+
+   Input
+   ses:         Datastructure previously initialized
+   value        The value to add to the measurement
+
+   Output
+   forecast         The forecasted value
+   confidence_band  The value +/- on which the value should fall is not an anomaly
+
+   Return code
+   0                Too early: we're still in the learning phase. Output values are zero.
+   1                Normal processing: forecast and confidence_band are meaningful
+*/
+int ndpi_ses_add_value(struct ndpi_ses_struct *ses, const u_int32_t _value, double *forecast, double *confidence_band) {
+  double value = (double)_value, error;
+  int rc;
+
+  if(ses->num_values == 0)
+    *forecast = value;
+  else
+    *forecast = (ses->params.alpha * (ses->last_value - ses->last_forecast)) + ses->last_forecast;
+
+  error  = value - *forecast;
+  ses->sum_square_error += error * error;
+
+  if(ses->num_values > 0) {
+    double sq = sqrt(ses->sum_square_error / (ses->num_values+1));
+
+    *confidence_band = ses->params.ro * sq;
+    rc = 1;
+  } else
+    *confidence_band = 0, rc = 0;
+
+  ses->num_values++, ses->last_value = value, ses->last_forecast = *forecast;
+
+#ifdef SES_DEBUG
+  printf("[num_values: %u][[error: %.3f][forecast: %.3f][sqe: %.3f][sq: %.3f][confidence_band: %.3f]\n",
+	   ses->num_values, error, *forecast, ses->sum_square_error, sq, *confidence_band);
+#endif
+
+  return(rc);
+}