jsPsych/docs/demos/js/webgazer/worker_scripts/mat.js

(function () {
  "use strict";

  self.webgazer = self.webgazer || {};
  self.webgazer.mat = self.webgazer.mat || {};

  /**
   * Transposes an mxn array
   * @param {Array.<Array.<Number>>} matrix - of 'M x N' dimensionality
   * @return {Array.<Array.<Number>>} transposed matrix
   */
  self.webgazer.mat.transpose = function (matrix) {
    var m = matrix.length;
    var n = matrix[0].length;
    var transposedMatrix = new Array(n);

    for (var i = 0; i < m; i++) {
      for (var j = 0; j < n; j++) {
        if (i === 0) transposedMatrix[j] = new Array(m);
        transposedMatrix[j][i] = matrix[i][j];
      }
    }

    return transposedMatrix;
  };

  /**
   * Get a sub-matrix of matrix
   * @param {Array.<Array.<Number>>} matrix - original matrix
   * @param {Array.<Number>} r - Array of row indices
   * @param {Number} j0 - Initial column index
   * @param {Number} j1 - Final column index
   * @returns {Array} The sub-matrix matrix(r(:),j0:j1)
   */
  self.webgazer.mat.getMatrix = function (matrix, r, j0, j1) {
    var X = new Array(r.length),
      m = j1 - j0 + 1;

    for (var i = 0; i < r.length; i++) {
      X[i] = new Array(m);
      for (var j = j0; j <= j1; j++) {
        X[i][j - j0] = matrix[r[i]][j];
      }
    }
    return X;
  };

  /**
   * Get a submatrix of matrix
   * @param {Array.<Array.<Number>>} matrix - original matrix
   * @param {Number} i0 - Initial row index
   * @param {Number} i1 - Final row index
   * @param {Number} j0 - Initial column index
   * @param {Number} j1 - Final column index
   * @return {Array} The sub-matrix matrix(i0:i1,j0:j1)
   */
  self.webgazer.mat.getSubMatrix = function (matrix, i0, i1, j0, j1) {
    var size = j1 - j0 + 1,
      X = new Array(i1 - i0 + 1);

    for (var i = i0; i <= i1; i++) {
      var subI = i - i0;

      X[subI] = new Array(size);

      for (var j = j0; j <= j1; j++) {
        X[subI][j - j0] = matrix[i][j];
      }
    }
    return X;
  };

  /**
   * Linear algebraic matrix multiplication, matrix1 * matrix2
   * @param {Array.<Array.<Number>>} matrix1
   * @param {Array.<Array.<Number>>} matrix2
   * @return {Array.<Array.<Number>>} Matrix product, matrix1 * matrix2
   */
  self.webgazer.mat.mult = function (matrix1, matrix2) {
    if (matrix2.length != matrix1[0].length) {
      console.log("Matrix inner dimensions must agree:");
    }

    var X = new Array(matrix1.length),
      Bcolj = new Array(matrix1[0].length);

    for (var j = 0; j < matrix2[0].length; j++) {
      for (var k = 0; k < matrix1[0].length; k++) {
        Bcolj[k] = matrix2[k][j];
      }
      for (var i = 0; i < matrix1.length; i++) {
        if (j === 0) X[i] = new Array(matrix2[0].length);

        var Arowi = matrix1[i];
        var s = 0;
        for (var k = 0; k < matrix1[0].length; k++) {
          s += Arowi[k] * Bcolj[k];
        }
        X[i][j] = s;
      }
    }
    return X;
  };

  /**
   * LUDecomposition to solve A*X = B, based on WEKA code
   * @param {Array.<Array.<Number>>} A - left matrix of equation to be solved
   * @param {Array.<Array.<Number>>} B - right matrix of equation to be solved
   * @return {Array.<Array.<Number>>} X so that L*U*X = B(piv,:)
   */
  self.webgazer.mat.LUDecomposition = function (A, B) {
    var LU = new Array(A.length);

    for (var i = 0; i < A.length; i++) {
      LU[i] = new Array(A[0].length);
      for (var j = 0; j < A[0].length; j++) {
        LU[i][j] = A[i][j];
      }
    }

    var m = A.length;
    var n = A[0].length;
    var piv = new Array(m);
    for (var i = 0; i < m; i++) {
      piv[i] = i;
    }
    var pivsign = 1;
    var LUrowi = new Array();
    var LUcolj = new Array(m);
    // Outer loop.
    for (var j = 0; j < n; j++) {
      // Make a copy of the j-th column to localize references.
      for (var i = 0; i < m; i++) {
        LUcolj[i] = LU[i][j];
      }
      // Apply previous transformations.
      for (var i = 0; i < m; i++) {
        LUrowi = LU[i];
        // Most of the time is spent in the following dot product.
        var kmax = Math.min(i, j);
        var s = 0;
        for (var k = 0; k < kmax; k++) {
          s += LUrowi[k] * LUcolj[k];
        }
        LUrowi[j] = LUcolj[i] -= s;
      }
      // Find pivot and exchange if necessary.
      var p = j;
      for (var i = j + 1; i < m; i++) {
        if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
          p = i;
        }
      }
      if (p != j) {
        for (var k = 0; k < n; k++) {
          var t = LU[p][k];
          LU[p][k] = LU[j][k];
          LU[j][k] = t;
        }
        var k = piv[p];
        piv[p] = piv[j];
        piv[j] = k;
        pivsign = -pivsign;
      }
      // Compute multipliers.
      if ((j < m) & (LU[j][j] != 0)) {
        for (var i = j + 1; i < m; i++) {
          LU[i][j] /= LU[j][j];
        }
      }
    }
    if (B.length != m) {
      console.log("Matrix row dimensions must agree.");
    }
    for (var j = 0; j < n; j++) {
      if (LU[j][j] === 0) {
        console.log("Matrix is singular.");
      }
    }
    var nx = B[0].length;
    var X = self.webgazer.mat.getMatrix(B, piv, 0, nx - 1);
    // Solve L*Y = B(piv,:)
    for (var k = 0; k < n; k++) {
      for (var i = k + 1; i < n; i++) {
        for (var j = 0; j < nx; j++) {
          X[i][j] -= X[k][j] * LU[i][k];
        }
      }
    }
    // Solve U*X = Y;
    for (var k = n - 1; k >= 0; k--) {
      for (var j = 0; j < nx; j++) {
        X[k][j] /= LU[k][k];
      }
      for (var i = 0; i < k; i++) {
        for (var j = 0; j < nx; j++) {
          X[i][j] -= X[k][j] * LU[i][k];
        }
      }
    }
    return X;
  };

  /**
   * Least squares solution of A*X = B, based on WEKA code
   * @param {Array.<Array.<Number>>} A - left side matrix to be solved
   * @param {Array.<Array.<Number>>} B - a matrix with as many rows as A and any number of columns.
   * @return {Array.<Array.<Number>>} X - that minimizes the two norms of QR*X-B.
   */
  self.webgazer.mat.QRDecomposition = function (A, B) {
    // Initialize.
    var QR = new Array(A.length);

    for (var i = 0; i < A.length; i++) {
      QR[i] = new Array(A[0].length);
      for (var j = 0; j < A[0].length; j++) {
        QR[i][j] = A[i][j];
      }
    }
    var m = A.length;
    var n = A[0].length;
    var Rdiag = new Array(n);
    var nrm;

    // Main loop.
    for (var k = 0; k < n; k++) {
      // Compute 2-norm of k-th column without under/overflow.
      nrm = 0;
      for (var i = k; i < m; i++) {
        nrm = Math.hypot(nrm, QR[i][k]);
      }
      if (nrm != 0) {
        // Form k-th Householder vector.
        if (QR[k][k] < 0) {
          nrm = -nrm;
        }
        for (var i = k; i < m; i++) {
          QR[i][k] /= nrm;
        }
        QR[k][k] += 1;

        // Apply transformation to remaining columns.
        for (var j = k + 1; j < n; j++) {
          var s = 0;
          for (var i = k; i < m; i++) {
            s += QR[i][k] * QR[i][j];
          }
          s = -s / QR[k][k];
          for (var i = k; i < m; i++) {
            QR[i][j] += s * QR[i][k];
          }
        }
      }
      Rdiag[k] = -nrm;
    }
    if (B.length != m) {
      console.log("Matrix row dimensions must agree.");
    }
    for (var j = 0; j < n; j++) {
      if (Rdiag[j] === 0) console.log("Matrix is rank deficient");
    }
    // Copy right hand side
    var nx = B[0].length;
    var X = new Array(B.length);
    for (var i = 0; i < B.length; i++) {
      X[i] = new Array(B[0].length);
    }
    for (var i = 0; i < B.length; i++) {
      for (var j = 0; j < B[0].length; j++) {
        X[i][j] = B[i][j];
      }
    }
    // Compute Y = transpose(Q)*B
    for (var k = 0; k < n; k++) {
      for (var j = 0; j < nx; j++) {
        var s = 0.0;
        for (var i = k; i < m; i++) {
          s += QR[i][k] * X[i][j];
        }
        s = -s / QR[k][k];
        for (var i = k; i < m; i++) {
          X[i][j] += s * QR[i][k];
        }
      }
    }
    // Solve R*X = Y;
    for (var k = n - 1; k >= 0; k--) {
      for (var j = 0; j < nx; j++) {
        X[k][j] /= Rdiag[k];
      }
      for (var i = 0; i < k; i++) {
        for (var j = 0; j < nx; j++) {
          X[i][j] -= X[k][j] * QR[i][k];
        }
      }
    }
    return self.webgazer.mat.getSubMatrix(X, 0, n - 1, 0, nx - 1);
  };
})();