new implementation of fitMaker that relies on DOT (see #49)

js/curve-fitting/model/FitMaker.js

@@ -1,10 +1,9 @@
 // Copyright 2015-2016, University of Colorado Boulder
 
-//TODO this needs documentation and cleanup
 /**
  * Fit maker for 'Curve Fitting' simulation.
  *
- * @author Andrey Zelenkov (Mlearner)
+ * @author Martin Veillette (Berea College)
  */
 define( function( require ) {
   'use strict';
@@ -13,112 +12,86 @@ define( function( require ) {
   var curveFitting = require( 'CURVE_FITTING/curveFitting' );
   var CurveFittingConstants = require( 'CURVE_FITTING/curve-fitting/CurveFittingConstants' );
   var inherit = require( 'PHET_CORE/inherit' );
+  var Matrix = require( 'DOT/Matrix' );
+
+  // constants
+  var solutionArrayLength = CurveFittingConstants.MAX_ORDER_OF_FIT + 1;
 
   /**
    * @constructor
    */
   function FitMaker() {
 
-    // @private max size of matrix
-    this.maxM = CurveFittingConstants.MAX_ORDER_OF_FIT + 1;
-    this.m = null;
-
     // @private
     this.solutionArray = [];
-
-    // @private matrix for further computing
-    this.matrix = [];
-    for ( var i = 0; i < this.maxM; i++ ) {
-      this.matrix[ i ] = [];
-    }
   }
 
   curveFitting.register( 'FitMaker', FitMaker );
 
   return inherit( Object, FitMaker, {
 
     /**
+     * Returns a numerical array of the coefficients of the polynomial
      * @param {Array.<Point>} points
      * @param {number} order
      * @returns {Array}
      * @public
      */
     getFit: function( points, order ) {
-      this.makeAugmentedMatrix( points, order );
-      this.reduceMatrix();
-      this.solveReducedMatrix();
+
+      // column matrix with the coefficients
+      var solutionMatrix = this.getSolutionMatrix( points, order );
+
+      // size of the column matrix
+      var numberOfRows = solutionMatrix.getRowDimension();
+
+      // the solutionArray has solutionArrayLength elements, regardless of the dimensions of the solution matrix
+      // the missing elements are padded with zeros
+      for ( var k = 0; k < solutionArrayLength; ++k ) {
+        this.solutionArray[ k ] = (k < numberOfRows) ? solutionMatrix.get( k, 0 ) : 0;
+      }
       return this.solutionArray;
     },
 
     /**
+     * Returns a column solution matrix, A, containing the coefficients of the polynomial
+     * The solution matrix is found by solving the equation, Y = X A
+     * where X is a square matrix, and Y is a column matrix.
+     *
+     * the number of rows of the solution matrix is given by the number of points, or
+     * the order +1, whichever is smaller.
+     *
+     * see http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html
+     *
      * @param {Array.<Point>} points
      * @param {number} order
+     * @returns {Matrix} the column matrix A
      * @private
      */
-    makeAugmentedMatrix: function( points, order ) {
+    getSolutionMatrix: function( points, order ) {
 
-      var numberOfPoints = points.length;
-      this.m = Math.min( order + 1, numberOfPoints );
+      // the rank of the matrix cannot be larger than the number of points
+      var m = Math.min( order + 1, points.length );
 
-      this.zeroMatrix();
+      var squareMatrix = new Matrix( m, m ); // matrix X
+      var columnMatrix = new Matrix( m, 1 ); // matrix Y
 
-      for ( var i = 0; i < numberOfPoints; ++i ) {
-        var delta = points[ i ].delta;
-        var deltaSquared = delta * delta;
-        var xPos = points[ i ].position.x;
-        var yPos = points[ i ].position.y;
-        for ( var j = 0; j < this.m; ++j ) {
-          for ( var k = 0; k < this.m; ++k ) {
-            this.matrix[ j ][ k ] = this.matrix[ j ][ k ] + Math.pow( xPos, j + k ) / deltaSquared;
-          }
-          this.matrix[ j ][ this.m ] = this.matrix[ j ][ this.m ] + Math.pow( xPos, j ) * yPos / deltaSquared;
-        }
-      }
-    },
+      // fill out the elements of the matrices Y and X
+      points.forEach( function( point ) {
+        var deltaSquared = point.delta * point.delta;
+        var x = point.position.x;
+        var y = point.position.y;
 
-    /**
-     * @private
-     */
-    reduceMatrix: function() {
-      for ( var i = 0; i < this.m - 1; ++i ) {
-        for ( var j = i + 1; j < this.m; ++j ) {
-          var l = this.matrix[ j ][ i ] / this.matrix[ i ][ i ];
-          for ( var k = 0; k < this.m + 1; ++k ) {
-            this.matrix[ j ][ k ] = this.matrix[ j ][ k ] - l * this.matrix[ i ][ k ];
+        for ( var j = 0; j < m; ++j ) {
+          for ( var k = 0; k < m; ++k ) {
+            squareMatrix.set( j, k, squareMatrix.get( j, k ) + Math.pow( x, j + k ) / deltaSquared );
           }
+          columnMatrix.set( j, 0, columnMatrix.get( j, 0 ) + Math.pow( x, j ) * y / deltaSquared );
         }
-      }
-      this.solveReducedMatrix();
-    },
+      } );
 
-    /**
-     * @private
-     */
-    solveReducedMatrix: function() {
-      var m = this.m;
-      for ( var i = m - 1; i >= 0; --i ) {
-        this.matrix[ i ][ m ] = this.matrix[ i ][ m ] / this.matrix[ i ][ i ];
-        this.matrix[ i ][ i ] = 1;
-        this.solutionArray[ i ] = this.matrix[ i ][ m ];
-        for ( var j = 0; j < i; ++j ) {
-          this.matrix[ j ][ m ] = this.matrix[ j ][ m ] - this.matrix[ j ][ i ] * this.solutionArray[ i ];
-          this.matrix[ j ][ i ] = 0;
-        }
-      }
-      for ( var k = m; k < this.maxM; ++k ) {
-        this.solutionArray[ k ] = 0;
-      }
-    },
-
-    /**
-     * @private
-     */
-    zeroMatrix: function() {
-      for ( var i = 0; i < this.maxM; ++i ) {
-        for ( var j = 0; j < this.maxM + 1; ++j ) {
-          this.matrix[ i ][ j ] = 0;
-        }
-      }
+      // the solution matrix, A, is X^-1 * Y
+      return squareMatrix.solve( columnMatrix );
     }
   } );
 } );