quantiles: add support for quantile estimation

Add estimation of quantiles using the Frugal-2U streaming algorithm
(https://arxiv.org/pdf/1407.1121v1.pdf). It does not need to save
previous samples and adapts to changes in the distribution.

Allow multiple estimates of the same quantile and select the median for
better stability.

diff --git a/quantiles.c b/quantiles.c
new file mode 100644 (file)
index 0000000..8b72a61
--- /dev/null
+++ b/quantiles.c
@@ -0,0 +1,201 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2022
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ **********************************************************************
+
+  =======================================================================
+
+  Estimation of quantiles using the Frugal-2U streaming algorithm
+  (https://arxiv.org/pdf/1407.1121v1.pdf)
+  */
+
+#include "config.h"
+
+#include "logging.h"
+#include "memory.h"
+#include "quantiles.h"
+#include "regress.h"
+#include "util.h"
+
+/* Maximum number of repeated estimates for stabilisation */
+#define MAX_REPEAT 64
+
+struct Quantile {
+  double est;
+  double step;
+  int sign;
+};
+
+struct QNT_Instance_Record {
+  struct Quantile *quants;
+  int n_quants;
+  int repeat;
+  int q;
+  int min_k;
+  double min_step;
+  int n_set;
+};
+
+/* ================================================== */
+
+QNT_Instance
+QNT_CreateInstance(int min_k, int max_k, int q, int repeat, 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)
+    assert(0);
+
+  inst = MallocNew(struct QNT_Instance_Record);
+  inst->n_quants = (max_k - min_k + 1) * repeat;
+  inst->quants = MallocArray(struct Quantile, inst->n_quants);
+  inst->repeat = repeat;
+  inst->q = q;
+  inst->min_k = min_k;
+  inst->min_step = min_step;
+
+  QNT_Reset(inst);
+
+  /* Seed the random number generator, which will not be isolated from
+     other instances and other random() users */
+  UTI_GetRandomBytes(&seed, sizeof (seed));
+  srandom(seed);
+
+  return inst;
+}
+
+/* ================================================== */
+
+void
+QNT_DestroyInstance(QNT_Instance inst)
+{
+  Free(inst->quants);
+  Free(inst);
+}
+
+/* ================================================== */
+
+void
+QNT_Reset(QNT_Instance inst)
+{
+  int i;
+
+  inst->n_set = 0;
+
+  for (i = 0; i < inst->n_quants; i++) {
+    inst->quants[i].est = 0.0;
+    inst->quants[i].step = inst->min_step;
+    inst->quants[i].sign = 1;
+  }
+}
+
+/* ================================================== */
+
+static void
+insert_initial_value(QNT_Instance inst, double value)
+{
+  int i, j, r = inst->repeat;
+
+  if (inst->n_set * r >= inst->n_quants)
+    assert(0);
+
+  /* Keep the initial estimates repeated and ordered */
+  for (i = inst->n_set; i > 0 && inst->quants[(i - 1) * r].est > value; i--) {
+    for (j = 0; j < r; j++)
+      inst->quants[i * r + j].est = inst->quants[(i - 1) * r].est;
+  }
+
+  for (j = 0; j < r; j++)
+    inst->quants[i * r + j].est = value;
+  inst->n_set++;
+
+  /* Duplicate the largest value in unset quantiles */
+  for (i = inst->n_set * r; i < inst->n_quants; i++)
+    inst->quants[i].est = inst->quants[i - 1].est;
+}
+
+/* ================================================== */
+
+static void
+update_estimate(struct Quantile *quantile, double value, double p, double rand,
+                double min_step)
+{
+  if (value > quantile->est && rand > (1.0 - p)) {
+    quantile->step += quantile->sign > 0 ? min_step : -min_step;
+    quantile->est += quantile->step > 0.0 ? fabs(quantile->step) : min_step;
+    if (quantile->est > value) {
+      quantile->step += value - quantile->est;
+      quantile->est = value;
+    }
+    if (quantile->sign < 0 && quantile->step > min_step)
+      quantile->step = min_step;
+    quantile->sign = 1;
+  } else if (value < quantile->est && rand > p) {
+    quantile->step += quantile->sign < 0 ? min_step : -min_step;
+    quantile->est -= quantile->step > 0.0 ? fabs(quantile->step) : min_step;
+    if (quantile->est < value) {
+      quantile->step += quantile->est - value;
+      quantile->est = value;
+    }
+    if (quantile->sign > 0 && quantile->step > min_step)
+      quantile->step = min_step;
+    quantile->sign = -1;
+  }
+}
+
+/* ================================================== */
+
+void
+QNT_Accumulate(QNT_Instance inst, double value)
+{
+  double p, rand;
+  int i;
+
+  /* Initialise the estimates with first received values */
+  if (inst->n_set * inst->repeat < inst->n_quants) {
+    insert_initial_value(inst, value);
+    return;
+  }
+
+  for (i = 0; i < inst->n_quants; i++) {
+    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);
+  }
+}
+
+/* ================================================== */
+
+double
+QNT_GetQuantile(QNT_Instance inst, int k)
+{
+  double estimates[MAX_REPEAT];
+  int i;
+
+  if (k < inst->min_k || k - inst->min_k >= inst->n_quants)
+    assert(0);
+
+  for (i = 0; i < inst->repeat; i++)
+    estimates[i] = inst->quants[(k - inst->min_k) * inst->repeat + i].est;
+
+  return RGR_FindMedian(estimates, inst->repeat);
+}

