diff options
Diffstat (limited to 'third_party/webkit/PerformanceTests/ARES-6/ml/index.js')
| -rw-r--r-- | third_party/webkit/PerformanceTests/ARES-6/ml/index.js | 6330 |
1 files changed, 6330 insertions, 0 deletions
diff --git a/third_party/webkit/PerformanceTests/ARES-6/ml/index.js b/third_party/webkit/PerformanceTests/ARES-6/ml/index.js new file mode 100644 index 0000000000..03871ba081 --- /dev/null +++ b/third_party/webkit/PerformanceTests/ARES-6/ml/index.js @@ -0,0 +1,6330 @@ +/* + * 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; +} |