lsq_solver.dart 5.35 KB
Newer Older
1 2 3 4
// Copyright 2015 The Chromium Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

5 6
import 'dart:math' as math;
import 'dart:typed_data';
7

8
// TODO(abarth): Consider using vector_math.
9
class _Vector {
10
  _Vector(int size) : _offset = 0, _length = size, _elements = new Float64List(size);
11

12
  _Vector.fromVOL(List<double> values, int offset, int length)
13
    : _offset = offset, _length = length, _elements = values;
14

Ian Hickson's avatar
Ian Hickson committed
15 16 17 18 19
  final int _offset;

  final int _length;

  final List<double> _elements;
20

Ian Hickson's avatar
Ian Hickson committed
21 22 23 24
  double operator [](int i) => _elements[i + _offset];
  void operator []=(int i, double value) {
    _elements[i + _offset] = value;
  }
25

Ian Hickson's avatar
Ian Hickson committed
26
  double operator *(_Vector a) {
27
    double result = 0.0;
Ian Hickson's avatar
Ian Hickson committed
28
    for (int i = 0; i < _length; i += 1)
29 30 31 32 33 34 35
      result += this[i] * a[i];
    return result;
  }

  double norm() => math.sqrt(this * this);
}

36
// TODO(abarth): Consider using vector_math.
37 38
class _Matrix {
  _Matrix(int rows, int cols)
39
  : _columns = cols,
40
    _elements = new Float64List(rows * cols);
41

Ian Hickson's avatar
Ian Hickson committed
42 43 44
  final int _columns;
  final List<double> _elements;

45
  double get(int row, int col) => _elements[row * _columns + col];
46
  void set(int row, int col, double value) {
47
    _elements[row * _columns + col] = value;
48 49
  }

50 51 52 53 54
  _Vector getRow(int row) => new _Vector.fromVOL(
    _elements,
    row * _columns,
    _columns
  );
55 56
}

57
/// An nth degree polynomial fit to a dataset.
58
class PolynomialFit {
59 60 61
  /// Creates a polynomial fit of the given degree.
  ///
  /// There are n + 1 coefficients in a fit of degree n.
62 63
  PolynomialFit(int degree) : coefficients = new Float64List(degree + 1);

64
  /// The polynomial coefficients of the fit.
65
  final List<double> coefficients;
66 67 68 69

  /// An indicator of the quality of the fit.
  ///
  /// Larger values indicate greater quality.
70 71 72
  double confidence;
}

73
/// Uses the least-squares algorithm to fit a polynomial to a set of data.
74
class LeastSquaresSolver {
75 76
  /// Creates a least-squares solver.
  ///
77
  /// The [x], [y], and [w] arguments must not be null.
78 79 80 81 82
  LeastSquaresSolver(this.x, this.y, this.w) {
    assert(x.length == y.length);
    assert(y.length == w.length);
  }

83
  /// The x-coordinates of each data point.
84
  final List<double> x;
85 86

  /// The y-coordinates of each data point.
87
  final List<double> y;
88 89

