gst/gstclock-linreg.c

   1 /* GStreamer
   2  * Copyright (C) 1999,2000 Erik Walthinsen <omega@cse.ogi.edu>
   3  *                    2000 Wim Taymans <wtay@chello.be>
   4  *                    2004 Wim Taymans <wim@fluendo.com>
   5  *                    2015 Jan Schmidt <jan@centricular.com>
   6  *
   7  * gstclock-linreg.c: Linear regression implementation, used in clock slaving
   8  *
   9  * This library is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Library General Public
  11  * License as published by the Free Software Foundation; either
  12  * version 2 of the License, or (at your option) any later version.
  13  *
  14  * This library is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Library General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Library General Public
  20  * License along with this library; if not, write to the
  21  * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
  22  * Boston, MA 02110-1301, USA.
  23  */
  24
  25 #include "gst_private.h"
  26 #include <time.h>
  27
  28 #include "gstclock.h"
  29 #include "gstinfo.h"
  30 #include "gstutils.h"
  31 #include "glib-compat-private.h"
  32
  33 /* Compute log2 of the passed 64-bit number by finding the highest set bit */
  34 static guint
  35 gst_log2 (GstClockTime in)
  36 {
  37   const guint64 b[] =
  38       { 0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000, 0xFFFFFFFF00000000LL };
  39   const guint64 S[] = { 1, 2, 4, 8, 16, 32 };
  40   int i;
  41
  42   guint count = 0;
  43   for (i = 5; i >= 0; i--) {
  44     if (in & b[i]) {
  45       in >>= S[i];
  46       count |= S[i];
  47     }
  48   }
  49
  50   return count;
  51 }
  52
  53 /* http://mathworld.wolfram.com/LeastSquaresFitting.html
  54  * with SLAVE_LOCK
  55  */
  56 gboolean
  57 _priv_gst_do_linear_regression (GstClockTime * times, guint n,
  58     GstClockTime * m_num, GstClockTime * m_denom, GstClockTime * b,
  59     GstClockTime * xbase, gdouble * r_squared)
  60 {
  61   GstClockTime *newx, *newy;
  62   GstClockTime xmin, ymin, xbar, ybar, xbar4, ybar4;
  63   GstClockTime xmax, ymax;
  64   GstClockTimeDiff sxx, sxy, syy;
  65   GstClockTime *x, *y;
  66   gint i, j;
  67   gint pshift = 0;
  68   gint max_bits;
  69
  70   xbar = ybar = sxx = syy = sxy = 0;
  71
  72   x = times;
  73   y = times + 2;
  74
  75   xmin = ymin = G_MAXUINT64;
  76   xmax = ymax = 0;
  77   for (i = j = 0; i < n; i++, j += 4) {
  78     xmin = MIN (xmin, x[j]);
  79     ymin = MIN (ymin, y[j]);
  80
  81     xmax = MAX (xmax, x[j]);
  82     ymax = MAX (ymax, y[j]);
  83   }
  84
  85   newx = times + 1;
  86   newy = times + 3;
  87
  88   /* strip off unnecessary bits of precision */
  89   for (i = j = 0; i < n; i++, j += 4) {
  90     newx[j] = x[j] - xmin;
  91     newy[j] = y[j] - ymin;
  92   }
  93
  94 #ifdef DEBUGGING_ENABLED
  95   GST_CAT_DEBUG (GST_CAT_CLOCK, "reduced numbers:");
  96   for (i = j = 0; i < n; i++, j += 4)
  97     GST_CAT_DEBUG (GST_CAT_CLOCK,
  98         "  %" G_GUINT64_FORMAT "  %" G_GUINT64_FORMAT, newx[j], newy[j]);
  99 #endif
 100
 101   /* have to do this precisely otherwise the results are pretty much useless.
 102    * should guarantee that none of these accumulators can overflow */
 103
 104   /* quantities on the order of 1e10 to 1e13 -> 30-35 bits;
 105    * window size a max of 2^10, so
 106    this addition could end up around 2^45 or so -- ample headroom */
 107   for (i = j = 0; i < n; i++, j += 4) {
 108     /* Just in case assumptions about headroom prove false, let's check */
 109     if ((newx[j] > 0 && G_MAXUINT64 - xbar <= newx[j]) ||
 110         (newy[j] > 0 && G_MAXUINT64 - ybar <= newy[j])) {
 111       GST_CAT_WARNING (GST_CAT_CLOCK,
 112           "Regression overflowed in clock slaving! xbar %"
 113           G_GUINT64_FORMAT " newx[j] %" G_GUINT64_FORMAT " ybar %"
 114           G_GUINT64_FORMAT " newy[j] %" G_GUINT64_FORMAT, xbar, newx[j], ybar,
 115           newy[j]);
 116       return FALSE;
 117     }
 118
 119     xbar += newx[j];
 120     ybar += newy[j];
 121   }
 122   xbar /= n;
 123   ybar /= n;
 124
 125   /* multiplying directly would give quantities on the order of 1e20-1e26 ->
 126    * 60 bits to 70 bits times the window size that's 80 which is too much.
 127    * Instead we (1) subtract off the xbar*ybar in the loop instead of after,
 128    * to avoid accumulation; (2) shift off some estimated number of bits from
 129    * each multiplicand to limit the expected ceiling. For strange
 130    * distributions of input values, things can still overflow, in which
 131    * case we drop precision and retry - at most a few times, in practice rarely
 132    */
 133
 134   /* Guess how many bits we might need for the usual distribution of input,
 135    * with a fallback loop that drops precision if things go pear-shaped */
 136   max_bits = gst_log2 (MAX (xmax - xmin, ymax - ymin)) * 7 / 8 + gst_log2 (n);
 137   if (max_bits > 64)
 138     pshift = max_bits - 64;
 139
 140   i = 0;
 141   do {
 142 #ifdef DEBUGGING_ENABLED
 143     GST_CAT_DEBUG (GST_CAT_CLOCK,
 144         "Restarting regression with precision shift %u", pshift);
 145 #endif
 146
 147     xbar4 = xbar >> pshift;
 148     ybar4 = ybar >> pshift;
 149     sxx = syy = sxy = 0;
 150     for (i = j = 0; i < n; i++, j += 4) {
 151       GstClockTime newx4, newy4;
 152       GstClockTimeDiff tmp;
 153
 154       newx4 = newx[j] >> pshift;
 155       newy4 = newy[j] >> pshift;
 156
 157       tmp = (newx4 + xbar4) * (newx4 - xbar4);
 158       if (G_UNLIKELY (tmp > 0 && sxx > 0 && (G_MAXINT64 - sxx <= tmp))) {
 159         do {
 160           /* Drop some precision and restart */
 161           pshift++;
 162           sxx /= 4;
 163           tmp /= 4;
 164         } while (G_MAXINT64 - sxx <= tmp);
 165         break;
 166       } else if (G_UNLIKELY (tmp < 0 && sxx < 0 && (G_MAXINT64 - sxx >= tmp))) {
 167         do {
 168           /* Drop some precision and restart */
 169           pshift++;
 170           sxx /= 4;
 171           tmp /= 4;
 172         } while (G_MININT64 - sxx >= tmp);
 173         break;
 174       }
 175       sxx += tmp;
 176
 177       tmp = newy4 * newy4 - ybar4 * ybar4;
 178       if (G_UNLIKELY (tmp > 0 && syy > 0 && (G_MAXINT64 - syy <= tmp))) {
 179         do {
 180           pshift++;
 181           syy /= 4;
 182           tmp /= 4;
 183         } while (G_MAXINT64 - syy <= tmp);
 184         break;
 185       } else if (G_UNLIKELY (tmp < 0 && syy < 0 && (G_MAXINT64 - syy >= tmp))) {
 186         do {
 187           pshift++;
 188           syy /= 4;
 189           tmp /= 4;
 190         } while (G_MININT64 - syy >= tmp);
 191         break;
 192       }
 193       syy += tmp;
 194
 195       tmp = newx4 * newy4 - xbar4 * ybar4;
 196       if (G_UNLIKELY (tmp > 0 && sxy > 0 && (G_MAXINT64 - sxy <= tmp))) {
 197         do {
 198           pshift++;
 199           sxy /= 4;
 200           tmp /= 4;
 201         } while (G_MAXINT64 - sxy <= tmp);
 202         break;
 203       } else if (G_UNLIKELY (tmp < 0 && sxy < 0 && (G_MININT64 - sxy >= tmp))) {
 204         do {
 205           pshift++;
 206           sxy /= 4;
 207           tmp /= 4;
 208         } while (G_MININT64 - sxy >= tmp);
 209         break;
 210       }
 211       sxy += tmp;
 212     }
 213   } while (i < n);
 214
 215   if (G_UNLIKELY (sxx == 0))
 216     goto invalid;
 217
 218   *m_num = sxy;
 219   *m_denom = sxx;
 220   *b = (ymin + ybar) - gst_util_uint64_scale (xbar, *m_num, *m_denom);
 221   /* Report base starting from the most recent observation */
 222   *xbase = xmax;
 223   *b += gst_util_uint64_scale (xmax - xmin, *m_num, *m_denom);
 224
 225   *r_squared = ((double) sxy * (double) sxy) / ((double) sxx * (double) syy);
 226
 227 #ifdef DEBUGGING_ENABLED
 228   GST_CAT_DEBUG (GST_CAT_CLOCK, "  m      = %g", ((double) *m_num) / *m_denom);
 229   GST_CAT_DEBUG (GST_CAT_CLOCK, "  b      = %" G_GUINT64_FORMAT, *b);
 230   GST_CAT_DEBUG (GST_CAT_CLOCK, "  xbase  = %" G_GUINT64_FORMAT, *xbase);
 231   GST_CAT_DEBUG (GST_CAT_CLOCK, "  r2     = %g", *r_squared);
 232 #endif
 233
 234   return TRUE;
 235
 236 invalid:
 237   {
 238     GST_CAT_DEBUG (GST_CAT_CLOCK, "sxx == 0, regression failed");
 239     return FALSE;
 240   }
 241 }