зеркало из https://github.com/mozilla/gecko-dev.git
6331 строка
207 KiB
JavaScript
6331 строка
207 KiB
JavaScript
/*
|
|
* 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
|
|
// http://jsperf.com/access-and-write-array-subclass
|
|
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;
|
|
|
|
// https://github.com/lutzroeder/Mapack/blob/master/Source/LuDecomposition.cs
|
|
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;
|
|
|
|
// https://github.com/lutzroeder/Mapack/blob/master/Source/SingularValueDecomposition.cs
|
|
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;
|
|
},
|
|
// https://github.com/accord-net/framework/blob/development/Sources/Accord.Math/Decompositions/SingularValueDecomposition.cs
|
|
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
|
|
* See https://en.wikipedia.org/wiki/Kronecker_product
|
|
* @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;
|
|
|
|
// https://github.com/lutzroeder/Mapack/blob/master/Source/CholeskyDecomposition.cs
|
|
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
|
|
};
|
|
|
|
// https://github.com/lutzroeder/Mapack/blob/master/Source/EigenvalueDecomposition.cs
|
|
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;
|
|
|
|
//https://github.com/lutzroeder/Mapack/blob/master/Source/QrDecomposition.cs
|
|
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;
|
|
}
|