Adding upstream version 4.0.upstream/4.0upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 704 insertions, 0 deletions

diff --git a/regress.c b/regress.c
new file mode 100644
index 0000000..e767e2f
--- /dev/null
+++ b/regress.c
@@ -0,0 +1,704 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Richard P. Curnow  1997-2003
+ * Copyright (C) Miroslav Lichvar  2011, 2016-2017
+ * 
+ * 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.
+ * 
+ **********************************************************************
+
+  =======================================================================
+
+  Regression algorithms.
+
+  */
+
+#include "config.h"
+
+#include "sysincl.h"
+
+#include "regress.h"
+#include "logging.h"
+#include "util.h"
+
+#define MAX_POINTS 64
+
+void
+RGR_WeightedRegression
+(double *x,                     /* independent variable */
+ double *y,                     /* measured data */
+ double *w,                     /* weightings (large => data
+                                   less reliable) */
+ 
+ int n,                         /* number of data points */
+
+ /* And now the results */
+
+ double *b0,                    /* estimated y axis intercept */
+ double *b1,                    /* estimated slope */
+ double *s2,                    /* estimated variance of data points */
+ 
+ double *sb0,                   /* estimated standard deviation of
+                                   intercept */
+ double *sb1                    /* estimated standard deviation of
+                                   slope */
+
+ /* Could add correlation stuff later if required */
+)
+{
+  double P, Q, U, V, W;
+  double diff;
+  double u, ui, aa;
+  int i;
+
+  assert(n >= 3);
+
+  W = U = 0;
+  for (i=0; i<n; i++) {
+    U += x[i]        / w[i];
+    W += 1.0         / w[i];
+  }
+
+  u = U / W;
+
+  /* Calculate statistics from data */
+  P = Q = V = 0.0;
+  for (i=0; i<n; i++) {
+    ui = x[i] - u;
+    P += y[i]        / w[i];
+    Q += y[i] * ui   / w[i];
+    V += ui   * ui   / w[i];
+  }
+
+  *b1 = Q / V;
+  *b0 = (P / W) - (*b1) * u;
+
+  *s2 = 0.0;
+  for (i=0; i<n; i++) {
+    diff = y[i] - *b0 - *b1*x[i];
+    *s2 += diff*diff / w[i];
+  }
+
+  *s2 /= (double)(n-2);
+
+  *sb1 = sqrt(*s2 / V);
+  aa = u * (*sb1);
+  *sb0 = sqrt(*s2 / W + aa * aa);
+
+  *s2 *= (n / W); /* Giving weighted average of variances */
+}
+
+/* ================================================== */
+/* Get the coefficient to multiply the standard deviation by, to get a
+   particular size of confidence interval (assuming a t-distribution) */
+  
+double
+RGR_GetTCoef(int dof)
+{
+  /* Assuming now the 99.95% quantile */
+  static const float coefs[] =
+  { 636.6, 31.6, 12.92, 8.61, 6.869,
+    5.959, 5.408, 5.041, 4.781, 4.587,
+    4.437, 4.318, 4.221, 4.140, 4.073,
+    4.015, 3.965, 3.922, 3.883, 3.850,
+    3.819, 3.792, 3.768, 3.745, 3.725,
+    3.707, 3.690, 3.674, 3.659, 3.646,
+    3.633, 3.622, 3.611, 3.601, 3.591,
+    3.582, 3.574, 3.566, 3.558, 3.551};
+
+  if (dof <= 40) {
+    return coefs[dof-1];
+  } else {
+    return 3.5; /* Until I can be bothered to do something better */
+  }
+}
+
+/* ================================================== */
+/* Get 90% quantile of chi-square distribution */
+
+double
+RGR_GetChi2Coef(int dof)
+{
+  static const float coefs[] = {
+    2.706, 4.605, 6.251, 7.779, 9.236, 10.645, 12.017, 13.362,
+    14.684, 15.987, 17.275, 18.549, 19.812, 21.064, 22.307, 23.542,
+    24.769, 25.989, 27.204, 28.412, 29.615, 30.813, 32.007, 33.196,
+    34.382, 35.563, 36.741, 37.916, 39.087, 40.256, 41.422, 42.585,
+    43.745, 44.903, 46.059, 47.212, 48.363, 49.513, 50.660, 51.805,
+    52.949, 54.090, 55.230, 56.369, 57.505, 58.641, 59.774, 60.907,
+    62.038, 63.167, 64.295, 65.422, 66.548, 67.673, 68.796, 69.919,
+    71.040, 72.160, 73.279, 74.397, 75.514, 76.630, 77.745, 78.860
+  };
+
+  if (dof <= 64) {
+    return coefs[dof-1];
+  } else {
+    return 1.2 * dof; /* Until I can be bothered to do something better */
+  }
+}
+
+/* ================================================== */
+/* Critical value for number of runs of residuals with same sign.
+   5% critical region for now. */
+
+static char critical_runs[] = {
+  0,  0,  0,  0,  0,  0,  0,  0,  2,  3,
+  3,  3,  4,  4,  5,  5,  5,  6,  6,  7,
+  7,  7,  8,  8,  9,  9,  9, 10, 10, 11,
+ 11, 11, 12, 12, 13, 13, 14, 14, 14, 15,
+ 15, 16, 16, 17, 17, 18, 18, 18, 19, 19,
+ 20, 20, 21, 21, 21, 22, 22, 23, 23, 24,
+ 24, 25, 25, 26, 26, 26, 27, 27, 28, 28,
+ 29, 29, 30, 30, 30, 31, 31, 32, 32, 33,
+ 33, 34, 34, 35, 35, 35, 36, 36, 37, 37,
+ 38, 38, 39, 39, 40, 40, 40, 41, 41, 42,
+ 42, 43, 43, 44, 44, 45, 45, 46, 46, 46,
+ 47, 47, 48, 48, 49, 49, 50, 50, 51, 51,
+ 52, 52, 52, 53, 53, 54, 54, 55, 55, 56
+};
+
+/* ================================================== */
+
+static int
+n_runs_from_residuals(double *resid, int n)
+{
+  int nruns;
+  int i;
+  
+  nruns = 1;
+  for (i=1; i<n; i++) {
+    if (((resid[i-1] < 0.0) && (resid[i] < 0.0)) ||
+        ((resid[i-1] > 0.0) && (resid[i] > 0.0))) {
+      /* Nothing to do */
+    } else {
+      nruns++;
+    }
+  }
+  
+  return nruns;
+}
+
+/* ================================================== */
+/* Return a boolean indicating whether we had enough points for
+   regression */
+
+int
+RGR_FindBestRegression 
+(double *x,                     /* independent variable */
+ double *y,                     /* measured data */
+ double *w,                     /* weightings (large => data
+                                   less reliable) */
+ 
+ int n,                         /* number of data points */
+ int m,                         /* number of extra samples in x and y arrays
+                                   (negative index) which can be used to
+                                   extend runs test */
+ int min_samples,               /* minimum number of samples to be kept after
+                                   changing the starting index to pass the runs
+                                   test */
+
+ /* And now the results */
+
+ double *b0,                    /* estimated y axis intercept */
+ double *b1,                    /* estimated slope */
+ double *s2,                    /* estimated variance of data points */
+ 
+ double *sb0,                   /* estimated standard deviation of
+                                   intercept */
+ double *sb1,                   /* estimated standard deviation of
+                                   slope */
+
+ int *new_start,                /* the new starting index to make the
+                                   residuals pass the two tests */
+ 
+ int *n_runs,                   /* number of runs amongst the residuals */
+
+ int *dof                       /* degrees of freedom in statistics (needed
+                                   to get confidence intervals later) */
+
+)
+{
+  double P, Q, U, V, W; /* total */
+  double resid[MAX_POINTS * REGRESS_RUNS_RATIO];
+  double ss;
+  double a, b, u, ui, aa;
+
+  int start, resid_start, nruns, npoints;
+  int i;
+
+  assert(n <= MAX_POINTS && m >= 0);
+  assert(n * REGRESS_RUNS_RATIO < sizeof (critical_runs) / sizeof (critical_runs[0]));
+
+  if (n < MIN_SAMPLES_FOR_REGRESS) {
+    return 0;
+  }
+
+  start = 0;
+  do {
+
+    W = U = 0;
+    for (i=start; i<n; i++) {
+      U += x[i]        / w[i];
+      W += 1.0         / w[i];
+    }
+
+    u = U / W;
+
+    P = Q = V = 0.0;
+    for (i=start; i<n; i++) {
+      ui = x[i] - u;
+      P += y[i]        / w[i];
+      Q += y[i] * ui   / w[i];
+      V += ui   * ui   / w[i];
+    }
+
+    b = Q / V;
+    a = (P / W) - (b * u);
+
+    /* Get residuals also for the extra samples before start */
+    resid_start = n - (n - start) * REGRESS_RUNS_RATIO;
+    if (resid_start < -m)
+      resid_start = -m;
+
+    for (i=resid_start; i<n; i++) {
+      resid[i - resid_start] = y[i] - a - b*x[i];
+    }
+
+    /* Count number of runs */
+    nruns = n_runs_from_residuals(resid, n - resid_start); 
+
+    if (nruns > critical_runs[n - resid_start] ||
+        n - start <= MIN_SAMPLES_FOR_REGRESS ||
+        n - start <= min_samples) {
+      if (start != resid_start) {
+        /* Ignore extra samples in returned nruns */
+        nruns = n_runs_from_residuals(resid + (start - resid_start), n - start);
+      }
+      break;
+    } else {
+      /* Try dropping one sample at a time until the runs test passes. */
+      ++start;
+    }
+
+  } while (1);
+
+  /* Work out statistics from full dataset */
+  *b1 = b;
+  *b0 = a;
+
+  ss = 0.0;
+  for (i=start; i<n; i++) {
+    ss += resid[i - resid_start]*resid[i - resid_start] / w[i];
+  }
+
+  npoints = n - start;
+  ss /= (double)(npoints - 2);
+  *sb1 = sqrt(ss / V);
+  aa = u * (*sb1);
+  *sb0 = sqrt((ss / W) + (aa * aa));
+  *s2 = ss * (double) npoints / W;
+
+  *new_start = start;
+  *dof = npoints - 2;
+  *n_runs = nruns;
+
+  return 1;
+
+}
+
+/* ================================================== */
+
+#define EXCH(a,b) temp=(a); (a)=(b); (b)=temp
+
+/* ================================================== */
+/* Find the index'th biggest element in the array x of n elements.
+   flags is an array where a 1 indicates that the corresponding entry
+   in x is known to be sorted into its correct position and a 0
+   indicates that the corresponding entry is not sorted.  However, if
+   flags[m] = flags[n] = 1 with m<n, then x[m] must be <= x[n] and for
+   all i with m<i<n, x[m] <= x[i] <= x[n].  In practice, this means
+   flags[] has to be the result of a previous call to this routine
+   with the same array x, and is used to remember which parts of the
+   x[] array we have already sorted.
+
+   The approach used is a cut-down quicksort, where we only bother to
+   keep sorting the partition that contains the index we are after.
+   The approach comes from Numerical Recipes in C (ISBN
+   0-521-43108-5). */
+
+static double
+find_ordered_entry_with_flags(double *x, int n, int index, char *flags)
+{
+  int u, v, l, r;
+  double temp;
+  double piv;
+  int pivind;
+
+  assert(index >= 0);
+
+  /* If this bit of the array is already sorted, simple! */
+  if (flags[index]) {
+    return x[index];
+  }
+  
+  /* Find subrange to look at */
+  u = v = index;
+  while (u > 0 && !flags[u]) u--;
+  if (flags[u]) u++;
+
+  while (v < (n-1) && !flags[v]) v++;
+  if (flags[v]) v--;
+
+  do {
+    if (v - u < 2) {
+      if (x[v] < x[u]) {
+        EXCH(x[v], x[u]);
+      }
+      flags[v] = flags[u] = 1;
+      return x[index];
+    } else { 
+      pivind = (u + v) >> 1;
+      EXCH(x[u], x[pivind]);
+      piv = x[u]; /* New value */
+      l = u + 1;
+      r = v;
+      do {
+        while (l < v && x[l] < piv) l++;
+        while (x[r] > piv) r--;
+        if (r <= l) break;
+        EXCH(x[l], x[r]);
+        l++;
+        r--;
+      } while (1);
+      EXCH(x[u], x[r]);
+      flags[r] = 1; /* Pivot now in correct place */
+      if (index == r) {
+        return x[r];
+      } else if (index < r) {
+        v = r - 1;
+      } else if (index > r) {
+        u = l;
+      }
+    }
+  } while (1);
+}
+
+/* ================================================== */
+
+#if 0
+/* Not used, but this is how it can be done */
+static double
+find_ordered_entry(double *x, int n, int index)
+{
+  char flags[MAX_POINTS];
+
+  memset(flags, 0, n * sizeof (flags[0]));
+  return find_ordered_entry_with_flags(x, n, index, flags);
+}
+#endif
+
+/* ================================================== */
+/* Find the median entry of an array x[] with n elements. */
+
+static double
+find_median(double *x, int n)
+{
+  int k;
+  char flags[MAX_POINTS];
+
+  memset(flags, 0, n * sizeof (flags[0]));
+  k = n>>1;
+  if (n&1) {
+    return find_ordered_entry_with_flags(x, n, k, flags);
+  } else {
+    return 0.5 * (find_ordered_entry_with_flags(x, n, k, flags) +
+                  find_ordered_entry_with_flags(x, n, k-1, flags));
+  }
+}
+
+/* ================================================== */
+
+double
+RGR_FindMedian(double *x, int n)
+{
+  double tmp[MAX_POINTS];
+
+  assert(n > 0 && n <= MAX_POINTS);
+  memcpy(tmp, x, n * sizeof (tmp[0]));
+
+  return find_median(tmp, n);
+}
+
+/* ================================================== */
+/* This function evaluates the equation
+
+   \sum_{i=0}^{n-1} x_i sign(y_i - a - b x_i)
+
+   and chooses the value of a that minimises the absolute value of the
+   result.  (See pp703-704 of Numerical Recipes in C). */
+
+static void
+eval_robust_residual
+(double *x,                     /* The independent points */
+ double *y,                     /* The dependent points */
+ int n,                         /* Number of points */
+ double b,                      /* Slope */
+ double *aa,                    /* Intercept giving smallest absolute
+                                   value for the above equation */
+ double *rr                     /* Corresponding value of equation */
+)
+{
+  int i;
+  double a, res, del;
+  double d[MAX_POINTS];
+
+  for (i=0; i<n; i++) {
+    d[i] = y[i] - b * x[i];
+  }
+  
+  a = find_median(d, n);
+
+  res = 0.0;
+  for (i=0; i<n; i++) {
+    del = y[i] - a - b * x[i];
+    if (del > 0.0) {
+      res += x[i];
+    } else if (del < 0.0) {
+      res -= x[i];
+    }
+  }
+
+  *aa = a;
+  *rr = res;
+}
+
+/* ================================================== */
+/* This routine performs a 'robust' regression, i.e. one which has low
+   susceptibility to outliers amongst the data.  If one thinks of a
+   normal (least squares) linear regression in 2D being analogous to
+   the arithmetic mean in 1D, this algorithm in 2D is roughly
+   analogous to the median in 1D.  This algorithm seems to work quite
+   well until the number of outliers is approximately half the number
+   of data points.
+
+   The return value is a status indicating whether there were enough
+   data points to run the routine or not. */
+
+int
+RGR_FindBestRobustRegression
+(double *x,                     /* The independent axis points */
+ double *y,                     /* The dependent axis points (which
+                                   may contain outliers). */
+ int n,                         /* The number of points */
+ double tol,                    /* The tolerance required in
+                                   determining the value of b1 */
+ double *b0,                    /* The estimated Y-axis intercept */
+ double *b1,                    /* The estimated slope */
+ int *n_runs,                   /* The number of runs of residuals */
+ int *best_start                /* The best starting index */
+)
+{
+  int i;
+  int start;
+  int n_points;
+  double a, b;
+  double P, U, V, W, X;
+  double resid, resids[MAX_POINTS];
+  double blo, bhi, bmid, rlo, rhi, rmid;
+  double s2, sb, incr;
+  double mx, dx, my, dy;
+  int nruns = 0;
+
+  assert(n <= MAX_POINTS);
+
+  if (n < 2) {
+    return 0;
+  } else if (n == 2) {
+    /* Just a straight line fit (we need this for the manual mode) */
+    *b1 = (y[1] - y[0]) / (x[1] - x[0]);
+    *b0 = y[0] - (*b1) * x[0];
+    *n_runs = 0;
+    *best_start = 0;
+    return 1;
+  }
+
+  /* else at least 3 points, apply normal algorithm */
+
+  start = 0;
+
+  /* Loop to strip oldest points that cause the regression residuals
+     to fail the number of runs test */
+  do {
+
+    n_points = n - start;
+
+    /* Use standard least squares regression to get starting estimate */
+
+    P = U = 0.0;
+    for (i=start; i<n; i++) {
+      P += y[i];
+      U += x[i];
+    }
+
+    W = (double) n_points;
+
+    my = P/W;
+    mx = U/W;
+
+    X = V = 0.0;
+    for (i=start; i<n; i++) {
+      dy = y[i] - my;
+      dx = x[i] - mx;
+      X += dy * dx;
+      V += dx * dx;
+    }
+
+    b = X / V;
+    a = my - b*mx;
+
+    s2 = 0.0;
+    for (i=start; i<n; i++) {
+      resid = y[i] - a - b * x[i];
+      s2 += resid * resid;
+    }
+
+    /* Need to expand range of b to get a root in the interval.
+       Estimate standard deviation of b and expand range about b based
+       on that. */
+    sb = sqrt(s2 * W/V);
+    incr = MAX(sb, tol);
+  
+    do {
+      incr *= 2.0;
+
+      /* Give up if the interval is too large */
+      if (incr > 100.0)
+        return 0;
+
+      blo = b - incr;
+      bhi = b + incr;
+
+      /* We don't want 'a' yet */
+      eval_robust_residual(x + start, y + start, n_points, blo, &a, &rlo);
+      eval_robust_residual(x + start, y + start, n_points, bhi, &a, &rhi);
+
+    } while (rlo * rhi >= 0.0); /* fn vals have same sign or one is zero,
+                                   i.e. root not in interval (rlo, rhi). */
+
+    /* OK, so the root for b lies in (blo, bhi). Start bisecting */
+    do {
+      bmid = 0.5 * (blo + bhi);
+      if (!(blo < bmid && bmid < bhi))
+        break;
+      eval_robust_residual(x + start, y + start, n_points, bmid, &a, &rmid);
+      if (rmid == 0.0) {
+        break;
+      } else if (rmid * rlo > 0.0) {
+        blo = bmid;
+        rlo = rmid;
+      } else if (rmid * rhi > 0.0) {
+        bhi = bmid;
+        rhi = rmid;
+      } else {
+        assert(0);
+      }
+    } while (bhi - blo > tol);
+
+    *b0 = a;
+    *b1 = bmid;
+
+    /* Number of runs test, but not if we're already down to the
+       minimum number of points */
+    if (n_points == MIN_SAMPLES_FOR_REGRESS) {
+      break;
+    }
+
+    for (i=start; i<n; i++) {
+      resids[i] = y[i] - a - bmid * x[i];
+    }
+
+    nruns = n_runs_from_residuals(resids + start, n_points);
+
+    if (nruns > critical_runs[n_points]) {
+      break;
+    } else {
+      start++;
+    }
+
+  } while (1);
+
+  *n_runs = nruns;
+  *best_start = start;
+
+  return 1;
+
+}
+
+/* ================================================== */
+/* This routine performs linear regression with two independent variables.
+   It returns non-zero status if there were enough data points and there
+   was a solution. */
+
+int
+RGR_MultipleRegress
+(double *x1,                    /* first independent variable */
+ double *x2,                    /* second independent variable */
+ double *y,                     /* measured data */
+
+ int n,                         /* number of data points */
+
+ /* The results */
+ double *b2                     /* estimated second slope */
+                                /* other values are not needed yet */
+)
+{
+  double Sx1, Sx2, Sx1x1, Sx1x2, Sx2x2, Sx1y, Sx2y, Sy;
+  double U, V, V1, V2, V3;
+  int i;
+
+  if (n < 4)
+    return 0;
+
+  Sx1 = Sx2 = Sx1x1 = Sx1x2 = Sx2x2 = Sx1y = Sx2y = Sy = 0.0;
+
+  for (i = 0; i < n; i++) {
+    Sx1 += x1[i];
+    Sx2 += x2[i];
+    Sx1x1 += x1[i] * x1[i];
+    Sx1x2 += x1[i] * x2[i];
+    Sx2x2 += x2[i] * x2[i];
+    Sx1y += x1[i] * y[i];
+    Sx2y += x2[i] * y[i];
+    Sy += y[i];
+  }
+
+  U = n * (Sx1x2 * Sx1y - Sx1x1 * Sx2y) +
+      Sx1 * Sx1 * Sx2y - Sx1 * Sx2 * Sx1y +
+      Sy * (Sx2 * Sx1x1 - Sx1 * Sx1x2);
+
+  V1 = n * (Sx1x2 * Sx1x2 - Sx1x1 * Sx2x2);
+  V2 = Sx1 * Sx1 * Sx2x2 + Sx2 * Sx2 * Sx1x1;
+  V3 = -2.0 * Sx1 * Sx2 * Sx1x2;
+  V = V1 + V2 + V3;
+
+  /* Check if there is a (numerically stable) solution */
+  if (fabs(V) * 1.0e10 <= -V1 + V2 + fabs(V3))
+    return 0;
+
+  *b2 = U / V;
+
+  return 1;
+}