Added double exponential smoothing implementation

4 files changed, 155 insertions, 2 deletions

diff --git a/example/ndpiReader.c b/example/ndpiReader.c
index 3e9c54f8d..e41cbb21a 100644
--- a/example/ndpiReader.c
+++ b/example/ndpiReader.c
@@ -3939,7 +3939,7 @@ void sesUnitTest() {
   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");
+  FILE *fd = fopen("/tmp/ses_result.csv", "w");
   
   assert(ndpi_ses_init(&ses, alpha, 0.05) == 0);
   
@@ -3974,6 +3974,77 @@ void sesUnitTest() {
 
 /* *********************************************** */
 
+void desUnitTest() {
+  struct ndpi_des_struct des;
+  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, beta = 0.5;
+  FILE *fd = fopen("/tmp/des_result.csv", "w");
+  
+  assert(ndpi_des_init(&des, alpha, beta, 0.05) == 0);
+  
+  if(trace) {
+    printf("\nDouble Exponential Smoothing [alpha: %.1f][beta: %.1f]\n", alpha, beta);
+
+    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_des_add_value(&des, 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);
+}
+
+/* *********************************************** */
+
 void hwUnitTest3() {
   struct ndpi_hw_struct hw;
   u_int num_learning_points = 3;
@@ -4063,7 +4134,8 @@ int original_main(int argc, char **argv) {
 #endif
 
     sesUnitTest();
-
+    desUnitTest();
+    
     /* Internal checks */    
     // binUnitTest();
     //hwUnitTest();
diff --git a/src/include/ndpi_api.h.in b/src/include/ndpi_api.h.in
index cb5e45931..118240916 100644
--- a/src/include/ndpi_api.h.in
+++ b/src/include/ndpi_api.h.in
@@ -1479,6 +1479,11 @@ extern "C" {
 
   /* ******************************* */
 
+  int ndpi_des_init(struct ndpi_des_struct *des, double alpha, double beta, float significance);
+  int ndpi_des_add_value(struct ndpi_des_struct *des, 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 de9a404c8..aac7b56db 100644
--- a/src/include/ndpi_typedefs.h
+++ b/src/include/ndpi_typedefs.h
@@ -1672,6 +1672,15 @@ struct ndpi_ses_struct {
   double sum_square_error, last_forecast, last_value;
 };
 
+struct ndpi_des_struct {
+  struct {
+    double alpha, beta, ro;
+  } params;
+
+  u_int32_t num_values;
+  double sum_square_error, last_forecast, last_trend, last_value;
+};
+
 /* **************************************** */
 
 /* Prototype used to define custom DGA detection function */
diff --git a/src/lib/ndpi_analyze.c b/src/lib/ndpi_analyze.c
index 19938b14f..3863dc910 100644
--- a/src/lib/ndpi_analyze.c
+++ b/src/lib/ndpi_analyze.c
@@ -1214,3 +1214,70 @@ int ndpi_ses_add_value(struct ndpi_ses_struct *ses, const u_int32_t _value, doub
 
   return(rc);
 }
+
+/* *********************************************************** */
+/* *********************************************************** */
+
+/*
+  Double Exponential Smoothing
+*/
+
+int ndpi_des_init(struct ndpi_des_struct *des, double alpha, double beta, float significance) {
+  memset(des, 0, sizeof(struct ndpi_des_struct));
+
+  des->params.alpha = alpha;
+
+  if((significance < 0) || (significance > 1)) significance = 0.05;
+  des->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
+   des:         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_des_add_value(struct ndpi_des_struct *des, const u_int32_t _value, double *forecast, double *confidence_band) {
+  double value = (double)_value, error;
+  int rc;
+
+  if(des->num_values == 0)
+    *forecast = value, des->last_trend = 0;
+  else {
+    *forecast = (des->params.alpha * value) + ((1 - des->params.alpha) * (des->last_forecast + des->last_trend));
+    des->last_trend = (des->params.beta * (*forecast - des->last_forecast)) + ((1 - des->params.beta) * des->last_trend);
+  }
+  
+  error  = value - *forecast;
+  des->sum_square_error += error * error;
+  
+  if(des->num_values > 0) {
+    double sq = sqrt(des->sum_square_error / (des->num_values+1));
+    
+    *confidence_band = des->params.ro * sq;
+    rc = 1;
+  } else
+    *confidence_band = 0, rc = 0;
+
+  des->num_values++, des->last_value = value, des->last_forecast = *forecast;
+  
+#ifdef DES_DEBUG
+  printf("[num_values: %u][[error: %.3f][forecast: %.3f][trend: %.3f[sqe: %.3f][sq: %.3f][confidence_band: %.3f]\n",
+	 des->num_values, error, *forecast, des->last_trend, des->sum_square_error, sq, *confidence_band);
+#endif
+  
+  return(rc);
+}