Source code
Revision control
Copy as Markdown
Other Tools
/*
* The MIT License (MIT)
*
* Copyright (c) 2015 ml.js
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
'use strict';
// ml-stat array.js
const MLStatArray = {};
{
function compareNumbers(a, b) {
return a - b;
}
/**
* Computes the sum of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.sum = function sum(values) {
var sum = 0;
for (var i = 0; i < values.length; i++) {
sum += values[i];
}
return sum;
};
/**
* Computes the maximum of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.max = function max(values) {
var max = values[0];
var l = values.length;
for (var i = 1; i < l; i++) {
if (values[i] > max) max = values[i];
}
return max;
};
/**
* Computes the minimum of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.min = function min(values) {
var min = values[0];
var l = values.length;
for (var i = 1; i < l; i++) {
if (values[i] < min) min = values[i];
}
return min;
};
/**
* Computes the min and max of the given values
* @param {Array} values
* @returns {{min: number, max: number}}
*/
MLStatArray.minMax = function minMax(values) {
var min = values[0];
var max = values[0];
var l = values.length;
for (var i = 1; i < l; i++) {
if (values[i] < min) min = values[i];
if (values[i] > max) max = values[i];
}
return {
min: min,
max: max
};
};
/**
* Computes the arithmetic mean of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.arithmeticMean = function arithmeticMean(values) {
var sum = 0;
var l = values.length;
for (var i = 0; i < l; i++) {
sum += values[i];
}
return sum / l;
};
/**
* {@link arithmeticMean}
*/
MLStatArray.mean = MLStatArray.arithmeticMean;
/**
* Computes the geometric mean of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.geometricMean = function geometricMean(values) {
var mul = 1;
var l = values.length;
for (var i = 0; i < l; i++) {
mul *= values[i];
}
return Math.pow(mul, 1 / l);
};
/**
* Computes the mean of the log of the given values
* If the return value is exponentiated, it gives the same result as the
* geometric mean.
* @param {Array} values
* @returns {number}
*/
MLStatArray.logMean = function logMean(values) {
var lnsum = 0;
var l = values.length;
for (var i = 0; i < l; i++) {
lnsum += Math.log(values[i]);
}
return lnsum / l;
};
/**
* Computes the weighted grand mean for a list of means and sample sizes
* @param {Array} means - Mean values for each set of samples
* @param {Array} samples - Number of original values for each set of samples
* @returns {number}
*/
MLStatArray.grandMean = function grandMean(means, samples) {
var sum = 0;
var n = 0;
var l = means.length;
for (var i = 0; i < l; i++) {
sum += samples[i] * means[i];
n += samples[i];
}
return sum / n;
};
/**
* Computes the truncated mean of the given values using a given percentage
* @param {Array} values
* @param {number} percent - The percentage of values to keep (range: [0,1])
* @param {boolean} [alreadySorted=false]
* @returns {number}
*/
MLStatArray.truncatedMean = function truncatedMean(values, percent, alreadySorted) {
if (alreadySorted === undefined) alreadySorted = false;
if (!alreadySorted) {
values = [].concat(values).sort(compareNumbers);
}
var l = values.length;
var k = Math.floor(l * percent);
var sum = 0;
for (var i = k; i < (l - k); i++) {
sum += values[i];
}
return sum / (l - 2 * k);
};
/**
* Computes the harmonic mean of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.harmonicMean = function harmonicMean(values) {
var sum = 0;
var l = values.length;
for (var i = 0; i < l; i++) {
if (values[i] === 0) {
throw new RangeError('value at index ' + i + 'is zero');
}
sum += 1 / values[i];
}
return l / sum;
};
/**
* Computes the contraharmonic mean of the given values
* @param {Array} values
* @returns {number}
*/
MLStatArray.contraHarmonicMean = function contraHarmonicMean(values) {
var r1 = 0;
var r2 = 0;
var l = values.length;
for (var i = 0; i < l; i++) {
r1 += values[i] * values[i];
r2 += values[i];
}
if (r2 < 0) {
throw new RangeError('sum of values is negative');
}
return r1 / r2;
};
/**
* Computes the median of the given values
* @param {Array} values
* @param {boolean} [alreadySorted=false]
* @returns {number}
*/
MLStatArray.median = function median(values, alreadySorted) {
if (alreadySorted === undefined) alreadySorted = false;
if (!alreadySorted) {
values = [].concat(values).sort(compareNumbers);
}
var l = values.length;
var half = Math.floor(l / 2);
if (l % 2 === 0) {
return (values[half - 1] + values[half]) * 0.5;
} else {
return values[half];
}
};
/**
* Computes the variance of the given values
* @param {Array} values
* @param {boolean} [unbiased=true] - if true, divide by (n-1); if false, divide by n.
* @returns {number}
*/
MLStatArray.variance = function variance(values, unbiased) {
if (unbiased === undefined) unbiased = true;
var theMean = MLStatArray.mean(values);
var theVariance = 0;
var l = values.length;
for (var i = 0; i < l; i++) {
var x = values[i] - theMean;
theVariance += x * x;
}
if (unbiased) {
return theVariance / (l - 1);
} else {
return theVariance / l;
}
};
/**
* Computes the standard deviation of the given values
* @param {Array} values
* @param {boolean} [unbiased=true] - if true, divide by (n-1); if false, divide by n.
* @returns {number}
*/
MLStatArray.standardDeviation = function standardDeviation(values, unbiased) {
return Math.sqrt(MLStatArray.variance(values, unbiased));
};
MLStatArray.standardError = function standardError(values) {
return MLStatArray.standardDeviation(values) / Math.sqrt(values.length);
};
/**
* IEEE Transactions on biomedical engineering, vol. 52, no. 1, january 2005, p. 76-
* Calculate the standard deviation via the Median of the absolute deviation
* The formula for the standard deviation only holds for Gaussian random variables.
* @returns {{mean: number, stdev: number}}
*/
MLStatArray.robustMeanAndStdev = function robustMeanAndStdev(y) {
var mean = 0, stdev = 0;
var length = y.length, i = 0;
for (i = 0; i < length; i++) {
mean += y[i];
}
mean /= length;
var averageDeviations = new Array(length);
for (i = 0; i < length; i++)
averageDeviations[i] = Math.abs(y[i] - mean);
averageDeviations.sort(compareNumbers);
if (length % 2 === 1) {
stdev = averageDeviations[(length - 1) / 2] / 0.6745;
} else {
stdev = 0.5 * (averageDeviations[length / 2] + averageDeviations[length / 2 - 1]) / 0.6745;
}
return {
mean: mean,
stdev: stdev
};
};
MLStatArray.quartiles = function quartiles(values, alreadySorted) {
if (typeof (alreadySorted) === 'undefined') alreadySorted = false;
if (!alreadySorted) {
values = [].concat(values).sort(compareNumbers);
}
var quart = values.length / 4;
var q1 = values[Math.ceil(quart) - 1];
var q2 = MLStatArray.median(values, true);
var q3 = values[Math.ceil(quart * 3) - 1];
return {q1: q1, q2: q2, q3: q3};
};
MLStatArray.pooledStandardDeviation = function pooledStandardDeviation(samples, unbiased) {
return Math.sqrt(MLStatArray.pooledVariance(samples, unbiased));
};
MLStatArray.pooledVariance = function pooledVariance(samples, unbiased) {
if (typeof (unbiased) === 'undefined') unbiased = true;
var sum = 0;
var length = 0, l = samples.length;
for (var i = 0; i < l; i++) {
var values = samples[i];
var vari = MLStatArray.variance(values);
sum += (values.length - 1) * vari;
if (unbiased)
length += values.length - 1;
else
length += values.length;
}
return sum / length;
};
MLStatArray.mode = function mode(values) {
var l = values.length,
itemCount = new Array(l),
i;
for (i = 0; i < l; i++) {
itemCount[i] = 0;
}
var itemArray = new Array(l);
var count = 0;
for (i = 0; i < l; i++) {
var index = itemArray.indexOf(values[i]);
if (index >= 0)
itemCount[index]++;
else {
itemArray[count] = values[i];
itemCount[count] = 1;
count++;
}
}
var maxValue = 0, maxIndex = 0;
for (i = 0; i < count; i++) {
if (itemCount[i] > maxValue) {
maxValue = itemCount[i];
maxIndex = i;
}
}
return itemArray[maxIndex];
};
MLStatArray.covariance = function covariance(vector1, vector2, unbiased) {
if (typeof (unbiased) === 'undefined') unbiased = true;
var mean1 = MLStatArray.mean(vector1);
var mean2 = MLStatArray.mean(vector2);
if (vector1.length !== vector2.length)
throw 'Vectors do not have the same dimensions';
var cov = 0, l = vector1.length;
for (var i = 0; i < l; i++) {
var x = vector1[i] - mean1;
var y = vector2[i] - mean2;
cov += x * y;
}
if (unbiased)
return cov / (l - 1);
else
return cov / l;
};
MLStatArray.skewness = function skewness(values, unbiased) {
if (typeof (unbiased) === 'undefined') unbiased = true;
var theMean = MLStatArray.mean(values);
var s2 = 0, s3 = 0, l = values.length;
for (var i = 0; i < l; i++) {
var dev = values[i] - theMean;
s2 += dev * dev;
s3 += dev * dev * dev;
}
var m2 = s2 / l;
var m3 = s3 / l;
var g = m3 / (Math.pow(m2, 3 / 2.0));
if (unbiased) {
var a = Math.sqrt(l * (l - 1));
var b = l - 2;
return (a / b) * g;
} else {
return g;
}
};
MLStatArray.kurtosis = function kurtosis(values, unbiased) {
if (typeof (unbiased) === 'undefined') unbiased = true;
var theMean = MLStatArray.mean(values);
var n = values.length, s2 = 0, s4 = 0;
for (var i = 0; i < n; i++) {
var dev = values[i] - theMean;
s2 += dev * dev;
s4 += dev * dev * dev * dev;
}
var m2 = s2 / n;
var m4 = s4 / n;
if (unbiased) {
var v = s2 / (n - 1);
var a = (n * (n + 1)) / ((n - 1) * (n - 2) * (n - 3));
var b = s4 / (v * v);
var c = ((n - 1) * (n - 1)) / ((n - 2) * (n - 3));
return a * b - 3 * c;
} else {
return m4 / (m2 * m2) - 3;
}
};
MLStatArray.entropy = function entropy(values, eps) {
if (typeof (eps) === 'undefined') eps = 0;
var sum = 0, l = values.length;
for (var i = 0; i < l; i++)
sum += values[i] * Math.log(values[i] + eps);
return -sum;
};
MLStatArray.weightedMean = function weightedMean(values, weights) {
var sum = 0, l = values.length;
for (var i = 0; i < l; i++)
sum += values[i] * weights[i];
return sum;
};
MLStatArray.weightedStandardDeviation = function weightedStandardDeviation(values, weights) {
return Math.sqrt(MLStatArray.weightedVariance(values, weights));
};
MLStatArray.weightedVariance = function weightedVariance(values, weights) {
var theMean = MLStatArray.weightedMean(values, weights);
var vari = 0, l = values.length;
var a = 0, b = 0;
for (var i = 0; i < l; i++) {
var z = values[i] - theMean;
var w = weights[i];
vari += w * (z * z);
b += w;
a += w * w;
}
return vari * (b / (b * b - a));
};
MLStatArray.center = function center(values, inPlace) {
if (typeof (inPlace) === 'undefined') inPlace = false;
var result = values;
if (!inPlace)
result = [].concat(values);
var theMean = MLStatArray.mean(result), l = result.length;
for (var i = 0; i < l; i++)
result[i] -= theMean;
};
MLStatArray.standardize = function standardize(values, standardDev, inPlace) {
if (typeof (standardDev) === 'undefined') standardDev = MLStatArray.standardDeviation(values);
if (typeof (inPlace) === 'undefined') inPlace = false;
var l = values.length;
var result = inPlace ? values : new Array(l);
for (var i = 0; i < l; i++)
result[i] = values[i] / standardDev;
return result;
};
MLStatArray.cumulativeSum = function cumulativeSum(array) {
var l = array.length;
var result = new Array(l);
result[0] = array[0];
for (var i = 1; i < l; i++)
result[i] = result[i - 1] + array[i];
return result;
};
}
// ml-stat matrix.js
const MLStatMatrix = {};
{
let arrayStat = MLStatArray;
function compareNumbers(a, b) {
return a - b;
}
MLStatMatrix.max = function max(matrix) {
var max = -Infinity;
for (var i = 0; i < matrix.length; i++) {
for (var j = 0; j < matrix[i].length; j++) {
if (matrix[i][j] > max) max = matrix[i][j];
}
}
return max;
};
MLStatMatrix.min = function min(matrix) {
var min = Infinity;
for (var i = 0; i < matrix.length; i++) {
for (var j = 0; j < matrix[i].length; j++) {
if (matrix[i][j] < min) min = matrix[i][j];
}
}
return min;
};
MLStatMatrix.minMax = function minMax(matrix) {
var min = Infinity;
var max = -Infinity;
for (var i = 0; i < matrix.length; i++) {
for (var j = 0; j < matrix[i].length; j++) {
if (matrix[i][j] < min) min = matrix[i][j];
if (matrix[i][j] > max) max = matrix[i][j];
}
}
return {
min:min,
max:max
};
};
MLStatMatrix.entropy = function entropy(matrix, eps) {
if (typeof (eps) === 'undefined') {
eps = 0;
}
var sum = 0,
l1 = matrix.length,
l2 = matrix[0].length;
for (var i = 0; i < l1; i++) {
for (var j = 0; j < l2; j++) {
sum += matrix[i][j] * Math.log(matrix[i][j] + eps);
}
}
return -sum;
};
MLStatMatrix.mean = function mean(matrix, dimension) {
if (typeof (dimension) === 'undefined') {
dimension = 0;
}
var rows = matrix.length,
cols = matrix[0].length,
theMean, N, i, j;
if (dimension === -1) {
theMean = [0];
N = rows * cols;
for (i = 0; i < rows; i++) {
for (j = 0; j < cols; j++) {
theMean[0] += matrix[i][j];
}
}
theMean[0] /= N;
} else if (dimension === 0) {
theMean = new Array(cols);
N = rows;
for (j = 0; j < cols; j++) {
theMean[j] = 0;
for (i = 0; i < rows; i++) {
theMean[j] += matrix[i][j];
}
theMean[j] /= N;
}
} else if (dimension === 1) {
theMean = new Array(rows);
N = cols;
for (j = 0; j < rows; j++) {
theMean[j] = 0;
for (i = 0; i < cols; i++) {
theMean[j] += matrix[j][i];
}
theMean[j] /= N;
}
} else {
throw new Error('Invalid dimension');
}
return theMean;
};
MLStatMatrix.sum = function sum(matrix, dimension) {
if (typeof (dimension) === 'undefined') {
dimension = 0;
}
var rows = matrix.length,
cols = matrix[0].length,
theSum, i, j;
if (dimension === -1) {
theSum = [0];
for (i = 0; i < rows; i++) {
for (j = 0; j < cols; j++) {
theSum[0] += matrix[i][j];
}
}
} else if (dimension === 0) {
theSum = new Array(cols);
for (j = 0; j < cols; j++) {
theSum[j] = 0;
for (i = 0; i < rows; i++) {
theSum[j] += matrix[i][j];
}
}
} else if (dimension === 1) {
theSum = new Array(rows);
for (j = 0; j < rows; j++) {
theSum[j] = 0;
for (i = 0; i < cols; i++) {
theSum[j] += matrix[j][i];
}
}
} else {
throw new Error('Invalid dimension');
}
return theSum;
};
MLStatMatrix.product = function product(matrix, dimension) {
if (typeof (dimension) === 'undefined') {
dimension = 0;
}
var rows = matrix.length,
cols = matrix[0].length,
theProduct, i, j;
if (dimension === -1) {
theProduct = [1];
for (i = 0; i < rows; i++) {
for (j = 0; j < cols; j++) {
theProduct[0] *= matrix[i][j];
}
}
} else if (dimension === 0) {
theProduct = new Array(cols);
for (j = 0; j < cols; j++) {
theProduct[j] = 1;
for (i = 0; i < rows; i++) {
theProduct[j] *= matrix[i][j];
}
}
} else if (dimension === 1) {
theProduct = new Array(rows);
for (j = 0; j < rows; j++) {
theProduct[j] = 1;
for (i = 0; i < cols; i++) {
theProduct[j] *= matrix[j][i];
}
}
} else {
throw new Error('Invalid dimension');
}
return theProduct;
};
MLStatMatrix.standardDeviation = function standardDeviation(matrix, means, unbiased) {
var vari = MLStatMatrix.variance(matrix, means, unbiased), l = vari.length;
for (var i = 0; i < l; i++) {
vari[i] = Math.sqrt(vari[i]);
}
return vari;
};
MLStatMatrix.variance = function variance(matrix, means, unbiased) {
if (typeof (unbiased) === 'undefined') {
unbiased = true;
}
means = means || MLStatMatrix.mean(matrix);
var rows = matrix.length;
if (rows === 0) return [];
var cols = matrix[0].length;
var vari = new Array(cols);
for (var j = 0; j < cols; j++) {
var sum1 = 0, sum2 = 0, x = 0;
for (var i = 0; i < rows; i++) {
x = matrix[i][j] - means[j];
sum1 += x;
sum2 += x * x;
}
if (unbiased) {
vari[j] = (sum2 - ((sum1 * sum1) / rows)) / (rows - 1);
} else {
vari[j] = (sum2 - ((sum1 * sum1) / rows)) / rows;
}
}
return vari;
};
MLStatMatrix.median = function median(matrix) {
var rows = matrix.length, cols = matrix[0].length;
var medians = new Array(cols);
for (var i = 0; i < cols; i++) {
var data = new Array(rows);
for (var j = 0; j < rows; j++) {
data[j] = matrix[j][i];
}
data.sort(compareNumbers);
var N = data.length;
if (N % 2 === 0) {
medians[i] = (data[N / 2] + data[(N / 2) - 1]) * 0.5;
} else {
medians[i] = data[Math.floor(N / 2)];
}
}
return medians;
};
MLStatMatrix.mode = function mode(matrix) {
var rows = matrix.length,
cols = matrix[0].length,
modes = new Array(cols),
i, j;
for (i = 0; i < cols; i++) {
var itemCount = new Array(rows);
for (var k = 0; k < rows; k++) {
itemCount[k] = 0;
}
var itemArray = new Array(rows);
var count = 0;
for (j = 0; j < rows; j++) {
var index = itemArray.indexOf(matrix[j][i]);
if (index >= 0) {
itemCount[index]++;
} else {
itemArray[count] = matrix[j][i];
itemCount[count] = 1;
count++;
}
}
var maxValue = 0, maxIndex = 0;
for (j = 0; j < count; j++) {
if (itemCount[j] > maxValue) {
maxValue = itemCount[j];
maxIndex = j;
}
}
modes[i] = itemArray[maxIndex];
}
return modes;
};
MLStatMatrix.skewness = function skewness(matrix, unbiased) {
if (typeof (unbiased) === 'undefined') unbiased = true;
var means = MLStatMatrix.mean(matrix);
var n = matrix.length, l = means.length;
var skew = new Array(l);
for (var j = 0; j < l; j++) {
var s2 = 0, s3 = 0;
for (var i = 0; i < n; i++) {
var dev = matrix[i][j] - means[j];
s2 += dev * dev;
s3 += dev * dev * dev;
}
var m2 = s2 / n;
var m3 = s3 / n;
var g = m3 / Math.pow(m2, 3 / 2);
if (unbiased) {
var a = Math.sqrt(n * (n - 1));
var b = n - 2;
skew[j] = (a / b) * g;
} else {
skew[j] = g;
}
}
return skew;
};
MLStatMatrix.kurtosis = function kurtosis(matrix, unbiased) {
if (typeof (unbiased) === 'undefined') unbiased = true;
var means = MLStatMatrix.mean(matrix);
var n = matrix.length, m = matrix[0].length;
var kurt = new Array(m);
for (var j = 0; j < m; j++) {
var s2 = 0, s4 = 0;
for (var i = 0; i < n; i++) {
var dev = matrix[i][j] - means[j];
s2 += dev * dev;
s4 += dev * dev * dev * dev;
}
var m2 = s2 / n;
var m4 = s4 / n;
if (unbiased) {
var v = s2 / (n - 1);
var a = (n * (n + 1)) / ((n - 1) * (n - 2) * (n - 3));
var b = s4 / (v * v);
var c = ((n - 1) * (n - 1)) / ((n - 2) * (n - 3));
kurt[j] = a * b - 3 * c;
} else {
kurt[j] = m4 / (m2 * m2) - 3;
}
}
return kurt;
};
MLStatMatrix.standardError = function standardError(matrix) {
var samples = matrix.length;
var standardDeviations = MLStatMatrix.standardDeviation(matrix);
var l = standardDeviations.length;
var standardErrors = new Array(l);
var sqrtN = Math.sqrt(samples);
for (var i = 0; i < l; i++) {
standardErrors[i] = standardDeviations[i] / sqrtN;
}
return standardErrors;
};
MLStatMatrix.covariance = function covariance(matrix, dimension) {
return MLStatMatrix.scatter(matrix, undefined, dimension);
};
MLStatMatrix.scatter = function scatter(matrix, divisor, dimension) {
if (typeof (dimension) === 'undefined') {
dimension = 0;
}
if (typeof (divisor) === 'undefined') {
if (dimension === 0) {
divisor = matrix.length - 1;
} else if (dimension === 1) {
divisor = matrix[0].length - 1;
}
}
var means = MLStatMatrix.mean(matrix, dimension);
var rows = matrix.length;
if (rows === 0) {
return [[]];
}
var cols = matrix[0].length,
cov, i, j, s, k;
if (dimension === 0) {
cov = new Array(cols);
for (i = 0; i < cols; i++) {
cov[i] = new Array(cols);
}
for (i = 0; i < cols; i++) {
for (j = i; j < cols; j++) {
s = 0;
for (k = 0; k < rows; k++) {
s += (matrix[k][j] - means[j]) * (matrix[k][i] - means[i]);
}
s /= divisor;
cov[i][j] = s;
cov[j][i] = s;
}
}
} else if (dimension === 1) {
cov = new Array(rows);
for (i = 0; i < rows; i++) {
cov[i] = new Array(rows);
}
for (i = 0; i < rows; i++) {
for (j = i; j < rows; j++) {
s = 0;
for (k = 0; k < cols; k++) {
s += (matrix[j][k] - means[j]) * (matrix[i][k] - means[i]);
}
s /= divisor;
cov[i][j] = s;
cov[j][i] = s;
}
}
} else {
throw new Error('Invalid dimension');
}
return cov;
};
MLStatMatrix.correlation = function correlation(matrix) {
var means = MLStatMatrix.mean(matrix),
standardDeviations = MLStatMatrix.standardDeviation(matrix, true, means),
scores = MLStatMatrix.zScores(matrix, means, standardDeviations),
rows = matrix.length,
cols = matrix[0].length,
i, j;
var cor = new Array(cols);
for (i = 0; i < cols; i++) {
cor[i] = new Array(cols);
}
for (i = 0; i < cols; i++) {
for (j = i; j < cols; j++) {
var c = 0;
for (var k = 0, l = scores.length; k < l; k++) {
c += scores[k][j] * scores[k][i];
}
c /= rows - 1;
cor[i][j] = c;
cor[j][i] = c;
}
}
return cor;
};
MLStatMatrix.zScores = function zScores(matrix, means, standardDeviations) {
means = means || MLStatMatrix.mean(matrix);
if (typeof (standardDeviations) === 'undefined') standardDeviations = MLStatMatrix.standardDeviation(matrix, true, means);
return MLStatMatrix.standardize(MLStatMatrix.center(matrix, means, false), standardDeviations, true);
};
MLStatMatrix.center = function center(matrix, means, inPlace) {
means = means || MLStatMatrix.mean(matrix);
var result = matrix,
l = matrix.length,
i, j, jj;
if (!inPlace) {
result = new Array(l);
for (i = 0; i < l; i++) {
result[i] = new Array(matrix[i].length);
}
}
for (i = 0; i < l; i++) {
var row = result[i];
for (j = 0, jj = row.length; j < jj; j++) {
row[j] = matrix[i][j] - means[j];
}
}
return result;
};
MLStatMatrix.standardize = function standardize(matrix, standardDeviations, inPlace) {
if (typeof (standardDeviations) === 'undefined') standardDeviations = MLStatMatrix.standardDeviation(matrix);
var result = matrix,
l = matrix.length,
i, j, jj;
if (!inPlace) {
result = new Array(l);
for (i = 0; i < l; i++) {
result[i] = new Array(matrix[i].length);
}
}
for (i = 0; i < l; i++) {
var resultRow = result[i];
var sourceRow = matrix[i];
for (j = 0, jj = resultRow.length; j < jj; j++) {
if (standardDeviations[j] !== 0 && !isNaN(standardDeviations[j])) {
resultRow[j] = sourceRow[j] / standardDeviations[j];
}
}
}
return result;
};
MLStatMatrix.weightedVariance = function weightedVariance(matrix, weights) {
var means = MLStatMatrix.mean(matrix);
var rows = matrix.length;
if (rows === 0) return [];
var cols = matrix[0].length;
var vari = new Array(cols);
for (var j = 0; j < cols; j++) {
var sum = 0;
var a = 0, b = 0;
for (var i = 0; i < rows; i++) {
var z = matrix[i][j] - means[j];
var w = weights[i];
sum += w * (z * z);
b += w;
a += w * w;
}
vari[j] = sum * (b / (b * b - a));
}
return vari;
};
MLStatMatrix.weightedMean = function weightedMean(matrix, weights, dimension) {
if (typeof (dimension) === 'undefined') {
dimension = 0;
}
var rows = matrix.length;
if (rows === 0) return [];
var cols = matrix[0].length,
means, i, ii, j, w, row;
if (dimension === 0) {
means = new Array(cols);
for (i = 0; i < cols; i++) {
means[i] = 0;
}
for (i = 0; i < rows; i++) {
row = matrix[i];
w = weights[i];
for (j = 0; j < cols; j++) {
means[j] += row[j] * w;
}
}
} else if (dimension === 1) {
means = new Array(rows);
for (i = 0; i < rows; i++) {
means[i] = 0;
}
for (j = 0; j < rows; j++) {
row = matrix[j];
w = weights[j];
for (i = 0; i < cols; i++) {
means[j] += row[i] * w;
}
}
} else {
throw new Error('Invalid dimension');
}
var weightSum = arrayStat.sum(weights);
if (weightSum !== 0) {
for (i = 0, ii = means.length; i < ii; i++) {
means[i] /= weightSum;
}
}
return means;
};
MLStatMatrix.weightedCovariance = function weightedCovariance(matrix, weights, means, dimension) {
dimension = dimension || 0;
means = means || MLStatMatrix.weightedMean(matrix, weights, dimension);
var s1 = 0, s2 = 0;
for (var i = 0, ii = weights.length; i < ii; i++) {
s1 += weights[i];
s2 += weights[i] * weights[i];
}
var factor = s1 / (s1 * s1 - s2);
return MLStatMatrix.weightedScatter(matrix, weights, means, factor, dimension);
};
MLStatMatrix.weightedScatter = function weightedScatter(matrix, weights, means, factor, dimension) {
dimension = dimension || 0;
means = means || MLStatMatrix.weightedMean(matrix, weights, dimension);
if (typeof (factor) === 'undefined') {
factor = 1;
}
var rows = matrix.length;
if (rows === 0) {
return [[]];
}
var cols = matrix[0].length,
cov, i, j, k, s;
if (dimension === 0) {
cov = new Array(cols);
for (i = 0; i < cols; i++) {
cov[i] = new Array(cols);
}
for (i = 0; i < cols; i++) {
for (j = i; j < cols; j++) {
s = 0;
for (k = 0; k < rows; k++) {
s += weights[k] * (matrix[k][j] - means[j]) * (matrix[k][i] - means[i]);
}
cov[i][j] = s * factor;
cov[j][i] = s * factor;
}
}
} else if (dimension === 1) {
cov = new Array(rows);
for (i = 0; i < rows; i++) {
cov[i] = new Array(rows);
}
for (i = 0; i < rows; i++) {
for (j = i; j < rows; j++) {
s = 0;
for (k = 0; k < cols; k++) {
s += weights[k] * (matrix[j][k] - means[j]) * (matrix[i][k] - means[i]);
}
cov[i][j] = s * factor;
cov[j][i] = s * factor;
}
}
} else {
throw new Error('Invalid dimension');
}
return cov;
};
}
// ml-stat index.js
const MLStat = {};
{
MLStat.array = MLStatArray;
MLStat.matrix = MLStatMatrix;
}
// ml-array-utils ArrayUtils.js
const MLArrayUtilsArrayUtils = {};
{
const Stat = MLStat.array;
/**
* Function that returns an array of points given 1D array as follows:
*
* [x1, y1, .. , x2, y2, ..]
*
* And receive the number of dimensions of each point.
* @param array
* @param dimensions
* @returns {Array} - Array of points.
*/
function coordArrayToPoints(array, dimensions) {
if(array.length % dimensions !== 0) {
throw new RangeError('Dimensions number must be accordance with the size of the array.');
}
var length = array.length / dimensions;
var pointsArr = new Array(length);
var k = 0;
for(var i = 0; i < array.length; i += dimensions) {
var point = new Array(dimensions);
for(var j = 0; j < dimensions; ++j) {
point[j] = array[i + j];
}
pointsArr[k] = point;
k++;
}
return pointsArr;
}
/**
* Function that given an array as follows:
* [x1, y1, .. , x2, y2, ..]
*
* Returns an array as follows:
* [[x1, x2, ..], [y1, y2, ..], [ .. ]]
*
* And receives the number of dimensions of each coordinate.
* @param array
* @param dimensions
* @returns {Array} - Matrix of coordinates
*/
function coordArrayToCoordMatrix(array, dimensions) {
if(array.length % dimensions !== 0) {
throw new RangeError('Dimensions number must be accordance with the size of the array.');
}
var coordinatesArray = new Array(dimensions);
var points = array.length / dimensions;
for (var i = 0; i < coordinatesArray.length; i++) {
coordinatesArray[i] = new Array(points);
}
for(i = 0; i < array.length; i += dimensions) {
for(var j = 0; j < dimensions; ++j) {
var currentPoint = Math.floor(i / dimensions);
coordinatesArray[j][currentPoint] = array[i + j];
}
}
return coordinatesArray;
}
/**
* Function that receives a coordinate matrix as follows:
* [[x1, x2, ..], [y1, y2, ..], [ .. ]]
*
* Returns an array of coordinates as follows:
* [x1, y1, .. , x2, y2, ..]
*
* @param coordMatrix
* @returns {Array}
*/
function coordMatrixToCoordArray(coordMatrix) {
var coodinatesArray = new Array(coordMatrix.length * coordMatrix[0].length);
var k = 0;
for(var i = 0; i < coordMatrix[0].length; ++i) {
for(var j = 0; j < coordMatrix.length; ++j) {
coodinatesArray[k] = coordMatrix[j][i];
++k;
}
}
return coodinatesArray;
}
/**
* Tranpose a matrix, this method is for coordMatrixToPoints and
* pointsToCoordMatrix, that because only transposing the matrix
* you can change your representation.
*
* @param matrix
* @returns {Array}
*/
function transpose(matrix) {
var resultMatrix = new Array(matrix[0].length);
for(var i = 0; i < resultMatrix.length; ++i) {
resultMatrix[i] = new Array(matrix.length);
}
for (i = 0; i < matrix.length; ++i) {
for(var j = 0; j < matrix[0].length; ++j) {
resultMatrix[j][i] = matrix[i][j];
}
}
return resultMatrix;
}
/**
* Function that transform an array of points into a coordinates array
* as follows:
* [x1, y1, .. , x2, y2, ..]
*
* @param points
* @returns {Array}
*/
function pointsToCoordArray(points) {
var coodinatesArray = new Array(points.length * points[0].length);
var k = 0;
for(var i = 0; i < points.length; ++i) {
for(var j = 0; j < points[0].length; ++j) {
coodinatesArray[k] = points[i][j];
++k;
}
}
return coodinatesArray;
}
/**
* Apply the dot product between the smaller vector and a subsets of the
* largest one.
*
* @param firstVector
* @param secondVector
* @returns {Array} each dot product of size of the difference between the
* larger and the smallest one.
*/
function applyDotProduct(firstVector, secondVector) {
var largestVector, smallestVector;
if(firstVector.length <= secondVector.length) {
smallestVector = firstVector;
largestVector = secondVector;
} else {
smallestVector = secondVector;
largestVector = firstVector;
}
var difference = largestVector.length - smallestVector.length + 1;
var dotProductApplied = new Array(difference);
for (var i = 0; i < difference; ++i) {
var sum = 0;
for (var j = 0; j < smallestVector.length; ++j) {
sum += smallestVector[j] * largestVector[i + j];
}
dotProductApplied[i] = sum;
}
return dotProductApplied;
}
/**
* To scale the input array between the specified min and max values. The operation is performed inplace
* if the options.inplace is specified. If only one of the min or max parameters is specified, then the scaling
* will multiply the input array by min/min(input) or max/max(input)
* @param input
* @param options
* @returns {*}
*/
function scale(input, options){
var y;
if(options.inPlace){
y = input;
}
else{
y = new Array(input.length);
}
const max = options.max;
const min = options.min;
if(typeof max === "number"){
if(typeof min === "number"){
var minMax = Stat.minMax(input);
var factor = (max - min)/(minMax.max-minMax.min);
for(var i=0;i< y.length;i++){
y[i]=(input[i]-minMax.min)*factor+min;
}
}
else{
var currentMin = Stat.max(input);
var factor = max/currentMin;
for(var i=0;i< y.length;i++){
y[i] = input[i]*factor;
}
}
}
else{
if(typeof min === "number"){
var currentMin = Stat.min(input);
var factor = min/currentMin;
for(var i=0;i< y.length;i++){
y[i] = input[i]*factor;
}
}
}
return y;
}
MLArrayUtilsArrayUtils.coordArrayToPoints = coordArrayToPoints;
MLArrayUtilsArrayUtils.coordArrayToCoordMatrix = coordArrayToCoordMatrix;
MLArrayUtilsArrayUtils.coordMatrixToCoordArray = coordMatrixToCoordArray;
MLArrayUtilsArrayUtils.coordMatrixToPoints = transpose;
MLArrayUtilsArrayUtils.pointsToCoordArray = pointsToCoordArray;
MLArrayUtilsArrayUtils.pointsToCoordMatrix = transpose;
MLArrayUtilsArrayUtils.applyDotProduct = applyDotProduct;
MLArrayUtilsArrayUtils.scale = scale;
}
// ml-array-utils getEquallySpaced.js
const MLArrayUtilsGetEquallySpaced = {};
{
/**
*
* Function that returns a Number array of equally spaced numberOfPoints
* containing a representation of intensities of the spectra arguments x
* and y.
*
* The options parameter contains an object in the following form:
* from: starting point
* to: last point
* numberOfPoints: number of points between from and to
* variant: "slot" or "smooth" - smooth is the default option
*
* The slot variant consist that each point in the new array is calculated
* averaging the existing points between the slot that belongs to the current
* value. The smooth variant is the same but takes the integral of the range
* of the slot and divide by the step size between two points in the new array.
*
* @param x - sorted increasing x values
* @param y
* @param options
* @returns {Array} new array with the equally spaced data.
*
*/
function getEquallySpacedData(x, y, options) {
if (x.length>1 && x[0]>x[1]) {
x=x.slice().reverse();
y=y.slice().reverse();
}
var xLength = x.length;
if(xLength !== y.length)
throw new RangeError("the x and y vector doesn't have the same size.");
if (options === undefined) options = {};
var from = options.from === undefined ? x[0] : options.from
if (isNaN(from) || !isFinite(from)) {
throw new RangeError("'From' value must be a number");
}
var to = options.to === undefined ? x[x.length - 1] : options.to;
if (isNaN(to) || !isFinite(to)) {
throw new RangeError("'To' value must be a number");
}
var reverse = from > to;
if(reverse) {
var temp = from;
from = to;
to = temp;
}
var numberOfPoints = options.numberOfPoints === undefined ? 100 : options.numberOfPoints;
if (isNaN(numberOfPoints) || !isFinite(numberOfPoints)) {
throw new RangeError("'Number of points' value must be a number");
}
if(numberOfPoints < 1)
throw new RangeError("the number of point must be higher than 1");
var algorithm = options.variant === "slot" ? "slot" : "smooth"; // default value: smooth
var output = algorithm === "slot" ? getEquallySpacedSlot(x, y, from, to, numberOfPoints) : getEquallySpacedSmooth(x, y, from, to, numberOfPoints);
return reverse ? output.reverse() : output;
}
/**
* function that retrieves the getEquallySpacedData with the variant "smooth"
*
* @param x
* @param y
* @param from - Initial point
* @param to - Final point
* @param numberOfPoints
* @returns {Array} - Array of y's equally spaced with the variant "smooth"
*/
function getEquallySpacedSmooth(x, y, from, to, numberOfPoints) {
var xLength = x.length;
var step = (to - from) / (numberOfPoints - 1);
var halfStep = step / 2;
var start = from - halfStep;
var output = new Array(numberOfPoints);
var initialOriginalStep = x[1] - x[0];
var lastOriginalStep = x[x.length - 1] - x[x.length - 2];
// Init main variables
var min = start;
var max = start + step;
var previousX = Number.MIN_VALUE;
var previousY = 0;
var nextX = x[0] - initialOriginalStep;
var nextY = 0;
var currentValue = 0;
var slope = 0;
var intercept = 0;
var sumAtMin = 0;
var sumAtMax = 0;
var i = 0; // index of input
var j = 0; // index of output
function getSlope(x0, y0, x1, y1) {
return (y1 - y0) / (x1 - x0);
}
main: while(true) {
while (nextX - max >= 0) {
// no overlap with original point, just consume current value
var add = integral(0, max - previousX, slope, previousY);
sumAtMax = currentValue + add;
output[j] = (sumAtMax - sumAtMin) / step;
j++;
if (j === numberOfPoints)
break main;
min = max;
max += step;
sumAtMin = sumAtMax;
}
if(previousX <= min && min <= nextX) {
add = integral(0, min - previousX, slope, previousY);
sumAtMin = currentValue + add;
}
currentValue += integral(previousX, nextX, slope, intercept);
previousX = nextX;
previousY = nextY;
if (i < xLength) {
nextX = x[i];
nextY = y[i];
i++;
} else if (i === xLength) {
nextX += lastOriginalStep;
nextY = 0;
}
// updating parameters
slope = getSlope(previousX, previousY, nextX, nextY);
intercept = -slope*previousX + previousY;
}
return output;
}
/**
* function that retrieves the getEquallySpacedData with the variant "slot"
*
* @param x
* @param y
* @param from - Initial point
* @param to - Final point
* @param numberOfPoints
* @returns {Array} - Array of y's equally spaced with the variant "slot"
*/
function getEquallySpacedSlot(x, y, from, to, numberOfPoints) {
var xLength = x.length;
var step = (to - from) / (numberOfPoints - 1);
var halfStep = step / 2;
var lastStep = x[x.length - 1] - x[x.length - 2];
var start = from - halfStep;
var output = new Array(numberOfPoints);
// Init main variables
var min = start;
var max = start + step;
var previousX = -Number.MAX_VALUE;
var previousY = 0;
var nextX = x[0];
var nextY = y[0];
var frontOutsideSpectra = 0;
var backOutsideSpectra = true;
var currentValue = 0;
// for slot algorithm
var currentPoints = 0;
var i = 1; // index of input
var j = 0; // index of output
main: while(true) {
if (previousX>=nextX) throw (new Error('x must be an increasing serie'));
while (previousX - max > 0) {
// no overlap with original point, just consume current value
if(backOutsideSpectra) {
currentPoints++;
backOutsideSpectra = false;
}
output[j] = currentPoints <= 0 ? 0 : currentValue / currentPoints;
j++;
if (j === numberOfPoints)
break main;
min = max;
max += step;
currentValue = 0;
currentPoints = 0;
}
if(previousX > min) {
currentValue += previousY;
currentPoints++;
}
if(previousX === -Number.MAX_VALUE || frontOutsideSpectra > 1)
currentPoints--;
previousX = nextX;
previousY = nextY;
if (i < xLength) {
nextX = x[i];
nextY = y[i];
i++;
} else {
nextX += lastStep;
nextY = 0;
frontOutsideSpectra++;
}
}
return output;
}
/**
* Function that calculates the integral of the line between two
* x-coordinates, given the slope and intercept of the line.
*
* @param x0
* @param x1
* @param slope
* @param intercept
* @returns {number} integral value.
*/
function integral(x0, x1, slope, intercept) {
return (0.5 * slope * x1 * x1 + intercept * x1) - (0.5 * slope * x0 * x0 + intercept * x0);
}
MLArrayUtilsGetEquallySpaced.getEquallySpacedData = getEquallySpacedData;
MLArrayUtilsGetEquallySpaced.integral = integral;
}
// ml-array-utils snv.js
const MLArrayUtilsSNV = {};
{
MLArrayUtilsSNV.SNV = SNV;
let Stat = MLStat.array;
/**
* Function that applies the standard normal variate (SNV) to an array of values.
*
* @param data - Array of values.
* @returns {Array} - applied the SNV.
*/
function SNV(data) {
var mean = Stat.mean(data);
var std = Stat.standardDeviation(data);
var result = data.slice();
for (var i = 0; i < data.length; i++) {
result[i] = (result[i] - mean) / std;
}
return result;
}
}
// ml-array-utils index.js
const MLArrayUtils = {};
{
MLArrayUtils.getEquallySpacedData = MLArrayUtilsGetEquallySpaced.getEquallySpacedData;
MLArrayUtils.SNV = MLArrayUtilsSNV.SNV;
}
// do this early so things can use it. This is from ml-matrix src/matrix.js
const MLMatrixMatrix = {};
// ml-matrix src/util.js
const MLMatrixUtil = {};
{
let exports = MLMatrixUtil;
let Matrix = MLMatrixMatrix;
/**
* @private
* Check that a row index is not out of bounds
* @param {Matrix} matrix
* @param {number} index
* @param {boolean} [outer]
*/
exports.checkRowIndex = function checkRowIndex(matrix, index, outer) {
var max = outer ? matrix.rows : matrix.rows - 1;
if (index < 0 || index > max) {
throw new RangeError('Row index out of range');
}
};
/**
* @private
* Check that a column index is not out of bounds
* @param {Matrix} matrix
* @param {number} index
* @param {boolean} [outer]
*/
exports.checkColumnIndex = function checkColumnIndex(matrix, index, outer) {
var max = outer ? matrix.columns : matrix.columns - 1;
if (index < 0 || index > max) {
throw new RangeError('Column index out of range');
}
};
/**
* @private
* Check that the provided vector is an array with the right length
* @param {Matrix} matrix
* @param {Array|Matrix} vector
* @return {Array}
* @throws {RangeError}
*/
exports.checkRowVector = function checkRowVector(matrix, vector) {
if (vector.to1DArray) {
vector = vector.to1DArray();
}
if (vector.length !== matrix.columns) {
throw new RangeError('vector size must be the same as the number of columns');
}
return vector;
};
/**
* @private
* Check that the provided vector is an array with the right length
* @param {Matrix} matrix
* @param {Array|Matrix} vector
* @return {Array}
* @throws {RangeError}
*/
exports.checkColumnVector = function checkColumnVector(matrix, vector) {
if (vector.to1DArray) {
vector = vector.to1DArray();
}
if (vector.length !== matrix.rows) {
throw new RangeError('vector size must be the same as the number of rows');
}
return vector;
};
exports.checkIndices = function checkIndices(matrix, rowIndices, columnIndices) {
var rowOut = rowIndices.some(r => {
return r < 0 || r >= matrix.rows;
});
var columnOut = columnIndices.some(c => {
return c < 0 || c >= matrix.columns;
});
if (rowOut || columnOut) {
throw new RangeError('Indices are out of range');
}
if (typeof rowIndices !== 'object' || typeof columnIndices !== 'object') {
throw new TypeError('Unexpected type for row/column indices');
}
if (!Array.isArray(rowIndices)) rowIndices = Array.from(rowIndices);
if (!Array.isArray(columnIndices)) rowIndices = Array.from(columnIndices);
return {
row: rowIndices,
column: columnIndices
};
};
exports.checkRange = function checkRange(matrix, startRow, endRow, startColumn, endColumn) {
if (arguments.length !== 5) throw new TypeError('Invalid argument type');
var notAllNumbers = Array.from(arguments).slice(1).some(function (arg) {
return typeof arg !== 'number';
});
if (notAllNumbers) throw new TypeError('Invalid argument type');
if (startRow > endRow || startColumn > endColumn || startRow < 0 || startRow >= matrix.rows || endRow < 0 || endRow >= matrix.rows || startColumn < 0 || startColumn >= matrix.columns || endColumn < 0 || endColumn >= matrix.columns) {
throw new RangeError('Submatrix indices are out of range');
}
};
exports.getRange = function getRange(from, to) {
var arr = new Array(to - from + 1);
for (var i = 0; i < arr.length; i++) {
arr[i] = from + i;
}
return arr;
};
exports.sumByRow = function sumByRow(matrix) {
var sum = Matrix.Matrix.zeros(matrix.rows, 1);
for (var i = 0; i < matrix.rows; ++i) {
for (var j = 0; j < matrix.columns; ++j) {
sum.set(i, 0, sum.get(i, 0) + matrix.get(i, j));
}
}
return sum;
};
exports.sumByColumn = function sumByColumn(matrix) {
var sum = Matrix.Matrix.zeros(1, matrix.columns);
for (var i = 0; i < matrix.rows; ++i) {
for (var j = 0; j < matrix.columns; ++j) {
sum.set(0, j, sum.get(0, j) + matrix.get(i, j));
}
}
return sum;
};
exports.sumAll = function sumAll(matrix) {
var v = 0;
for (var i = 0; i < matrix.rows; i++) {
for (var j = 0; j < matrix.columns; j++) {
v += matrix.get(i, j);
}
}
return v;
};
}
// ml-matrix symbolsspecies.js
if (!Symbol.species) {
Symbol.species = Symbol.for('@@species');
}
// ml-matrix src/dc/util.js
const MLMatrixDCUtil = {};
{
let exports = MLMatrixDCUtil;
exports.hypotenuse = function hypotenuse(a, b) {
var r;
if (Math.abs(a) > Math.abs(b)) {
r = b / a;
return Math.abs(a) * Math.sqrt(1 + r * r);
}
if (b !== 0) {
r = a / b;
return Math.abs(b) * Math.sqrt(1 + r * r);
}
return 0;
};
// For use in the decomposition algorithms. With big matrices, access time is
// too long on elements from array subclass
// todo check when it is fixed in v8
exports.getEmpty2DArray = function (rows, columns) {
var array = new Array(rows);
for (var i = 0; i < rows; i++) {
array[i] = new Array(columns);
}
return array;
};
exports.getFilled2DArray = function (rows, columns, value) {
var array = new Array(rows);
for (var i = 0; i < rows; i++) {
array[i] = new Array(columns);
for (var j = 0; j < columns; j++) {
array[i][j] = value;
}
}
return array;
};
}
// ml-matrix src/dc/lu.js
let MLMatrixDCLU = {};
{
let Matrix = MLMatrixMatrix;
function LuDecomposition(matrix) {
if (!(this instanceof LuDecomposition)) {
return new LuDecomposition(matrix);
}
matrix = Matrix.Matrix.checkMatrix(matrix);
var lu = matrix.clone(),
rows = lu.rows,
columns = lu.columns,
pivotVector = new Array(rows),
pivotSign = 1,
i, j, k, p, s, t, v,
LUrowi, LUcolj, kmax;
for (i = 0; i < rows; i++) {
pivotVector[i] = i;
}
LUcolj = new Array(rows);
for (j = 0; j < columns; j++) {
for (i = 0; i < rows; i++) {
LUcolj[i] = lu[i][j];
}
for (i = 0; i < rows; i++) {
LUrowi = lu[i];
kmax = Math.min(i, j);
s = 0;
for (k = 0; k < kmax; k++) {
s += LUrowi[k] * LUcolj[k];
}
LUrowi[j] = LUcolj[i] -= s;
}
p = j;
for (i = j + 1; i < rows; i++) {
if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) {
p = i;
}
}
if (p !== j) {
for (k = 0; k < columns; k++) {
t = lu[p][k];
lu[p][k] = lu[j][k];
lu[j][k] = t;
}
v = pivotVector[p];
pivotVector[p] = pivotVector[j];
pivotVector[j] = v;
pivotSign = -pivotSign;
}
if (j < rows && lu[j][j] !== 0) {
for (i = j + 1; i < rows; i++) {
lu[i][j] /= lu[j][j];
}
}
}
this.LU = lu;
this.pivotVector = pivotVector;
this.pivotSign = pivotSign;
}
LuDecomposition.prototype = {
isSingular: function () {
var data = this.LU,
col = data.columns;
for (var j = 0; j < col; j++) {
if (data[j][j] === 0) {
return true;
}
}
return false;
},
get determinant() {
var data = this.LU;
if (!data.isSquare()) {
throw new Error('Matrix must be square');
}
var determinant = this.pivotSign, col = data.columns;
for (var j = 0; j < col; j++) {
determinant *= data[j][j];
}
return determinant;
},
get lowerTriangularMatrix() {
var data = this.LU,
rows = data.rows,
columns = data.columns,
X = new Matrix.Matrix(rows, columns);
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
if (i > j) {
X[i][j] = data[i][j];
} else if (i === j) {
X[i][j] = 1;
} else {
X[i][j] = 0;
}
}
}
return X;
},
get upperTriangularMatrix() {
var data = this.LU,
rows = data.rows,
columns = data.columns,
X = new Matrix.Matrix(rows, columns);
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
if (i <= j) {
X[i][j] = data[i][j];
} else {
X[i][j] = 0;
}
}
}
return X;
},
get pivotPermutationVector() {
return this.pivotVector.slice();
},
solve: function (value) {
value = Matrix.Matrix.checkMatrix(value);
var lu = this.LU,
rows = lu.rows;
if (rows !== value.rows) {
throw new Error('Invalid matrix dimensions');
}
if (this.isSingular()) {
throw new Error('LU matrix is singular');
}
var count = value.columns;
var X = value.subMatrixRow(this.pivotVector, 0, count - 1);
var columns = lu.columns;
var i, j, k;
for (k = 0; k < columns; k++) {
for (i = k + 1; i < columns; i++) {
for (j = 0; j < count; j++) {
X[i][j] -= X[k][j] * lu[i][k];
}
}
}
for (k = columns - 1; k >= 0; k--) {
for (j = 0; j < count; j++) {
X[k][j] /= lu[k][k];
}
for (i = 0; i < k; i++) {
for (j = 0; j < count; j++) {
X[i][j] -= X[k][j] * lu[i][k];
}
}
}
return X;
}
};
MLMatrixDCLU = LuDecomposition;
}
// ml-matrix src/dc/svd.js
let MLMatrixDCSVD = {};
{
let Matrix = MLMatrixMatrix;
let util = MLMatrixDCUtil;
let hypotenuse = util.hypotenuse;
let getFilled2DArray = util.getFilled2DArray;
function SingularValueDecomposition(value, options) {
if (!(this instanceof SingularValueDecomposition)) {
return new SingularValueDecomposition(value, options);
}
value = Matrix.Matrix.checkMatrix(value);
options = options || {};
var m = value.rows,
n = value.columns,
nu = Math.min(m, n);
var wantu = true, wantv = true;
if (options.computeLeftSingularVectors === false) wantu = false;
if (options.computeRightSingularVectors === false) wantv = false;
var autoTranspose = options.autoTranspose === true;
var swapped = false;
var a;
if (m < n) {
if (!autoTranspose) {
a = value.clone();
// eslint-disable-next-line no-console
console.warn('Computing SVD on a matrix with more columns than rows. Consider enabling autoTranspose');
} else {
a = value.transpose();
m = a.rows;
n = a.columns;
swapped = true;
var aux = wantu;
wantu = wantv;
wantv = aux;
}
} else {
a = value.clone();
}
var s = new Array(Math.min(m + 1, n)),
U = getFilled2DArray(m, nu, 0),
V = getFilled2DArray(n, n, 0),
e = new Array(n),
work = new Array(m);
var nct = Math.min(m - 1, n);
var nrt = Math.max(0, Math.min(n - 2, m));
var i, j, k, p, t, ks, f, cs, sn, max, kase,
scale, sp, spm1, epm1, sk, ek, b, c, shift, g;
for (k = 0, max = Math.max(nct, nrt); k < max; k++) {
if (k < nct) {
s[k] = 0;
for (i = k; i < m; i++) {
s[k] = hypotenuse(s[k], a[i][k]);
}
if (s[k] !== 0) {
if (a[k][k] < 0) {
s[k] = -s[k];
}
for (i = k; i < m; i++) {
a[i][k] /= s[k];
}
a[k][k] += 1;
}
s[k] = -s[k];
}
for (j = k + 1; j < n; j++) {
if ((k < nct) && (s[k] !== 0)) {
t = 0;
for (i = k; i < m; i++) {
t += a[i][k] * a[i][j];
}
t = -t / a[k][k];
for (i = k; i < m; i++) {
a[i][j] += t * a[i][k];
}
}
e[j] = a[k][j];
}
if (wantu && (k < nct)) {
for (i = k; i < m; i++) {
U[i][k] = a[i][k];
}
}
if (k < nrt) {
e[k] = 0;
for (i = k + 1; i < n; i++) {
e[k] = hypotenuse(e[k], e[i]);
}
if (e[k] !== 0) {
if (e[k + 1] < 0) {
e[k] = 0 - e[k];
}
for (i = k + 1; i < n; i++) {
e[i] /= e[k];
}
e[k + 1] += 1;
}
e[k] = -e[k];
if ((k + 1 < m) && (e[k] !== 0)) {
for (i = k + 1; i < m; i++) {
work[i] = 0;
}
for (j = k + 1; j < n; j++) {
for (i = k + 1; i < m; i++) {
work[i] += e[j] * a[i][j];
}
}
for (j = k + 1; j < n; j++) {
t = -e[j] / e[k + 1];
for (i = k + 1; i < m; i++) {
a[i][j] += t * work[i];
}
}
}
if (wantv) {
for (i = k + 1; i < n; i++) {
V[i][k] = e[i];
}
}
}
}
p = Math.min(n, m + 1);
if (nct < n) {
s[nct] = a[nct][nct];
}
if (m < p) {
s[p - 1] = 0;
}
if (nrt + 1 < p) {
e[nrt] = a[nrt][p - 1];
}
e[p - 1] = 0;
if (wantu) {
for (j = nct; j < nu; j++) {
for (i = 0; i < m; i++) {
U[i][j] = 0;
}
U[j][j] = 1;
}
for (k = nct - 1; k >= 0; k--) {
if (s[k] !== 0) {
for (j = k + 1; j < nu; j++) {
t = 0;
for (i = k; i < m; i++) {
t += U[i][k] * U[i][j];
}
t = -t / U[k][k];
for (i = k; i < m; i++) {
U[i][j] += t * U[i][k];
}
}
for (i = k; i < m; i++) {
U[i][k] = -U[i][k];
}
U[k][k] = 1 + U[k][k];
for (i = 0; i < k - 1; i++) {
U[i][k] = 0;
}
} else {
for (i = 0; i < m; i++) {
U[i][k] = 0;
}
U[k][k] = 1;
}
}
}
if (wantv) {
for (k = n - 1; k >= 0; k--) {
if ((k < nrt) && (e[k] !== 0)) {
for (j = k + 1; j < n; j++) {
t = 0;
for (i = k + 1; i < n; i++) {
t += V[i][k] * V[i][j];
}
t = -t / V[k + 1][k];
for (i = k + 1; i < n; i++) {
V[i][j] += t * V[i][k];
}
}
}
for (i = 0; i < n; i++) {
V[i][k] = 0;
}
V[k][k] = 1;
}
}
var pp = p - 1,
iter = 0,
eps = Math.pow(2, -52);
while (p > 0) {
for (k = p - 2; k >= -1; k--) {
if (k === -1) {
break;
}
if (Math.abs(e[k]) <= eps * (Math.abs(s[k]) + Math.abs(s[k + 1]))) {
e[k] = 0;
break;
}
}
if (k === p - 2) {
kase = 4;
} else {
for (ks = p - 1; ks >= k; ks--) {
if (ks === k) {
break;
}
t = (ks !== p ? Math.abs(e[ks]) : 0) + (ks !== k + 1 ? Math.abs(e[ks - 1]) : 0);
if (Math.abs(s[ks]) <= eps * t) {
s[ks] = 0;
break;
}
}
if (ks === k) {
kase = 3;
} else if (ks === p - 1) {
kase = 1;
} else {
kase = 2;
k = ks;
}
}
k++;
switch (kase) {
case 1: {
f = e[p - 2];
e[p - 2] = 0;
for (j = p - 2; j >= k; j--) {
t = hypotenuse(s[j], f);
cs = s[j] / t;
sn = f / t;
s[j] = t;
if (j !== k) {
f = -sn * e[j - 1];
e[j - 1] = cs * e[j - 1];
}
if (wantv) {
for (i = 0; i < n; i++) {
t = cs * V[i][j] + sn * V[i][p - 1];
V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
V[i][j] = t;
}
}
}
break;
}
case 2 : {
f = e[k - 1];
e[k - 1] = 0;
for (j = k; j < p; j++) {
t = hypotenuse(s[j], f);
cs = s[j] / t;
sn = f / t;
s[j] = t;
f = -sn * e[j];
e[j] = cs * e[j];
if (wantu) {
for (i = 0; i < m; i++) {
t = cs * U[i][j] + sn * U[i][k - 1];
U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
U[i][j] = t;
}
}
}
break;
}
case 3 : {
scale = Math.max(Math.max(Math.max(Math.max(Math.abs(s[p - 1]), Math.abs(s[p - 2])), Math.abs(e[p - 2])), Math.abs(s[k])), Math.abs(e[k]));
sp = s[p - 1] / scale;
spm1 = s[p - 2] / scale;
epm1 = e[p - 2] / scale;
sk = s[k] / scale;
ek = e[k] / scale;
b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2;
c = (sp * epm1) * (sp * epm1);
shift = 0;
if ((b !== 0) || (c !== 0)) {
shift = Math.sqrt(b * b + c);
if (b < 0) {
shift = -shift;
}
shift = c / (b + shift);
}
f = (sk + sp) * (sk - sp) + shift;
g = sk * ek;
for (j = k; j < p - 1; j++) {
t = hypotenuse(f, g);
cs = f / t;
sn = g / t;
if (j !== k) {
e[j - 1] = t;
}
f = cs * s[j] + sn * e[j];
e[j] = cs * e[j] - sn * s[j];
g = sn * s[j + 1];
s[j + 1] = cs * s[j + 1];
if (wantv) {
for (i = 0; i < n; i++) {
t = cs * V[i][j] + sn * V[i][j + 1];
V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
V[i][j] = t;
}
}
t = hypotenuse(f, g);
cs = f / t;
sn = g / t;
s[j] = t;
f = cs * e[j] + sn * s[j + 1];
s[j + 1] = -sn * e[j] + cs * s[j + 1];
g = sn * e[j + 1];
e[j + 1] = cs * e[j + 1];
if (wantu && (j < m - 1)) {
for (i = 0; i < m; i++) {
t = cs * U[i][j] + sn * U[i][j + 1];
U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
U[i][j] = t;
}
}
}
e[p - 2] = f;
iter = iter + 1;
break;
}
case 4: {
if (s[k] <= 0) {
s[k] = (s[k] < 0 ? -s[k] : 0);
if (wantv) {
for (i = 0; i <= pp; i++) {
V[i][k] = -V[i][k];
}
}
}
while (k < pp) {
if (s[k] >= s[k + 1]) {
break;
}
t = s[k];
s[k] = s[k + 1];
s[k + 1] = t;
if (wantv && (k < n - 1)) {
for (i = 0; i < n; i++) {
t = V[i][k + 1];
V[i][k + 1] = V[i][k];
V[i][k] = t;
}
}
if (wantu && (k < m - 1)) {
for (i = 0; i < m; i++) {
t = U[i][k + 1];
U[i][k + 1] = U[i][k];
U[i][k] = t;
}
}
k++;
}
iter = 0;
p--;
break;
}
// no default
}
}
if (swapped) {
var tmp = V;
V = U;
U = tmp;
}
this.m = m;
this.n = n;
this.s = s;
this.U = U;
this.V = V;
}
SingularValueDecomposition.prototype = {
get condition() {
return this.s[0] / this.s[Math.min(this.m, this.n) - 1];
},
get norm2() {
return this.s[0];
},
get rank() {
var eps = Math.pow(2, -52),
tol = Math.max(this.m, this.n) * this.s[0] * eps,
r = 0,
s = this.s;
for (var i = 0, ii = s.length; i < ii; i++) {
if (s[i] > tol) {
r++;
}
}
return r;
},
get diagonal() {
return this.s;
},
get threshold() {
return (Math.pow(2, -52) / 2) * Math.max(this.m, this.n) * this.s[0];
},
get leftSingularVectors() {
if (!Matrix.Matrix.isMatrix(this.U)) {
this.U = new Matrix.Matrix(this.U);
}
return this.U;
},
get rightSingularVectors() {
if (!Matrix.Matrix.isMatrix(this.V)) {
this.V = new Matrix.Matrix(this.V);
}
return this.V;
},
get diagonalMatrix() {
return Matrix.Matrix.diag(this.s);
},
solve: function (value) {
var Y = value,
e = this.threshold,
scols = this.s.length,
Ls = Matrix.Matrix.zeros(scols, scols),
i;
for (i = 0; i < scols; i++) {
if (Math.abs(this.s[i]) <= e) {
Ls[i][i] = 0;
} else {
Ls[i][i] = 1 / this.s[i];
}
}
var U = this.U;
var V = this.rightSingularVectors;
var VL = V.mmul(Ls),
vrows = V.rows,
urows = U.length,
VLU = Matrix.Matrix.zeros(vrows, urows),
j, k, sum;
for (i = 0; i < vrows; i++) {
for (j = 0; j < urows; j++) {
sum = 0;
for (k = 0; k < scols; k++) {
sum += VL[i][k] * U[j][k];
}
VLU[i][j] = sum;
}
}
return VLU.mmul(Y);
},
solveForDiagonal: function (value) {
return this.solve(Matrix.Matrix.diag(value));
},
inverse: function () {
var V = this.V;
var e = this.threshold,
vrows = V.length,
vcols = V[0].length,
X = new Matrix.Matrix(vrows, this.s.length),
i, j;
for (i = 0; i < vrows; i++) {
for (j = 0; j < vcols; j++) {
if (Math.abs(this.s[j]) > e) {
X[i][j] = V[i][j] / this.s[j];
} else {
X[i][j] = 0;
}
}
}
var U = this.U;
var urows = U.length,
ucols = U[0].length,
Y = new Matrix.Matrix(vrows, urows),
k, sum;
for (i = 0; i < vrows; i++) {
for (j = 0; j < urows; j++) {
sum = 0;
for (k = 0; k < ucols; k++) {
sum += X[i][k] * U[j][k];
}
Y[i][j] = sum;
}
}
return Y;
}
};
MLMatrixDCSVD = SingularValueDecomposition;
}
// ml-matrix src/abstractMatrix.js
let MLMatrixAbstractMatrix;
{
let LuDecomposition = MLMatrixDCLU;
let SvDecomposition = MLMatrixDCSVD;
let arrayUtils = MLArrayUtils;
let util = MLMatrixUtil;
MLMatrixAbstractMatrix = function abstractMatrix(superCtor) {
if (superCtor === undefined) superCtor = Object;
/**
* Real matrix
* @class Matrix
* @param {number|Array|Matrix} nRows - Number of rows of the new matrix,
* 2D array containing the data or Matrix instance to clone
* @param {number} [nColumns] - Number of columns of the new matrix
*/
class Matrix extends superCtor {
static get [Symbol.species]() {
return this;
}
/**
* Constructs a Matrix with the chosen dimensions from a 1D array
* @param {number} newRows - Number of rows
* @param {number} newColumns - Number of columns
* @param {Array} newData - A 1D array containing data for the matrix
* @return {Matrix} - The new matrix
*/
static from1DArray(newRows, newColumns, newData) {
var length = newRows * newColumns;
if (length !== newData.length) {
throw new RangeError('Data length does not match given dimensions');
}
var newMatrix = new this(newRows, newColumns);
for (var row = 0; row < newRows; row++) {
for (var column = 0; column < newColumns; column++) {
newMatrix.set(row, column, newData[row * newColumns + column]);
}
}
return newMatrix;
}
/**
* Creates a row vector, a matrix with only one row.
* @param {Array} newData - A 1D array containing data for the vector
* @return {Matrix} - The new matrix
*/
static rowVector(newData) {
var vector = new this(1, newData.length);
for (var i = 0; i < newData.length; i++) {
vector.set(0, i, newData[i]);
}
return vector;
}
/**
* Creates a column vector, a matrix with only one column.
* @param {Array} newData - A 1D array containing data for the vector
* @return {Matrix} - The new matrix
*/
static columnVector(newData) {
var vector = new this(newData.length, 1);
for (var i = 0; i < newData.length; i++) {
vector.set(i, 0, newData[i]);
}
return vector;
}
/**
* Creates an empty matrix with the given dimensions. Values will be undefined. Same as using new Matrix(rows, columns).
* @param {number} rows - Number of rows
* @param {number} columns - Number of columns
* @return {Matrix} - The new matrix
*/
static empty(rows, columns) {
return new this(rows, columns);
}
/**
* Creates a matrix with the given dimensions. Values will be set to zero.
* @param {number} rows - Number of rows
* @param {number} columns - Number of columns
* @return {Matrix} - The new matrix
*/
static zeros(rows, columns) {
return this.empty(rows, columns).fill(0);
}
/**
* Creates a matrix with the given dimensions. Values will be set to one.
* @param {number} rows - Number of rows
* @param {number} columns - Number of columns
* @return {Matrix} - The new matrix
*/
static ones(rows, columns) {
return this.empty(rows, columns).fill(1);
}
/**
* Creates a matrix with the given dimensions. Values will be randomly set.
* @param {number} rows - Number of rows
* @param {number} columns - Number of columns
* @param {function} [rng=Math.random] - Random number generator
* @return {Matrix} The new matrix
*/
static rand(rows, columns, rng) {
if (rng === undefined) rng = Math.random;
var matrix = this.empty(rows, columns);
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
matrix.set(i, j, rng());
}
}
return matrix;
}
/**
* Creates a matrix with the given dimensions. Values will be random integers.
* @param {number} rows - Number of rows
* @param {number} columns - Number of columns
* @param {number} [maxValue=1000] - Maximum value
* @param {function} [rng=Math.random] - Random number generator
* @return {Matrix} The new matrix
*/
static randInt(rows, columns, maxValue, rng) {
if (maxValue === undefined) maxValue = 1000;
if (rng === undefined) rng = Math.random;
var matrix = this.empty(rows, columns);
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
var value = Math.floor(rng() * maxValue);
matrix.set(i, j, value);
}
}
return matrix;
}
/**
* Creates an identity matrix with the given dimension. Values of the diagonal will be 1 and others will be 0.
* @param {number} rows - Number of rows
* @param {number} [columns=rows] - Number of columns
* @param {number} [value=1] - Value to fill the diagonal with
* @return {Matrix} - The new identity matrix
*/
static eye(rows, columns, value) {
if (columns === undefined) columns = rows;
if (value === undefined) value = 1;
var min = Math.min(rows, columns);
var matrix = this.zeros(rows, columns);
for (var i = 0; i < min; i++) {
matrix.set(i, i, value);
}
return matrix;
}
/**
* Creates a diagonal matrix based on the given array.
* @param {Array} data - Array containing the data for the diagonal
* @param {number} [rows] - Number of rows (Default: data.length)
* @param {number} [columns] - Number of columns (Default: rows)
* @return {Matrix} - The new diagonal matrix
*/
static diag(data, rows, columns) {
var l = data.length;
if (rows === undefined) rows = l;
if (columns === undefined) columns = rows;
var min = Math.min(l, rows, columns);
var matrix = this.zeros(rows, columns);
for (var i = 0; i < min; i++) {
matrix.set(i, i, data[i]);
}
return matrix;
}
/**
* Returns a matrix whose elements are the minimum between matrix1 and matrix2
* @param {Matrix} matrix1
* @param {Matrix} matrix2
* @return {Matrix}
*/
static min(matrix1, matrix2) {
matrix1 = this.checkMatrix(matrix1);
matrix2 = this.checkMatrix(matrix2);
var rows = matrix1.rows;
var columns = matrix1.columns;
var result = new this(rows, columns);
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
result.set(i, j, Math.min(matrix1.get(i, j), matrix2.get(i, j)));
}
}
return result;
}
/**
* Returns a matrix whose elements are the maximum between matrix1 and matrix2
* @param {Matrix} matrix1
* @param {Matrix} matrix2
* @return {Matrix}
*/
static max(matrix1, matrix2) {
matrix1 = this.checkMatrix(matrix1);
matrix2 = this.checkMatrix(matrix2);
var rows = matrix1.rows;
var columns = matrix1.columns;
var result = new this(rows, columns);
for (var i = 0; i < rows; i++) {
for (var j = 0; j < columns; j++) {
result.set(i, j, Math.max(matrix1.get(i, j), matrix2.get(i, j)));
}
}
return result;
}
/**
* Check that the provided value is a Matrix and tries to instantiate one if not
* @param {*} value - The value to check
* @return {Matrix}
*/
static checkMatrix(value) {
return Matrix.isMatrix(value) ? value : new this(value);
}
/**
* Returns true if the argument is a Matrix, false otherwise
* @param {*} value - The value to check
* @return {boolean}
*/
static isMatrix(value) {
return (value != null) && (value.klass === 'Matrix');
}
/**
* @prop {number} size - The number of elements in the matrix.
*/
get size() {
return this.rows * this.columns;
}
/**
* Applies a callback for each element of the matrix. The function is called in the matrix (this) context.
* @param {function} callback - Function that will be called with two parameters : i (row) and j (column)
* @return {Matrix} this
*/
apply(callback) {
if (typeof callback !== 'function') {
throw new TypeError('callback must be a function');
}
var ii = this.rows;
var jj = this.columns;
for (var i = 0; i < ii; i++) {
for (var j = 0; j < jj; j++) {
callback.call(this, i, j);
}
}
return this;
}
/**
* Returns a new 1D array filled row by row with the matrix values
* @return {Array}
*/
to1DArray() {
var array = new Array(this.size);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
array[i * this.columns + j] = this.get(i, j);
}
}
return array;
}
/**
* Returns a 2D array containing a copy of the data
* @return {Array}
*/
to2DArray() {
var copy = new Array(this.rows);
for (var i = 0; i < this.rows; i++) {
copy[i] = new Array(this.columns);
for (var j = 0; j < this.columns; j++) {
copy[i][j] = this.get(i, j);
}
}
return copy;
}
/**
* @return {boolean} true if the matrix has one row
*/
isRowVector() {
return this.rows === 1;
}
/**
* @return {boolean} true if the matrix has one column
*/
isColumnVector() {
return this.columns === 1;
}
/**
* @return {boolean} true if the matrix has one row or one column
*/
isVector() {
return (this.rows === 1) || (this.columns === 1);
}
/**
* @return {boolean} true if the matrix has the same number of rows and columns
*/
isSquare() {
return this.rows === this.columns;
}
/**
* @return {boolean} true if the matrix is square and has the same values on both sides of the diagonal
*/
isSymmetric() {
if (this.isSquare()) {
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j <= i; j++) {
if (this.get(i, j) !== this.get(j, i)) {
return false;
}
}
}
return true;
}
return false;
}
/**
* Sets a given element of the matrix. mat.set(3,4,1) is equivalent to mat[3][4]=1
* @abstract
* @param {number} rowIndex - Index of the row
* @param {number} columnIndex - Index of the column
* @param {number} value - The new value for the element
* @return {Matrix} this
*/
set(rowIndex, columnIndex, value) { // eslint-disable-line no-unused-vars
throw new Error('set method is unimplemented');
}
/**
* Returns the given element of the matrix. mat.get(3,4) is equivalent to matrix[3][4]
* @abstract
* @param {number} rowIndex - Index of the row
* @param {number} columnIndex - Index of the column
* @return {number}
*/
get(rowIndex, columnIndex) { // eslint-disable-line no-unused-vars
throw new Error('get method is unimplemented');
}
/**
* Creates a new matrix that is a repetition of the current matrix. New matrix has rowRep times the number of
* rows of the matrix, and colRep times the number of columns of the matrix
* @param {number} rowRep - Number of times the rows should be repeated
* @param {number} colRep - Number of times the columns should be re
* @return {Matrix}
* @example
* var matrix = new Matrix([[1,2]]);
* matrix.repeat(2); // [[1,2],[1,2]]
*/
repeat(rowRep, colRep) {
rowRep = rowRep || 1;
colRep = colRep || 1;
var matrix = new this.constructor[Symbol.species](this.rows * rowRep, this.columns * colRep);
for (var i = 0; i < rowRep; i++) {
for (var j = 0; j < colRep; j++) {
matrix.setSubMatrix(this, this.rows * i, this.columns * j);
}
}
return matrix;
}
/**
* Fills the matrix with a given value. All elements will be set to this value.
* @param {number} value - New value
* @return {Matrix} this
*/
fill(value) {
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, value);
}
}
return this;
}
/**
* Negates the matrix. All elements will be multiplied by (-1)
* @return {Matrix} this
*/
neg() {
return this.mulS(-1);
}
/**
* Returns a new array from the given row index
* @param {number} index - Row index
* @return {Array}
*/
getRow(index) {
util.checkRowIndex(this, index);
var row = new Array(this.columns);
for (var i = 0; i < this.columns; i++) {
row[i] = this.get(index, i);
}
return row;
}
/**
* Returns a new row vector from the given row index
* @param {number} index - Row index
* @return {Matrix}
*/
getRowVector(index) {
return this.constructor.rowVector(this.getRow(index));
}
/**
* Sets a row at the given index
* @param {number} index - Row index
* @param {Array|Matrix} array - Array or vector
* @return {Matrix} this
*/
setRow(index, array) {
util.checkRowIndex(this, index);
array = util.checkRowVector(this, array);
for (var i = 0; i < this.columns; i++) {
this.set(index, i, array[i]);
}
return this;
}
/**
* Swaps two rows
* @param {number} row1 - First row index
* @param {number} row2 - Second row index
* @return {Matrix} this
*/
swapRows(row1, row2) {
util.checkRowIndex(this, row1);
util.checkRowIndex(this, row2);
for (var i = 0; i < this.columns; i++) {
var temp = this.get(row1, i);
this.set(row1, i, this.get(row2, i));
this.set(row2, i, temp);
}
return this;
}
/**
* Returns a new array from the given column index
* @param {number} index - Column index
* @return {Array}
*/
getColumn(index) {
util.checkColumnIndex(this, index);
var column = new Array(this.rows);
for (var i = 0; i < this.rows; i++) {
column[i] = this.get(i, index);
}
return column;
}
/**
* Returns a new column vector from the given column index
* @param {number} index - Column index
* @return {Matrix}
*/
getColumnVector(index) {
return this.constructor.columnVector(this.getColumn(index));
}
/**
* Sets a column at the given index
* @param {number} index - Column index
* @param {Array|Matrix} array - Array or vector
* @return {Matrix} this
*/
setColumn(index, array) {
util.checkColumnIndex(this, index);
array = util.checkColumnVector(this, array);
for (var i = 0; i < this.rows; i++) {
this.set(i, index, array[i]);
}
return this;
}
/**
* Swaps two columns
* @param {number} column1 - First column index
* @param {number} column2 - Second column index
* @return {Matrix} this
*/
swapColumns(column1, column2) {
util.checkColumnIndex(this, column1);
util.checkColumnIndex(this, column2);
for (var i = 0; i < this.rows; i++) {
var temp = this.get(i, column1);
this.set(i, column1, this.get(i, column2));
this.set(i, column2, temp);
}
return this;
}
/**
* Adds the values of a vector to each row
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
addRowVector(vector) {
vector = util.checkRowVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) + vector[j]);
}
}
return this;
}
/**
* Subtracts the values of a vector from each row
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
subRowVector(vector) {
vector = util.checkRowVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) - vector[j]);
}
}
return this;
}
/**
* Multiplies the values of a vector with each row
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
mulRowVector(vector) {
vector = util.checkRowVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) * vector[j]);
}
}
return this;
}
/**
* Divides the values of each row by those of a vector
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
divRowVector(vector) {
vector = util.checkRowVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) / vector[j]);
}
}
return this;
}
/**
* Adds the values of a vector to each column
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
addColumnVector(vector) {
vector = util.checkColumnVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) + vector[i]);
}
}
return this;
}
/**
* Subtracts the values of a vector from each column
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
subColumnVector(vector) {
vector = util.checkColumnVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) - vector[i]);
}
}
return this;
}
/**
* Multiplies the values of a vector with each column
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
mulColumnVector(vector) {
vector = util.checkColumnVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) * vector[i]);
}
}
return this;
}
/**
* Divides the values of each column by those of a vector
* @param {Array|Matrix} vector - Array or vector
* @return {Matrix} this
*/
divColumnVector(vector) {
vector = util.checkColumnVector(this, vector);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) / vector[i]);
}
}
return this;
}
/**
* Multiplies the values of a row with a scalar
* @param {number} index - Row index
* @param {number} value
* @return {Matrix} this
*/
mulRow(index, value) {
util.checkRowIndex(this, index);
for (var i = 0; i < this.columns; i++) {
this.set(index, i, this.get(index, i) * value);
}
return this;
}
/**
* Multiplies the values of a column with a scalar
* @param {number} index - Column index
* @param {number} value
* @return {Matrix} this
*/
mulColumn(index, value) {
util.checkColumnIndex(this, index);
for (var i = 0; i < this.rows; i++) {
this.set(i, index, this.get(i, index) * value);
}
return this;
}
/**
* Returns the maximum value of the matrix
* @return {number}
*/
max() {
var v = this.get(0, 0);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
if (this.get(i, j) > v) {
v = this.get(i, j);
}
}
}
return v;
}
/**
* Returns the index of the maximum value
* @return {Array}
*/
maxIndex() {
var v = this.get(0, 0);
var idx = [0, 0];
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
if (this.get(i, j) > v) {
v = this.get(i, j);
idx[0] = i;
idx[1] = j;
}
}
}
return idx;
}
/**
* Returns the minimum value of the matrix
* @return {number}
*/
min() {
var v = this.get(0, 0);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
if (this.get(i, j) < v) {
v = this.get(i, j);
}
}
}
return v;
}
/**
* Returns the index of the minimum value
* @return {Array}
*/
minIndex() {
var v = this.get(0, 0);
var idx = [0, 0];
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
if (this.get(i, j) < v) {
v = this.get(i, j);
idx[0] = i;
idx[1] = j;
}
}
}
return idx;
}
/**
* Returns the maximum value of one row
* @param {number} row - Row index
* @return {number}
*/
maxRow(row) {
util.checkRowIndex(this, row);
var v = this.get(row, 0);
for (var i = 1; i < this.columns; i++) {
if (this.get(row, i) > v) {
v = this.get(row, i);
}
}
return v;
}
/**
* Returns the index of the maximum value of one row
* @param {number} row - Row index
* @return {Array}
*/
maxRowIndex(row) {
util.checkRowIndex(this, row);
var v = this.get(row, 0);
var idx = [row, 0];
for (var i = 1; i < this.columns; i++) {
if (this.get(row, i) > v) {
v = this.get(row, i);
idx[1] = i;
}
}
return idx;
}
/**
* Returns the minimum value of one row
* @param {number} row - Row index
* @return {number}
*/
minRow(row) {
util.checkRowIndex(this, row);
var v = this.get(row, 0);
for (var i = 1; i < this.columns; i++) {
if (this.get(row, i) < v) {
v = this.get(row, i);
}
}
return v;
}
/**
* Returns the index of the maximum value of one row
* @param {number} row - Row index
* @return {Array}
*/
minRowIndex(row) {
util.checkRowIndex(this, row);
var v = this.get(row, 0);
var idx = [row, 0];
for (var i = 1; i < this.columns; i++) {
if (this.get(row, i) < v) {
v = this.get(row, i);
idx[1] = i;
}
}
return idx;
}
/**
* Returns the maximum value of one column
* @param {number} column - Column index
* @return {number}
*/
maxColumn(column) {
util.checkColumnIndex(this, column);
var v = this.get(0, column);
for (var i = 1; i < this.rows; i++) {
if (this.get(i, column) > v) {
v = this.get(i, column);
}
}
return v;
}
/**
* Returns the index of the maximum value of one column
* @param {number} column - Column index
* @return {Array}
*/
maxColumnIndex(column) {
util.checkColumnIndex(this, column);
var v = this.get(0, column);
var idx = [0, column];
for (var i = 1; i < this.rows; i++) {
if (this.get(i, column) > v) {
v = this.get(i, column);
idx[0] = i;
}
}
return idx;
}
/**
* Returns the minimum value of one column
* @param {number} column - Column index
* @return {number}
*/
minColumn(column) {
util.checkColumnIndex(this, column);
var v = this.get(0, column);
for (var i = 1; i < this.rows; i++) {
if (this.get(i, column) < v) {
v = this.get(i, column);
}
}
return v;
}
/**
* Returns the index of the minimum value of one column
* @param {number} column - Column index
* @return {Array}
*/
minColumnIndex(column) {
util.checkColumnIndex(this, column);
var v = this.get(0, column);
var idx = [0, column];
for (var i = 1; i < this.rows; i++) {
if (this.get(i, column) < v) {
v = this.get(i, column);
idx[0] = i;
}
}
return idx;
}
/**
* Returns an array containing the diagonal values of the matrix
* @return {Array}
*/
diag() {
var min = Math.min(this.rows, this.columns);
var diag = new Array(min);
for (var i = 0; i < min; i++) {
diag[i] = this.get(i, i);
}
return diag;
}
/**
* Returns the sum by the argument given, if no argument given,
* it returns the sum of all elements of the matrix.
* @param {string} by - sum by 'row' or 'column'.
* @return {Matrix|number}
*/
sum(by) {
switch (by) {
case 'row':
return util.sumByRow(this);
case 'column':
return util.sumByColumn(this);
default:
return util.sumAll(this);
}
}
/**
* Returns the mean of all elements of the matrix
* @return {number}
*/
mean() {
return this.sum() / this.size;
}
/**
* Returns the product of all elements of the matrix
* @return {number}
*/
prod() {
var prod = 1;
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
prod *= this.get(i, j);
}
}
return prod;
}
/**
* Computes the cumulative sum of the matrix elements (in place, row by row)
* @return {Matrix} this
*/
cumulativeSum() {
var sum = 0;
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
sum += this.get(i, j);
this.set(i, j, sum);
}
}
return this;
}
/**
* Computes the dot (scalar) product between the matrix and another
* @param {Matrix} vector2 vector
* @return {number}
*/
dot(vector2) {
if (Matrix.isMatrix(vector2)) vector2 = vector2.to1DArray();
var vector1 = this.to1DArray();
if (vector1.length !== vector2.length) {
throw new RangeError('vectors do not have the same size');
}
var dot = 0;
for (var i = 0; i < vector1.length; i++) {
dot += vector1[i] * vector2[i];
}
return dot;
}
/**
* Returns the matrix product between this and other
* @param {Matrix} other
* @return {Matrix}
*/
mmul(other) {
other = this.constructor.checkMatrix(other);
if (this.columns !== other.rows) {
// eslint-disable-next-line no-console
console.warn('Number of columns of left matrix are not equal to number of rows of right matrix.');
}
var m = this.rows;
var n = this.columns;
var p = other.columns;
var result = new this.constructor[Symbol.species](m, p);
var Bcolj = new Array(n);
for (var j = 0; j < p; j++) {
for (var k = 0; k < n; k++) {
Bcolj[k] = other.get(k, j);
}
for (var i = 0; i < m; i++) {
var s = 0;
for (k = 0; k < n; k++) {
s += this.get(i, k) * Bcolj[k];
}
result.set(i, j, s);
}
}
return result;
}
strassen2x2(other) {
var result = new this.constructor[Symbol.species](2, 2);
const a11 = this.get(0, 0);
const b11 = other.get(0, 0);
const a12 = this.get(0, 1);
const b12 = other.get(0, 1);
const a21 = this.get(1, 0);
const b21 = other.get(1, 0);
const a22 = this.get(1, 1);
const b22 = other.get(1, 1);
// Compute intermediate values.
const m1 = (a11 + a22) * (b11 + b22);
const m2 = (a21 + a22) * b11;
const m3 = a11 * (b12 - b22);
const m4 = a22 * (b21 - b11);
const m5 = (a11 + a12) * b22;
const m6 = (a21 - a11) * (b11 + b12);
const m7 = (a12 - a22) * (b21 + b22);
// Combine intermediate values into the output.
const c00 = m1 + m4 - m5 + m7;
const c01 = m3 + m5;
const c10 = m2 + m4;
const c11 = m1 - m2 + m3 + m6;
result.set(0, 0, c00);
result.set(0, 1, c01);
result.set(1, 0, c10);
result.set(1, 1, c11);
return result;
}
strassen3x3(other) {
var result = new this.constructor[Symbol.species](3, 3);
const a00 = this.get(0, 0);
const a01 = this.get(0, 1);
const a02 = this.get(0, 2);
const a10 = this.get(1, 0);
const a11 = this.get(1, 1);
const a12 = this.get(1, 2);
const a20 = this.get(2, 0);
const a21 = this.get(2, 1);
const a22 = this.get(2, 2);
const b00 = other.get(0, 0);
const b01 = other.get(0, 1);
const b02 = other.get(0, 2);
const b10 = other.get(1, 0);
const b11 = other.get(1, 1);
const b12 = other.get(1, 2);
const b20 = other.get(2, 0);
const b21 = other.get(2, 1);
const b22 = other.get(2, 2);
const m1 = (a00 + a01 + a02 - a10 - a11 - a21 - a22) * b11;
const m2 = (a00 - a10) * (-b01 + b11);
const m3 = a11 * (-b00 + b01 + b10 - b11 - b12 - b20 + b22);
const m4 = (-a00 + a10 + a11) * (b00 - b01 + b11);
const m5 = (a10 + a11) * (-b00 + b01);
const m6 = a00 * b00;
const m7 = (-a00 + a20 + a21) * (b00 - b02 + b12);
const m8 = (-a00 + a20) * (b02 - b12);
const m9 = (a20 + a21) * (-b00 + b02);
const m10 = (a00 + a01 + a02 - a11 - a12 - a20 - a21) * b12;
const m11 = a21 * (-b00 + b02 + b10 - b11 - b12 - b20 + b21);
const m12 = (-a02 + a21 + a22) * (b11 + b20 - b21);
const m13 = (a02 - a22) * (b11 - b21);
const m14 = a02 * b20;
const m15 = (a21 + a22) * (-b20 + b21);
const m16 = (-a02 + a11 + a12) * (b12 + b20 - b22);
const m17 = (a02 - a12) * (b12 - b22);
const m18 = (a11 + a12) * (-b20 + b22);
const m19 = a01 * b10;
const m20 = a12 * b21;
const m21 = a10 * b02;
const m22 = a20 * b01;
const m23 = a22 * b22;
const c00 = m6 + m14 + m19;
const c01 = m1 + m4 + m5 + m6 + m12 + m14 + m15;
const c02 = m6 + m7 + m9 + m10 + m14 + m16 + m18;
const c10 = m2 + m3 + m4 + m6 + m14 + m16 + m17;
const c11 = m2 + m4 + m5 + m6 + m20;
const c12 = m14 + m16 + m17 + m18 + m21;
const c20 = m6 + m7 + m8 + m11 + m12 + m13 + m14;
const c21 = m12 + m13 + m14 + m15 + m22;
const c22 = m6 + m7 + m8 + m9 + m23;
result.set(0, 0, c00);
result.set(0, 1, c01);
result.set(0, 2, c02);
result.set(1, 0, c10);
result.set(1, 1, c11);
result.set(1, 2, c12);
result.set(2, 0, c20);
result.set(2, 1, c21);
result.set(2, 2, c22);
return result;
}
/**
* Returns the matrix product between x and y. More efficient than mmul(other) only when we multiply squared matrix and when the size of the matrix is > 1000.
* @param {Matrix} y
* @return {Matrix}
*/
mmulStrassen(y) {
var x = this.clone();
var r1 = x.rows;
var c1 = x.columns;
var r2 = y.rows;
var c2 = y.columns;
if (c1 !== r2) {
// eslint-disable-next-line no-console
console.warn(`Multiplying ${r1} x ${c1} and ${r2} x ${c2} matrix: dimensions do not match.`);
}
// Put a matrix into the top left of a matrix of zeros.
// `rows` and `cols` are the dimensions of the output matrix.
function embed(mat, rows, cols) {
var r = mat.rows;
var c = mat.columns;
if ((r === rows) && (c === cols)) {
return mat;
} else {
var resultat = Matrix.zeros(rows, cols);
resultat = resultat.setSubMatrix(mat, 0, 0);
return resultat;
}
}
// Make sure both matrices are the same size.
// This is exclusively for simplicity:
// this algorithm can be implemented with matrices of different sizes.
var r = Math.max(r1, r2);
var c = Math.max(c1, c2);
x = embed(x, r, c);
y = embed(y, r, c);
// Our recursive multiplication function.
function blockMult(a, b, rows, cols) {
// For small matrices, resort to naive multiplication.
if (rows <= 512 || cols <= 512) {
return a.mmul(b); // a is equivalent to this
}
// Apply dynamic padding.
if ((rows % 2 === 1) && (cols % 2 === 1)) {
a = embed(a, rows + 1, cols + 1);
b = embed(b, rows + 1, cols + 1);
} else if (rows % 2 === 1) {
a = embed(a, rows + 1, cols);
b = embed(b, rows + 1, cols);
} else if (cols % 2 === 1) {
a = embed(a, rows, cols + 1);
b = embed(b, rows, cols + 1);
}
var halfRows = parseInt(a.rows / 2);
var halfCols = parseInt(a.columns / 2);
// Subdivide input matrices.
var a11 = a.subMatrix(0, halfRows - 1, 0, halfCols - 1);
var b11 = b.subMatrix(0, halfRows - 1, 0, halfCols - 1);
var a12 = a.subMatrix(0, halfRows - 1, halfCols, a.columns - 1);
var b12 = b.subMatrix(0, halfRows - 1, halfCols, b.columns - 1);
var a21 = a.subMatrix(halfRows, a.rows - 1, 0, halfCols - 1);
var b21 = b.subMatrix(halfRows, b.rows - 1, 0, halfCols - 1);
var a22 = a.subMatrix(halfRows, a.rows - 1, halfCols, a.columns - 1);
var b22 = b.subMatrix(halfRows, b.rows - 1, halfCols, b.columns - 1);
// Compute intermediate values.
var m1 = blockMult(Matrix.add(a11, a22), Matrix.add(b11, b22), halfRows, halfCols);
var m2 = blockMult(Matrix.add(a21, a22), b11, halfRows, halfCols);
var m3 = blockMult(a11, Matrix.sub(b12, b22), halfRows, halfCols);
var m4 = blockMult(a22, Matrix.sub(b21, b11), halfRows, halfCols);
var m5 = blockMult(Matrix.add(a11, a12), b22, halfRows, halfCols);
var m6 = blockMult(Matrix.sub(a21, a11), Matrix.add(b11, b12), halfRows, halfCols);
var m7 = blockMult(Matrix.sub(a12, a22), Matrix.add(b21, b22), halfRows, halfCols);
// Combine intermediate values into the output.
var c11 = Matrix.add(m1, m4);
c11.sub(m5);
c11.add(m7);
var c12 = Matrix.add(m3, m5);
var c21 = Matrix.add(m2, m4);
var c22 = Matrix.sub(m1, m2);
c22.add(m3);
c22.add(m6);
//Crop output to the desired size (undo dynamic padding).
var resultat = Matrix.zeros(2 * c11.rows, 2 * c11.columns);
resultat = resultat.setSubMatrix(c11, 0, 0);
resultat = resultat.setSubMatrix(c12, c11.rows, 0);
resultat = resultat.setSubMatrix(c21, 0, c11.columns);
resultat = resultat.setSubMatrix(c22, c11.rows, c11.columns);
return resultat.subMatrix(0, rows - 1, 0, cols - 1);
}
return blockMult(x, y, r, c);
}
/**
* Returns a row-by-row scaled matrix
* @param {number} [min=0] - Minimum scaled value
* @param {number} [max=1] - Maximum scaled value
* @return {Matrix} - The scaled matrix
*/
scaleRows(min, max) {
min = min === undefined ? 0 : min;
max = max === undefined ? 1 : max;
if (min >= max) {
throw new RangeError('min should be strictly smaller than max');
}
var newMatrix = this.constructor.empty(this.rows, this.columns);
for (var i = 0; i < this.rows; i++) {
var scaled = arrayUtils.scale(this.getRow(i), {min, max});
newMatrix.setRow(i, scaled);
}
return newMatrix;
}
/**
* Returns a new column-by-column scaled matrix
* @param {number} [min=0] - Minimum scaled value
* @param {number} [max=1] - Maximum scaled value
* @return {Matrix} - The new scaled matrix
* @example
* var matrix = new Matrix([[1,2],[-1,0]]);
* var scaledMatrix = matrix.scaleColumns(); // [[1,1],[0,0]]
*/
scaleColumns(min, max) {
min = min === undefined ? 0 : min;
max = max === undefined ? 1 : max;
if (min >= max) {
throw new RangeError('min should be strictly smaller than max');
}
var newMatrix = this.constructor.empty(this.rows, this.columns);
for (var i = 0; i < this.columns; i++) {
var scaled = arrayUtils.scale(this.getColumn(i), {
min: min,
max: max
});
newMatrix.setColumn(i, scaled);
}
return newMatrix;
}
/**
* Returns the Kronecker product (also known as tensor product) between this and other
* @param {Matrix} other
* @return {Matrix}
*/
kroneckerProduct(other) {
other = this.constructor.checkMatrix(other);
var m = this.rows;
var n = this.columns;
var p = other.rows;
var q = other.columns;
var result = new this.constructor[Symbol.species](m * p, n * q);
for (var i = 0; i < m; i++) {
for (var j = 0; j < n; j++) {
for (var k = 0; k < p; k++) {
for (var l = 0; l < q; l++) {
result[p * i + k][q * j + l] = this.get(i, j) * other.get(k, l);
}
}
}
}
return result;
}
/**
* Transposes the matrix and returns a new one containing the result
* @return {Matrix}
*/
transpose() {
var result = new this.constructor[Symbol.species](this.columns, this.rows);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
result.set(j, i, this.get(i, j));
}
}
return result;
}
/**
* Sorts the rows (in place)
* @param {function} compareFunction - usual Array.prototype.sort comparison function
* @return {Matrix} this
*/
sortRows(compareFunction) {
if (compareFunction === undefined) compareFunction = compareNumbers;
for (var i = 0; i < this.rows; i++) {
this.setRow(i, this.getRow(i).sort(compareFunction));
}
return this;
}
/**
* Sorts the columns (in place)
* @param {function} compareFunction - usual Array.prototype.sort comparison function
* @return {Matrix} this
*/
sortColumns(compareFunction) {
if (compareFunction === undefined) compareFunction = compareNumbers;
for (var i = 0; i < this.columns; i++) {
this.setColumn(i, this.getColumn(i).sort(compareFunction));
}
return this;
}
/**
* Returns a subset of the matrix
* @param {number} startRow - First row index
* @param {number} endRow - Last row index
* @param {number} startColumn - First column index
* @param {number} endColumn - Last column index
* @return {Matrix}
*/
subMatrix(startRow, endRow, startColumn, endColumn) {
util.checkRange(this, startRow, endRow, startColumn, endColumn);
var newMatrix = new this.constructor[Symbol.species](endRow - startRow + 1, endColumn - startColumn + 1);
for (var i = startRow; i <= endRow; i++) {
for (var j = startColumn; j <= endColumn; j++) {
newMatrix[i - startRow][j - startColumn] = this.get(i, j);
}
}
return newMatrix;
}
/**
* Returns a subset of the matrix based on an array of row indices
* @param {Array} indices - Array containing the row indices
* @param {number} [startColumn = 0] - First column index
* @param {number} [endColumn = this.columns-1] - Last column index
* @return {Matrix}
*/
subMatrixRow(indices, startColumn, endColumn) {
if (startColumn === undefined) startColumn = 0;
if (endColumn === undefined) endColumn = this.columns - 1;
if ((startColumn > endColumn) || (startColumn < 0) || (startColumn >= this.columns) || (endColumn < 0) || (endColumn >= this.columns)) {
throw new RangeError('Argument out of range');
}
var newMatrix = new this.constructor[Symbol.species](indices.length, endColumn - startColumn + 1);
for (var i = 0; i < indices.length; i++) {
for (var j = startColumn; j <= endColumn; j++) {
if (indices[i] < 0 || indices[i] >= this.rows) {
throw new RangeError('Row index out of range: ' + indices[i]);
}
newMatrix.set(i, j - startColumn, this.get(indices[i], j));
}
}
return newMatrix;
}
/**
* Returns a subset of the matrix based on an array of column indices
* @param {Array} indices - Array containing the column indices
* @param {number} [startRow = 0] - First row index
* @param {number} [endRow = this.rows-1] - Last row index
* @return {Matrix}
*/
subMatrixColumn(indices, startRow, endRow) {
if (startRow === undefined) startRow = 0;
if (endRow === undefined) endRow = this.rows - 1;
if ((startRow > endRow) || (startRow < 0) || (startRow >= this.rows) || (endRow < 0) || (endRow >= this.rows)) {
throw new RangeError('Argument out of range');
}
var newMatrix = new this.constructor[Symbol.species](endRow - startRow + 1, indices.length);
for (var i = 0; i < indices.length; i++) {
for (var j = startRow; j <= endRow; j++) {
if (indices[i] < 0 || indices[i] >= this.columns) {
throw new RangeError('Column index out of range: ' + indices[i]);
}
newMatrix.set(j - startRow, i, this.get(j, indices[i]));
}
}
return newMatrix;
}
/**
* Set a part of the matrix to the given sub-matrix
* @param {Matrix|Array< Array >} matrix - The source matrix from which to extract values.
* @param {number} startRow - The index of the first row to set
* @param {number} startColumn - The index of the first column to set
* @return {Matrix}
*/
setSubMatrix(matrix, startRow, startColumn) {
matrix = this.constructor.checkMatrix(matrix);
var endRow = startRow + matrix.rows - 1;
var endColumn = startColumn + matrix.columns - 1;
util.checkRange(this, startRow, endRow, startColumn, endColumn);
for (var i = 0; i < matrix.rows; i++) {
for (var j = 0; j < matrix.columns; j++) {
this[startRow + i][startColumn + j] = matrix.get(i, j);
}
}
return this;
}
/**
* Return a new matrix based on a selection of rows and columns
* @param {Array<number>} rowIndices - The row indices to select. Order matters and an index can be more than once.
* @param {Array<number>} columnIndices - The column indices to select. Order matters and an index can be use more than once.
* @return {Matrix} The new matrix
*/
selection(rowIndices, columnIndices) {
var indices = util.checkIndices(this, rowIndices, columnIndices);
var newMatrix = new this.constructor[Symbol.species](rowIndices.length, columnIndices.length);
for (var i = 0; i < indices.row.length; i++) {
var rowIndex = indices.row[i];
for (var j = 0; j < indices.column.length; j++) {
var columnIndex = indices.column[j];
newMatrix[i][j] = this.get(rowIndex, columnIndex);
}
}
return newMatrix;
}
/**
* Returns the trace of the matrix (sum of the diagonal elements)
* @return {number}
*/
trace() {
var min = Math.min(this.rows, this.columns);
var trace = 0;
for (var i = 0; i < min; i++) {
trace += this.get(i, i);
}
return trace;
}
/*
Matrix views
*/
/**
* Returns a view of the transposition of the matrix
* @return {MatrixTransposeView}
*/
transposeView() {
return new MLMatrixTransposeView(this);
}
/**
* Returns a view of the row vector with the given index
* @param {number} row - row index of the vector
* @return {MatrixRowView}
*/
rowView(row) {
util.checkRowIndex(this, row);
return new MLMatrixRowView(this, row);
}
/**
* Returns a view of the column vector with the given index
* @param {number} column - column index of the vector
* @return {MatrixColumnView}
*/
columnView(column) {
util.checkColumnIndex(this, column);
return new MLMatrixColumnView(this, column);
}
/**
* Returns a view of the matrix flipped in the row axis
* @return {MatrixFlipRowView}
*/
flipRowView() {
return new MLMatrixFlipRowView(this);
}
/**
* Returns a view of the matrix flipped in the column axis
* @return {MatrixFlipColumnView}
*/
flipColumnView() {
return new MLMatrixFlipColumnView(this);
}
/**
* Returns a view of a submatrix giving the index boundaries
* @param {number} startRow - first row index of the submatrix
* @param {number} endRow - last row index of the submatrix
* @param {number} startColumn - first column index of the submatrix
* @param {number} endColumn - last column index of the submatrix
* @return {MatrixSubView}
*/
subMatrixView(startRow, endRow, startColumn, endColumn) {
return new MLMatrixSubView(this, startRow, endRow, startColumn, endColumn);
}
/**
* Returns a view of the cross of the row indices and the column indices
* @example
* // resulting vector is [[2], [2]]
* var matrix = new Matrix([[1,2,3], [4,5,6]]).selectionView([0, 0], [1])
* @param {Array<number>} rowIndices
* @param {Array<number>} columnIndices
* @return {MatrixSelectionView}
*/
selectionView(rowIndices, columnIndices) {
return new MLMatrixSelectionView(this, rowIndices, columnIndices);
}
/**
* Calculates and returns the determinant of a matrix as a Number
* @example
* new Matrix([[1,2,3], [4,5,6]]).det()
* @return {number}
*/
det() {
if (this.isSquare()) {
var a, b, c, d;
if (this.columns === 2) {
// 2 x 2 matrix
a = this.get(0, 0);
b = this.get(0, 1);
c = this.get(1, 0);
d = this.get(1, 1);
return a * d - (b * c);
} else if (this.columns === 3) {
// 3 x 3 matrix
var subMatrix0, subMatrix1, subMatrix2;
subMatrix0 = this.selectionView([1, 2], [1, 2]);
subMatrix1 = this.selectionView([1, 2], [0, 2]);
subMatrix2 = this.selectionView([1, 2], [0, 1]);
a = this.get(0, 0);
b = this.get(0, 1);
c = this.get(0, 2);
return a * subMatrix0.det() - b * subMatrix1.det() + c * subMatrix2.det();
} else {
// general purpose determinant using the LU decomposition
return new LuDecomposition(this).determinant;
}
} else {
throw Error('Determinant can only be calculated for a square matrix.');
}
}
/**
* Returns inverse of a matrix if it exists or the pseudoinverse
* @param {number} threshold - threshold for taking inverse of singular values (default = 1e-15)
* @return {Matrix} the (pseudo)inverted matrix.
*/
pseudoInverse(threshold) {
if (threshold === undefined) threshold = Number.EPSILON;
var svdSolution = new SvDecomposition(this, {autoTranspose: true});
var U = svdSolution.leftSingularVectors;
var V = svdSolution.rightSingularVectors;
var s = svdSolution.diagonal;
for (var i = 0; i < s.length; i++) {
if (Math.abs(s[i]) > threshold) {
s[i] = 1.0 / s[i];
} else {
s[i] = 0.0;
}
}
// convert list to diagonal
s = this.constructor[Symbol.species].diag(s);
return V.mmul(s.mmul(U.transposeView()));
}
}
Matrix.prototype.klass = 'Matrix';
/**
* @private
* Check that two matrices have the same dimensions
* @param {Matrix} matrix
* @param {Matrix} otherMatrix
*/
function checkDimensions(matrix, otherMatrix) { // eslint-disable-line no-unused-vars
if (matrix.rows !== otherMatrix.rows ||
matrix.columns !== otherMatrix.columns) {
throw new RangeError('Matrices dimensions must be equal');
}
}
function compareNumbers(a, b) {
return a - b;
}
/*
Synonyms
*/
Matrix.random = Matrix.rand;
Matrix.diagonal = Matrix.diag;
Matrix.prototype.diagonal = Matrix.prototype.diag;
Matrix.identity = Matrix.eye;
Matrix.prototype.negate = Matrix.prototype.neg;
Matrix.prototype.tensorProduct = Matrix.prototype.kroneckerProduct;
Matrix.prototype.determinant = Matrix.prototype.det;
/*
Add dynamically instance and static methods for mathematical operations
*/
var inplaceOperator = `
(function %name%(value) {
if (typeof value === 'number') return this.%name%S(value);
return this.%name%M(value);
})
`;
var inplaceOperatorScalar = `
(function %name%S(value) {
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) %op% value);
}
}
return this;
})
`;
var inplaceOperatorMatrix = `
(function %name%M(matrix) {
matrix = this.constructor.checkMatrix(matrix);
checkDimensions(this, matrix);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, this.get(i, j) %op% matrix.get(i, j));
}
}
return this;
})
`;
var staticOperator = `
(function %name%(matrix, value) {
var newMatrix = new this[Symbol.species](matrix);
return newMatrix.%name%(value);
})
`;
var inplaceMethod = `
(function %name%() {
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, %method%(this.get(i, j)));
}
}
return this;
})
`;
var staticMethod = `
(function %name%(matrix) {
var newMatrix = new this[Symbol.species](matrix);
return newMatrix.%name%();
})
`;
var inplaceMethodWithArgs = `
(function %name%(%args%) {
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, %method%(this.get(i, j), %args%));
}
}
return this;
})
`;
var staticMethodWithArgs = `
(function %name%(matrix, %args%) {
var newMatrix = new this[Symbol.species](matrix);
return newMatrix.%name%(%args%);
})
`;
var inplaceMethodWithOneArgScalar = `
(function %name%S(value) {
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, %method%(this.get(i, j), value));
}
}
return this;
})
`;
var inplaceMethodWithOneArgMatrix = `
(function %name%M(matrix) {
matrix = this.constructor.checkMatrix(matrix);
checkDimensions(this, matrix);
for (var i = 0; i < this.rows; i++) {
for (var j = 0; j < this.columns; j++) {
this.set(i, j, %method%(this.get(i, j), matrix.get(i, j)));
}
}
return this;
})
`;
var inplaceMethodWithOneArg = `
(function %name%(value) {
if (typeof value === 'number') return this.%name%S(value);
return this.%name%M(value);
})
`;
var staticMethodWithOneArg = staticMethodWithArgs;
var operators = [
// Arithmetic operators
['+', 'add'],
['-', 'sub', 'subtract'],
['*', 'mul', 'multiply'],
['/', 'div', 'divide'],
['%', 'mod', 'modulus'],
// Bitwise operators
['&', 'and'],
['|', 'or'],
['^', 'xor'],
['<<', 'leftShift'],
['>>', 'signPropagatingRightShift'],
['>>>', 'rightShift', 'zeroFillRightShift']
];
var i;
for (var operator of operators) {
var inplaceOp = eval(fillTemplateFunction(inplaceOperator, {name: operator[1], op: operator[0]}));
var inplaceOpS = eval(fillTemplateFunction(inplaceOperatorScalar, {name: operator[1] + 'S', op: operator[0]}));
var inplaceOpM = eval(fillTemplateFunction(inplaceOperatorMatrix, {name: operator[1] + 'M', op: operator[0]}));
var staticOp = eval(fillTemplateFunction(staticOperator, {name: operator[1]}));
for (i = 1; i < operator.length; i++) {
Matrix.prototype[operator[i]] = inplaceOp;
Matrix.prototype[operator[i] + 'S'] = inplaceOpS;
Matrix.prototype[operator[i] + 'M'] = inplaceOpM;
Matrix[operator[i]] = staticOp;
}
}
var methods = [
['~', 'not']
];
[
'abs', 'acos', 'acosh', 'asin', 'asinh', 'atan', 'atanh', 'cbrt', 'ceil',
'clz32', 'cos', 'cosh', 'exp', 'expm1', 'floor', 'fround', 'log', 'log1p',
'log10', 'log2', 'round', 'sign', 'sin', 'sinh', 'sqrt', 'tan', 'tanh', 'trunc'
].forEach(function (mathMethod) {
methods.push(['Math.' + mathMethod, mathMethod]);
});
for (var method of methods) {
var inplaceMeth = eval(fillTemplateFunction(inplaceMethod, {name: method[1], method: method[0]}));
var staticMeth = eval(fillTemplateFunction(staticMethod, {name: method[1]}));
for (i = 1; i < method.length; i++) {
Matrix.prototype[method[i]] = inplaceMeth;
Matrix[method[i]] = staticMeth;
}
}
var methodsWithArgs = [
['Math.pow', 1, 'pow']
];
for (var methodWithArg of methodsWithArgs) {
var args = 'arg0';
for (i = 1; i < methodWithArg[1]; i++) {
args += `, arg${i}`;
}
if (methodWithArg[1] !== 1) {
var inplaceMethWithArgs = eval(fillTemplateFunction(inplaceMethodWithArgs, {
name: methodWithArg[2],
method: methodWithArg[0],
args: args
}));
var staticMethWithArgs = eval(fillTemplateFunction(staticMethodWithArgs, {name: methodWithArg[2], args: args}));
for (i = 2; i < methodWithArg.length; i++) {
Matrix.prototype[methodWithArg[i]] = inplaceMethWithArgs;
Matrix[methodWithArg[i]] = staticMethWithArgs;
}
} else {
var tmplVar = {
name: methodWithArg[2],
args: args,
method: methodWithArg[0]
};
var inplaceMethod2 = eval(fillTemplateFunction(inplaceMethodWithOneArg, tmplVar));
var inplaceMethodS = eval(fillTemplateFunction(inplaceMethodWithOneArgScalar, tmplVar));
var inplaceMethodM = eval(fillTemplateFunction(inplaceMethodWithOneArgMatrix, tmplVar));
var staticMethod2 = eval(fillTemplateFunction(staticMethodWithOneArg, tmplVar));
for (i = 2; i < methodWithArg.length; i++) {
Matrix.prototype[methodWithArg[i]] = inplaceMethod2;
Matrix.prototype[methodWithArg[i] + 'M'] = inplaceMethodM;
Matrix.prototype[methodWithArg[i] + 'S'] = inplaceMethodS;
Matrix[methodWithArg[i]] = staticMethod2;
}
}
}
function fillTemplateFunction(template, values) {
for (var value in values) {
template = template.replace(new RegExp('%' + value + '%', 'g'), values[value]);
}
return template;
}
return Matrix;
}
}
// ml-matrix src/views/base
let MLMatrixBaseView;
{
let abstractMatrix = MLMatrixAbstractMatrix;
let Matrix = MLMatrixMatrix;
class BaseView extends abstractMatrix() {
constructor(matrix, rows, columns) {
super();
this.matrix = matrix;
this.rows = rows;
this.columns = columns;
}
static get [Symbol.species]() {
return Matrix.Matrix;
}
}
MLMatrixBaseView = BaseView;
}
// ml-matrix src/views/column.js
let MLMatrixColumnView;
{
let BaseView = MLMatrixBaseView;
class MatrixColumnView extends BaseView {
constructor(matrix, column) {
super(matrix, matrix.rows, 1);
this.column = column;
}
set(rowIndex, columnIndex, value) {
this.matrix.set(rowIndex, this.column, value);
return this;
}
get(rowIndex) {
return this.matrix.get(rowIndex, this.column);
}
}
MLMatrixColumnView = MatrixColumnView;
}
// ml-matrix src/views/flipColumn.js
let MLMatrixFlipColumnView;
{
let BaseView = MLMatrixBaseView
class MatrixFlipColumnView extends BaseView {
constructor(matrix) {
super(matrix, matrix.rows, matrix.columns);
}
set(rowIndex, columnIndex, value) {
this.matrix.set(rowIndex, this.columns - columnIndex - 1, value);
return this;
}
get(rowIndex, columnIndex) {
return this.matrix.get(rowIndex, this.columns - columnIndex - 1);
}
}
MLMatrixFlipColumnView = MatrixFlipColumnView;
}
// ml-matrix src/views/flipRow.js
let MLMatrixFlipRowView;
{
let BaseView = MLMatrixBaseView
class MatrixFlipRowView extends BaseView {
constructor(matrix) {
super(matrix, matrix.rows, matrix.columns);
}
set(rowIndex, columnIndex, value) {
this.matrix.set(this.rows - rowIndex - 1, columnIndex, value);
return this;
}
get(rowIndex, columnIndex) {
return this.matrix.get(this.rows - rowIndex - 1, columnIndex);
}
}
MLMatrixFlipRowView = MatrixFlipRowView;
}
// ml-matrix src/views/row.js
let MLMatrixRowView;
{
let BaseView = MLMatrixBaseView;
class MatrixRowView extends BaseView {
constructor(matrix, row) {
super(matrix, 1, matrix.columns);
this.row = row;
}
set(rowIndex, columnIndex, value) {
this.matrix.set(this.row, columnIndex, value);
return this;
}
get(rowIndex, columnIndex) {
return this.matrix.get(this.row, columnIndex);
}
}
MLMatrixRowView = MatrixRowView;
}
// ml-matrix src/views/selection.js
let MLMatrixSelectionView;
{
let BaseView = MLMatrixBaseView;
let util = MLMatrixUtil;
class MatrixSelectionView extends BaseView {
constructor(matrix, rowIndices, columnIndices) {
var indices = util.checkIndices(matrix, rowIndices, columnIndices);
super(matrix, indices.row.length, indices.column.length);
this.rowIndices = indices.row;
this.columnIndices = indices.column;
}
set(rowIndex, columnIndex, value) {
this.matrix.set(this.rowIndices[rowIndex], this.columnIndices[columnIndex], value);
return this;
}
get(rowIndex, columnIndex) {
return this.matrix.get(this.rowIndices[rowIndex], this.columnIndices[columnIndex]);
}
}
MLMatrixSelectionView = MatrixSelectionView;
}
// ml-matrix src/views/sub.js
let MLMatrixSubView;
{
let BaseView = MLMatrixBaseView;
let util = MLMatrixUtil;
class MatrixSubView extends BaseView {
constructor(matrix, startRow, endRow, startColumn, endColumn) {
util.checkRange(matrix, startRow, endRow, startColumn, endColumn);
super(matrix, endRow - startRow + 1, endColumn - startColumn + 1);
this.startRow = startRow;
this.startColumn = startColumn;
}
set(rowIndex, columnIndex, value) {
this.matrix.set(this.startRow + rowIndex, this.startColumn + columnIndex, value);
return this;
}
get(rowIndex, columnIndex) {
return this.matrix.get(this.startRow + rowIndex, this.startColumn + columnIndex);
}
}
MLMatrixSubView = MatrixSubView;
}
// ml-matrix src/views/transpose.js
let MLMatrixTransposeView;
{
let BaseView = MLMatrixBaseView;
class MatrixTransposeView extends BaseView {
constructor(matrix) {
super(matrix, matrix.columns, matrix.rows);
}
set(rowIndex, columnIndex, value) {
this.matrix.set(columnIndex, rowIndex, value);
return this;
}
get(rowIndex, columnIndex) {
return this.matrix.get(columnIndex, rowIndex);
}
}
MLMatrixTransposeView = MatrixTransposeView;
}
// mlmatrix src/matrix.js
{
let abstractMatrix = MLMatrixAbstractMatrix;
let util = MLMatrixUtil;
class Matrix extends abstractMatrix(Array) {
constructor(nRows, nColumns) {
var i;
if (arguments.length === 1 && typeof nRows === 'number') {
return new Array(nRows);
}
if (Matrix.isMatrix(nRows)) {
return nRows.clone();
} else if (Number.isInteger(nRows) && nRows > 0) { // Create an empty matrix
super(nRows);
if (Number.isInteger(nColumns) && nColumns > 0) {
for (i = 0; i < nRows; i++) {
this[i] = new Array(nColumns);
}
} else {
throw new TypeError('nColumns must be a positive integer');
}
} else if (Array.isArray(nRows)) { // Copy the values from the 2D array
const matrix = nRows;
nRows = matrix.length;
nColumns = matrix[0].length;
if (typeof nColumns !== 'number' || nColumns === 0) {
throw new TypeError('Data must be a 2D array with at least one element');
}
super(nRows);
for (i = 0; i < nRows; i++) {
if (matrix[i].length !== nColumns) {
throw new RangeError('Inconsistent array dimensions');
}
this[i] = [].concat(matrix[i]);
}
} else {
throw new TypeError('First argument must be a positive number or an array');
}
this.rows = nRows;
this.columns = nColumns;
return this;
}
set(rowIndex, columnIndex, value) {
this[rowIndex][columnIndex] = value;
return this;
}
get(rowIndex, columnIndex) {
return this[rowIndex][columnIndex];
}
/**
* Creates an exact and independent copy of the matrix
* @return {Matrix}
*/
clone() {
var newMatrix = new this.constructor[Symbol.species](this.rows, this.columns);
for (var row = 0; row < this.rows; row++) {
for (var column = 0; column < this.columns; column++) {
newMatrix.set(row, column, this.get(row, column));
}
}
return newMatrix;
}
/**
* Removes a row from the given index
* @param {number} index - Row index
* @return {Matrix} this
*/
removeRow(index) {
util.checkRowIndex(this, index);
if (this.rows === 1) {
throw new RangeError('A matrix cannot have less than one row');
}
this.splice(index, 1);
this.rows -= 1;
return this;
}
/**
* Adds a row at the given index
* @param {number} [index = this.rows] - Row index
* @param {Array|Matrix} array - Array or vector
* @return {Matrix} this
*/
addRow(index, array) {
if (array === undefined) {
array = index;
index = this.rows;
}
util.checkRowIndex(this, index, true);
array = util.checkRowVector(this, array, true);
this.splice(index, 0, array);
this.rows += 1;
return this;
}
/**
* Removes a column from the given index
* @param {number} index - Column index
* @return {Matrix} this
*/
removeColumn(index) {
util.checkColumnIndex(this, index);
if (this.columns === 1) {
throw new RangeError('A matrix cannot have less than one column');
}
for (var i = 0; i < this.rows; i++) {
this[i].splice(index, 1);
}
this.columns -= 1;
return this;
}
/**
* Adds a column at the given index
* @param {number} [index = this.columns] - Column index
* @param {Array|Matrix} array - Array or vector
* @return {Matrix} this
*/
addColumn(index, array) {
if (typeof array === 'undefined') {
array = index;
index = this.columns;
}
util.checkColumnIndex(this, index, true);
array = util.checkColumnVector(this, array);
for (var i = 0; i < this.rows; i++) {
this[i].splice(index, 0, array[i]);
}
this.columns += 1;
return this;
}
}
MLMatrixMatrix.Matrix = Matrix;
Matrix.abstractMatrix = abstractMatrix;
}
// ml-matrix src/dc/cholesky.js
let MLMatrixDCCholesky = {};
{
let Matrix = MLMatrixMatrix.Matrix;
function CholeskyDecomposition(value) {
if (!(this instanceof CholeskyDecomposition)) {
return new CholeskyDecomposition(value);
}
value = Matrix.checkMatrix(value);
if (!value.isSymmetric()) {
throw new Error('Matrix is not symmetric');
}
var a = value,
dimension = a.rows,
l = new Matrix(dimension, dimension),
positiveDefinite = true,
i, j, k;
for (j = 0; j < dimension; j++) {
var Lrowj = l[j];
var d = 0;
for (k = 0; k < j; k++) {
var Lrowk = l[k];
var s = 0;
for (i = 0; i < k; i++) {
s += Lrowk[i] * Lrowj[i];
}
Lrowj[k] = s = (a[j][k] - s) / l[k][k];
d = d + s * s;
}
d = a[j][j] - d;
positiveDefinite &= (d > 0);
l[j][j] = Math.sqrt(Math.max(d, 0));
for (k = j + 1; k < dimension; k++) {
l[j][k] = 0;
}
}
if (!positiveDefinite) {
throw new Error('Matrix is not positive definite');
}
this.L = l;
}
CholeskyDecomposition.prototype = {
get lowerTriangularMatrix() {
return this.L;
},
solve: function (value) {
value = Matrix.checkMatrix(value);
var l = this.L,
dimension = l.rows;
if (value.rows !== dimension) {
throw new Error('Matrix dimensions do not match');
}
var count = value.columns,
B = value.clone(),
i, j, k;
for (k = 0; k < dimension; k++) {
for (j = 0; j < count; j++) {
for (i = 0; i < k; i++) {
B[k][j] -= B[i][j] * l[k][i];
}
B[k][j] /= l[k][k];
}
}
for (k = dimension - 1; k >= 0; k--) {
for (j = 0; j < count; j++) {
for (i = k + 1; i < dimension; i++) {
B[k][j] -= B[i][j] * l[i][k];
}
B[k][j] /= l[k][k];
}
}
return B;
}
};
MLMatrixDCCholesky = CholeskyDecomposition;
}
// ml-matrix src/dc/evd.js
let MLMatrixDCEVD;
{
const Matrix = MLMatrixMatrix.Matrix;
const util = MLMatrixDCUtil;
const hypotenuse = util.hypotenuse;
const getFilled2DArray = util.getFilled2DArray;
const defaultOptions = {
assumeSymmetric: false
};
function EigenvalueDecomposition(matrix, options) {
options = Object.assign({}, defaultOptions, options);
if (!(this instanceof EigenvalueDecomposition)) {
return new EigenvalueDecomposition(matrix, options);
}
matrix = Matrix.checkMatrix(matrix);
if (!matrix.isSquare()) {
throw new Error('Matrix is not a square matrix');
}
var n = matrix.columns,
V = getFilled2DArray(n, n, 0),
d = new Array(n),
e = new Array(n),
value = matrix,
i, j;
var isSymmetric = false;
if (options.assumeSymmetric) {
isSymmetric = true;
} else {
isSymmetric = matrix.isSymmetric();
}
if (isSymmetric) {
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
V[i][j] = value.get(i, j);
}
}
tred2(n, e, d, V);
tql2(n, e, d, V);
} else {
var H = getFilled2DArray(n, n, 0),
ort = new Array(n);
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
H[i][j] = value.get(i, j);
}
}
orthes(n, H, ort, V);
hqr2(n, e, d, V, H);
}
this.n = n;
this.e = e;
this.d = d;
this.V = V;
}
EigenvalueDecomposition.prototype = {
get realEigenvalues() {
return this.d;
},
get imaginaryEigenvalues() {
return this.e;
},
get eigenvectorMatrix() {
if (!Matrix.isMatrix(this.V)) {
this.V = new Matrix(this.V);
}
return this.V;
},
get diagonalMatrix() {
var n = this.n,
e = this.e,
d = this.d,
X = new Matrix(n, n),
i, j;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
X[i][j] = 0;
}
X[i][i] = d[i];
if (e[i] > 0) {
X[i][i + 1] = e[i];
} else if (e[i] < 0) {
X[i][i - 1] = e[i];
}
}
return X;
}
};
function tred2(n, e, d, V) {
var f, g, h, i, j, k,
hh, scale;
for (j = 0; j < n; j++) {
d[j] = V[n - 1][j];
}
for (i = n - 1; i > 0; i--) {
scale = 0;
h = 0;
for (k = 0; k < i; k++) {
scale = scale + Math.abs(d[k]);
}
if (scale === 0) {
e[i] = d[i - 1];
for (j = 0; j < i; j++) {
d[j] = V[i - 1][j];
V[i][j] = 0;
V[j][i] = 0;
}
} else {
for (k = 0; k < i; k++) {
d[k] /= scale;
h += d[k] * d[k];
}
f = d[i - 1];
g = Math.sqrt(h);
if (f > 0) {
g = -g;
}
e[i] = scale * g;
h = h - f * g;
d[i - 1] = f - g;
for (j = 0; j < i; j++) {
e[j] = 0;
}
for (j = 0; j < i; j++) {
f = d[j];
V[j][i] = f;
g = e[j] + V[j][j] * f;
for (k = j + 1; k <= i - 1; k++) {
g += V[k][j] * d[k];
e[k] += V[k][j] * f;
}
e[j] = g;
}
f = 0;
for (j = 0; j < i; j++) {
e[j] /= h;
f += e[j] * d[j];
}
hh = f / (h + h);
for (j = 0; j < i; j++) {
e[j] -= hh * d[j];
}
for (j = 0; j < i; j++) {
f = d[j];
g = e[j];
for (k = j; k <= i - 1; k++) {
V[k][j] -= (f * e[k] + g * d[k]);
}
d[j] = V[i - 1][j];
V[i][j] = 0;
}
}
d[i] = h;
}
for (i = 0; i < n - 1; i++) {
V[n - 1][i] = V[i][i];
V[i][i] = 1;
h = d[i + 1];
if (h !== 0) {
for (k = 0; k <= i; k++) {
d[k] = V[k][i + 1] / h;
}
for (j = 0; j <= i; j++) {
g = 0;
for (k = 0; k <= i; k++) {
g += V[k][i + 1] * V[k][j];
}
for (k = 0; k <= i; k++) {
V[k][j] -= g * d[k];
}
}
}
for (k = 0; k <= i; k++) {
V[k][i + 1] = 0;
}
}
for (j = 0; j < n; j++) {
d[j] = V[n - 1][j];
V[n - 1][j] = 0;
}
V[n - 1][n - 1] = 1;
e[0] = 0;
}
function tql2(n, e, d, V) {
var g, h, i, j, k, l, m, p, r,
dl1, c, c2, c3, el1, s, s2,
iter;
for (i = 1; i < n; i++) {
e[i - 1] = e[i];
}
e[n - 1] = 0;
var f = 0,
tst1 = 0,
eps = Math.pow(2, -52);
for (l = 0; l < n; l++) {
tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
m = l;
while (m < n) {
if (Math.abs(e[m]) <= eps * tst1) {
break;
}
m++;
}
if (m > l) {
iter = 0;
do {
iter = iter + 1;
g = d[l];
p = (d[l + 1] - g) / (2 * e[l]);
r = hypotenuse(p, 1);
if (p < 0) {
r = -r;
}
d[l] = e[l] / (p + r);
d[l + 1] = e[l] * (p + r);
dl1 = d[l + 1];
h = g - d[l];
for (i = l + 2; i < n; i++) {
d[i] -= h;
}
f = f + h;
p = d[m];
c = 1;
c2 = c;
c3 = c;
el1 = e[l + 1];
s = 0;
s2 = 0;
for (i = m - 1; i >= l; i--) {
c3 = c2;
c2 = c;
s2 = s;
g = c * e[i];
h = c * p;
r = hypotenuse(p, e[i]);
e[i + 1] = s * r;
s = e[i] / r;
c = p / r;
p = c * d[i] - s * g;
d[i + 1] = h + s * (c * g + s * d[i]);
for (k = 0; k < n; k++) {
h = V[k][i + 1];
V[k][i + 1] = s * V[k][i] + c * h;
V[k][i] = c * V[k][i] - s * h;
}
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
d[l] = c * p;
}
while (Math.abs(e[l]) > eps * tst1);
}
d[l] = d[l] + f;
e[l] = 0;
}
for (i = 0; i < n - 1; i++) {
k = i;
p = d[i];
for (j = i + 1; j < n; j++) {
if (d[j] < p) {
k = j;
p = d[j];
}
}
if (k !== i) {
d[k] = d[i];
d[i] = p;
for (j = 0; j < n; j++) {
p = V[j][i];
V[j][i] = V[j][k];
V[j][k] = p;
}
}
}
}
function orthes(n, H, ort, V) {
var low = 0,
high = n - 1,
f, g, h, i, j, m,
scale;
for (m = low + 1; m <= high - 1; m++) {
scale = 0;
for (i = m; i <= high; i++) {
scale = scale + Math.abs(H[i][m - 1]);
}
if (scale !== 0) {
h = 0;
for (i = high; i >= m; i--) {
ort[i] = H[i][m - 1] / scale;
h += ort[i] * ort[i];
}
g = Math.sqrt(h);
if (ort[m] > 0) {
g = -g;
}
h = h - ort[m] * g;
ort[m] = ort[m] - g;
for (j = m; j < n; j++) {
f = 0;
for (i = high; i >= m; i--) {
f += ort[i] * H[i][j];
}
f = f / h;
for (i = m; i <= high; i++) {
H[i][j] -= f * ort[i];
}
}
for (i = 0; i <= high; i++) {
f = 0;
for (j = high; j >= m; j--) {
f += ort[j] * H[i][j];
}
f = f / h;
for (j = m; j <= high; j++) {
H[i][j] -= f * ort[j];
}
}
ort[m] = scale * ort[m];
H[m][m - 1] = scale * g;
}
}
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
V[i][j] = (i === j ? 1 : 0);
}
}
for (m = high - 1; m >= low + 1; m--) {
if (H[m][m - 1] !== 0) {
for (i = m + 1; i <= high; i++) {
ort[i] = H[i][m - 1];
}
for (j = m; j <= high; j++) {
g = 0;
for (i = m; i <= high; i++) {
g += ort[i] * V[i][j];
}
g = (g / ort[m]) / H[m][m - 1];
for (i = m; i <= high; i++) {
V[i][j] += g * ort[i];
}
}
}
}
}
function hqr2(nn, e, d, V, H) {
var n = nn - 1,
low = 0,
high = nn - 1,
eps = Math.pow(2, -52),
exshift = 0,
norm = 0,
p = 0,
q = 0,
r = 0,
s = 0,
z = 0,
iter = 0,
i, j, k, l, m, t, w, x, y,
ra, sa, vr, vi,
notlast, cdivres;
for (i = 0; i < nn; i++) {
if (i < low || i > high) {
d[i] = H[i][i];
e[i] = 0;
}
for (j = Math.max(i - 1, 0); j < nn; j++) {
norm = norm + Math.abs(H[i][j]);
}
}
while (n >= low) {
l = n;
while (l > low) {
s = Math.abs(H[l - 1][l - 1]) + Math.abs(H[l][l]);
if (s === 0) {
s = norm;
}
if (Math.abs(H[l][l - 1]) < eps * s) {
break;
}
l--;
}
if (l === n) {
H[n][n] = H[n][n] + exshift;
d[n] = H[n][n];
e[n] = 0;
n--;
iter = 0;
} else if (l === n - 1) {
w = H[n][n - 1] * H[n - 1][n];
p = (H[n - 1][n - 1] - H[n][n]) / 2;
q = p * p + w;
z = Math.sqrt(Math.abs(q));
H[n][n] = H[n][n] + exshift;
H[n - 1][n - 1] = H[n - 1][n - 1] + exshift;
x = H[n][n];
if (q >= 0) {
z = (p >= 0) ? (p + z) : (p - z);
d[n - 1] = x + z;
d[n] = d[n - 1];
if (z !== 0) {
d[n] = x - w / z;
}
e[n - 1] = 0;
e[n] = 0;
x = H[n][n - 1];
s = Math.abs(x) + Math.abs(z);
p = x / s;
q = z / s;
r = Math.sqrt(p * p + q * q);
p = p / r;
q = q / r;
for (j = n - 1; j < nn; j++) {
z = H[n - 1][j];
H[n - 1][j] = q * z + p * H[n][j];
H[n][j] = q * H[n][j] - p * z;
}
for (i = 0; i <= n; i++) {
z = H[i][n - 1];
H[i][n - 1] = q * z + p * H[i][n];
H[i][n] = q * H[i][n] - p * z;
}
for (i = low; i <= high; i++) {
z = V[i][n - 1];
V[i][n - 1] = q * z + p * V[i][n];
V[i][n] = q * V[i][n] - p * z;
}
} else {
d[n - 1] = x + p;
d[n] = x + p;
e[n - 1] = z;
e[n] = -z;
}
n = n - 2;
iter = 0;
} else {
x = H[n][n];
y = 0;
w = 0;
if (l < n) {
y = H[n - 1][n - 1];
w = H[n][n - 1] * H[n - 1][n];
}
if (iter === 10) {
exshift += x;
for (i = low; i <= n; i++) {
H[i][i] -= x;
}
s = Math.abs(H[n][n - 1]) + Math.abs(H[n - 1][n - 2]);
x = y = 0.75 * s;
w = -0.4375 * s * s;
}
if (iter === 30) {
s = (y - x) / 2;
s = s * s + w;
if (s > 0) {
s = Math.sqrt(s);
if (y < x) {
s = -s;
}
s = x - w / ((y - x) / 2 + s);
for (i = low; i <= n; i++) {
H[i][i] -= s;
}
exshift += s;
x = y = w = 0.964;
}
}
iter = iter + 1;
m = n - 2;
while (m >= l) {
z = H[m][m];
r = x - z;
s = y - z;
p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
q = H[m + 1][m + 1] - z - r - s;
r = H[m + 2][m + 1];
s = Math.abs(p) + Math.abs(q) + Math.abs(r);
p = p / s;
q = q / s;
r = r / s;
if (m === l) {
break;
}
if (Math.abs(H[m][m - 1]) * (Math.abs(q) + Math.abs(r)) < eps * (Math.abs(p) * (Math.abs(H[m - 1][m - 1]) + Math.abs(z) + Math.abs(H[m + 1][m + 1])))) {
break;
}
m--;
}
for (i = m + 2; i <= n; i++) {
H[i][i - 2] = 0;
if (i > m + 2) {
H[i][i - 3] = 0;
}
}
for (k = m; k <= n - 1; k++) {
notlast = (k !== n - 1);
if (k !== m) {
p = H[k][k - 1];
q = H[k + 1][k - 1];
r = (notlast ? H[k + 2][k - 1] : 0);
x = Math.abs(p) + Math.abs(q) + Math.abs(r);
if (x !== 0) {
p = p / x;
q = q / x;
r = r / x;
}
}
if (x === 0) {
break;
}
s = Math.sqrt(p * p + q * q + r * r);
if (p < 0) {
s = -s;
}
if (s !== 0) {
if (k !== m) {
H[k][k - 1] = -s * x;
} else if (l !== m) {
H[k][k - 1] = -H[k][k - 1];
}
p = p + s;
x = p / s;
y = q / s;
z = r / s;
q = q / p;
r = r / p;
for (j = k; j < nn; j++) {
p = H[k][j] + q * H[k + 1][j];
if (notlast) {
p = p + r * H[k + 2][j];
H[k + 2][j] = H[k + 2][j] - p * z;
}
H[k][j] = H[k][j] - p * x;
H[k + 1][j] = H[k + 1][j] - p * y;
}
for (i = 0; i <= Math.min(n, k + 3); i++) {
p = x * H[i][k] + y * H[i][k + 1];
if (notlast) {
p = p + z * H[i][k + 2];
H[i][k + 2] = H[i][k + 2] - p * r;
}
H[i][k] = H[i][k] - p;
H[i][k + 1] = H[i][k + 1] - p * q;
}
for (i = low; i <= high; i++) {
p = x * V[i][k] + y * V[i][k + 1];
if (notlast) {
p = p + z * V[i][k + 2];
V[i][k + 2] = V[i][k + 2] - p * r;
}
V[i][k] = V[i][k] - p;
V[i][k + 1] = V[i][k + 1] - p * q;
}
}
}
}
}
if (norm === 0) {
return;
}
for (n = nn - 1; n >= 0; n--) {
p = d[n];
q = e[n];
if (q === 0) {
l = n;
H[n][n] = 1;
for (i = n - 1; i >= 0; i--) {
w = H[i][i] - p;
r = 0;
for (j = l; j <= n; j++) {
r = r + H[i][j] * H[j][n];
}
if (e[i] < 0) {
z = w;
s = r;
} else {
l = i;
if (e[i] === 0) {
H[i][n] = (w !== 0) ? (-r / w) : (-r / (eps * norm));
} else {
x = H[i][i + 1];
y = H[i + 1][i];
q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
t = (x * s - z * r) / q;
H[i][n] = t;
H[i + 1][n] = (Math.abs(x) > Math.abs(z)) ? ((-r - w * t) / x) : ((-s - y * t) / z);
}
t = Math.abs(H[i][n]);
if ((eps * t) * t > 1) {
for (j = i; j <= n; j++) {
H[j][n] = H[j][n] / t;
}
}
}
}
} else if (q < 0) {
l = n - 1;
if (Math.abs(H[n][n - 1]) > Math.abs(H[n - 1][n])) {
H[n - 1][n - 1] = q / H[n][n - 1];
H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1];
} else {
cdivres = cdiv(0, -H[n - 1][n], H[n - 1][n - 1] - p, q);
H[n - 1][n - 1] = cdivres[0];
H[n - 1][n] = cdivres[1];
}
H[n][n - 1] = 0;
H[n][n] = 1;
for (i = n - 2; i >= 0; i--) {
ra = 0;
sa = 0;
for (j = l; j <= n; j++) {
ra = ra + H[i][j] * H[j][n - 1];
sa = sa + H[i][j] * H[j][n];
}
w = H[i][i] - p;
if (e[i] < 0) {
z = w;
r = ra;
s = sa;
} else {
l = i;
if (e[i] === 0) {
cdivres = cdiv(-ra, -sa, w, q);
H[i][n - 1] = cdivres[0];
H[i][n] = cdivres[1];
} else {
x = H[i][i + 1];
y = H[i + 1][i];
vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
vi = (d[i] - p) * 2 * q;
if (vr === 0 && vi === 0) {
vr = eps * norm * (Math.abs(w) + Math.abs(q) + Math.abs(x) + Math.abs(y) + Math.abs(z));
}
cdivres = cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi);
H[i][n - 1] = cdivres[0];
H[i][n] = cdivres[1];
if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x;
H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x;
} else {
cdivres = cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q);
H[i + 1][n - 1] = cdivres[0];
H[i + 1][n] = cdivres[1];
}
}
t = Math.max(Math.abs(H[i][n - 1]), Math.abs(H[i][n]));
if ((eps * t) * t > 1) {
for (j = i; j <= n; j++) {
H[j][n - 1] = H[j][n - 1] / t;
H[j][n] = H[j][n] / t;
}
}
}
}
}
}
for (i = 0; i < nn; i++) {
if (i < low || i > high) {
for (j = i; j < nn; j++) {
V[i][j] = H[i][j];
}
}
}
for (j = nn - 1; j >= low; j--) {
for (i = low; i <= high; i++) {
z = 0;
for (k = low; k <= Math.min(j, high); k++) {
z = z + V[i][k] * H[k][j];
}
V[i][j] = z;
}
}
}
function cdiv(xr, xi, yr, yi) {
var r, d;
if (Math.abs(yr) > Math.abs(yi)) {
r = yi / yr;
d = yr + r * yi;
return [(xr + r * xi) / d, (xi - r * xr) / d];
} else {
r = yr / yi;
d = yi + r * yr;
return [(r * xr + xi) / d, (r * xi - xr) / d];
}
}
MLMatrixDCEVD = EigenvalueDecomposition;
}
// ml-matrix src/dc/qr.js
let MLMatrixDCQR;
{
let Matrix = MLMatrixMatrix.Matrix;
let hypotenuse = MLMatrixDCUtil.hypotenuse;
function QrDecomposition(value) {
if (!(this instanceof QrDecomposition)) {
return new QrDecomposition(value);
}
value = Matrix.checkMatrix(value);
var qr = value.clone(),
m = value.rows,
n = value.columns,
rdiag = new Array(n),
i, j, k, s;
for (k = 0; k < n; k++) {
var nrm = 0;
for (i = k; i < m; i++) {
nrm = hypotenuse(nrm, qr[i][k]);
}
if (nrm !== 0) {
if (qr[k][k] < 0) {
nrm = -nrm;
}
for (i = k; i < m; i++) {
qr[i][k] /= nrm;
}
qr[k][k] += 1;
for (j = k + 1; j < n; j++) {
s = 0;
for (i = k; i < m; i++) {
s += qr[i][k] * qr[i][j];
}
s = -s / qr[k][k];
for (i = k; i < m; i++) {
qr[i][j] += s * qr[i][k];
}
}
}
rdiag[k] = -nrm;
}
this.QR = qr;
this.Rdiag = rdiag;
}
QrDecomposition.prototype = {
solve: function (value) {
value = Matrix.checkMatrix(value);
var qr = this.QR,
m = qr.rows;
if (value.rows !== m) {
throw new Error('Matrix row dimensions must agree');
}
if (!this.isFullRank()) {
throw new Error('Matrix is rank deficient');
}
var count = value.columns;
var X = value.clone();
var n = qr.columns;
var i, j, k, s;
for (k = 0; k < n; k++) {
for (j = 0; j < count; j++) {
s = 0;
for (i = k; i < m; i++) {
s += qr[i][k] * X[i][j];
}
s = -s / qr[k][k];
for (i = k; i < m; i++) {
X[i][j] += s * qr[i][k];
}
}
}
for (k = n - 1; k >= 0; k--) {
for (j = 0; j < count; j++) {
X[k][j] /= this.Rdiag[k];
}
for (i = 0; i < k; i++) {
for (j = 0; j < count; j++) {
X[i][j] -= X[k][j] * qr[i][k];
}
}
}
return X.subMatrix(0, n - 1, 0, count - 1);
},
isFullRank: function () {
var columns = this.QR.columns;
for (var i = 0; i < columns; i++) {
if (this.Rdiag[i] === 0) {
return false;
}
}
return true;
},
get upperTriangularMatrix() {
var qr = this.QR,
n = qr.columns,
X = new Matrix(n, n),
i, j;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
if (i < j) {
X[i][j] = qr[i][j];
} else if (i === j) {
X[i][j] = this.Rdiag[i];
} else {
X[i][j] = 0;
}
}
}
return X;
},
get orthogonalMatrix() {
var qr = this.QR,
rows = qr.rows,
columns = qr.columns,
X = new Matrix(rows, columns),
i, j, k, s;
for (k = columns - 1; k >= 0; k--) {
for (i = 0; i < rows; i++) {
X[i][k] = 0;
}
X[k][k] = 1;
for (j = k; j < columns; j++) {
if (qr[k][k] !== 0) {
s = 0;
for (i = k; i < rows; i++) {
s += qr[i][k] * X[i][j];
}
s = -s / qr[k][k];
for (i = k; i < rows; i++) {
X[i][j] += s * qr[i][k];
}
}
}
}
return X;
}
};
MLMatrixDCQR = QrDecomposition;
}
// ml-matric src/decompositions.js
let MLMatrixDecompositions = {};
{
let Matrix = MLMatrixMatrix.Matrix;
let SingularValueDecomposition = MLMatrixDCSVD;
let EigenvalueDecomposition = MLMatrixDCEVD;
let LuDecomposition = MLMatrixDCLU;
let QrDecomposition = MLMatrixDCQR
let CholeskyDecomposition = MLMatrixDCCholesky;
function inverse(matrix) {
matrix = Matrix.checkMatrix(matrix);
return solve(matrix, Matrix.eye(matrix.rows));
}
/**
* Returns the inverse
* @memberOf Matrix
* @static
* @param {Matrix} matrix
* @return {Matrix} matrix
* @alias inv
*/
Matrix.inverse = Matrix.inv = inverse;
/**
* Returns the inverse
* @memberOf Matrix
* @static
* @param {Matrix} matrix
* @return {Matrix} matrix
* @alias inv
*/
Matrix.prototype.inverse = Matrix.prototype.inv = function () {
return inverse(this);
};
function solve(leftHandSide, rightHandSide) {
leftHandSide = Matrix.checkMatrix(leftHandSide);
rightHandSide = Matrix.checkMatrix(rightHandSide);
return leftHandSide.isSquare() ? new LuDecomposition(leftHandSide).solve(rightHandSide) : new QrDecomposition(leftHandSide).solve(rightHandSide);
}
Matrix.solve = solve;
Matrix.prototype.solve = function (other) {
return solve(this, other);
};
MLMatrixDecompositions = {
SingularValueDecomposition: SingularValueDecomposition,
SVD: SingularValueDecomposition,
EigenvalueDecomposition: EigenvalueDecomposition,
EVD: EigenvalueDecomposition,
LuDecomposition: LuDecomposition,
LU: LuDecomposition,
QrDecomposition: QrDecomposition,
QR: QrDecomposition,
CholeskyDecomposition: CholeskyDecomposition,
CHO: CholeskyDecomposition,
inverse: inverse,
solve: solve
};
}
// ml-matrix src/index.js
let MLMatrix = {};
{
MLMatrix = MLMatrixMatrix.Matrix;
MLMatrix.Decompositions = MLMatrix.DC = MLMatrixDecompositions;
}
// feedforward-neural-networks utils.js
let FeedforwardNeuralNetworksUtils;
{
let Matrix = MLMatrix;
/**
* @private
* Retrieves the sum at each row of the given matrix.
* @param {Matrix} matrix
* @return {Matrix}
*/
function sumRow(matrix) {
var sum = Matrix.zeros(matrix.rows, 1);
for (var i = 0; i < matrix.rows; ++i) {
for (var j = 0; j < matrix.columns; ++j) {
sum[i][0] += matrix[i][j];
}
}
return sum;
}
/**
* @private
* Retrieves the sum at each column of the given matrix.
* @param {Matrix} matrix
* @return {Matrix}
*/
function sumCol(matrix) {
var sum = Matrix.zeros(1, matrix.columns);
for (var i = 0; i < matrix.rows; ++i) {
for (var j = 0; j < matrix.columns; ++j) {
sum[0][j] += matrix[i][j];
}
}
return sum;
}
/**
* @private
* Method that given an array of labels(predictions), returns two dictionaries, one to transform from labels to
* numbers and other in the reverse way
* @param {Array} array
* @return {object}
*/
function dictOutputs(array) {
var inputs = {}, outputs = {}, l = array.length, index = 0;
for (var i = 0; i < l; i += 1) {
if (inputs[array[i]] === undefined) {
inputs[array[i]] = index;
outputs[index] = array[i];
index++;
}
}
return {
inputs: inputs,
outputs: outputs
};
}
FeedforwardNeuralNetworksUtils = {
dictOutputs: dictOutputs,
sumCol: sumCol,
sumRow: sumRow
};
}
// feedforward-neural-networks activationFunctions.js
let FeedforwardNeuralNetworksActivationFunctions;
{
function logistic(val) {
return 1 / (1 + Math.exp(-val));
}
function expELU(val, param) {
return val < 0 ? param * (Math.exp(val) - 1) : val;
}
function softExponential(val, param) {
if (param < 0) {
return -Math.log(1 - param * (val + param)) / param;
}
if (param > 0) {
return ((Math.exp(param * val) - 1) / param) + param;
}
return val;
}
function softExponentialPrime(val, param) {
if (param < 0) {
return 1 / (1 - param * (param + val));
} else {
return Math.exp(param * val);
}
}
const ACTIVATION_FUNCTIONS = {
'tanh': {
activation: Math.tanh,
derivate: val => 1 - (val * val)
},
'identity': {
activation: val => val,
derivate: () => 1
},
'logistic': {
activation: logistic,
derivate: val => logistic(val) * (1 - logistic(val))
},
'arctan': {
activation: Math.atan,
derivate: val => 1 / (val * val + 1)
},
'softsign': {
activation: val => val / (1 + Math.abs(val)),
derivate: val => 1 / ((1 + Math.abs(val)) * (1 + Math.abs(val)))
},
'relu': {
activation: val => val < 0 ? 0 : val,
derivate: val => val < 0 ? 0 : 1
},
'softplus': {
activation: val => Math.log(1 + Math.exp(val)),
derivate: val => 1 / (1 + Math.exp(-val))
},
'bent': {
activation: val => ((Math.sqrt(val * val + 1) - 1) / 2) + val,
derivate: val => (val / (2 * Math.sqrt(val * val + 1))) + 1
},
'sinusoid': {
activation: Math.sin,
derivate: Math.cos
},
'sinc': {
activation: val => val === 0 ? 1 : Math.sin(val) / val,
derivate: val => val === 0 ? 0 : (Math.cos(val) / val) - (Math.sin(val) / (val * val))
},
'gaussian': {
activation: val => Math.exp(-(val * val)),
derivate: val => -2 * val * Math.exp(-(val * val))
},
'parametric-relu': {
activation: (val, param) => val < 0 ? param * val : val,
derivate: (val, param) => val < 0 ? param : 1
},
'exponential-elu': {
activation: expELU,
derivate: (val, param) => val < 0 ? expELU(val, param) + param : 1
},
'soft-exponential': {
activation: softExponential,
derivate: softExponentialPrime
}
};
FeedforwardNeuralNetworksActivationFunctions = ACTIVATION_FUNCTIONS;
}
// feedforward-neural-networks Layer.js
let FeedforwardNeuralNetworksLayer;
{
let Matrix = MLMatrix;
let Utils = FeedforwardNeuralNetworksUtils;
const ACTIVATION_FUNCTIONS = FeedforwardNeuralNetworksActivationFunctions;
class Layer {
/**
* @private
* Create a new layer with the given options
* @param {object} options
* @param {number} [options.inputSize] - Number of conections that enter the neurons.
* @param {number} [options.outputSize] - Number of conections that leave the neurons.
* @param {number} [options.regularization] - Regularization parameter.
* @param {number} [options.epsilon] - Learning rate parameter.
* @param {string} [options.activation] - Activation function parameter from the FeedForwardNeuralNetwork class.
* @param {number} [options.activationParam] - Activation parameter if needed.
*/
constructor(options) {
this.inputSize = options.inputSize;
this.outputSize = options.outputSize;
this.regularization = options.regularization;
this.epsilon = options.epsilon;
this.activation = options.activation;
this.activationParam = options.activationParam;
var selectedFunction = ACTIVATION_FUNCTIONS[options.activation];
var params = selectedFunction.activation.length;
var actFunction = params > 1 ? val => selectedFunction.activation(val, options.activationParam) : selectedFunction.activation;
var derFunction = params > 1 ? val => selectedFunction.derivate(val, options.activationParam) : selectedFunction.derivate;
this.activationFunction = function (i, j) {
this[i][j] = actFunction(this[i][j]);
};
this.derivate = function (i, j) {
this[i][j] = derFunction(this[i][j]);
};
if (options.model) {
// load model
this.W = Matrix.checkMatrix(options.W);
this.b = Matrix.checkMatrix(options.b);
} else {
// default constructor
this.W = Matrix.rand(this.inputSize, this.outputSize);
this.b = Matrix.zeros(1, this.outputSize);
this.W.apply(function (i, j) {
this[i][j] /= Math.sqrt(options.inputSize);
});
}
}
/**
* @private
* propagate the given input through the current layer.
* @param {Matrix} X - input.
* @return {Matrix} output at the current layer.
*/
forward(X) {
var z = X.mmul(this.W).addRowVector(this.b);
z.apply(this.activationFunction);
this.a = z.clone();
return z;
}
/**
* @private
* apply backpropagation algorithm at the current layer
* @param {Matrix} delta - delta values estimated at the following layer.
* @param {Matrix} a - 'a' values from the following layer.
* @return {Matrix} the new delta values for the next layer.
*/
backpropagation(delta, a) {
this.dW = a.transposeView().mmul(delta);
this.db = Utils.sumCol(delta);
var aCopy = a.clone();
return delta.mmul(this.W.transposeView()).mul(aCopy.apply(this.derivate));
}
/**
* @private
* Function that updates the weights at the current layer with the derivatives.
*/
update() {
this.dW.add(this.W.clone().mul(this.regularization));
this.W.add(this.dW.mul(-this.epsilon));
this.b.add(this.db.mul(-this.epsilon));
}
/**
* @private
* Export the current layer to JSON.
* @return {object} model
*/
toJSON() {
return {
model: 'Layer',
inputSize: this.inputSize,
outputSize: this.outputSize,
regularization: this.regularization,
epsilon: this.epsilon,
activation: this.activation,
W: this.W,
b: this.b
};
}
/**
* @private
* Creates a new Layer with the given model.
* @param {object} model
* @return {Layer}
*/
static load(model) {
if (model.model !== 'Layer') {
throw new RangeError('the current model is not a Layer model');
}
return new Layer(model);
}
}
FeedforwardNeuralNetworksLayer = Layer;
}
// feedforward-neural-networks OutputLayer.js
let FeedforwardNeuralNetworksOutputLayer;
{
let Layer = FeedforwardNeuralNetworksLayer;
class OutputLayer extends Layer {
constructor(options) {
super(options);
this.activationFunction = function (i, j) {
this[i][j] = Math.exp(this[i][j]);
};
}
static load(model) {
if (model.model !== 'Layer') {
throw new RangeError('the current model is not a Layer model');
}
return new OutputLayer(model);
}
}
FeedforwardNeuralNetworksOutputLayer = OutputLayer;
}
// feedforward-neural-networks FeedForwardNeuralNetwork.js
let FeedforwardNeuralNetwork;
{
const Matrix = MLMatrix;
const Layer = FeedforwardNeuralNetworksLayer;
const OutputLayer = FeedforwardNeuralNetworksOutputLayer;
const Utils = FeedforwardNeuralNetworksUtils;
const ACTIVATION_FUNCTIONS = FeedforwardNeuralNetworksActivationFunctions;
class FeedForwardNeuralNetworks {
/**
* Create a new Feedforword neural network model.
* @param {object} options
* @param {Array} [options.hiddenLayers=[10]] - Array that contains the sizes of the hidden layers.
* @oaram {number} [options.iterations=50] - Number of iterations at the training step.
* @param {number} [options.learningRate=0.01] - Learning rate of the neural net (also known as epsilon).
* @poram {number} [options.regularization=0.01] - Regularization parameter af the neural net.
* @poram {string} [options.activation='tanh'] - activation function to be used. (options: 'tanh'(default),
* 'identity', 'logistic', 'arctan', 'softsign', 'relu', 'softplus', 'bent', 'sinusoid', 'sinc', 'gaussian').
* (single-parametric options: 'parametric-relu', 'exponential-relu', 'soft-exponential').
* @param {number} [options.activationParam=1] - if the selected activation function needs a parameter.
*/
constructor(options) {
options = options || {};
if (options.model) {
// load network
this.hiddenLayers = options.hiddenLayers;
this.iterations = options.iterations;
this.learningRate = options.learningRate;
this.regularization = options.regularization;
this.dicts = options.dicts;
this.activation = options.activation;
this.activationParam = options.activationParam;
this.model = new Array(options.layers.length);
for (var i = 0; i < this.model.length - 1; ++i) {
this.model[i] = Layer.load(options.layers[i]);
}
this.model[this.model.length - 1] = OutputLayer.load(options.layers[this.model.length - 1]);
} else {
// default constructor
this.hiddenLayers = options.hiddenLayers === undefined ? [10] : options.hiddenLayers;
this.iterations = options.iterations === undefined ? 50 : options.iterations;
this.learningRate = options.learningRate === undefined ? 0.01 : options.learningRate;
//this.momentum = options.momentum === undefined ? 0.1 : options.momentum;
this.regularization = options.regularization === undefined ? 0.01 : options.regularization;
this.activation = options.activation === undefined ? 'tanh' : options.activation;
this.activationParam = options.activationParam === undefined ? 1 : options.activationParam;
if (!(this.activation in Object.keys(ACTIVATION_FUNCTIONS))) {
this.activation = 'tanh';
}
}
}
/**
* @private
* Function that build and initialize the neural net.
* @param {number} inputSize - total of features to fit.
* @param {number} outputSize - total of labels of the prediction set.
*/
buildNetwork(inputSize, outputSize) {
var size = 2 + (this.hiddenLayers.length - 1);
this.model = new Array(size);
// input layer
this.model[0] = new Layer({
inputSize: inputSize,
outputSize: this.hiddenLayers[0],
activation: this.activation,
activationParam: this.activationParam,
regularization: this.regularization,
epsilon: this.learningRate
});
// hidden layers
for (var i = 1; i < this.hiddenLayers.length; ++i) {
this.model[i] = new Layer({
inputSize: this.hiddenLayers[i - 1],
outputSize: this.hiddenLayers[i],
activation: this.activation,
activationParam: this.activationParam,
regularization: this.regularization,
epsilon: this.learningRate
});
}
// output layer
this.model[size - 1] = new OutputLayer({
inputSize: this.hiddenLayers[this.hiddenLayers.length - 1],
outputSize: outputSize,
activation: this.activation,
activationParam: this.activationParam,
regularization: this.regularization,
epsilon: this.learningRate
});
}
/**
* Train the neural net with the given features and labels.
* @param {Matrix|Array} features
* @param {Matrix|Array} labels
*/
train(features, labels) {
features = Matrix.checkMatrix(features);
this.dicts = Utils.dictOutputs(labels);
var inputSize = features.columns;
var outputSize = Object.keys(this.dicts.inputs).length;
this.buildNetwork(inputSize, outputSize);
for (var i = 0; i < this.iterations; ++i) {
var probabilities = this.propagate(features);
this.backpropagation(features, labels, probabilities);
}
}
/**
* @private
* Propagate the input(training set) and retrives the probabilities of each class.
* @param {Matrix} X
* @return {Matrix} probabilities of each class.
*/
propagate(X) {
var input = X;
for (var i = 0; i < this.model.length; ++i) {
//console.log(i);
input = this.model[i].forward(input);
}
// get probabilities
return input.divColumnVector(Utils.sumRow(input));
}
/**
* @private
* Function that applies the backpropagation algorithm on each layer of the network
* in order to fit the features and labels.
* @param {Matrix} features
* @param {Array} labels
* @param {Matrix} probabilities - probabilities of each class of the feature set.
*/
backpropagation(features, labels, probabilities) {
for (var i = 0; i < probabilities.length; ++i) {
probabilities[i][this.dicts.inputs[labels[i]]] -= 1;
}
// remember, the last delta doesn't matter
var delta = probabilities;
for (i = this.model.length - 1; i >= 0; --i) {
var a = i > 0 ? this.model[i - 1].a : features;
delta = this.model[i].backpropagation(delta, a);
}
for (i = 0; i < this.model.length; ++i) {
this.model[i].update();
}
}
/**
* Predict the output given the feature set.
* @param {Array|Matrix} features
* @return {Array}
*/
predict(features) {
features = Matrix.checkMatrix(features);
var outputs = new Array(features.rows);
var probabilities = this.propagate(features);
for (var i = 0; i < features.rows; ++i) {
outputs[i] = this.dicts.outputs[probabilities.maxRowIndex(i)[1]];
}
return outputs;
}
/**
* Export the current model to JSOM.
* @return {object} model
*/
toJSON() {
var model = {
model: 'FNN',
hiddenLayers: this.hiddenLayers,
iterations: this.iterations,
learningRate: this.learningRate,
regularization: this.regularization,
activation: this.activation,
activationParam: this.activationParam,
dicts: this.dicts,
layers: new Array(this.model.length)
};
for (var i = 0; i < this.model.length; ++i) {
model.layers[i] = this.model[i].toJSON();
}
return model;
}
/**
* Load a Feedforward Neural Network with the current model.
* @param {object} model
* @return {FeedForwardNeuralNetworks}
*/
static load(model) {
if (model.model !== 'FNN') {
throw new RangeError('the current model is not a feed forward network');
}
return new FeedForwardNeuralNetworks(model);
}
}
FeedforwardNeuralNetwork = FeedForwardNeuralNetworks;
}