diff --git a/quantiles.h b/quantiles.h
new file mode 100644 (file)
index 0000000..590105e
--- /dev/null
+++ b/quantiles.h
@@ -0,0 +1,40 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2022
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ **********************************************************************
+
+  =======================================================================
+
+  Header file for estimation of quantiles.
+
+  */
+
+#ifndef GOT_QUANTILES_H
+#define GOT_QUANTILES_H
+
+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 void QNT_DestroyInstance(QNT_Instance inst);
+
+extern void QNT_Reset(QNT_Instance inst);
+extern void QNT_Accumulate(QNT_Instance inst, double value);
+extern double QNT_GetQuantile(QNT_Instance inst, int k);
+
+#endif

diff --git a/test/unit/quantiles.c b/test/unit/quantiles.c
new file mode 100644 (file)
index 0000000..fc6364b
--- /dev/null
+++ b/test/unit/quantiles.c
@@ -0,0 +1,66 @@
+/*
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2022
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ **********************************************************************
+ */
+
+#include <local.h>
+#include "test.h"
+
+#include <quantiles.c>
+
+void
+test_unit(void)
+{
+  int i, j, k, min_k, max_k, q, r, in_order, out_order;
+  QNT_Instance inst;
+  double x;
+
+  in_order = out_order = 0;
+
+  for (i = 0; i < 100; i++) {
+    r = random() % 10 + 1;
+    q = random() % 20 + 2;
+    do {
+      min_k = random() % (q - 1) + 1;
+      max_k = random() % (q - 1) + 1;
+    } while (min_k > max_k);
+
+    inst = QNT_CreateInstance(min_k, max_k, q, r, 1e-9);
+
+    for (j = 0; j < 3000; j++) {
+      x = TST_GetRandomDouble(0.0, 2e-6);
+      QNT_Accumulate(inst, x);
+      for (k = min_k; k < max_k; k++)
+        if (j < max_k - min_k) {
+          TEST_CHECK(QNT_GetQuantile(inst, k) <= QNT_GetQuantile(inst, k + 1));
+        } else if (j > 1000) {
+          if (QNT_GetQuantile(inst, k) <= QNT_GetQuantile(inst, k + 1))
+            in_order++;
+          else
+            out_order++;
+        }
+    }
+
+    QNT_Reset(inst);
+    TEST_CHECK(inst->n_set == 0);
+
+    QNT_DestroyInstance(inst);
+  }
+
+  TEST_CHECK(in_order > 100 * out_order);
+}