curve fitting and assorted fixes

lib/curve-fitter.js

@@ -0,0 +1,162 @@
+var invert = require('./matrix-invert.js');
+var matrices = [];
+
+var binomialCoefficients = [[1],[1,1]];
+
+function binomial(n,k) {
+  if (n===0) return 1;
+  var lut = binomialCoefficients;
+  while(n >= lut.length) {
+    var s = lut.length;
+    var nextRow = [1];
+    for(var i=1,prev=s-1; i<=s; i++) {
+      nextRow[i] = lut[prev][i-1] + lut[prev][i];
+    }
+    nextRow[s] = 1;
+    lut.push(nextRow);
+  }
+  return lut[n][k];
+}
+
+
+function dist(p1,p2) {
+  var dx = p1.x - p2.x, dy = p1.y - p2.y;
+  return Math.sqrt(dx*dx + dy*dy);
+}
+
+function transpose(M) {
+  var Mt = [];
+  M.forEach(row => Mt.push([]));
+  M.forEach((row,r) => row.forEach((v,c) => Mt[c][r] = v));
+  return Mt;
+}
+
+function row(M,i) {
+  return M[i];
+}
+
+function col(M,i) {
+  var col = [];
+  for(var r=0, l=M.length; r<l; r++) {
+    col.push(M[r][i]);
+  }
+  return col;
+}
+
+function multiply(M1, M2) {
+  // prep
+  var M = [];
+  var dims = [M1.length, M1[0].length, M2.length, M2[0].length];
+  // work
+  for (var r=0, c; r<dims[0]; r++) {
+    M[r] = [];
+    var _row = row(M1, r);
+    for (c=0; c<dims[3]; c++) {
+      var _col = col(M2,c);
+      var reducer = (a,v,i) => a + _col[i]*v;
+      M[r][c] = _row.reduce(reducer, 0);
+    }
+  }
+  return M;
+}
+
+function getValueColumn(P, prop) {
+  var col = [];
+  P.forEach(v => col.push([v[prop]]));
+  return col;
+}
+
+function computeBasisMatrix(n) {
+  /*
+    We can form any basis matrix using a generative approach:
+
+     - it's an M = (n x n) matrix
+     - it's a lower triangular matrix: all the entries above the main diagonal are zero
+     - the main diagonal consists of the binomial coefficients for n
+     - all entries are symmetric about the antidiagonal.
+
+    What's more, if we number rows and columns starting at 0, then
+    the value at position M[r,c], with row=r and column=c, can be
+    expressed as:
+
+      M[r,c] = (r choose c) * M[r,r] * S,
+
+      where S = 1 if r+c is even, or -1 otherwise
+
+    That is: the values in column c are directly computed off of the
+    binomial coefficients on the main diagonal, through multiplication
+    by a binomial based on matrix position, with the sign of the value
+    also determined by matrix position. This is actually very easy to
+    write out in code:
+  */
+
+  // form the square matrix, and set it to all zeroes
+  var M = [], i = n;
+  while (i--) { M[i] = "0".repeat(n).split('').map(v => parseInt(v)); }
+
+  // populate the main diagonal
+  var k = n - 1;
+  for (i=0; i<n; i++) { M[i][i] = binomial(k, i); }
+
+  // compute the remaining values
+  for (var c=0, r; c<n; c++) {
+    for (r=c+1; r<n; r++) {
+      var sign = (r+c)%2 ? -1 : 1;
+      var value = binomial(r, c) * M[r][r];
+      M[r][c] = sign * value; }}
+
+  return M;
+}
+
+function computeTimeValues(P, n) {
+  n = n || P.length;
+  var D = [0];
+  for(var i = 1; i<n; i++) {
+    D[i] = D[i-1] + dist(P[i-1], P[i]);
+  }
+  var S = [0], len = D[n-1];
+  D.forEach((v,i) => { S[i] = v/len; });
+  return S;
+}
+
+function raiseRowPower(row, i) {
+  return row.map(v => Math.pow(v,i));
+}
+
+function formTMatrix(S, n) {
+  n = n || S.length;
+  var Tp = [];
+  // it's easier to generate the transposed matrix:
+  for(var i=0; i<n; i++) Tp.push( raiseRowPower(S, i));
+  return {
+    Tt: Tp,
+    T: transpose(Tp) // and then transpose "again" to get the real matrix
+  };
+}
+
+function computeBestFit(P, M, S, n) {
+  n = n || P.length;
+  var tm = formTMatrix(S, n),
+      T = tm.T,
+      Tt = tm.Tt,
+      M1 = invert(M),
+      TtT1 = invert(multiply(Tt,T)),
+      step1 = multiply(TtT1, Tt),
+      step2 = multiply(M1, step1),
+      X = getValueColumn(P,'x'),
+      Cx = multiply(step2, X),
+      Y = getValueColumn(P,'y'),
+      Cy = multiply(step2, Y);
+  return { x: Cx, y: Cy };
+}
+
+function fit(points) {
+  var n = points.length,
+      P = Array.from(points),
+      M = computeBasisMatrix(n),
+      S = computeTimeValues(P, n),
+      C = computeBestFit(P, M, S, n);
+  return { n, P, M, S, C };
+}
+
+module.exports = window.makeFit = fit;

lib/latex-loader.js

