/* GStreamer * Copyright (C) 1999,2000 Erik Walthinsen * 2000 Wim Taymans * 2004 Wim Taymans * 2015 Jan Schmidt * * gstclock-linreg.c: Linear regression implementation, used in clock slaving * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library 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 * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301, USA. */ #include "gst_private.h" #include #include "gstclock.h" #include "gstinfo.h" #include "gstutils.h" #include "glib-compat-private.h" /* Compute log2 of the passed 64-bit number by finding the highest set bit */ static guint gst_log2 (GstClockTime in) { const guint64 b[] = { 0x2, 0xC, 0xF0, 0xFF00, 0xFFFF0000, 0xFFFFFFFF00000000LL }; const guint64 S[] = { 1, 2, 4, 8, 16, 32 }; int i; guint count = 0; for (i = 5; i >= 0; i--) { if (in & b[i]) { in >>= S[i]; count |= S[i]; } } return count; } /* http://mathworld.wolfram.com/LeastSquaresFitting.html * with SLAVE_LOCK */ gboolean _priv_gst_do_linear_regression (GstClockTime * times, guint n, GstClockTime * m_num, GstClockTime * m_denom, GstClockTime * b, GstClockTime * xbase, gdouble * r_squared) { GstClockTime *newx, *newy; GstClockTime xmin, ymin, xbar, ybar, xbar4, ybar4; GstClockTime xmax, ymax; GstClockTimeDiff sxx, sxy, syy; GstClockTime *x, *y; gint i, j; gint pshift = 0; gint max_bits; xbar = ybar = sxx = syy = sxy = 0; x = times; y = times + 2; xmin = ymin = G_MAXUINT64; xmax = ymax = 0; for (i = j = 0; i < n; i++, j += 4) { xmin = MIN (xmin, x[j]); ymin = MIN (ymin, y[j]); xmax = MAX (xmax, x[j]); ymax = MAX (ymax, y[j]); } newx = times + 1; newy = times + 3; /* strip off unnecessary bits of precision */ for (i = j = 0; i < n; i++, j += 4) { newx[j] = x[j] - xmin; newy[j] = y[j] - ymin; } #ifdef DEBUGGING_ENABLED GST_CAT_DEBUG (GST_CAT_CLOCK, "reduced numbers:"); for (i = j = 0; i < n; i++, j += 4) GST_CAT_DEBUG (GST_CAT_CLOCK, " %" G_GUINT64_FORMAT " %" G_GUINT64_FORMAT, newx[j], newy[j]); #endif /* have to do this precisely otherwise the results are pretty much useless. * should guarantee that none of these accumulators can overflow */ /* quantities on the order of 1e10 to 1e13 -> 30-35 bits; * window size a max of 2^10, so this addition could end up around 2^45 or so -- ample headroom */ for (i = j = 0; i < n; i++, j += 4) { /* Just in case assumptions about headroom prove false, let's check */ if ((newx[j] > 0 && G_MAXUINT64 - xbar <= newx[j]) || (newy[j] > 0 && G_MAXUINT64 - ybar <= newy[j])) { GST_CAT_WARNING (GST_CAT_CLOCK, "Regression overflowed in clock slaving! xbar %" G_GUINT64_FORMAT " newx[j] %" G_GUINT64_FORMAT " ybar %" G_GUINT64_FORMAT " newy[j] %" G_GUINT64_FORMAT, xbar, newx[j], ybar, newy[j]); return FALSE; } xbar += newx[j]; ybar += newy[j]; } xbar /= n; ybar /= n; /* multiplying directly would give quantities on the order of 1e20-1e26 -> * 60 bits to 70 bits times the window size that's 80 which is too much. * Instead we (1) subtract off the xbar*ybar in the loop instead of after, * to avoid accumulation; (2) shift off some estimated number of bits from * each multiplicand to limit the expected ceiling. For strange * distributions of input values, things can still overflow, in which * case we drop precision and retry - at most a few times, in practice rarely */ /* Guess how many bits we might need for the usual distribution of input, * with a fallback loop that drops precision if things go pear-shaped */ max_bits = gst_log2 (MAX (xmax - xmin, ymax - ymin)) * 7 / 8 + gst_log2 (n); if (max_bits > 64) pshift = max_bits - 64; i = 0; do { #ifdef DEBUGGING_ENABLED GST_CAT_DEBUG (GST_CAT_CLOCK, "Restarting regression with precision shift %u", pshift); #endif xbar4 = xbar >> pshift; ybar4 = ybar >> pshift; sxx = syy = sxy = 0; for (i = j = 0; i < n; i++, j += 4) { GstClockTime newx4, newy4; GstClockTimeDiff tmp; newx4 = newx[j] >> pshift; newy4 = newy[j] >> pshift; tmp = (newx4 + xbar4) * (newx4 - xbar4); if (G_UNLIKELY (tmp > 0 && sxx > 0 && (G_MAXINT64 - sxx <= tmp))) { do { /* Drop some precision and restart */ pshift++; sxx /= 4; tmp /= 4; } while (G_MAXINT64 - sxx <= tmp); break; } else if (G_UNLIKELY (tmp < 0 && sxx < 0 && (G_MAXINT64 - sxx >= tmp))) { do { /* Drop some precision and restart */ pshift++; sxx /= 4; tmp /= 4; } while (G_MININT64 - sxx >= tmp); break; } sxx += tmp; tmp = newy4 * newy4 - ybar4 * ybar4; if (G_UNLIKELY (tmp > 0 && syy > 0 && (G_MAXINT64 - syy <= tmp))) { do { pshift++; syy /= 4; tmp /= 4; } while (G_MAXINT64 - syy <= tmp); break; } else if (G_UNLIKELY (tmp < 0 && syy < 0 && (G_MAXINT64 - syy >= tmp))) { do { pshift++; syy /= 4; tmp /= 4; } while (G_MININT64 - syy >= tmp); break; } syy += tmp; tmp = newx4 * newy4 - xbar4 * ybar4; if (G_UNLIKELY (tmp > 0 && sxy > 0 && (G_MAXINT64 - sxy <= tmp))) { do { pshift++; sxy /= 4; tmp /= 4; } while (G_MAXINT64 - sxy <= tmp); break; } else if (G_UNLIKELY (tmp < 0 && sxy < 0 && (G_MININT64 - sxy >= tmp))) { do { pshift++; sxy /= 4; tmp /= 4; } while (G_MININT64 - sxy >= tmp); break; } sxy += tmp; } } while (i < n); if (G_UNLIKELY (sxx == 0)) goto invalid; *m_num = sxy; *m_denom = sxx; *b = (ymin + ybar) - gst_util_uint64_scale (xbar, *m_num, *m_denom); /* Report base starting from the most recent observation */ *xbase = xmax; *b += gst_util_uint64_scale (xmax - xmin, *m_num, *m_denom); *r_squared = ((double) sxy * (double) sxy) / ((double) sxx * (double) syy); #ifdef DEBUGGING_ENABLED GST_CAT_DEBUG (GST_CAT_CLOCK, " m = %g", ((double) *m_num) / *m_denom); GST_CAT_DEBUG (GST_CAT_CLOCK, " b = %" G_GUINT64_FORMAT, *b); GST_CAT_DEBUG (GST_CAT_CLOCK, " xbase = %" G_GUINT64_FORMAT, *xbase); GST_CAT_DEBUG (GST_CAT_CLOCK, " r2 = %g", *r_squared); #endif return TRUE; invalid: { GST_CAT_DEBUG (GST_CAT_CLOCK, "sxx == 0, regression failed"); return FALSE; } }