quantiles: add parameter to limit negative step

Add a new parameter to limit the negative value of the step state
variable. It's set as a maximum delay in number of updates before the
actual step applied to the quantile estimate starts growing from the
minimum step when the input value is consistently larger or smaller than
the estimate.

This prevents the algorithm from effectively becoming the slower 1U
variant if the quantile estimate is stable most of the time.

Set it to 100 updates for the NTP delay and 1000 updates for the hwclock
delay. An option could be added later to make it configurable.

diff --git a/hwclock.c b/hwclock.c
index e2f72e7aeba99b76ec36945f05a2714cf73060a6..06515e7b1c2a35a3a35ac0dc16433e596d5c20e2 100644 (file)
--- a/hwclock.c
+++ b/hwclock.c
@@ -49,6 +49,7 @@
 #define DELAY_QUANT_MAX_K 2
 #define DELAY_QUANT_Q 10
 #define DELAY_QUANT_REPEAT 7
+#define DELAY_QUANT_LARGE_STEP_DELAY 1000
 #define DELAY_QUANT_MIN_STEP 1.0e-9
 
 struct HCL_Instance_Record {
@@ -127,6 +128,7 @@ HCL_CreateInstance(int min_samples, int max_samples, double min_separation, doub
   clock->precision = precision;
   clock->delay_quants = QNT_CreateInstance(DELAY_QUANT_MIN_K, DELAY_QUANT_MAX_K,
                                            DELAY_QUANT_Q, DELAY_QUANT_REPEAT,
+                                           DELAY_QUANT_LARGE_STEP_DELAY,
                                            DELAY_QUANT_MIN_STEP);
 
   LCL_AddParameterChangeHandler(handle_slew, clock);

diff --git a/ntp_core.c b/ntp_core.c
index e802e1832379a30c9fda350688d769123c7ac0ab..19b70a9b0c876b35edbffd03fb5dae22936d2c0a 100644 (file)
--- a/ntp_core.c
+++ b/ntp_core.c
@@ -286,6 +286,7 @@ static ARR_Instance broadcasts;
 
 /* Parameters for the peer delay quantile */
 #define DELAY_QUANT_Q 100
+#define DELAY_QUANT_LARGE_STEP_DELAY 100
 #define DELAY_QUANT_REPEAT 7
 
 /* Minimum and maximum allowed poll interval */
@@ -690,6 +691,7 @@ NCR_CreateInstance(NTP_Remote_Address *remote_addr, NTP_Source_Type type,
   if (params->max_delay_quant > 0.0) {
     int k = round(CLAMP(0.05, params->max_delay_quant, 0.95) * DELAY_QUANT_Q);
     result->delay_quant = QNT_CreateInstance(k, k, DELAY_QUANT_Q, DELAY_QUANT_REPEAT,
+                                             DELAY_QUANT_LARGE_STEP_DELAY,
                                              LCL_GetSysPrecisionAsQuantum() / 2.0);
   } else {
     result->delay_quant = NULL;

diff --git a/quantiles.c b/quantiles.c
index f4783d4d5ab9b2533821476fd704485ab071c05c..3f0800dd40e54cf360b9b83f6cd93bd551964bf1 100644 (file)
--- a/quantiles.c
+++ b/quantiles.c
@@ -49,19 +49,21 @@ struct QNT_Instance_Record {
   int q;
   int min_k;
   double min_step;
+  double neg_step_limit;
   int n_set;
 };
 
 /* ================================================== */
 
 QNT_Instance
-QNT_CreateInstance(int min_k, int max_k, int q, int repeat, double min_step)
+QNT_CreateInstance(int min_k, int max_k, int q, int repeat,
+                   int large_step_delay, double min_step)
 {
   QNT_Instance inst;
   long seed;
 
   if (q < 2 || min_k > max_k || min_k < 1 || max_k >= q ||
-      repeat < 1 || repeat > MAX_REPEAT || min_step <= 0.0)
+      repeat < 1 || repeat > MAX_REPEAT || min_step <= 0.0 || large_step_delay < 0)
     assert(0);
 
   inst = MallocNew(struct QNT_Instance_Record);
@@ -71,6 +73,7 @@ QNT_CreateInstance(int min_k, int max_k, int q, int repeat, double min_step)
   inst->q = q;
   inst->min_k = min_k;
   inst->min_step = min_step;
+  inst->neg_step_limit = -large_step_delay * min_step;
 
   QNT_Reset(inst);
 
@@ -136,7 +139,7 @@ insert_initial_value(QNT_Instance inst, double value)
 
 static void
 update_estimate(struct Quantile *quantile, double value, double p, double rand,
-                double min_step)
+                double min_step, double neg_step_limit)
 {
   if (value >= quantile->est) {
     if (rand < (1.0 - p))
@@ -163,6 +166,9 @@ update_estimate(struct Quantile *quantile, double value, double p, double rand,
       quantile->step = min_step;
     quantile->sign = -1;
   }
+
+  if (quantile->step < neg_step_limit)
+    quantile->step = neg_step_limit;
 }
 
 /* ================================================== */
@@ -183,7 +189,7 @@ QNT_Accumulate(QNT_Instance inst, double value)
     p = (double)(i / inst->repeat + inst->min_k) / inst->q;
     rand = (double)random() / ((1U << 31) - 1);
 
-    update_estimate(&inst->quants[i], value, p, rand, inst->min_step);
+    update_estimate(&inst->quants[i], value, p, rand, inst->min_step, inst->neg_step_limit);
   }
 }
 

diff --git a/quantiles.h b/quantiles.h
index b9665d508ff74df9e5b00a9e1249f0b96bef5dfa..bd335dedd44c869acf6ebb4da951f30f4999168e 100644 (file)
--- a/quantiles.h
+++ b/quantiles.h
@@ -30,7 +30,8 @@
 
 typedef struct QNT_Instance_Record *QNT_Instance;
 
-extern QNT_Instance QNT_CreateInstance(int min_k, int max_k, int q, int repeat, double min_step);
+extern QNT_Instance QNT_CreateInstance(int min_k, int max_k, int q, int repeat,
+                                       int large_step_delay, double min_step);
 extern void QNT_DestroyInstance(QNT_Instance inst);
 
 extern void QNT_Reset(QNT_Instance inst);

diff --git a/test/unit/quantiles.c b/test/unit/quantiles.c
index e08801e18aec9c82abf36468d6516bb3b4db2483..0462ec7d996616874154921b4a0ace592678d643 100644 (file)
--- a/test/unit/quantiles.c
+++ b/test/unit/quantiles.c
@@ -40,7 +40,7 @@ test_unit(void)
       max_k = random() % (q - 1) + 1;
     } while (min_k > max_k);
 
-    inst = QNT_CreateInstance(min_k, max_k, q, r, 1e-9);
+    inst = QNT_CreateInstance(min_k, max_k, q, r, 900, 1e-9);
 
     TEST_CHECK(min_k == QNT_GetMinK(inst));
     TEST_CHECK(max_k == QNT_GetMaxK(inst));
@@ -76,11 +76,11 @@ test_unit(void)
         TEST_CHECK(inst->quants[k].est > x - 0.4e-9);
         TEST_CHECK(inst->quants[k].est < x2 + 0.4e-9);
         TEST_CHECK(inst->quants[k].step < -15e-9);
-        TEST_CHECK(inst->quants[k].step > -1000e-9);
+        TEST_CHECK(inst->quants[k].step > -901e-9);
         if (min_k * 2 == q && k < inst->repeat) {
           if (x == x2) {
             TEST_CHECK(inst->quants[k].step < -750e-9);
-            TEST_CHECK(inst->quants[k].step > -1000e-9);
+            TEST_CHECK(inst->quants[k].step > -901e-9);
           } else {
             TEST_CHECK(inst->quants[k].step < -350e-9);
             TEST_CHECK(inst->quants[k].step > -600e-9);