  /// The weight to use for each data point.
90 91
  final List<double> w;

92
  /// Fits a polynomial of the given degree to the data points.
93
  PolynomialFit solve(int degree) {
Florian Loitsch's avatar
Florian Loitsch committed
94
    if (degree > x.length) // Not enough data to fit a curve.
95 96 97 98
      return null;

    PolynomialFit result = new PolynomialFit(degree);

Florian Loitsch's avatar
Florian Loitsch committed
99
    // Shorthands for the purpose of notation equivalence to original C++ code.
100 101 102 103
    final int m = x.length;
    final int n = degree + 1;

    // Expand the X vector to a matrix A, pre-multiplied by the weights.
104
    _Matrix a = new _Matrix(n, m);
Ian Hickson's avatar
Ian Hickson committed
105
    for (int h = 0; h < m; h += 1) {
106
      a.set(0, h, w[h]);
Ian Hickson's avatar
Ian Hickson committed
107
      for (int i = 1; i < n; i += 1)
108 109 110 111 112 113
        a.set(i, h, a.get(i - 1, h) * x[h]);
    }

    // Apply the Gram-Schmidt process to A to obtain its QR decomposition.

    // Orthonormal basis, column-major ordVectorer.
114
    _Matrix q = new _Matrix(n, m);
115
    // Upper triangular matrix, row-major order.
116
    _Matrix r = new _Matrix(n, n);
Ian Hickson's avatar
Ian Hickson committed
117 118
    for (int j = 0; j < n; j += 1) {
      for (int h = 0; h < m; h += 1)
119
        q.set(j, h, a.get(j, h));
Ian Hickson's avatar
Ian Hickson committed
120
      for (int i = 0; i < j; i += 1) {
121
        double dot = q.getRow(j) * q.getRow(i);
Ian Hickson's avatar
Ian Hickson committed
122
        for (int h = 0; h < m; h += 1)
123 124 125 126 127
          q.set(j, h, q.get(j, h) - dot * q.get(i, h));
      }

      double norm = q.getRow(j).norm();
      if (norm < 0.000001) {
Florian Loitsch's avatar
Florian Loitsch committed
128
        // Vectors are linearly dependent or zero so no solution.
129 130 131
        return null;
      }

132
      double inverseNorm = 1.0 / norm;
Ian Hickson's avatar
Ian Hickson committed
133
      for (int h = 0; h < m; h += 1)
134
        q.set(j, h, q.get(j, h) * inverseNorm);
Ian Hickson's avatar
Ian Hickson committed
135
      for (int i = 0; i < n; i += 1)
136
        r.set(j, i, i < j ? 0.0 : q.getRow(j) * a.getRow(i));
137 138 139 140
    }

    // 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.
141
    _Vector wy = new _Vector(m);
Ian Hickson's avatar
Ian Hickson committed
142
    for (int h = 0; h < m; h += 1)
143
      wy[h] = y[h] * w[h];
Ian Hickson's avatar
Ian Hickson committed
144
    for (int i = n - 1; i >= 0; i -= 1) {
145
      result.coefficients[i] = q.getRow(i) * wy;
Ian Hickson's avatar
Ian Hickson committed
146
      for (int j = n - 1; j > i; j -= 1)
147 148
        result.coefficients[i] -= r.get(i, j) * result.coefficients[j];
      result.coefficients[i] /= r.get(i, i);
149 150
    }

151
    // Calculate the coefficient of determination (confidence) as:
Ian Hickson's avatar
Ian Hickson committed
152 153
    //   1 - (sumSquaredError / sumSquaredTotal)
    // ...where sumSquaredError is the residual sum of squares (variance of the
154 155 156
    // error), and sumSquaredTotal is the total sum of squares (variance of the
    // data) where each has been weighted.
    double yMean = 0.0;
Ian Hickson's avatar
Ian Hickson committed
157
    for (int h = 0; h < m; h += 1)
158 159
      yMean += y[h];
    yMean /= m;
160

161 162
    double sumSquaredError = 0.0;
    double sumSquaredTotal = 0.0;
Ian Hickson's avatar
Ian Hickson committed
163
    for (int h = 0; h < m; h += 1) {
164
      double term = 1.0;
Ian Hickson's avatar
Ian Hickson committed
165 166
      double err = y[h] - result.coefficients[0];
      for (int i = 1; i < n; i += 1) {
167
        term *= x[h];
168
        err -= term * result.coefficients[i];
169
      }
170
      sumSquaredError += w[h] * w[h] * err * err;
Ian Hickson's avatar
Ian Hickson committed
171
      final double v = y[h] - yMean;
172
      sumSquaredTotal += w[h] * w[h] * v * v;
173 174
    }

Ian Hickson's avatar
Ian Hickson committed
175 176
    result.confidence = sumSquaredTotal <= 0.000001 ? 1.0 :
                          1.0 - (sumSquaredError / sumSquaredTotal);
177 178 179 180 181

    return result;
  }

}