854 lines
21 KiB
JavaScript
854 lines
21 KiB
JavaScript
// Copyright (c) 2011, Chris Umbel, James Coglan
|
|
// Matrix class - depends on Vector.
|
|
|
|
var fs = require('fs');
|
|
var Sylvester = require('./sylvester');
|
|
var Vector = require('./vector');
|
|
|
|
// augment a matrix M with identity rows/cols
|
|
function identSize(M, m, n, k) {
|
|
var e = M.elements;
|
|
var i = k - 1;
|
|
|
|
while(i--) {
|
|
var row = [];
|
|
|
|
for(var j = 0; j < n; j++)
|
|
row.push(j == i ? 1 : 0);
|
|
|
|
e.unshift(row);
|
|
}
|
|
|
|
for(var i = k - 1; i < m; i++) {
|
|
while(e[i].length < n)
|
|
e[i].unshift(0);
|
|
}
|
|
|
|
return $M(e);
|
|
}
|
|
|
|
function pca(X) {
|
|
var Sigma = X.transpose().x(X).x(1 / X.rows());
|
|
var svd = Sigma.svd();
|
|
return {U: svd.U, S: svd.S};
|
|
}
|
|
|
|
function Matrix() {}
|
|
Matrix.prototype = {
|
|
pcaProject: function(k, U) {
|
|
var U = U || pca(this).U;
|
|
var Ureduce= U.slice(1, U.rows(), 1, k);
|
|
return {Z: this.x(Ureduce), U: U};
|
|
},
|
|
|
|
pcaRecover: function(U) {
|
|
var k = this.cols();
|
|
var Ureduce = U.slice(1, U.rows(), 1, k);
|
|
return this.x(Ureduce.transpose());
|
|
},
|
|
|
|
triu: function(k) {
|
|
if(!k)
|
|
k = 0;
|
|
|
|
return this.map(function(x, i, j) {
|
|
return j - i >= k ? x : 0;
|
|
});
|
|
},
|
|
|
|
svd: function() {
|
|
var A = this;
|
|
var U = Matrix.I(A.rows());
|
|
var S = A.transpose();
|
|
var V = Matrix.I(A.cols());
|
|
var err = Number.MAX_VALUE;
|
|
var i = 0;
|
|
var maxLoop = 100;
|
|
|
|
while(err > 2.2737e-13 && i < maxLoop) {
|
|
var qr = S.transpose().qr();
|
|
S = qr.R;
|
|
U = U.x(qr.Q);
|
|
qr = S.transpose().qr();
|
|
V = V.x(qr.Q);
|
|
S = qr.R;
|
|
|
|
var e = S.triu(1).unroll().norm();
|
|
var f = S.diagonal().norm();
|
|
|
|
if(f == 0)
|
|
f = 1;
|
|
|
|
err = e / f;
|
|
|
|
i++;
|
|
}
|
|
|
|
var ss = S.diagonal();
|
|
var s = [];
|
|
|
|
for(var i = 1; i <= ss.cols(); i++) {
|
|
var ssn = ss.e(i);
|
|
s.push(Math.abs(ssn));
|
|
|
|
if(ssn < 0) {
|
|
for(var j = 0; j < U.rows(); j++) {
|
|
U.elements[j][i - 1] = -(U.elements[j][i - 1]);
|
|
}
|
|
}
|
|
}
|
|
|
|
return {U: U, S: $V(s).toDiagonalMatrix(), V: V};
|
|
},
|
|
|
|
unroll: function() {
|
|
var v = [];
|
|
|
|
for(var i = 1; i <= this.cols(); i++) {
|
|
for(var j = 1; j <= this.rows(); j++) {
|
|
v.push(this.e(j, i));
|
|
}
|
|
}
|
|
|
|
return $V(v);
|
|
},
|
|
|
|
qr: function() {
|
|
var m = this.rows();
|
|
var n = this.cols();
|
|
var Q = Matrix.I(m);
|
|
var A = this;
|
|
|
|
for(var k = 1; k < Math.min(m, n); k++) {
|
|
var ak = A.slice(k, 0, k, k).col(1);
|
|
var oneZero = [1];
|
|
|
|
while(oneZero.length <= m - k)
|
|
oneZero.push(0);
|
|
|
|
oneZero = $V(oneZero);
|
|
var vk = ak.add(oneZero.x(ak.norm() * Math.sign(ak.e(1))));
|
|
var Vk = $M(vk);
|
|
var Hk = Matrix.I(m - k + 1).subtract(Vk.x(2).x(Vk.transpose()).div(Vk.transpose().x(Vk).e(1, 1)));
|
|
var Qk = identSize(Hk, m, n, k);
|
|
A = Qk.x(A);
|
|
Q = Q.x(Qk);
|
|
}
|
|
|
|
return {Q: Q, R: A};
|
|
},
|
|
|
|
|
|
slice: function(startRow, endRow, startCol, endCol) {
|
|
var x = [];
|
|
|
|
if(endRow == 0)
|
|
endRow = this.rows();
|
|
|
|
if(endCol == 0)
|
|
endCol = this.cols();
|
|
|
|
for(i = startRow; i <= endRow; i++) {
|
|
var row = [];
|
|
|
|
for(j = startCol; j <= endCol; j++) {
|
|
row.push(this.e(i, j));
|
|
}
|
|
|
|
x.push(row);
|
|
}
|
|
|
|
return $M(x);
|
|
},
|
|
|
|
// Returns element (i,j) of the matrix
|
|
e: function(i,j) {
|
|
if (i < 1 || i > this.elements.length || j < 1 || j > this.elements[0].length) { return null; }
|
|
return this.elements[i - 1][j - 1];
|
|
},
|
|
|
|
// Returns row k of the matrix as a vector
|
|
row: function(i) {
|
|
if (i > this.elements.length) { return null; }
|
|
return $V(this.elements[i - 1]);
|
|
},
|
|
|
|
// Returns column k of the matrix as a vector
|
|
col: function(j) {
|
|
if (j > this.elements[0].length) { return null; }
|
|
var col = [], n = this.elements.length;
|
|
for (var i = 0; i < n; i++) { col.push(this.elements[i][j - 1]); }
|
|
return $V(col);
|
|
},
|
|
|
|
// Returns the number of rows/columns the matrix has
|
|
dimensions: function() {
|
|
return {rows: this.elements.length, cols: this.elements[0].length};
|
|
},
|
|
|
|
// Returns the number of rows in the matrix
|
|
rows: function() {
|
|
return this.elements.length;
|
|
},
|
|
|
|
// Returns the number of columns in the matrix
|
|
cols: function() {
|
|
return this.elements[0].length;
|
|
},
|
|
|
|
// Returns true iff the matrix is equal to the argument. You can supply
|
|
// a vector as the argument, in which case the receiver must be a
|
|
// one-column matrix equal to the vector.
|
|
eql: function(matrix) {
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined') { M = Matrix.create(M).elements; }
|
|
if (this.elements.length != M.length ||
|
|
this.elements[0].length != M[0].length) { return false; }
|
|
var i = this.elements.length, nj = this.elements[0].length, j;
|
|
while (i--) { j = nj;
|
|
while (j--) {
|
|
if (Math.abs(this.elements[i][j] - M[i][j]) > Sylvester.precision) { return false; }
|
|
}
|
|
}
|
|
return true;
|
|
},
|
|
|
|
// Returns a copy of the matrix
|
|
dup: function() {
|
|
return Matrix.create(this.elements);
|
|
},
|
|
|
|
// Maps the matrix to another matrix (of the same dimensions) according to the given function
|
|
map: function(fn) {
|
|
var els = [], i = this.elements.length, nj = this.elements[0].length, j;
|
|
while (i--) { j = nj;
|
|
els[i] = [];
|
|
while (j--) {
|
|
els[i][j] = fn(this.elements[i][j], i + 1, j + 1);
|
|
}
|
|
}
|
|
return Matrix.create(els);
|
|
},
|
|
|
|
// Returns true iff the argument has the same dimensions as the matrix
|
|
isSameSizeAs: function(matrix) {
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined') { M = Matrix.create(M).elements; }
|
|
return (this.elements.length == M.length &&
|
|
this.elements[0].length == M[0].length);
|
|
},
|
|
|
|
// Returns the result of adding the argument to the matrix
|
|
add: function(matrix) {
|
|
if(typeof(matrix) == 'number') {
|
|
return this.map(function(x, i, j) { return x + matrix});
|
|
} else {
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined') { M = Matrix.create(M).elements; }
|
|
if (!this.isSameSizeAs(M)) { return null; }
|
|
return this.map(function(x, i, j) { return x + M[i - 1][j - 1]; });
|
|
}
|
|
},
|
|
|
|
// Returns the result of subtracting the argument from the matrix
|
|
subtract: function(matrix) {
|
|
if(typeof(matrix) == 'number') {
|
|
return this.map(function(x, i, j) { return x - matrix});
|
|
} else {
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined') { M = Matrix.create(M).elements; }
|
|
if (!this.isSameSizeAs(M)) { return null; }
|
|
return this.map(function(x, i, j) { return x - M[i - 1][j - 1]; });
|
|
}
|
|
},
|
|
|
|
// Returns true iff the matrix can multiply the argument from the left
|
|
canMultiplyFromLeft: function(matrix) {
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined') { M = Matrix.create(M).elements; }
|
|
// this.columns should equal matrix.rows
|
|
return (this.elements[0].length == M.length);
|
|
},
|
|
|
|
// Returns the result of a multiplication-style operation the matrix from the right by the argument.
|
|
// If the argument is a scalar then just operate on all the elements. If the argument is
|
|
// a vector, a vector is returned, which saves you having to remember calling
|
|
// col(1) on the result.
|
|
mulOp: function(matrix, op) {
|
|
if (!matrix.elements) {
|
|
return this.map(function(x) { return op(x, matrix); });
|
|
}
|
|
|
|
var returnVector = matrix.modulus ? true : false;
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined')
|
|
M = Matrix.create(M).elements;
|
|
if (!this.canMultiplyFromLeft(M))
|
|
return null;
|
|
var e = this.elements, rowThis, rowElem, elements = [],
|
|
sum, m = e.length, n = M[0].length, o = e[0].length, i = m, j, k;
|
|
|
|
while (i--) {
|
|
rowElem = [];
|
|
rowThis = e[i];
|
|
j = n;
|
|
|
|
while (j--) {
|
|
sum = 0;
|
|
k = o;
|
|
|
|
while (k--) {
|
|
sum += op(rowThis[k], M[k][j]);
|
|
}
|
|
|
|
rowElem[j] = sum;
|
|
}
|
|
|
|
elements[i] = rowElem;
|
|
}
|
|
|
|
var M = Matrix.create(elements);
|
|
return returnVector ? M.col(1) : M;
|
|
},
|
|
|
|
// Returns the result of dividing the matrix from the right by the argument.
|
|
// If the argument is a scalar then just divide all the elements. If the argument is
|
|
// a vector, a vector is returned, which saves you having to remember calling
|
|
// col(1) on the result.
|
|
div: function(matrix) {
|
|
return this.mulOp(matrix, function(x, y) { return x / y});
|
|
},
|
|
|
|
// Returns the result of multiplying the matrix from the right by the argument.
|
|
// If the argument is a scalar then just multiply all the elements. If the argument is
|
|
// a vector, a vector is returned, which saves you having to remember calling
|
|
// col(1) on the result.
|
|
multiply: function(matrix) {
|
|
return this.mulOp(matrix, function(x, y) { return x * y});
|
|
},
|
|
|
|
x: function(matrix) { return this.multiply(matrix); },
|
|
|
|
elementMultiply: function(v) {
|
|
return this.map(function(k, i, j) {
|
|
return v.e(i, j) * k;
|
|
});
|
|
},
|
|
|
|
sum: function() {
|
|
var sum = 0;
|
|
|
|
this.map(function(x) { sum += x;});
|
|
|
|
return sum;
|
|
},
|
|
|
|
// Returns a Vector of each colum averaged.
|
|
mean: function() {
|
|
var dim = this.dimensions();
|
|
var r = [];
|
|
for (var i = 1; i <= dim.cols; i++) {
|
|
r.push(this.col(i).sum() / dim.rows);
|
|
}
|
|
return $V(r);
|
|
},
|
|
|
|
column: function(n) {
|
|
return this.col(n);
|
|
},
|
|
|
|
log: function() {
|
|
return this.map(function(x) { return Math.log(x); });
|
|
},
|
|
|
|
// Returns a submatrix taken from the matrix
|
|
// Argument order is: start row, start col, nrows, ncols
|
|
// Element selection wraps if the required index is outside the matrix's bounds, so you could
|
|
// use this to perform row/column cycling or copy-augmenting.
|
|
minor: function(a, b, c, d) {
|
|
var elements = [], ni = c, i, nj, j;
|
|
var rows = this.elements.length, cols = this.elements[0].length;
|
|
while (ni--) {
|
|
i = c - ni - 1;
|
|
elements[i] = [];
|
|
nj = d;
|
|
while (nj--) {
|
|
j = d - nj - 1;
|
|
elements[i][j] = this.elements[(a + i - 1) % rows][(b + j - 1) % cols];
|
|
}
|
|
}
|
|
return Matrix.create(elements);
|
|
},
|
|
|
|
// Returns the transpose of the matrix
|
|
transpose: function() {
|
|
var rows = this.elements.length, i, cols = this.elements[0].length, j;
|
|
var elements = [], i = cols;
|
|
while (i--) {
|
|
j = rows;
|
|
elements[i] = [];
|
|
while (j--) {
|
|
elements[i][j] = this.elements[j][i];
|
|
}
|
|
}
|
|
return Matrix.create(elements);
|
|
},
|
|
|
|
// Returns true iff the matrix is square
|
|
isSquare: function() {
|
|
return (this.elements.length == this.elements[0].length);
|
|
},
|
|
|
|
// Returns the (absolute) largest element of the matrix
|
|
max: function() {
|
|
var m = 0, i = this.elements.length, nj = this.elements[0].length, j;
|
|
while (i--) {
|
|
j = nj;
|
|
while (j--) {
|
|
if (Math.abs(this.elements[i][j]) > Math.abs(m)) { m = this.elements[i][j]; }
|
|
}
|
|
}
|
|
return m;
|
|
},
|
|
|
|
// Returns the indeces of the first match found by reading row-by-row from left to right
|
|
indexOf: function(x) {
|
|
var index = null, ni = this.elements.length, i, nj = this.elements[0].length, j;
|
|
for (i = 0; i < ni; i++) {
|
|
for (j = 0; j < nj; j++) {
|
|
if (this.elements[i][j] == x) { return {i: i + 1, j: j + 1}; }
|
|
}
|
|
}
|
|
return null;
|
|
},
|
|
|
|
// If the matrix is square, returns the diagonal elements as a vector.
|
|
// Otherwise, returns null.
|
|
diagonal: function() {
|
|
if (!this.isSquare) { return null; }
|
|
var els = [], n = this.elements.length;
|
|
for (var i = 0; i < n; i++) {
|
|
els.push(this.elements[i][i]);
|
|
}
|
|
return $V(els);
|
|
},
|
|
|
|
// Make the matrix upper (right) triangular by Gaussian elimination.
|
|
// This method only adds multiples of rows to other rows. No rows are
|
|
// scaled up or switched, and the determinant is preserved.
|
|
toRightTriangular: function() {
|
|
var M = this.dup(), els;
|
|
var n = this.elements.length, i, j, np = this.elements[0].length, p;
|
|
for (i = 0; i < n; i++) {
|
|
if (M.elements[i][i] == 0) {
|
|
for (j = i + 1; j < n; j++) {
|
|
if (M.elements[j][i] != 0) {
|
|
els = [];
|
|
for (p = 0; p < np; p++) { els.push(M.elements[i][p] + M.elements[j][p]); }
|
|
M.elements[i] = els;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
if (M.elements[i][i] != 0) {
|
|
for (j = i + 1; j < n; j++) {
|
|
var multiplier = M.elements[j][i] / M.elements[i][i];
|
|
els = [];
|
|
for (p = 0; p < np; p++) {
|
|
// Elements with column numbers up to an including the number
|
|
// of the row that we're subtracting can safely be set straight to
|
|
// zero, since that's the point of this routine and it avoids having
|
|
// to loop over and correct rounding errors later
|
|
els.push(p <= i ? 0 : M.elements[j][p] - M.elements[i][p] * multiplier);
|
|
}
|
|
M.elements[j] = els;
|
|
}
|
|
}
|
|
}
|
|
return M;
|
|
},
|
|
|
|
toUpperTriangular: function() { return this.toRightTriangular(); },
|
|
|
|
// Returns the determinant for square matrices
|
|
determinant: function() {
|
|
if (!this.isSquare()) { return null; }
|
|
if (this.cols == 1 && this.rows == 1) { return this.row(1); }
|
|
if (this.cols == 0 && this.rows == 0) { return 1; }
|
|
var M = this.toRightTriangular();
|
|
var det = M.elements[0][0], n = M.elements.length;
|
|
for (var i = 1; i < n; i++) {
|
|
det = det * M.elements[i][i];
|
|
}
|
|
return det;
|
|
},
|
|
det: function() { return this.determinant(); },
|
|
|
|
// Returns true iff the matrix is singular
|
|
isSingular: function() {
|
|
return (this.isSquare() && this.determinant() === 0);
|
|
},
|
|
|
|
// Returns the trace for square matrices
|
|
trace: function() {
|
|
if (!this.isSquare()) { return null; }
|
|
var tr = this.elements[0][0], n = this.elements.length;
|
|
for (var i = 1; i < n; i++) {
|
|
tr += this.elements[i][i];
|
|
}
|
|
return tr;
|
|
},
|
|
|
|
tr: function() { return this.trace(); },
|
|
|
|
// Returns the rank of the matrix
|
|
rank: function() {
|
|
var M = this.toRightTriangular(), rank = 0;
|
|
var i = this.elements.length, nj = this.elements[0].length, j;
|
|
while (i--) {
|
|
j = nj;
|
|
while (j--) {
|
|
if (Math.abs(M.elements[i][j]) > Sylvester.precision) { rank++; break; }
|
|
}
|
|
}
|
|
return rank;
|
|
},
|
|
|
|
rk: function() { return this.rank(); },
|
|
|
|
// Returns the result of attaching the given argument to the right-hand side of the matrix
|
|
augment: function(matrix) {
|
|
var M = matrix.elements || matrix;
|
|
if (typeof(M[0][0]) == 'undefined') { M = Matrix.create(M).elements; }
|
|
var T = this.dup(), cols = T.elements[0].length;
|
|
var i = T.elements.length, nj = M[0].length, j;
|
|
if (i != M.length) { return null; }
|
|
while (i--) {
|
|
j = nj;
|
|
while (j--) {
|
|
T.elements[i][cols + j] = M[i][j];
|
|
}
|
|
}
|
|
return T;
|
|
},
|
|
|
|
// Returns the inverse (if one exists) using Gauss-Jordan
|
|
inverse: function() {
|
|
if (!this.isSquare() || this.isSingular()) { return null; }
|
|
var n = this.elements.length, i = n, j;
|
|
var M = this.augment(Matrix.I(n)).toRightTriangular();
|
|
var np = M.elements[0].length, p, els, divisor;
|
|
var inverse_elements = [], new_element;
|
|
// Matrix is non-singular so there will be no zeros on the diagonal
|
|
// Cycle through rows from last to first
|
|
while (i--) {
|
|
// First, normalise diagonal elements to 1
|
|
els = [];
|
|
inverse_elements[i] = [];
|
|
divisor = M.elements[i][i];
|
|
for (p = 0; p < np; p++) {
|
|
new_element = M.elements[i][p] / divisor;
|
|
els.push(new_element);
|
|
// Shuffle off the current row of the right hand side into the results
|
|
// array as it will not be modified by later runs through this loop
|
|
if (p >= n) { inverse_elements[i].push(new_element); }
|
|
}
|
|
M.elements[i] = els;
|
|
// Then, subtract this row from those above it to
|
|
// give the identity matrix on the left hand side
|
|
j = i;
|
|
while (j--) {
|
|
els = [];
|
|
for (p = 0; p < np; p++) {
|
|
els.push(M.elements[j][p] - M.elements[i][p] * M.elements[j][i]);
|
|
}
|
|
M.elements[j] = els;
|
|
}
|
|
}
|
|
return Matrix.create(inverse_elements);
|
|
},
|
|
|
|
inv: function() { return this.inverse(); },
|
|
|
|
// Returns the result of rounding all the elements
|
|
round: function() {
|
|
return this.map(function(x) { return Math.round(x); });
|
|
},
|
|
|
|
// Returns a copy of the matrix with elements set to the given value if they
|
|
// differ from it by less than Sylvester.precision
|
|
snapTo: function(x) {
|
|
return this.map(function(p) {
|
|
return (Math.abs(p - x) <= Sylvester.precision) ? x : p;
|
|
});
|
|
},
|
|
|
|
// Returns a string representation of the matrix
|
|
inspect: function() {
|
|
var matrix_rows = [];
|
|
var n = this.elements.length;
|
|
for (var i = 0; i < n; i++) {
|
|
matrix_rows.push($V(this.elements[i]).inspect());
|
|
}
|
|
return matrix_rows.join('\n');
|
|
},
|
|
|
|
// Returns a array representation of the matrix
|
|
toArray: function() {
|
|
var matrix_rows = [];
|
|
var n = this.elements.length;
|
|
for (var i = 0; i < n; i++) {
|
|
matrix_rows.push(this.elements[i]);
|
|
}
|
|
return matrix_rows;
|
|
},
|
|
|
|
|
|
// Set the matrix's elements from an array. If the argument passed
|
|
// is a vector, the resulting matrix will be a single column.
|
|
setElements: function(els) {
|
|
var i, j, elements = els.elements || els;
|
|
if (typeof(elements[0][0]) != 'undefined') {
|
|
i = elements.length;
|
|
this.elements = [];
|
|
while (i--) {
|
|
j = elements[i].length;
|
|
this.elements[i] = [];
|
|
while (j--) {
|
|
this.elements[i][j] = elements[i][j];
|
|
}
|
|
}
|
|
return this;
|
|
}
|
|
var n = elements.length;
|
|
this.elements = [];
|
|
for (i = 0; i < n; i++) {
|
|
this.elements.push([elements[i]]);
|
|
}
|
|
return this;
|
|
},
|
|
|
|
maxColumnIndexes: function() {
|
|
var maxes = [];
|
|
|
|
for(var i = 1; i <= this.rows(); i++) {
|
|
var max = null;
|
|
var maxIndex = -1;
|
|
|
|
for(var j = 1; j <= this.cols(); j++) {
|
|
if(max === null || this.e(i, j) > max) {
|
|
max = this.e(i, j);
|
|
maxIndex = j;
|
|
}
|
|
}
|
|
|
|
maxes.push(maxIndex);
|
|
}
|
|
|
|
return $V(maxes);
|
|
},
|
|
|
|
maxColumns: function() {
|
|
var maxes = [];
|
|
|
|
for(var i = 1; i <= this.rows(); i++) {
|
|
var max = null;
|
|
|
|
for(var j = 1; j <= this.cols(); j++) {
|
|
if(max === null || this.e(i, j) > max) {
|
|
max = this.e(i, j);
|
|
}
|
|
}
|
|
|
|
maxes.push(max);
|
|
}
|
|
|
|
return $V(maxes);
|
|
},
|
|
|
|
minColumnIndexes: function() {
|
|
var mins = [];
|
|
|
|
for(var i = 1; i <= this.rows(); i++) {
|
|
var min = null;
|
|
var minIndex = -1;
|
|
|
|
for(var j = 1; j <= this.cols(); j++) {
|
|
if(min === null || this.e(i, j) < min) {
|
|
min = this.e(i, j);
|
|
minIndex = j;
|
|
}
|
|
}
|
|
|
|
mins.push(minIndex);
|
|
}
|
|
|
|
return $V(mins);
|
|
},
|
|
|
|
minColumns: function() {
|
|
var mins = [];
|
|
|
|
for(var i = 1; i <= this.rows(); i++) {
|
|
var min = null;
|
|
|
|
for(var j = 1; j <= this.cols(); j++) {
|
|
if(min === null || this.e(i, j) < min) {
|
|
min = this.e(i, j);
|
|
}
|
|
}
|
|
|
|
mins.push(min);
|
|
}
|
|
|
|
return $V(mins);
|
|
}
|
|
};
|
|
|
|
// Constructor function
|
|
Matrix.create = function(elements) {
|
|
var M = new Matrix();
|
|
return M.setElements(elements);
|
|
};
|
|
|
|
// Identity matrix of size n
|
|
Matrix.I = function(n) {
|
|
var els = [], i = n, j;
|
|
while (i--) {
|
|
j = n;
|
|
els[i] = [];
|
|
while (j--) {
|
|
els[i][j] = (i == j) ? 1 : 0;
|
|
}
|
|
}
|
|
return Matrix.create(els);
|
|
};
|
|
|
|
Matrix.loadFile = function(file) {
|
|
var contents = fs.readFileSync(file, 'utf-8');
|
|
var matrix = [];
|
|
|
|
var rowArray = contents.split('\n');
|
|
for (var i = 0; i < rowArray.length; i++) {
|
|
var d = rowArray[i].split(',');
|
|
if (d.length > 1) {
|
|
matrix.push(d);
|
|
}
|
|
}
|
|
|
|
var M = new Matrix();
|
|
return M.setElements(matrix);
|
|
};
|
|
|
|
// Diagonal matrix - all off-diagonal elements are zero
|
|
Matrix.Diagonal = function(elements) {
|
|
var i = elements.length;
|
|
var M = Matrix.I(i);
|
|
while (i--) {
|
|
M.elements[i][i] = elements[i];
|
|
}
|
|
return M;
|
|
};
|
|
|
|
// Rotation matrix about some axis. If no axis is
|
|
// supplied, assume we're after a 2D transform
|
|
Matrix.Rotation = function(theta, a) {
|
|
if (!a) {
|
|
return Matrix.create([
|
|
[Math.cos(theta), -Math.sin(theta)],
|
|
[Math.sin(theta), Math.cos(theta)]
|
|
]);
|
|
}
|
|
var axis = a.dup();
|
|
if (axis.elements.length != 3) { return null; }
|
|
var mod = axis.modulus();
|
|
var x = axis.elements[0] / mod, y = axis.elements[1] / mod, z = axis.elements[2] / mod;
|
|
var s = Math.sin(theta), c = Math.cos(theta), t = 1 - c;
|
|
// Formula derived here: http://www.gamedev.net/reference/articles/article1199.asp
|
|
// That proof rotates the co-ordinate system so theta
|
|
// becomes -theta and sin becomes -sin here.
|
|
return Matrix.create([
|
|
[t * x * x + c, t * x * y - s * z, t * x * z + s * y],
|
|
[t * x * y + s * z, t * y * y + c, t * y * z - s * x],
|
|
[t * x * z - s * y, t * y * z + s * x, t * z * z + c]
|
|
]);
|
|
};
|
|
|
|
// Special case rotations
|
|
Matrix.RotationX = function(t) {
|
|
var c = Math.cos(t), s = Math.sin(t);
|
|
return Matrix.create([
|
|
[1, 0, 0],
|
|
[0, c, -s],
|
|
[0, s, c]
|
|
]);
|
|
};
|
|
|
|
Matrix.RotationY = function(t) {
|
|
var c = Math.cos(t), s = Math.sin(t);
|
|
return Matrix.create([
|
|
[c, 0, s],
|
|
[0, 1, 0],
|
|
[-s, 0, c]
|
|
]);
|
|
};
|
|
|
|
Matrix.RotationZ = function(t) {
|
|
var c = Math.cos(t), s = Math.sin(t);
|
|
return Matrix.create([
|
|
[c, -s, 0],
|
|
[s, c, 0],
|
|
[0, 0, 1]
|
|
]);
|
|
};
|
|
|
|
// Random matrix of n rows, m columns
|
|
Matrix.Random = function(n, m) {
|
|
if (arguments.length === 1) m = n;
|
|
return Matrix.Zero(n, m).map(
|
|
function() { return Math.random(); }
|
|
);
|
|
};
|
|
|
|
Matrix.Fill = function(n, m, v) {
|
|
if (arguments.length === 2) {
|
|
v = m;
|
|
m = n;
|
|
}
|
|
|
|
var els = [], i = n, j;
|
|
|
|
while (i--) {
|
|
j = m;
|
|
els[i] = [];
|
|
|
|
while (j--) {
|
|
els[i][j] = v;
|
|
}
|
|
}
|
|
|
|
return Matrix.create(els);
|
|
};
|
|
|
|
// Matrix filled with zeros
|
|
Matrix.Zero = function(n, m) {
|
|
return Matrix.Fill(n, m, 0);
|
|
};
|
|
|
|
// Matrix filled with zeros
|
|
Matrix.Zeros = function(n, m) {
|
|
return Matrix.Zero(n, m);
|
|
};
|
|
|
|
// Matrix filled with ones
|
|
Matrix.One = function(n, m) {
|
|
return Matrix.Fill(n, m, 1);
|
|
};
|
|
|
|
// Matrix filled with ones
|
|
Matrix.Ones = function(n, m) {
|
|
return Matrix.One(n, m);
|
|
};
|
|
|
|
module.exports = Matrix;
|