using System; using System.Collections.Generic; using System.Linq; using UIWidgets.foundation; namespace UIWidgets.gestures { class _Vector { internal _Vector(int size) { this._offset = 0; this._length = size; this._elements = Enumerable.Repeat(0.0, size).ToList(); } private _Vector(List values, int offset, int length) { this._offset = offset; this._length = length; this._elements = values; } internal static _Vector fromVOL(List values, int offset, int length) { return new _Vector(values, offset, length); } readonly int _offset; readonly int _length; readonly List _elements; public double this[int i] { get { return this._elements[i + this._offset]; } set { this._elements[i + this._offset] = value; } } public static double operator *(_Vector a, _Vector b) { double result = 0.0; for (int i = 0; i < a._length; i += 1) result += a[i] * b[i]; return result; } public double norm() { return Math.Sqrt(this * this); } } class _Matrix { internal _Matrix(int rows, int cols) { this._columns = cols; this._elements = Enumerable.Repeat(0.0, rows * cols).ToList(); } readonly int _columns; readonly List _elements; public double this[int row, int col] { get { return this._elements[row * this._columns + col]; } set { this._elements[row * this._columns + col] = value; } } public _Vector getRow(int row) { return _Vector.fromVOL( this._elements, row * this._columns, this._columns ); } } public class PolynomialFit { public PolynomialFit(int degree) { this.coefficients = Enumerable.Repeat(0.0, degree + 1).ToList(); } public readonly List coefficients; public double confidence; } public class LeastSquaresSolver { public LeastSquaresSolver(List x, List y, List w) { D.assert(x != null && y != null && w != null); D.assert(x.Count == y.Count); D.assert(y.Count == w.Count); this.x = x; this.y = y; this.w = w; } public readonly List x; public readonly List y; public readonly List w; /// Fits a polynomial of the given degree to the data points. public PolynomialFit solve(int degree) { if (degree > this.x.Count) { // Not enough data to fit a curve. return null; } PolynomialFit result = new PolynomialFit(degree); // Shorthands for the purpose of notation equivalence to original C++ code. int m = x.Count; int n = degree + 1; // Expand the X vector to a matrix A, pre-multiplied by the weights. _Matrix a = new _Matrix(n, m); for (int h = 0; h < m; h += 1) { a[0, h] = this.w[h]; for (int i = 1; i < n; i += 1) a[i, h] = a[i - 1, h] * this.x[h]; } // Apply the Gram-Schmidt process to A to obtain its QR decomposition. // Orthonormal basis, column-major ordVectorer. _Matrix q = new _Matrix(n, m); // Upper triangular matrix, row-major order. _Matrix r = new _Matrix(n, n); for (int j = 0; j < n; j += 1) { for (int h = 0; h < m; h += 1) q[j, h] = a[j, h]; for (int i = 0; i < j; i += 1) { double dot = q.getRow(j) * q.getRow(i); for (int h = 0; h < m; h += 1) q[j, h] = q[j, h] - dot * q[i, h]; } double norm = q.getRow(j).norm(); if (norm < 0.000001) { // Vectors are linearly dependent or zero so no solution. return null; } double inverseNorm = 1.0 / norm; for (int h = 0; h < m; h += 1) q[j, h] = q[j, h] * inverseNorm; for (int i = 0; i < n; i += 1) r[j, i] = i < j ? 0.0 : q.getRow(j) * a.getRow(i); } // Solve R B = Qt W Y to find B. This is easy because R is upper triangular. // We just work from bottom-right to top-left calculating B's coefficients. _Vector wy = new _Vector(m); for (int h = 0; h < m; h += 1) wy[h] = y[h] * w[h]; for (int i = n - 1; i >= 0; i -= 1) { result.coefficients[i] = q.getRow(i) * wy; for (int j = n - 1; j > i; j -= 1) result.coefficients[i] -= r[i, j] * result.coefficients[j]; result.coefficients[i] /= r[i, i]; } // Calculate the coefficient of determination (confidence) as: // 1 - (sumSquaredError / sumSquaredTotal) // ...where sumSquaredError is the residual sum of squares (variance of the // error), and sumSquaredTotal is the total sum of squares (variance of the // data) where each has been weighted. double yMean = 0.0; for (int h = 0; h < m; h += 1) yMean += y[h]; yMean /= m; double sumSquaredError = 0.0; double sumSquaredTotal = 0.0; for (int h = 0; h < m; h += 1) { double term = 1.0; double err = y[h] - result.coefficients[0]; for (int i = 1; i < n; i += 1) { term *= x[h]; err -= term * result.coefficients[i]; } sumSquaredError += w[h] * w[h] * err * err; double v = y[h] - yMean; sumSquaredTotal += w[h] * w[h] * v * v; } result.confidence = sumSquaredTotal <= 0.000001 ? 1.0 : 1.0 - (sumSquaredError / sumSquaredTotal); return result; } } }