// Generated by CoffeeScript 1.10.0
/**
* API for dealing with matrices as simple 2D arrays.
* @module arrays
*/
/*
* Internal reference to module.exports.
* @private
*/
var combine, diagonalProduct, isEven, map, solve, t;
t = this;
/**
* @param {!Integer} size
* @return {Array.<Array.<Number>>} The identity 2D array of given size.
*/
this.createIdentity = function(size) {
var i, idnetity, j, m, n, ref, ref1;
if (size === 0) {
return [[]];
}
idnetity = new Array(size);
for (i = m = 0, ref = size; 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
idnetity[i] = new Array(size);
for (j = n = 0, ref1 = size; 0 <= ref1 ? n < ref1 : n > ref1; j = 0 <= ref1 ? ++n : --n) {
idnetity[i][j] = i === j ? 1 : 0;
}
}
return idnetity;
};
/**
* @param {!Integer} numOfRows
* @param {Integer} [numOfCols=numOfRows]
* @return {Array.<Array.<Number>>} The blank (zeros) 2D array of given size.
*/
this.createBlank = function(numOfRows, numOfCols) {
var blank, i, j, m, n, ref, ref1;
if (numOfCols == null) {
numOfCols = numOfRows;
}
if (numOfRows === 0 || numOfCols === 0) {
return [[]];
}
blank = new Array(numOfRows);
for (i = m = 0, ref = numOfRows; 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
blank[i] = new Array(numOfCols);
for (j = n = 0, ref1 = numOfCols; 0 <= ref1 ? n < ref1 : n > ref1; j = 0 <= ref1 ? ++n : --n) {
blank[i][j] = 0;
}
}
return blank;
};
/**
* @param {!Array.<Array.<Number>>} arrays - 2D array to copy.
* @return {Array.<Array.<Number>>}
*/
this.copy = function(arrays) {
var e, i, j, len, len1, m, n, newMatrix, row;
newMatrix = new Array(t.getNumOfRows(arrays));
for (i = m = 0, len = arrays.length; m < len; i = ++m) {
row = arrays[i];
newMatrix[i] = new Array(t.getNumOfColumns(arrays));
for (j = n = 0, len1 = row.length; n < len1; j = ++n) {
e = row[j];
newMatrix[i][j] = e;
}
}
return newMatrix;
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Boolean}
*/
this.isEmpty = function(arrays) {
if (arrays.length === 0 || arrays[0].length === 0) {
return true;
} else {
return false;
}
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Boolean}
*/
this.isSquare = function(arrays) {
if (t.isEmpty(arrays)) {
return false;
} else {
return t.getNumOfRows(arrays) === t.getNumOfColumns(arrays);
}
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Integer}
*/
this.getNumOfColumns = function(arrays) {
if (arrays.length === 0) {
return 0;
} else {
return arrays[0].length;
}
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Integer}
*/
this.getNumOfRows = function(arrays) {
if (t.isEmpty(arrays)) {
return 0;
} else if (arrays[0].length > 0) {
return arrays.length;
}
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Array.<Integer>} e.g. [2, 3] for a 2x3 matrix
*/
this.getDimensions = function(arrays) {
return [t.getNumOfRows(arrays), t.getNumOfColumns(arrays)];
};
/**
* Get a string representation of the 2D array size.
* @param {!Array.<Array.<Number>>} arrays
* @return {String}
*/
this.getSize = function(arrays) {
var numOfCols, numOfRows, ref;
ref = t.getDimensions(arrays), numOfRows = ref[0], numOfCols = ref[1];
return numOfRows + "x" + numOfCols;
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Boolean}
*/
this.isLowerTriangular = function(arrays) {
var i, j, m, n, numOfColumns, numOfRows, ref, ref1, ref2, ref3;
ref = t.getDimensions(arrays), numOfRows = ref[0], numOfColumns = ref[1];
if (numOfRows === 0 || numOfColumns === 0) {
return false;
}
if (numOfColumns === 1) {
return true;
}
for (i = m = 0, ref1 = numOfRows; 0 <= ref1 ? m < ref1 : m > ref1; i = 0 <= ref1 ? ++m : --m) {
for (j = n = ref2 = i + 1, ref3 = numOfColumns; ref2 <= ref3 ? n < ref3 : n > ref3; j = ref2 <= ref3 ? ++n : --n) {
if (arrays[i][j] !== 0) {
return false;
}
}
}
return true;
};
/**
* @param {!Array.<Array.<Number>>} arrays
* @return {Boolean}
*/
this.isUpperTriangular = function(arrays) {
var i, j, m, n, numOfColumns, numOfRows, ref, ref1, ref2, rowToCountTo;
ref = t.getDimensions(arrays), numOfRows = ref[0], numOfColumns = ref[1];
if (numOfRows === 0 || numOfColumns === 0) {
return false;
}
if (numOfRows === 1) {
return true;
}
for (i = m = 1, ref1 = numOfRows; 1 <= ref1 ? m < ref1 : m > ref1; i = 1 <= ref1 ? ++m : --m) {
rowToCountTo = Math.min(i - 1, numOfColumns - 1);
for (j = n = 0, ref2 = rowToCountTo; 0 <= ref2 ? n <= ref2 : n >= ref2; j = 0 <= ref2 ? ++n : --n) {
if (arrays[i][j] !== 0) {
return false;
}
}
}
return true;
};
/*
* Create a new 2D array from a combination of two.
* @private
* @param {!Array.<Array.<Number>>} a1
* @param {!Array.<Array.<Number>>} a2
* @param {function(number, number): number} f
* @return {Array.<Array.<Number>>}
*/
combine = function(a1, a2, f) {
var i, j, m, n, r, ref, ref1;
r = new Array(t.getNumOfRows(a1));
for (i = m = 0, ref = t.getNumOfRows(a1); 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
r[i] = new Array(t.getNumOfColumns(a1));
for (j = n = 0, ref1 = t.getNumOfColumns(a1); 0 <= ref1 ? n < ref1 : n > ref1; j = 0 <= ref1 ? ++n : --n) {
r[i][j] = f(a1[i][j], a2[i][j]);
}
}
return r;
};
/**
* Add two 2D arrays and return the result.
* @param {Array.<Array.<Number>>} a1
* @param {Array.<Array.<Number>>} a2
* @return {Array.<Array.<Number>>}
*/
this.add = function(a1, a2) {
return combine(a1, a2, function(n1, n2) {
return n1 + n2;
});
};
/**
* Subtract a2 from a1 and return the result.
* @param {Array.<Array.<Number>>} a1
* @param {Array.<Array.<Number>>} a2
* @return {Array.<Array.<Number>>}
*/
this.subtract = function(a1, a2) {
return combine(a1, a2, function(n1, n2) {
return n1 - n2;
});
};
/**
* Multiply two 2D arrays and return the result.
* @param {(Array.<Array.<Number>>|Number)} a1
* @param {(Array.<Array.<Number>>|Number)} a2
* @return {Array.<Array.<Number>>}
*/
this.multiply = function(a1, a2) {
var i, j, k, m, n, o, r, ref, ref1, ref2;
if (typeof a1 === 'number' && typeof a2 === 'number') {
return [[a1 * a2]];
}
if (typeof a1 === 'number') {
return map((function(e) {
return e * a1;
}), a2);
}
if (typeof a2 === 'number') {
return map((function(e) {
return e * a2;
}), a1);
}
r = new Array(t.getNumOfRows(a1));
for (i = m = 0, ref = t.getNumOfRows(a1); 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
r[i] = new Array(t.getNumOfColumns(a2));
for (j = n = 0, ref1 = t.getNumOfColumns(a2); 0 <= ref1 ? n < ref1 : n > ref1; j = 0 <= ref1 ? ++n : --n) {
r[i][j] = 0;
for (k = o = 0, ref2 = t.getNumOfColumns(a1); 0 <= ref2 ? o < ref2 : o > ref2; k = 0 <= ref2 ? ++o : --o) {
r[i][j] += a1[i][k] * a2[k][j];
}
}
}
return r;
};
/**
* Divide a1 by a2 and return the result.
* @param {(Array.<Array.<Number>>|Number)} a1
* @param {(Array.<Array.<Number>>|Number)} a2
* @return {Array.<Array.<Number>>}
*/
this.divide = function(a1, a2) {
return t.multiply(a1, t.invert(a2));
};
/**
* @typedef module:arrays.LUP
* @type {object}
* @property {Array.<Array.<Number>>} l - Lower triangular matrix.
* @property {Array.<Array.<Number>>} u - Upper triangular matrix.
* @property {Array.<Array.<Number>>} p - Permutation matrix.
*/
/**
* Decompose this 2D array into lower and upper triangular 2D arrays.
* @param {Array.<Array.<Number>>} arrays
* @throws {SingularMatrixException}
* @return {LUP} Lower and upper triangular 2D arrays, and a permutation 2D array.
*/
this.decompose = function(arrays) {
var e, first, i, j, k, l, m, maxFirstElement, n, numOfSwaps, o, p, pivotRowIndex, q, r, ref, ref1, ref10, ref11, ref12, ref2, ref3, ref4, ref5, ref6, ref7, ref8, ref9, s, size, u, v;
if (t.isEmpty(arrays)) {
return {
l: [[]],
u: [[]],
p: [[]]
};
}
if (!t.isSquare(arrays)) {
throw {
name: 'NotImplementedException',
message: 'LU Decomposition not implemented for non-square 2D arrays.'
};
}
size = t.getNumOfRows(arrays);
l = t.createIdentity(size);
u = t.copy(arrays);
p = t.createIdentity(size);
numOfSwaps = 0;
/*
Gaussian elimination w / partial pivoting
*/
for (i = m = 0, ref = size - 1; 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
pivotRowIndex = i;
maxFirstElement = u[i][i];
for (r = n = ref1 = i, ref2 = size; ref1 <= ref2 ? n < ref2 : n > ref2; r = ref1 <= ref2 ? ++n : --n) {
first = Math.abs(u[r][i]);
if (first > maxFirstElement) {
maxFirstElement = first;
pivotRowIndex = r;
}
}
if (maxFirstElement === 0) {
throw {
name: 'SingularMatrixException',
message: 'Singular matrix',
cause: arrays
};
}
if (pivotRowIndex !== i) {
numOfSwaps++;
for (e = o = ref3 = i, ref4 = size; ref3 <= ref4 ? o < ref4 : o > ref4; e = ref3 <= ref4 ? ++o : --o) {
ref5 = [u[pivotRowIndex][e], u[i][e]], u[i][e] = ref5[0], u[pivotRowIndex][e] = ref5[1];
}
for (e = q = 0, ref6 = i; 0 <= ref6 ? q < ref6 : q > ref6; e = 0 <= ref6 ? ++q : --q) {
ref7 = [l[pivotRowIndex][e], l[i][e]], l[i][e] = ref7[0], l[pivotRowIndex][e] = ref7[1];
}
ref8 = [p[pivotRowIndex], p[i]], p[i] = ref8[0], p[pivotRowIndex] = ref8[1];
}
for (j = s = ref9 = i + 1, ref10 = size; ref9 <= ref10 ? s < ref10 : s > ref10; j = ref9 <= ref10 ? ++s : --s) {
l[j][i] = u[j][i] / u[i][i];
for (k = v = ref11 = i, ref12 = size; ref11 <= ref12 ? v < ref12 : v > ref12; k = ref11 <= ref12 ? ++v : --v) {
u[j][k] -= l[j][i] * u[i][k];
}
}
}
return {
l: l,
u: u,
p: p,
numOfSwaps: numOfSwaps
};
};
/**
* Solve the matrix equation Ax = b for x with matrix size N.
* @param {Array.<Array.<Number>>} A - NxN
* @param {Array.<Array.<Number>>} b - Nx1
* @return {Array.<Array.<Number>>} Nx1
*/
this.solve = function(A, b) {
var l, p, ref, u;
ref = t.decompose(A), l = ref.l, u = ref.u, p = ref.p;
return solve({
l: l,
u: u,
p: p
}, b);
};
/**
* @private
* @param {LUP} lup
* @param {Array.<Array.<Number>>} b
* @return {Array.<Array.<Number>>}
*/
solve = function(arg, b) {
var colIndex, l, m, n, o, p, pb, q, ref, ref1, ref2, ref3, ref4, rowIndex, size, u, x, y;
l = arg.l, u = arg.u, p = arg.p;
size = t.getNumOfRows(b);
pb = t.multiply(p, b);
y = t.createBlank(size, 1);
/*
TODO optimise by treating x (and y) as a single array, then copy into a 1 column matrix?
*/
for (rowIndex = m = 0, ref = size; 0 <= ref ? m < ref : m > ref; rowIndex = 0 <= ref ? ++m : --m) {
y[rowIndex][0] = pb[rowIndex][0];
for (colIndex = n = 0, ref1 = rowIndex; 0 <= ref1 ? n < ref1 : n > ref1; colIndex = 0 <= ref1 ? ++n : --n) {
y[rowIndex][0] -= l[rowIndex][colIndex] * y[colIndex][0];
}
y[rowIndex][0] /= l[rowIndex][rowIndex];
}
x = t.createBlank(size, 1);
for (rowIndex = o = ref2 = size - 1; ref2 <= 0 ? o <= 0 : o >= 0; rowIndex = ref2 <= 0 ? ++o : --o) {
x[rowIndex][0] = y[rowIndex][0];
for (colIndex = q = ref3 = rowIndex + 1, ref4 = size; ref3 <= ref4 ? q < ref4 : q > ref4; colIndex = ref3 <= ref4 ? ++q : --q) {
x[rowIndex][0] -= u[rowIndex][colIndex] * x[colIndex][0];
}
x[rowIndex][0] /= u[rowIndex][rowIndex];
}
return x;
};
/**
* Return a transposed version of the given 2D array.
* @param {Array.<Array.<Number>>} arrays
* @return {Array.<Array.<Number>>}
*/
this.transpose = function(arrays) {
var i, j, m, n, ref, ref1, trans;
trans = new Array();
for (i = m = 0, ref = t.getNumOfColumns(arrays); 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
trans[i] = new Array();
for (j = n = 0, ref1 = t.getNumOfRows(arrays); 0 <= ref1 ? n < ref1 : n > ref1; j = 0 <= ref1 ? ++n : --n) {
trans[i][j] = arrays[j][i];
}
}
return trans;
};
/**
* Return an inverted version of the given 2D array.
* @param {Array.<Array.<Number>>} arrays
* @return {Array.<Array.<Number>>}
*/
this.invert = function(arrays) {
var columnOfIdentity, i, identity, inversion, j, l, m, n, p, ref, ref1, ref2, size, solution, u;
size = t.getNumOfRows(arrays);
identity = t.createIdentity(size);
ref = t.decompose(arrays), l = ref.l, u = ref.u, p = ref.p;
inversion = t.createBlank(size);
for (i = m = 0, ref1 = size; 0 <= ref1 ? m < ref1 : m > ref1; i = 0 <= ref1 ? ++m : --m) {
columnOfIdentity = t.transpose([identity[i]]);
solution = solve({
l: l,
u: u,
p: p
}, columnOfIdentity);
for (j = n = 0, ref2 = size; 0 <= ref2 ? n < ref2 : n > ref2; j = 0 <= ref2 ? ++n : --n) {
inversion[j][i] = solution[j][0];
}
}
return inversion;
};
/**
* Is it possible to invert the given 2D array?
* @param {Array.<Array.<Number>>} arrays
* @return {Boolean}
*/
this.isInvertible = function(arrays) {
return t.isSquare(arrays) && t.determinant(arrays) !== 0;
};
map = function(f, arrays) {
var i, j, m, n, r, ref, ref1;
r = new Array(t.getNumOfRows(arrays));
for (i = m = 0, ref = t.getNumOfRows(arrays); 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
r[i] = new Array(t.getNumOfColumns(arrays));
for (j = n = 0, ref1 = t.getNumOfColumns(arrays); 0 <= ref1 ? n < ref1 : n > ref1; j = 0 <= ref1 ? ++n : --n) {
r[i][j] = f(arrays[i][j]);
}
}
return r;
};
/**
* Return the determinant of the given 2D array.
* @param {Array.<Array.<Number>>} arrays - Square matrix
* @return {Number}
*/
this.determinant = function(arrays) {
var l, numOfSwaps, ref, u;
ref = t.decompose(arrays), l = ref.l, u = ref.u, numOfSwaps = ref.numOfSwaps;
return diagonalProduct(l) * diagonalProduct(u) * (isEven(numOfSwaps) ? 1 : -1);
};
diagonalProduct = function(triangularArrays) {
var i, m, r, ref, size;
size = t.getNumOfRows(triangularArrays);
r = 1;
for (i = m = 0, ref = size; 0 <= ref ? m < ref : m > ref; i = 0 <= ref ? ++m : --m) {
r *= triangularArrays[i][i];
}
return r;
};
isEven = function(number) {
return number % 2 === 0;
};