Added timeseries forecasting support implementing Holt-Winters with confidence interval

New API calls added
- ndpi_hw_init()
- ndpi_hw_add_value()
- ndpi_hw_free()

4 files changed, 286 insertions, 69 deletions

diff --git a/example/ndpiReader.c b/example/ndpiReader.c
index 5ece4e414..758d59ec9 100644
--- a/example/ndpiReader.c
+++ b/example/ndpiReader.c
@@ -3733,6 +3733,40 @@ void hashUnitTest() {
 
 /* *********************************************** */
 
+void hwUnitTest() {
+  struct ndpi_hw_struct hw;
+  double v[] = { 10, 14, 8, 25, 16, 22, 14, 35, 15, 27, 218, 40, 28, 40, 25, 65 };
+  u_int i, j, num = sizeof(v) / sizeof(double);
+  u_int num_learning_points = 2;
+  u_int8_t trace = 0;
+  
+  for(j=0; j<2; j++) {
+    assert(ndpi_hw_init(&hw, num_learning_points, j /* 0=multiplicative, 1=additive */, 0.9, 0.9, 0.1, 0.05) == 0);
+
+    if(trace)
+      printf("\nHolt-Winters %s method\n", (j == 0) ? "multiplicative" : "additive");
+    
+    for(i=0; i<num; i++) {
+      double prediction, confidence_band;
+      double lower, upper;
+      int rc = ndpi_hw_add_value(&hw, v[i], &prediction, &confidence_band);
+      
+      lower = prediction - confidence_band, upper = prediction + confidence_band;
+
+      if(trace)
+	printf("%2u)\t%.3f\t%.3f\t%.3f\t%.3f\t %s [%.3f]\n", i, v[i], prediction, lower, upper,
+	       ((rc == 0) || ((v[i] >= lower) && (v[i] <= upper))) ? "OK" : "ANOMALY",
+	       confidence_band);
+    }
+    
+    ndpi_hw_free(&hw);
+  }
+  
+  exit(0);
+}
+
+/* *********************************************** */
+
 /**
    @brief MAIN FUNCTION
 **/
@@ -3753,8 +3787,9 @@ int orginal_main(int argc, char **argv) {
 
     if(ndpi_info_mod == NULL) return -1;
 
-    /* Internal checks */
+    /* Internal checks */    
     // binUnitTest();
+    hwUnitTest();
     rsiUnitTest();
     hashUnitTest();
     dgaUnitTest();
diff --git a/src/include/ndpi_api.h.in b/src/include/ndpi_api.h.in
index 9055546fc..5e2082bbd 100644
--- a/src/include/ndpi_api.h.in
+++ b/src/include/ndpi_api.h.in
@@ -359,7 +359,7 @@ extern "C" {
 				 struct ndpi_flow_struct *flow,
 				 NDPI_SELECTION_BITMASK_PROTOCOL_SIZE *ndpi_selection_packet);
 
-  
+
   /**
    * Query the pointer to the layer 4 packet
    *
@@ -941,7 +941,7 @@ extern "C" {
 
   const char* ndpi_http_method2str(ndpi_http_method m);
   ndpi_http_method ndpi_http_str2method(const char* method, u_int16_t method_len);
-  
+
   /* ptree (trie) API */
   ndpi_ptree_t* ndpi_ptree_create(void);
   int ndpi_ptree_insert(ndpi_ptree_t *tree, const ndpi_ip_addr_t *addr, u_int8_t bits, u_int64_t user_data);
@@ -961,7 +961,7 @@ extern "C" {
 
   /**
    * Initialize a serializer handle (allocated by the caller).
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param fmt The serialization format (ndpi_serialization_format_json, ndpi_serialization_format_tlv, ndpi_serialization_format_csv)
    * @return 0 on success, a negative number otherwise
    */
@@ -969,28 +969,28 @@ extern "C" {
 
   /**
    * Initialize a serializer handle. Same as ndpi_init_serializer, but with some low-level settings.
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param fmt The serialization format (ndpi_serialization_format_json, ndpi_serialization_format_tlv, ndpi_serialization_format_csv)
    * @param buffer_size The initial internal buffer_size
-   * @return 0 on success, a negative number otherwise 
+   * @return 0 on success, a negative number otherwise
    */
   int ndpi_init_serializer_ll(ndpi_serializer *serializer, ndpi_serialization_format fmt, u_int32_t buffer_size);
 
   /**
    * Release all allocated data structure.
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    */
   void ndpi_term_serializer(ndpi_serializer *serializer);
 
   /**
    * Reset the serializer (cleanup the internal buffer to start a new serialization)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    */
   void ndpi_reset_serializer(ndpi_serializer *serializer);
 
   /**
    * Serialize a 32-bit unsigned int key and a 32-bit unsigned int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -998,8 +998,8 @@ extern "C" {
   int ndpi_serialize_uint32_uint32(ndpi_serializer *serializer, u_int32_t key, u_int32_t value);
 
   /**
-   * Serialize a 32-bit unsigned int key and a 64-bit unsigned int value 
-   * @param serializer The serializer handle 
+   * Serialize a 32-bit unsigned int key and a 64-bit unsigned int value
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1008,7 +1008,7 @@ extern "C" {
 
   /**
    * Serialize a 32-bit unsigned int key and a 32-bit signed int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1017,7 +1017,7 @@ extern "C" {
 
   /**
    * Serialize a 32-bit unsigned int key and a 64-bit signed int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1026,7 +1026,7 @@ extern "C" {
 
   /**
    * Serialize a 32-bit unsigned int key and a float value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @param format The float value format
@@ -1036,7 +1036,7 @@ extern "C" {
 
   /**
    * Serialize a 32-bit unsigned int key and a string value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1045,7 +1045,7 @@ extern "C" {
 
   /**
    * Serialize a 32-bit unsigned int key and a boolean value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1054,7 +1054,7 @@ extern "C" {
 
   /**
    * Serialize an unterminated string key and a 32-bit signed int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @param value The field value
@@ -1064,7 +1064,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a 32-bit signed int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1073,7 +1073,7 @@ extern "C" {
 
   /**
    * Serialize an unterminated string key and a 64-bit signed int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @param value The field value
@@ -1083,7 +1083,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a 64-bit signed int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1092,7 +1092,7 @@ extern "C" {
 
   /**
    * Serialize an unterminated string key and a 32-bit unsigned int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @param value The field value
@@ -1102,7 +1102,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a 32-bit unsigned int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1111,7 +1111,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a float value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @param format The float format
@@ -1121,7 +1121,7 @@ extern "C" {
 
   /**
    * Serialize an unterminated string key and a 64-bit unsigned int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @param value The field value
@@ -1131,7 +1131,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a 64-bit unsigned int value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1140,7 +1140,7 @@ extern "C" {
 
   /**
    * Serialize an unterminated string key and an unterminated string value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @param value The field value
@@ -1151,7 +1151,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a string value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1160,7 +1160,7 @@ extern "C" {
 
   /**
    * Serialize a string key and an unterminated string value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @param vlen The value length
@@ -1170,7 +1170,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a raw value (this is a string which is added to the JSON without any quote or escaping)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @param vlen The value length
@@ -1180,7 +1180,7 @@ extern "C" {
 
   /**
    * Serialize an unterminated string key and a float value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @param value The field value
@@ -1190,8 +1190,8 @@ extern "C" {
   int ndpi_serialize_binary_float(ndpi_serializer *_serializer, const char *key, u_int16_t klen, float value, const char *format /* e.f. "%.2f" */);
 
   /**
-   * Serialize a string key and a a float value 
-   * @param serializer The serializer handle 
+   * Serialize a string key and a a float value
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @param format The float format
@@ -1201,7 +1201,7 @@ extern "C" {
 
   /**
    * Serialize a string key and a boolean value
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1210,7 +1210,7 @@ extern "C" {
 
   /**
    * Serialize a raw record in an array (this is a low-level function and its use is not recommended)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param record The record value
    * @param record_len The record length
    * @return 0 on success, a negative number otherwise
@@ -1220,15 +1220,15 @@ extern "C" {
   /**
    * Serialize an End-Of-Record (the current object becomes is terminated and added to an array,
    * and a new object is created where the next items will be added)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @return 0 on success, a negative number otherwise
    */
   int ndpi_serialize_end_of_record(ndpi_serializer *serializer);
 
   /**
-   * Serialize the start of a list with an unterminated string key, where the next serialized items 
+   * Serialize the start of a list with an unterminated string key, where the next serialized items
    * will be added (note: keys for the new items are ignored)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @return 0 on success, a negative number otherwise
@@ -1238,7 +1238,7 @@ extern "C" {
   /**
    * Serialize the start of a list, where the next serialized items will be added (note: keys for
    * the new items are ignored)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1247,14 +1247,14 @@ extern "C" {
 
   /**
    * Serialize the end of a list
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @return 0 on success, a negative number otherwise
    */
   int ndpi_serialize_end_of_list(ndpi_serializer *serializer);
 
   /**
-   * Serialize the start of a block with an unterminated string key 
-   * @param serializer The serializer handle 
+   * Serialize the start of a block with an unterminated string key
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param klen The key length
    * @return 0 on success, a negative number otherwise
@@ -1262,16 +1262,16 @@ extern "C" {
   int ndpi_serialize_start_of_block_binary(ndpi_serializer *_serializer, const char *key, u_int16_t klen);
 
   /**
-   * Serialize the start of a block with a string key 
-   * @param serializer The serializer handle 
+   * Serialize the start of a block with a string key
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @return 0 on success, a negative number otherwise
    */
   int ndpi_serialize_start_of_block(ndpi_serializer *serializer, const char *key);
 
   /**
-   * Serialize the end of a block 
-   * @param serializer The serializer handle 
+   * Serialize the end of a block
+   * @param serializer The serializer handle
    * @param key The field name or ID
    * @param value The field value
    * @return 0 on success, a negative number otherwise
@@ -1280,7 +1280,7 @@ extern "C" {
 
   /**
    * Return the serialized buffer
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param buffer_len The buffer length (out)
    * @return The buffer
    */
@@ -1288,7 +1288,7 @@ extern "C" {
 
   /**
    * Return the current serialized buffer length
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @return The buffer length
    */
   u_int32_t ndpi_serializer_get_buffer_len(ndpi_serializer *serializer);
@@ -1302,29 +1302,29 @@ extern "C" {
 
   /**
    * Change the serializer buffer length
-   * @param serializer The serializer handle 
-   * @param l The new buffer length 
+   * @param serializer The serializer handle
+   * @param l The new buffer length
    * @return 0 on success, a negative number otherwise
    */
   int ndpi_serializer_set_buffer_len(ndpi_serializer *serializer, u_int32_t l);
 
   /**
    * Return the configured serialization format
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @return The serialization format
    */
   ndpi_serialization_format ndpi_serializer_get_format(ndpi_serializer *serializer);
 
   /**
    * Set the CSV separator
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param separator The separator
    */
   void ndpi_serializer_set_csv_separator(ndpi_serializer *serializer, char separator);
 
   /**
    * Return the header automatically built from keys (CSV only)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    * @param buffer_len The buffer length (out)
    * @return The header
    */
@@ -1332,13 +1332,13 @@ extern "C" {
 
   /**
    * Create a snapshot of the internal buffer for later rollback (ndpi_serializer_rollback_snapshot)
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    */
   void ndpi_serializer_create_snapshot(ndpi_serializer *serializer);
 
   /**
    * Rollback to the latest snapshot
-   * @param serializer The serializer handle 
+   * @param serializer The serializer handle
    */
   void ndpi_serializer_rollback_snapshot(ndpi_serializer *serializer);
 
@@ -1394,7 +1394,14 @@ extern "C" {
   int ndpi_alloc_rsi(struct ndpi_rsi_struct *s, u_int16_t num_learning_values);
   void  ndpi_free_rsi(struct ndpi_rsi_struct *s);
   float ndpi_rsi_add_value(struct ndpi_rsi_struct *s, const u_int32_t value);
-  
+
+  /* ******************************* */
+
+  int  ndpi_hw_init(struct ndpi_hw_struct *hw, u_int16_t num_periods, u_int8_t additive_seeasonal,
+		    double alpha, double beta, double gamma, float significance);
+  void ndpi_hw_free(struct ndpi_hw_struct *hw);
+  int  ndpi_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t value, double *forecast,  double *confidence_band);
+
   /* ******************************* */
 
   const char* ndpi_data_ratio2str(float ratio);
@@ -1417,7 +1424,7 @@ extern "C" {
   int ndpi_hll_init(struct ndpi_hll *hll, u_int8_t bits);
   void ndpi_hll_destroy(struct ndpi_hll *hll);
   void ndpi_hll_reset(struct ndpi_hll *hll);
-  
+
   /* Add values */
   void ndpi_hll_add(struct ndpi_hll *hll, const char *data, size_t data_len);
   void ndpi_hll_add_number(struct ndpi_hll *hll, u_int32_t value) ;
@@ -1442,7 +1449,7 @@ extern "C" {
 			struct ndpi_bin *centroids);
 
   /* ******************************* */
-  
+
   u_int32_t ndpi_quick_16_byte_hash(u_int8_t *in_16_bytes_long);
 
   /* ******************************* */
@@ -1451,7 +1458,7 @@ extern "C" {
   void ndpi_hash_free(ndpi_str_hash *h);
   int ndpi_hash_find_entry(ndpi_str_hash *h, char *key, u_int key_len, u_int8_t *value);
   int ndpi_hash_add_entry(ndpi_str_hash *h, char *key, u_int8_t key_len, u_int8_t value);
-  
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/include/ndpi_typedefs.h b/src/include/ndpi_typedefs.h
index 94556f6a3..7656cb603 100644
--- a/src/include/ndpi_typedefs.h
+++ b/src/include/ndpi_typedefs.h
@@ -1599,4 +1599,23 @@ typedef struct {
 
 /* **************************************** */
 
+#define HW_HISTORY_LEN   4
+
+struct ndpi_hw_struct {
+  struct {
+    u_int8_t use_hw_additive_seasonal;
+    double alpha, beta, gamma, ro;
+    u_int16_t num_season_periods; /* num of values of a season */
+  } params;
+
+  u_int32_t num_values;
+  double    u, v, sum_square_error;
+  
+  /* These two values need to store the signal history */
+  u_int32_t *y;
+  double    *s;
+};
+
+/* **************************************** */
+
 #endif /* __NDPI_TYPEDEFS_H__ */
diff --git a/src/lib/ndpi_analyze.c b/src/lib/ndpi_analyze.c
index 1ee6e99b2..38f74bc0b 100644
--- a/src/lib/ndpi_analyze.c
+++ b/src/lib/ndpi_analyze.c
@@ -105,7 +105,7 @@ void ndpi_data_add_value(struct ndpi_analyze_struct *s, const u_int32_t value) {
 
   /*
     Optimized stddev calculation
-    
+
     https://www.khanacademy.org/math/probability/data-distributions-a1/summarizing-spread-distributions/a/calculating-standard-deviation-step-by-step
     https://math.stackexchange.com/questions/683297/how-to-calculate-standard-deviation-without-detailed-historical-data
     http://mathcentral.uregina.ca/QQ/database/QQ.09.02/carlos1.html
@@ -140,7 +140,7 @@ u_int32_t ndpi_data_max(struct ndpi_analyze_struct *s) { return(s->max_val); }
 
 /* Compute the variance on all values */
 float ndpi_data_variance(struct ndpi_analyze_struct *s) {
-  return(s->num_data_entries ? ((float)s->stddev.sum_square_total - ((float)s->sum_total * (float)s->sum_total / (float)s->num_data_entries)) / (float)s->num_data_entries : 0.0);  
+  return(s->num_data_entries ? ((float)s->stddev.sum_square_total - ((float)s->sum_total * (float)s->sum_total / (float)s->num_data_entries)) / (float)s->num_data_entries : 0.0);
 }
 
 /* ********************************************************************************* */
@@ -468,7 +468,7 @@ void ndpi_normalize_bin(struct ndpi_bin *b) {
 	b->u.bins8[i] = (b->u.bins8[i]*100) / tot;
     }
     break;
-    
+
   case ndpi_bin_family16:
     for(i=0; i<b->num_bins; i++) tot += b->u.bins16[i];
 
@@ -588,7 +588,7 @@ float ndpi_bin_similarity(struct ndpi_bin *b1, struct ndpi_bin *b2, u_int8_t nor
       u_int32_t a = ndpi_get_bin_value(b1, i);
       u_int32_t b = ndpi_get_bin_value(b2, i);
       u_int32_t diff = (a > b) ? (a - b) : (b - a);
-      
+
       if(a != b) sum += pow(diff, 2);
 
       // printf("[a: %u][b: %u][sum: %u]\n", a, b, sum);
@@ -721,7 +721,7 @@ int ndpi_cluster_bins(struct ndpi_bin *bins, u_int16_t num_bins,
 
 	if(j == cluster_ids[i])
 	  current_similarity = similarity;
-	
+
 	if(verbose)
 	  printf("Bin %u / centroid %u [similarity: %f]\n", i, j, similarity);
 
@@ -743,7 +743,7 @@ int ndpi_cluster_bins(struct ndpi_bin *bins, u_int16_t num_bins,
 	*/
 	cluster_id = cluster_ids[i];
       }
-      
+
       bin_score[i] = best_similarity;
 
       if(cluster_ids[i] != cluster_id) {
@@ -818,7 +818,7 @@ int ndpi_cluster_bins(struct ndpi_bin *bins, u_int16_t num_bins,
 
 /* ********************************************************************************* */
 
-/* 
+/*
    RSI (Relative Strength Index)
 
    RSI = 100 − [ 100/ (1 + (Average gain/Average loss)) ]
@@ -856,12 +856,12 @@ void ndpi_free_rsi(struct ndpi_rsi_struct *s) {
 */
 float ndpi_rsi_add_value(struct ndpi_rsi_struct *s, const u_int32_t value) {
   float relative_strength;
-  
+
   if(!s->empty) {
     u_int32_t val;
 
     s->total_gains -= s->gains[s->next_index], s->total_losses -= s->losses[s->next_index];
-    
+
     if(value > s->last_value) {
       val = value - s->last_value;
       s->gains[s->next_index] = val, s->losses[s->next_index] = 0;
@@ -882,10 +882,10 @@ float ndpi_rsi_add_value(struct ndpi_rsi_struct *s, const u_int32_t value) {
     printf("[value: %u][total_gains: %u][total_losses: %u][cur_idx: %u]\n", value, s->total_gains, s->total_losses, s->next_index);
 #endif
   }
-  
+
   s->last_value = value, s->next_index = (s->next_index + 1) % s->num_values, s->empty = 0;
   if(s->next_index == 0) s->rsi_ready = 1; /* We have completed one round */
-  
+
   if(!s->rsi_ready)
     return(0);
   else if(s->total_losses == 0) /* Avoid division by zero (**) */
@@ -898,3 +898,159 @@ float ndpi_rsi_add_value(struct ndpi_rsi_struct *s, const u_int32_t value) {
     return(100. - (100. / (1. + relative_strength)));
   }
 }
+
+/* *********************************************************** */
+
+/* https://www.johndcook.com/blog/cpp_phi_inverse/ */
+
+static double ndpi_rational_approximation(double t) {
+  // Abramowitz and Stegun formula 26.2.23.
+  // The absolute value of the error should be less than 4.5 e-4.
+  double c[] = { 2.515517, 0.802853, 0.010328 };
+  double d[] = { 1.432788, 0.189269, 0.001308 };
+
+  return(t - ((c[2]*t + c[1])*t + c[0]) / (((d[2]*t + d[1])*t + d[0])*t + 1.0));
+}
+
+static double ndpi_normal_cdf_inverse(double p) {
+  if(p <= 0.0 || p >= 1.0)
+    return(0); /* Invalid argument: valid range 0 < X < 1 */
+
+  if(p < 0.5) {
+    // F^-1(p) = - G^-1(p)
+    return -ndpi_rational_approximation( sqrt(-2.0*log(p)) );
+  } else {
+    // F^-1(p) = G^-1(1-p)
+    return ndpi_rational_approximation( sqrt(-2.0*log(1-p)) );
+  }
+}
+
+double ndpi_avg_inline(u_int32_t *v, u_int num) {
+  double avg = 0;
+  u_int i;
+
+  for(i=0; i<num; i++)
+    avg += v[i];
+
+  return(avg / (u_int32_t)num);
+}
+
+/* *********************************************************** */
+
+/*
+  Initializes Holt-Winters with Confidence Interval
+
+   Input
+   hw:          Datastructure to initialize and that needs tobe freed with ndpi_hw_free()
+   num_periods  Number of observations of a season, or in ML-parlance the number of points that are required to make the forecast
+   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.
+
+   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
+
+   NOTE (See https://otexts.com/fpp2/holt-winters.html)
+   The additive method is preferred when the seasonal variations are roughly constant through the series,
+   while the multiplicative method is preferred when the seasonal variations are changing proportional to the level of the series.
+
+   For learning more about timeseries forecasting see
+   https://www.real-statistics.com/time-series-analysis/basic-time-series-forecasting/
+ */
+
+int ndpi_hw_init(struct ndpi_hw_struct *hw,
+		 u_int16_t num_periods, u_int8_t additive_seeasonal,
+		 double alpha, double beta, double gamma, float significance) {
+  memset(hw, 0, sizeof(struct ndpi_hw_struct));
+
+  hw->params.num_season_periods = num_periods;
+  hw->params.alpha      = alpha;
+  hw->params.beta       = beta;
+  hw->params.gamma      = gamma;
+  hw->params.use_hw_additive_seasonal = additive_seeasonal;
+
+  if((significance < 0) || (significance > 1)) significance = 0.05;
+  hw->params.ro         = ndpi_normal_cdf_inverse(1 - (significance / 2.));
+
+  if((hw->y = (u_int32_t*)ndpi_calloc(hw->params.num_season_periods, sizeof(u_int32_t))) == NULL)
+    return(-1);
+
+  if((hw->s = (double*)ndpi_calloc(hw->params.num_season_periods, sizeof(double))) == NULL) {
+    free(hw->y);
+    return(-1);
+  }
+
+  return(0);
+}
+
+/* *********************************************************** */
+
+/* Frees the memory allocated by ndpi_hw_init() */
+void ndpi_hw_free(struct ndpi_hw_struct *hw) {
+  if(hw->y) ndpi_free(hw->y);
+  if(hw->s) ndpi_free(hw->s);
+}
+
+/* *********************************************************** */
+
+/*
+   Returns the forecast and the band (forecast +/- band are the upper and lower values)
+
+   Input
+   hw:          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_hw_add_value(struct ndpi_hw_struct *hw, const u_int32_t _value, double *forecast,  double *confidence_band) {
+  if(hw->num_values < hw->params.num_season_periods) {
+    hw->y[hw->num_values++] = _value;
+
+    if(hw->num_values == hw->params.num_season_periods) {
+      double avg = ndpi_avg_inline(hw->y, hw->params.num_season_periods);
+      u_int i;
+
+      for(i=0; i<hw->params.num_season_periods; i++)
+	hw->s[i] = hw->y[i] / avg;
+
+      hw->u = _value / hw->s[hw->params.num_season_periods-1];
+      hw->v = 0;
+    }
+
+    *forecast = 0;
+    *confidence_band = 0;
+
+    return(0); /* Too early still training... */
+  } else {
+    u_int idx     = hw->num_values % hw->params.num_season_periods;
+    double prev_u = hw->u;
+    double prev_v = hw->v;
+    double prev_s = hw->s[idx];
+    double value  = (double)_value;;
+    double error;
+
+    hw->u      = ((hw->params.alpha * value) / prev_s) + (1 - hw->params.alpha) * (hw->u +  hw->v);
+    hw->v      = (hw->params.beta * (hw->u - prev_u))  + ((1 - hw->params.beta) * hw->v);
+    hw->s[idx] = (hw->params.gamma * (value / hw->u))  + ((1 - hw->params.gamma) * prev_s);
+
+    hw->num_values++, idx = (idx + 1) % hw->params.num_season_periods;
+
+    if(hw->params.use_hw_additive_seasonal)
+      *forecast = (prev_u + prev_v) + prev_s;
+    else
+      *forecast = (prev_u + prev_v) * prev_s;
+
+    error                 = value - *forecast;
+    hw->sum_square_error += error * error;
+    *confidence_band      = hw->params.ro * (sqrt(hw->sum_square_error / (hw->num_values - hw->params.num_season_periods)));
+
+    return(1); /* We're in business: forecast is meaningful now */
+  }
+}