Add ML to ARES6
https://bugs.webkit.org/show_bug.cgi?id=171206
Rubber stamped by Filip Pizlo.
This patch adds a new test to ARES6 called ML. ML is an implementation of
a feedforward neural network: https://github.com/mljs/feedforward-neural-networks.
It makes heavy use of classes, and does non-trivial matrix math using the
ml-matrix library: https://github.com/mljs/matrix
* ARES-6/about.html:
* ARES-6/cli.js:
* ARES-6/glue.js:
* ARES-6/index.html:
* ARES-6/ml: Added.
* ARES-6/ml/benchmark.js: Added.
* ARES-6/ml/index.js: Added.
* ARES-6/ml_benchmark.js: Added.
git-svn-id: http://svn.webkit.org/repository/webkit/trunk@215704 268f45cc-cd09-0410-ab3c-d52691b4dbfc
diff --git a/PerformanceTests/ARES-6/about.html b/PerformanceTests/ARES-6/about.html
index eec06f9..b92d849 100644
--- a/PerformanceTests/ARES-6/about.html
+++ b/PerformanceTests/ARES-6/about.html
@@ -16,7 +16,7 @@
<h2>About the Tests</h2>
- <p>ARES-6 measures the execution time of JavaScript’s newest features, including symbols, for-of, arrow functions, Map/Set/WeakMap, let/const, classes, proxies, string interpolation, destructuring, default arguments, spread, tail calls, and generators. ARES-6 is comprised of three sub-tests: Air, Basic, and Babylon.</p>
+ <p>ARES-6 measures the execution time of JavaScript’s newest features, including symbols, for-of, arrow functions, Map/Set/WeakMap, let/const, classes, proxies, string interpolation, destructuring, default arguments, spread, tail calls, and generators. ARES-6 is comprised of four sub-tests: Air, Basic, Babylon, and ML.</p>
<p>Air is an ES2015 port of the <a href="https://webkit.org/blog/5852/introducing-the-b3-jit-compiler/">WebKit B3 JIT</a>’s <a href="https://trac.webkit.org/changeset/201783"><code>Air::allocateStack</code> phase</a>. This code is a heavy user of Map, Set, classes, spread, and for-of. The benchmark runs <code>allocateStack</code> on hot function bodies from other popular JavaScript benchmarks: <code>executeIteration</code> from <a href="https://developers.google.com/octane/">Octane</a>/Gameboy, <code>gaussianBlur</code> from <a href="http://krakenbenchmark.mozilla.org">Kraken</a>, and <code>scanIdentifier</code> from Octane/Typescript. Air also runs <code>allocateStack</code> on a hot function from Air itself. <a href="https://trac.webkit.org/browser/trunk/PerformanceTests/ARES-6/Air?rev=211697">Browse the source.</a></p>
@@ -29,9 +29,16 @@
does non trivial string processing, and creates non-trivial object graphs.
</p>
- <p>ARES-6 rewards browsers that start up quickly and run smoothly. It's not enough to just measure the total running time of a workload. Browsers may perform differently for the same workload depending on how many times it has run. Garbage collection runs periodically, making some iterations take longer than others. Code that runs repeatedly gets optimized by the browser, so the first iteration of any workload is more expensive than the rest. ARES-6 runs each test for 200 iterations, and reports the execution time of the first iteration, the average of the 4 worst iterations, and the overall geometric mean. Each of these values is given an equal weight when computing the ovarall time. ARES-6 equally rewards fast start-up, low jank, and sophisticated adaptive optimizations for long-running code.</p>
+ <p>
+ <a href="https://github.com/mljs/feedforward-neural-networks">ML</a> is an implementation of a <a href="https://en.wikipedia.org/wiki/Feedforward_neural_network"> feedforward neural network</a>.
+ The benchmark trains several networks using different <a href="https://en.wikipedia.org/wiki/Activation_function">activation functions</a>
+ and several sample data sets. ML makes heavy use of classes. It relies on the <a href="https://github.com/mljs/matrix">ml-matrix</a> library
+ and does non-trivial matrix math.
+ </p>
+
+ <p>ARES-6 rewards browsers that start up quickly and run smoothly. It's not enough to just measure the total running time of a workload. Browsers may perform differently for the same workload depending on how many times it has run. Garbage collection runs periodically, making some iterations take longer than others. Code that runs repeatedly gets optimized by the browser, so the first iteration of any workload is more expensive than the rest. ARES-6 runs Air, Basic, and Babylon for 200 iterations, and ML for 60 iterations, and reports the execution time of the first iteration, the average of the 4 worst iterations, and the overall geometric mean. Each of these values is given an equal weight when computing the ovarall time. ARES-6 equally rewards fast start-up, low jank, and sophisticated adaptive optimizations for long-running code.</p>
- <p>Each ARES-6 sample has 200 iterations of Air, Basic, and Babylon. ARES-6 runs 8 samples, and reports the average with 95% confidence intervals. Each sample runs in a fresh <code>iframe</code> to simulate some of the effects of page navigation.</p>
+ <p>Each ARES-6 sample has 200 iterations of Air, Basic, and Babylon, and 60 iterations of ML. ARES-6 runs 6 samples, and reports the average with 95% confidence intervals. Each sample runs in a fresh <code>iframe</code> to simulate some of the effects of page navigation.</p>
<p><a href="index.html" class="button return">Return to Testing</a></p>
diff --git a/PerformanceTests/ARES-6/cli.js b/PerformanceTests/ARES-6/cli.js
index 25b6c67..ffd137d 100644
--- a/PerformanceTests/ARES-6/cli.js
+++ b/PerformanceTests/ARES-6/cli.js
@@ -25,7 +25,7 @@
const isInBrowser = false;
-function makeBenchmarkRunner(sources, benchmarkConstructor) {
+function makeBenchmarkRunner(sources, benchmarkConstructor, numIterations = 200) {
let source = "'use strict';"
for (let file of sources) {
source += readFile(file);
@@ -33,7 +33,7 @@
source += `
this.results = [];
var benchmark = new ${benchmarkConstructor}();
- var numIterations = 200;
+ var numIterations = ${numIterations};
for (var i = 0; i < numIterations; ++i) {
var before = currentTime();
benchmark.runIteration();
@@ -54,6 +54,7 @@
load("air_benchmark.js");
load("basic_benchmark.js");
load("babylon_benchmark.js");
+load("ml_benchmark.js");
load("glue.js");
-driver.start(8);
+driver.start(6);
diff --git a/PerformanceTests/ARES-6/glue.js b/PerformanceTests/ARES-6/glue.js
index 75b9f59..a5388ff 100644
--- a/PerformanceTests/ARES-6/glue.js
+++ b/PerformanceTests/ARES-6/glue.js
@@ -28,7 +28,7 @@
isInBrowser ? document.getElementById("status") : null,
isInBrowser ? document.getElementById("trigger") : null,
function() {
- driver.start(8)
+ driver.start(6)
},
isInBrowser ? document.getElementById("magic") : null,
isInBrowser ? document.getElementById("Geomean") : null,
@@ -41,4 +41,5 @@
driver.addBenchmark(AirBenchmarkRunner);
driver.addBenchmark(BasicBenchmarkRunner);
driver.addBenchmark(BabylonBenchmarkRunner);
+driver.addBenchmark(MLBenchmarkRunner);
driver.readyTrigger();
diff --git a/PerformanceTests/ARES-6/index.html b/PerformanceTests/ARES-6/index.html
index 230657b..e067b25 100644
--- a/PerformanceTests/ARES-6/index.html
+++ b/PerformanceTests/ARES-6/index.html
@@ -119,6 +119,33 @@
</span>
</div>
</div>
+
+ <div class="ML test">
+ <div id="MLMessage" class="indicator">.</div>
+ <h2>ML</h2>
+
+ <div class="score">
+ <label>First Iteration</label>
+ <span id="MLFirstIteration">
+ <span class="value">∅</span><span class="units">ms</span>
+ </span>
+ </div>
+
+ <div class="score">
+ <label>Worst 4 Iterations</label>
+ <span id="MLAverageWorstCase">
+ <span class="value">∅</span><span class="units">ms</span>
+ </span>
+ </div>
+
+ <div class="score">
+ <label>Average</label>
+ <span id="MLSteadyState">
+ <span class="value">∅</span><span class="units">ms</span>
+ </span>
+ </div>
+ </div>
+
</div>
</main>
@@ -126,6 +153,7 @@
<script src="air_benchmark.js"></script>
<script src="basic_benchmark.js"></script>
<script src="babylon_benchmark.js"></script>
+ <script src="ml_benchmark.js"></script>
<script src="glue.js"></script>
</body>
diff --git a/PerformanceTests/ARES-6/ml/benchmark.js b/PerformanceTests/ARES-6/ml/benchmark.js
new file mode 100644
index 0000000..bd6ccd3
--- /dev/null
+++ b/PerformanceTests/ARES-6/ml/benchmark.js
@@ -0,0 +1,171 @@
+/*
+ * Copyright (C) 2017 Apple Inc. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY APPLE INC. ``AS IS'' AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL APPLE INC. OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
+ * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+"use strict";
+
+let currentTime;
+if (this.performance && performance.now)
+ currentTime = function() { return performance.now() };
+else if (this.preciseTime)
+ currentTime = function() { return preciseTime() * 1000; };
+else
+ currentTime = function() { return +new Date(); };
+
+class MLBenchmark {
+ constructor() { }
+
+ runIteration()
+ {
+ let Matrix = MLMatrix;
+ let ACTIVATION_FUNCTIONS = FeedforwardNeuralNetworksActivationFunctions;
+
+ function run() {
+
+ let it = (name, f) => {
+ f();
+ };
+
+ function assert(b) {
+ if (!b)
+ throw new Error("Bad");
+ }
+
+ var functions = Object.keys(ACTIVATION_FUNCTIONS);
+
+ it('Training the neural network with XOR operator', function () {
+ var trainingSet = new Matrix([[0, 0], [0, 1], [1, 0], [1, 1]]);
+ var predictions = [false, true, true, false];
+
+ for (var i = 0; i < functions.length; ++i) {
+ var options = {
+ hiddenLayers: [4],
+ iterations: 40,
+ learningRate: 0.3,
+ activation: functions[i]
+ };
+ var xorNN = new FeedforwardNeuralNetwork(options);
+
+ xorNN.train(trainingSet, predictions);
+ var results = xorNN.predict(trainingSet);
+ }
+ });
+
+ it('Training the neural network with AND operator', function () {
+ var trainingSet = [[0, 0], [0, 1], [1, 0], [1, 1]];
+ var predictions = [[1, 0], [1, 0], [1, 0], [0, 1]];
+
+ for (var i = 0; i < functions.length; ++i) {
+ var options = {
+ hiddenLayers: [3],
+ iterations: 75,
+ learningRate: 0.3,
+ activation: functions[i]
+ };
+ var andNN = new FeedforwardNeuralNetwork(options);
+ andNN.train(trainingSet, predictions);
+
+ var results = andNN.predict(trainingSet);
+ }
+ });
+
+ it('Export and import', function () {
+ var trainingSet = [[0, 0], [0, 1], [1, 0], [1, 1]];
+ var predictions = [0, 1, 1, 1];
+
+ for (var i = 0; i < functions.length; ++i) {
+ var options = {
+ hiddenLayers: [4],
+ iterations: 40,
+ learningRate: 0.3,
+ activation: functions[i]
+ };
+ var orNN = new FeedforwardNeuralNetwork(options);
+ orNN.train(trainingSet, predictions);
+
+ var model = JSON.parse(JSON.stringify(orNN));
+ var networkNN = FeedforwardNeuralNetwork.load(model);
+
+ var results = networkNN.predict(trainingSet);
+ }
+ });
+
+ it('Multiclass clasification', function () {
+ var trainingSet = [[0, 0], [0, 1], [1, 0], [1, 1]];
+ var predictions = [2, 0, 1, 0];
+
+ for (var i = 0; i < functions.length; ++i) {
+ var options = {
+ hiddenLayers: [4],
+ iterations: 40,
+ learningRate: 0.5,
+ activation: functions[i]
+ };
+ var nn = new FeedforwardNeuralNetwork(options);
+ nn.train(trainingSet, predictions);
+
+ var result = nn.predict(trainingSet);
+ }
+ });
+
+ it('Big case', function () {
+ var trainingSet = [[1, 1], [1, 2], [2, 1], [2, 2], [3, 1], [1, 3], [1, 4], [4, 1],
+ [6, 1], [6, 2], [6, 3], [6, 4], [6, 5], [5, 5], [4, 5], [3, 5]];
+ var predictions = [[1, 0], [1, 0], [1, 0], [1, 0], [1, 0], [1, 0], [1, 0], [1, 0],
+ [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1]];
+ for (var i = 0; i < functions.length; ++i) {
+ var options = {
+ hiddenLayers: [20],
+ iterations: 60,
+ learningRate: 0.01,
+ activation: functions[i]
+ };
+ var nn = new FeedforwardNeuralNetwork(options);
+ nn.train(trainingSet, predictions);
+
+ var result = nn.predict([[5, 4]]);
+
+ assert(result[0][0] < result[0][1]);
+ }
+ });
+ }
+
+ run();
+ }
+}
+
+function runBenchmark()
+{
+ const numIterations = 60;
+
+ let before = currentTime();
+
+ let benchmark = new Benchmark();
+
+ for (let iteration = 0; iteration < numIterations; ++iteration)
+ benchmark.runIteration();
+
+ let after = currentTime();
+ return after - before;
+}
diff --git a/PerformanceTests/ARES-6/ml/index.js b/PerformanceTests/ARES-6/ml/index.js
new file mode 100644
index 0000000..03871ba
--- /dev/null
+++ b/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;
+}
diff --git a/PerformanceTests/ARES-6/ml_benchmark.js b/PerformanceTests/ARES-6/ml_benchmark.js
new file mode 100644
index 0000000..fa4ee17
--- /dev/null
+++ b/PerformanceTests/ARES-6/ml_benchmark.js
@@ -0,0 +1,68 @@
+/*
+ * Copyright (C) 2017 Apple Inc. All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY APPLE INC. ``AS IS'' AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL APPLE INC. OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
+ * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+"use strict";
+
+const MLBenchmarkCode = String.raw`
+<script src="ml/index.js"></script>
+<script src="ml/benchmark.js"></script>
+<script>
+var results = [];
+var benchmark = new MLBenchmark();
+var numIterations = 60;
+for (var i = 0; i < numIterations; ++i) {
+ var before = currentTime();
+ benchmark.runIteration();
+ var after = currentTime();
+ results.push(after - before);
+}
+reportResult(results);
+</script>`;
+
+
+let runMLBenchmark = null;
+if (!isInBrowser) {
+ let sources = [
+ "ml/index.js"
+ , "ml/benchmark.js"
+ ];
+
+ runMLBenchmark = makeBenchmarkRunner(sources, "MLBenchmark", 60);
+}
+
+const MLBenchmarkRunner = {
+ name: "ML",
+ code: MLBenchmarkCode,
+ run: runMLBenchmark,
+ cells: {}
+};
+
+if (isInBrowser) {
+ MLBenchmarkRunner.cells = {
+ firstIteration: document.getElementById("MLFirstIteration"),
+ averageWorstCase: document.getElementById("MLAverageWorstCase"),
+ steadyState: document.getElementById("MLSteadyState"),
+ message: document.getElementById("MLMessage")
+ };
+}
diff --git a/PerformanceTests/ChangeLog b/PerformanceTests/ChangeLog
index ec46bf3..e33386c 100644
--- a/PerformanceTests/ChangeLog
+++ b/PerformanceTests/ChangeLog
@@ -1,3 +1,24 @@
+2017-04-24 Saam Barati <sbarati@apple.com>
+
+ Add ML to ARES6
+ https://bugs.webkit.org/show_bug.cgi?id=171206
+
+ Rubber stamped by Filip Pizlo.
+
+ This patch adds a new test to ARES6 called ML. ML is an implementation of
+ a feedforward neural network: https://github.com/mljs/feedforward-neural-networks.
+ It makes heavy use of classes, and does non-trivial matrix math using the
+ ml-matrix library: https://github.com/mljs/matrix
+
+ * ARES-6/about.html:
+ * ARES-6/cli.js:
+ * ARES-6/glue.js:
+ * ARES-6/index.html:
+ * ARES-6/ml: Added.
+ * ARES-6/ml/benchmark.js: Added.
+ * ARES-6/ml/index.js: Added.
+ * ARES-6/ml_benchmark.js: Added.
+
2017-04-21 Zalan Bujtas <zalan@apple.com>
Simple line layout: Add performance test to measure mid-word line breaking with long text.