@@ -23,9 +23,8 @@ var options = {
   ],
 
   /**
-   * We look for MathJax/KaTeX style data, and make sure
-   * it is escaped properly so that JSX conversion will
-   * still work.
+   * We look for MathJax/KaTeX style data, and generate a static
+   * SVG file for it, returning the <img> code to load that file.
    */
   process: function escapeBlockLaTeX(latex) {
     latex = cleanup(latex);

lib/matrix-invert.js

@@ -0,0 +1,109 @@
+// Copied from http://blog.acipo.com/matrix-inversion-in-javascript/
+
+// Returns the inverse of matrix `M`.
+module.exports = function matrix_invert(M) {
+  // I use Guassian Elimination to calculate the inverse:
+  // (1) 'augment' the matrix (left) by the identity (on the right)
+  // (2) Turn the matrix on the left into the identity by elemetry row ops
+  // (3) The matrix on the right is the inverse (was the identity matrix)
+  // There are 3 elemtary row ops: (I combine b and c in my code)
+  // (a) Swap 2 rows
+  // (b) Multiply a row by a scalar
+  // (c) Add 2 rows
+
+  //if the matrix isn't square: exit (error)
+  if (M.length !== M[0].length) {
+    return;
+  }
+
+  //create the identity matrix (I), and a copy (C) of the original
+  var i = 0,
+      ii = 0,
+      j = 0,
+      dim = M.length,
+      e = 0,
+      t = 0;
+  var I = [],
+      C = [];
+  for (i = 0; i < dim; i += 1) {
+    // Create the row
+    I[I.length] = [];
+    C[C.length] = [];
+    for (j = 0; j < dim; j += 1) {
+      //if we're on the diagonal, put a 1 (for identity)
+      if (i == j) {
+        I[i][j] = 1;
+      } else {
+        I[i][j] = 0;
+      }
+
+      // Also, make the copy of the original
+      C[i][j] = M[i][j];
+    }
+  }
+
+  // Perform elementary row operations
+  for (i = 0; i < dim; i += 1) {
+    // get the element e on the diagonal
+    e = C[i][i];
+
+    // if we have a 0 on the diagonal (we'll need to swap with a lower row)
+    if (e == 0) {
+      //look through every row below the i'th row
+      for (ii = i + 1; ii < dim; ii += 1) {
+        //if the ii'th row has a non-0 in the i'th col
+        if (C[ii][i] != 0) {
+          //it would make the diagonal have a non-0 so swap it
+          for (j = 0; j < dim; j++) {
+            e = C[i][j]; //temp store i'th row
+            C[i][j] = C[ii][j]; //replace i'th row by ii'th
+            C[ii][j] = e; //repace ii'th by temp
+            e = I[i][j]; //temp store i'th row
+            I[i][j] = I[ii][j]; //replace i'th row by ii'th
+            I[ii][j] = e; //repace ii'th by temp
+          }
+          //don't bother checking other rows since we've swapped
+          break;
+        }
+      }
+      //get the new diagonal
+      e = C[i][i];
+      //if it's still 0, not invertable (error)
+      if (e == 0) {
+        return;
+      }
+    }
+
+    // Scale this row down by e (so we have a 1 on the diagonal)
+    for (j = 0; j < dim; j++) {
+      C[i][j] = C[i][j] / e; //apply to original matrix
+      I[i][j] = I[i][j] / e; //apply to identity
+    }
+
+    // Subtract this row (scaled appropriately for each row) from ALL of
+    // the other rows so that there will be 0's in this column in the
+    // rows above and below this one
+    for (ii = 0; ii < dim; ii++) {
+      // Only apply to other rows (we want a 1 on the diagonal)
+      if (ii == i) {
+        continue;
+      }
+
+      // We want to change this element to 0
+      e = C[ii][i];
+
+      // Subtract (the row above(or below) scaled by e) from (the
+      // current row) but start at the i'th column and assume all the
+      // stuff left of diagonal is 0 (which it should be if we made this
+      // algorithm correctly)
+      for (j = 0; j < dim; j++) {
+        C[ii][j] -= e * C[i][j]; //apply to original matrix
+        I[ii][j] -= e * I[i][j]; //apply to identity
+      }
+    }
+  }
+
+  //we've done all operations, C should be the identity
+  //matrix I should be the inverse:
+  return I;
+};

lib/site/handlers.js

@@ -2356,6 +2356,64 @@ return {
       api.drawCurve(curve);
     }
   }
+};
+ }())
+  },
+  "curvefitting": {
+    handler: (function() { var fit = require('../../../lib/curve-fitter.js');
+
+return {
+  setup: function(api) {
+    this.api = api;
+    this.reset();
+  },
+
+  reset: function() {
+    this.points = [];
+    this.curveset = false;
+    let api = this.api;
+    if (api) {
+      api.setCurve(false);
+      api.reset();
+      api.redraw();
+    }
+  },
+
+  draw: function(api, curve) {
+    api.setPanelCount(1);
+    api.reset();
+    api.setColor('lightgrey');
+    api.drawGrid(10,10);
+
+    api.setColor('black');
+    if (!this.curveset && this.points.length > 2) {
+      let bestFitData = fit(this.points),
+          x = bestFitData.C.x,
+          y = bestFitData.C.y,
+          bpoints = [];
+      x.forEach((r,i) => {
+        bpoints.push({
+          x: r[0],
+          y: y[i][0]
+        });
+      });
+      curve = new api.Bezier(bpoints);
+      api.setCurve(curve);
+      this.curveset = true;
+    }
+
+    if (curve) {
+      api.drawCurve(curve);
+      api.drawSkeleton(curve);
+    }
+    api.drawPoints(this.points);
+  },
+
+  onClick: function(evt, api) {
+    this.curveset = false;
+    this.points.push({x: api.mx, y: api.my });
+    api.redraw();
+  }
 };
  }())
   },