diff --git a/LICENSE b/LICENSE index 471d2b0..5f4a9a2 100644 --- a/LICENSE +++ b/LICENSE @@ -186,7 +186,7 @@ same "printed page" as the copyright notice for easier identification within third-party archives. - Copyright 2019 Google Inc. + Copyright 2019 Google LLC. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. diff --git a/jest.config.js b/jest.config.js index 6984c29..5f7c264 100644 --- a/jest.config.js +++ b/jest.config.js @@ -1,4 +1,4 @@ -/* Copyright 2019 Google Inc. All Rights Reserved. +/* Copyright 2019 Google LLC. All Rights Reserved. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -17,4 +17,3 @@ module.exports = { preset: 'ts-jest', testEnvironment: 'node', }; - diff --git a/lib/umap-js.d.ts b/lib/umap-js.d.ts deleted file mode 100644 index 3d52026..0000000 --- a/lib/umap-js.d.ts +++ /dev/null @@ -1,25 +0,0 @@ -// Generated by dts-bundle v0.7.3 - -declare module 'umap' { - export type DistanceFn = (x: Vector, y: Vector) => number; - export type EpochCallback = (epoch: number) => boolean | void; - export type Vector = number[]; - export type Vectors = Vector[]; - export interface UMAPParameters { - nComponents?: number; - nEpochs?: number; - nNeighbors?: number; - random?: () => number; - } - export class UMAP { - constructor(params?: UMAPParameters); - fit(X: Vectors): number[][]; - fitAsync(X: Vectors, callback?: (epochNumber: number) => void | boolean): Promise; - initializeFit(X: Vectors, knnIndices?: number[][], knnDistances?: number[][]): number; - step(): number; - getEmbedding(): number[][]; - } - export function euclidean(x: Vector, y: Vector): number; - export function cosine(x: Vector, y: Vector): number; -} - diff --git a/lib/umap-js.js b/lib/umap-js.js index b78b3b6..6430a44 100644 --- a/lib/umap-js.js +++ b/lib/umap-js.js @@ -81,7 +81,7 @@ /******/ /******/ /******/ // Load entry module and return exports -/******/ return __webpack_require__(__webpack_require__.s = 1); +/******/ return __webpack_require__(__webpack_require__.s = 2); /******/ }) /************************************************************************/ /******/ ([ @@ -90,6 +90,32 @@ "use strict"; + +const toString = Object.prototype.toString; + +function isAnyArray(object) { + return toString.call(object).endsWith('Array]'); +} + +module.exports = isAnyArray; + + +/***/ }), +/* 1 */ +/***/ (function(module, exports, __webpack_require__) { + +"use strict"; + +var __values = (this && this.__values) || function (o) { + var m = typeof Symbol === "function" && o[Symbol.iterator], i = 0; + if (m) return m.call(o); + return { + next: function () { + if (o && i >= o.length) o = void 0; + return { value: o && o[i++], done: !o }; + } + }; +}; Object.defineProperty(exports, "__esModule", { value: true }); function generateGaussian(mean, std, rng) { if (rng === void 0) { rng = Math.random; } @@ -123,10 +149,20 @@ function tauRand(random) { } exports.tauRand = tauRand; function norm(vec) { + var e_1, _a; var result = 0; - for (var _i = 0, vec_1 = vec; _i < vec_1.length; _i++) { - var item = vec_1[_i]; - result += Math.pow(item, 2); + try { + for (var vec_1 = __values(vec), vec_1_1 = vec_1.next(); !vec_1_1.done; vec_1_1 = vec_1.next()) { + var item = vec_1_1.value; + result += Math.pow(item, 2); + } + } + catch (e_1_1) { e_1 = { error: e_1_1 }; } + finally { + try { + if (vec_1_1 && !vec_1_1.done && (_a = vec_1.return)) _a.call(vec_1); + } + finally { if (e_1) throw e_1.error; } } return Math.sqrt(result); } @@ -155,6 +191,12 @@ function ones(n) { return filled(n, 1); } exports.ones = ones; +function linear(a, b, len) { + return empty(len).map(function (_, i) { + return a + i * ((b - a) / (len - 1)); + }); +} +exports.linear = linear; function sum(input) { return input.reduce(function (sum, val) { return sum + val; }); } @@ -184,18 +226,18 @@ exports.max2d = max2d; /***/ }), -/* 1 */ +/* 2 */ /***/ (function(module, exports, __webpack_require__) { "use strict"; Object.defineProperty(exports, "__esModule", { value: true }); -var umap_1 = __webpack_require__(2); +var umap_1 = __webpack_require__(3); window.UMAP = umap_1.UMAP; /***/ }), -/* 2 */ +/* 3 */ /***/ (function(module, exports, __webpack_require__) { "use strict"; @@ -235,28 +277,56 @@ var __generator = (this && this.__generator) || function (thisArg, body) { if (op[0] & 5) throw op[1]; return { value: op[0] ? op[1] : void 0, done: true }; } }; +var __read = (this && this.__read) || function (o, n) { + var m = typeof Symbol === "function" && o[Symbol.iterator]; + if (!m) return o; + var i = m.call(o), r, ar = [], e; + try { + while ((n === void 0 || n-- > 0) && !(r = i.next()).done) ar.push(r.value); + } + catch (error) { e = { error: error }; } + finally { + try { + if (r && !r.done && (m = i["return"])) m.call(i); + } + finally { if (e) throw e.error; } + } + return ar; +}; +var __spread = (this && this.__spread) || function () { + for (var ar = [], i = 0; i < arguments.length; i++) ar = ar.concat(__read(arguments[i])); + return ar; +}; Object.defineProperty(exports, "__esModule", { value: true }); -var matrix = __webpack_require__(3); -var nnDescent = __webpack_require__(4); -var tree = __webpack_require__(6); -var utils = __webpack_require__(0); +var matrix = __webpack_require__(4); +var nnDescent = __webpack_require__(5); +var tree = __webpack_require__(7); +var utils = __webpack_require__(1); +var LM = __webpack_require__(8); var SMOOTH_K_TOLERANCE = 1e-5; var MIN_K_DIST_SCALE = 1e-3; var UMAP = (function () { function UMAP(params) { if (params === void 0) { params = {}; } - this.nNeighbors = 15; + this.minDist = 0.1; this.nComponents = 2; this.nEpochs = 0; + this.nNeighbors = 15; this.random = Math.random; + this.spread = 1.0; + this.targetMetric = "categorical"; + this.targetWeight = 0.5; + this.targetNNeighbors = this.nNeighbors; this.distanceFn = euclidean; this.isInitialized = false; this.embedding = []; this.optimizationState = new OptimizationState(); + this.minDist = params.minDist || this.minDist; this.nComponents = params.nComponents || this.nComponents; this.nEpochs = params.nEpochs || this.nEpochs; this.nNeighbors = params.nNeighbors || this.nNeighbors; this.random = params.random || this.random; + this.spread = params.spread || this.spread; } UMAP.prototype.fit = function (X) { this.initializeFit(X); @@ -278,21 +348,29 @@ var UMAP = (function () { }); }); }; - UMAP.prototype.initializeFit = function (X, knnIndices, knnDistances) { - if (this.data === X && this.isInitialized) { + UMAP.prototype.setSupervisedProjection = function (Y, params) { + if (params === void 0) { params = {}; } + this.Y = Y; + this.targetMetric = params.targetMetric || this.targetMetric; + this.targetWeight = params.targetWeight || this.targetWeight; + this.targetNNeighbors = params.targetNNeighbors || this.targetNNeighbors; + }; + UMAP.prototype.setPrecomputedKNN = function (knnIndices, knnDistances) { + this.knnIndices = knnIndices; + this.knnDistances = knnDistances; + }; + UMAP.prototype.initializeFit = function (X) { + if (this.X === X && this.isInitialized) { return this.getNEpochs(); } - this.data = X; - if (knnIndices && knnDistances) { - this.knnIndices = knnIndices; - this.knnDistances = knnDistances; - } - else { + this.X = X; + if (!this.knnIndices && !this.knnDistances) { var knnResults = this.nearestNeighbors(X); this.knnIndices = knnResults.knnIndices; this.knnDistances = knnResults.knnDistances; } - this.graph = this.fuzzySimplicialSet(X); + this.graph = this.fuzzySimplicialSet(X, this.nNeighbors); + this.processGraphForSupervisedProjection(); var _a = this.initializeSimplicialSetEmbedding(), head = _a.head, tail = _a.tail, epochsPerSample = _a.epochsPerSample; this.optimizationState.head = head; this.optimizationState.tail = tail; @@ -300,6 +378,19 @@ var UMAP = (function () { this.isInitialized = true; return this.getNEpochs(); }; + UMAP.prototype.processGraphForSupervisedProjection = function () { + var _a = this, Y = _a.Y, X = _a.X; + if (Y) { + if (Y.length !== X.length) { + throw new Error('Length of X and y must be equal'); + } + if (this.targetMetric === "categorical") { + var lt = this.targetWeight < 1.0; + var farDist = lt ? 2.5 * (1.0 / (1.0 - this.targetWeight)) : 1.0e12; + this.graph = this.categoricalSimplicialSetIntersection(this.graph, Y, farDist); + } + } + }; UMAP.prototype.step = function () { var _a = this.optimizationState, currentEpoch = _a.currentEpoch, isInitialized = _a.isInitialized; if (!isInitialized) { @@ -327,22 +418,28 @@ var UMAP = (function () { var _b = metricNNDescent(X, leafArray, nNeighbors, nIters), indices = _b.indices, weights = _b.weights; return { knnIndices: indices, knnDistances: weights }; }; - UMAP.prototype.fuzzySimplicialSet = function (X, localConnectivity, setOpMixRatio) { + UMAP.prototype.fuzzySimplicialSet = function (X, nNeighbors, localConnectivity, setOpMixRatio) { if (localConnectivity === void 0) { localConnectivity = 1.0; } if (setOpMixRatio === void 0) { setOpMixRatio = 1.0; } - var _a = this, nNeighbors = _a.nNeighbors, _b = _a.knnIndices, knnIndices = _b === void 0 ? [] : _b, _c = _a.knnDistances, knnDistances = _c === void 0 ? [] : _c; + var _a = this, _b = _a.knnIndices, knnIndices = _b === void 0 ? [] : _b, _c = _a.knnDistances, knnDistances = _c === void 0 ? [] : _c; var _d = this.smoothKNNDistance(knnDistances, nNeighbors, localConnectivity), sigmas = _d.sigmas, rhos = _d.rhos; var _e = this.computeMembershipStrengths(knnIndices, knnDistances, sigmas, rhos), rows = _e.rows, cols = _e.cols, vals = _e.vals; var size = [X.length, X.length]; var sparseMatrix = new matrix.SparseMatrix(rows, cols, vals, size); var transpose = matrix.transpose(sparseMatrix); - var prodMatrix = matrix.dotMultiply(sparseMatrix, transpose); + var prodMatrix = matrix.pairwiseMultiply(sparseMatrix, transpose); var a = matrix.subtract(matrix.add(sparseMatrix, transpose), prodMatrix); var b = matrix.multiplyScalar(a, setOpMixRatio); var c = matrix.multiplyScalar(prodMatrix, 1.0 - setOpMixRatio); var result = matrix.add(b, c); return result; }; + UMAP.prototype.categoricalSimplicialSetIntersection = function (simplicialSet, target, farDist, unknownDist) { + if (unknownDist === void 0) { unknownDist = 1.0; } + var intersection = fastIntersection(simplicialSet, target, unknownDist, farDist); + intersection = matrix.eliminateZeros(intersection); + return resetLocalConnectivity(intersection); + }; UMAP.prototype.smoothKNNDistance = function (distances, k, localConnectivity, nIter, bandwidth) { if (localConnectivity === void 0) { localConnectivity = 1.0; } if (nIter === void 0) { nIter = 64; } @@ -505,14 +602,13 @@ var UMAP = (function () { var negativeSampleRate = 5; var nEpochs = this.getNEpochs(); var nVertices = this.graph.nCols; - var a = 1.5769434603113077; - var b = 0.8950608779109733; + var _b = findABParams(this.spread, this.minDist), a = _b.a, b = _b.b; var dim = headEmbedding[0].length; var moveOther = headEmbedding.length === tailEmbedding.length; var alpha = initialAlpha; var epochsPerNegativeSample = epochsPerSample.map(function (e) { return e / negativeSampleRate; }); - var epochOfNextNegativeSample = epochsPerNegativeSample.slice(); - var epochOfNextSample = epochsPerSample.slice(); + var epochOfNextNegativeSample = __spread(epochsPerNegativeSample); + var epochOfNextSample = __spread(epochsPerSample); Object.assign(this.optimizationState, { isInitialized: true, headEmbedding: headEmbedding, @@ -709,23 +805,107 @@ function rDist(x, y) { } return result; } +function findABParams(spread, minDist) { + var curve = function (_a) { + var _b = __read(_a, 2), a = _b[0], b = _b[1]; + return function (x) { + return 1.0 / (1.0 + a * Math.pow(x, (2 * b))); + }; + }; + var xv = utils + .linear(0, spread * 3, 300) + .map(function (val) { return (val < minDist ? 1.0 : val); }); + var yv = utils.zeros(xv.length).map(function (val, index) { + var gte = xv[index] >= minDist; + return gte ? Math.exp(-(xv[index] - minDist) / spread) : val; + }); + var initialValues = [0.5, 0.5]; + var data = { x: xv, y: yv }; + var options = { + damping: 1.5, + initialValues: initialValues, + gradientDifference: 10e-2, + maxIterations: 100, + errorTolerance: 10e-3, + }; + var parameterValues = LM(data, curve, options).parameterValues; + var _a = __read(parameterValues, 2), a = _a[0], b = _a[1]; + return { a: a, b: b }; +} +exports.findABParams = findABParams; +function fastIntersection(graph, target, unknownDist, farDist) { + if (unknownDist === void 0) { unknownDist = 1.0; } + if (farDist === void 0) { farDist = 5.0; } + return graph.map(function (value, row, col) { + if (target[row] === -1 || target[col] === -1) { + return value * Math.exp(-unknownDist); + } + else if (target[row] !== target[col]) { + return value * Math.exp(-farDist); + } + else { + return value; + } + }); +} +exports.fastIntersection = fastIntersection; +function resetLocalConnectivity(simplicialSet) { + simplicialSet = matrix.normalize(simplicialSet, "max"); + var transpose = matrix.transpose(simplicialSet); + var prodMatrix = matrix.pairwiseMultiply(transpose, simplicialSet); + simplicialSet = matrix.add(simplicialSet, matrix.subtract(transpose, prodMatrix)); + return matrix.eliminateZeros(simplicialSet); +} +exports.resetLocalConnectivity = resetLocalConnectivity; /***/ }), -/* 3 */ +/* 4 */ /***/ (function(module, exports, __webpack_require__) { "use strict"; +var __read = (this && this.__read) || function (o, n) { + var m = typeof Symbol === "function" && o[Symbol.iterator]; + if (!m) return o; + var i = m.call(o), r, ar = [], e; + try { + while ((n === void 0 || n-- > 0) && !(r = i.next()).done) ar.push(r.value); + } + catch (error) { e = { error: error }; } + finally { + try { + if (r && !r.done && (m = i["return"])) m.call(i); + } + finally { if (e) throw e.error; } + } + return ar; +}; +var __spread = (this && this.__spread) || function () { + for (var ar = [], i = 0; i < arguments.length; i++) ar = ar.concat(__read(arguments[i])); + return ar; +}; +var __values = (this && this.__values) || function (o) { + var m = typeof Symbol === "function" && o[Symbol.iterator], i = 0; + if (m) return m.call(o); + return { + next: function () { + if (o && i >= o.length) o = void 0; + return { value: o && o[i++], done: !o }; + } + }; +}; Object.defineProperty(exports, "__esModule", { value: true }); +var _a; +var utils = __webpack_require__(1); var SparseMatrix = (function () { function SparseMatrix(rows, cols, values, dims) { this.entries = new Map(); this.nRows = 0; this.nCols = 0; - this.rows = rows.slice(); - this.cols = cols.slice(); - this.values = values.slice(); + this.rows = __spread(rows); + this.cols = __spread(cols); + this.values = __spread(values); for (var i = 0; i < values.length; i++) { var key = this.makeKey(this.rows[i], this.cols[i]); this.entries.set(key, i); @@ -768,14 +948,17 @@ var SparseMatrix = (function () { return defaultValue; } }; + SparseMatrix.prototype.getDims = function () { + return [this.nRows, this.nCols]; + }; SparseMatrix.prototype.getRows = function () { - return this.rows.slice(); + return __spread(this.rows); }; SparseMatrix.prototype.getCols = function () { - return this.cols.slice(); + return __spread(this.cols); }; SparseMatrix.prototype.getValues = function () { - return this.values.slice(); + return __spread(this.values); }; SparseMatrix.prototype.forEach = function (fn) { for (var i = 0; i < this.values.length; i++) { @@ -792,10 +975,9 @@ var SparseMatrix = (function () { }; SparseMatrix.prototype.toArray = function () { var _this = this; - var rows = new Array(this.nRows).slice(); + var rows = utils.empty(this.nRows); var output = rows.map(function () { - var cols = new Array(_this.nCols).slice(); - return cols.map(function () { return 0; }); + return utils.zeros(_this.nCols); }); for (var i = 0; i < this.values.length; i++) { output[this.rows[i]][this.cols[i]] = this.values[i]; @@ -819,7 +1001,7 @@ function transpose(matrix) { } exports.transpose = transpose; function identity(size) { - var rows = size[0]; + var _a = __read(size, 1), rows = _a[0]; var matrix = new SparseMatrix([], [], [], size); for (var i = 0; i < rows; i++) { matrix.set(i, i, 1); @@ -827,10 +1009,10 @@ function identity(size) { return matrix; } exports.identity = identity; -function dotMultiply(a, b) { +function pairwiseMultiply(a, b) { return elementWise(a, b, function (x, y) { return x * y; }); } -exports.dotMultiply = dotMultiply; +exports.pairwiseMultiply = pairwiseMultiply; function add(a, b) { return elementWise(a, b, function (x, y) { return x + y; }); } @@ -845,6 +1027,81 @@ function multiplyScalar(a, scalar) { }); } exports.multiplyScalar = multiplyScalar; +function eliminateZeros(m) { + var zeroIndices = new Set(); + var values = m.getValues(); + var rows = m.getRows(); + var cols = m.getCols(); + for (var i = 0; i < values.length; i++) { + if (values[i] === 0) { + zeroIndices.add(i); + } + } + var removeByZeroIndex = function (_, index) { return !zeroIndices.has(index); }; + var nextValues = values.filter(removeByZeroIndex); + var nextRows = rows.filter(removeByZeroIndex); + var nextCols = cols.filter(removeByZeroIndex); + return new SparseMatrix(nextRows, nextCols, nextValues, m.getDims()); +} +exports.eliminateZeros = eliminateZeros; +function normalize(m, normType) { + if (normType === void 0) { normType = "l2"; } + var e_1, _a; + var normFn = normFns[normType]; + var colsByRow = new Map(); + m.forEach(function (_, row, col) { + var cols = colsByRow.get(row) || []; + cols.push(col); + colsByRow.set(row, cols); + }); + var nextMatrix = new SparseMatrix([], [], [], m.getDims()); + var _loop_1 = function (row) { + var cols = colsByRow.get(row).sort(); + var vals = cols.map(function (col) { return m.get(row, col); }); + var norm = normFn(vals); + for (var i = 0; i < norm.length; i++) { + nextMatrix.set(row, cols[i], norm[i]); + } + }; + try { + for (var _b = __values(colsByRow.keys()), _c = _b.next(); !_c.done; _c = _b.next()) { + var row = _c.value; + _loop_1(row); + } + } + catch (e_1_1) { e_1 = { error: e_1_1 }; } + finally { + try { + if (_c && !_c.done && (_a = _b.return)) _a.call(_b); + } + finally { if (e_1) throw e_1.error; } + } + return nextMatrix; +} +exports.normalize = normalize; +var normFns = (_a = {}, + _a["max"] = function (xs) { + var max = -Infinity; + for (var i = 0; i < xs.length; i++) { + max = xs[i] > max ? xs[i] : max; + } + return xs.map(function (x) { return x / max; }); + }, + _a["l1"] = function (xs) { + var sum = 0; + for (var i = 0; i < xs.length; i++) { + sum += xs[i]; + } + return xs.map(function (x) { return x / sum; }); + }, + _a["l2"] = function (xs) { + var sum = 0; + for (var i = 0; i < xs.length; i++) { + sum += Math.pow(xs[i], 2); + } + return xs.map(function (x) { return Math.sqrt(Math.pow(x, 2) / sum); }); + }, + _a); function elementWise(a, b, op) { var visited = new Set(); var rows = []; @@ -883,14 +1140,14 @@ function elementWise(a, b, op) { /***/ }), -/* 4 */ +/* 5 */ /***/ (function(module, exports, __webpack_require__) { "use strict"; Object.defineProperty(exports, "__esModule", { value: true }); -var heap = __webpack_require__(5); -var utils = __webpack_require__(0); +var heap = __webpack_require__(6); +var utils = __webpack_require__(1); function makeNNDescent(distanceFn, random) { return function nNDescent(data, leafArray, nNeighbors, nIters, maxCandidates, delta, rho, rpTreeInit) { if (nIters === void 0) { nIters = 10; } @@ -959,13 +1216,13 @@ exports.makeNNDescent = makeNNDescent; /***/ }), -/* 5 */ +/* 6 */ /***/ (function(module, exports, __webpack_require__) { "use strict"; Object.defineProperty(exports, "__esModule", { value: true }); -var utils = __webpack_require__(0); +var utils = __webpack_require__(1); function makeHeap(nPoints, size) { var makeArrays = function (fillValue) { return utils.empty(nPoints).map(function () { @@ -1128,13 +1385,43 @@ function siftDown(heap1, heap2, ceiling, elt) { /***/ }), -/* 6 */ +/* 7 */ /***/ (function(module, exports, __webpack_require__) { "use strict"; +var __read = (this && this.__read) || function (o, n) { + var m = typeof Symbol === "function" && o[Symbol.iterator]; + if (!m) return o; + var i = m.call(o), r, ar = [], e; + try { + while ((n === void 0 || n-- > 0) && !(r = i.next()).done) ar.push(r.value); + } + catch (error) { e = { error: error }; } + finally { + try { + if (r && !r.done && (m = i["return"])) m.call(i); + } + finally { if (e) throw e.error; } + } + return ar; +}; +var __spread = (this && this.__spread) || function () { + for (var ar = [], i = 0; i < arguments.length; i++) ar = ar.concat(__read(arguments[i])); + return ar; +}; +var __values = (this && this.__values) || function (o) { + var m = typeof Symbol === "function" && o[Symbol.iterator], i = 0; + if (m) return m.call(o); + return { + next: function () { + if (o && i >= o.length) o = void 0; + return { value: o && o[i++], done: !o }; + } + }; +}; Object.defineProperty(exports, "__esModule", { value: true }); -var utils = __webpack_require__(0); +var utils = __webpack_require__(1); var FlatTree = (function () { function FlatTree(hyperplanes, offsets, children, indices) { this.hyperplanes = hyperplanes; @@ -1255,7 +1542,7 @@ function recursiveFlatten(tree, hyperplanes, offsets, children, indices, nodeNum var _a; if (tree.isLeaf) { children[nodeNum][0] = -leafNum; - (_a = indices[leafNum]).splice.apply(_a, [0, tree.indices.length].concat(tree.indices)); + (_a = indices[leafNum]).splice.apply(_a, __spread([0, tree.indices.length], tree.indices)); leafNum += 1; return { nodeNum: nodeNum, leafNum: leafNum }; } @@ -1289,11 +1576,21 @@ function numLeaves(tree) { } } function makeLeafArray(rpForest) { + var e_1, _a; if (rpForest.length > 0) { var output = []; - for (var _i = 0, rpForest_1 = rpForest; _i < rpForest_1.length; _i++) { - var tree = rpForest_1[_i]; - output.push.apply(output, tree.indices); + try { + for (var rpForest_1 = __values(rpForest), rpForest_1_1 = rpForest_1.next(); !rpForest_1_1.done; rpForest_1_1 = rpForest_1.next()) { + var tree = rpForest_1_1.value; + output.push.apply(output, __spread(tree.indices)); + } + } + catch (e_1_1) { e_1 = { error: e_1_1 }; } + finally { + try { + if (rpForest_1_1 && !rpForest_1_1.done && (_a = rpForest_1.return)) _a.call(rpForest_1); + } + finally { if (e_1) throw e_1.error; } } return output; } @@ -1304,5 +1601,4878 @@ function makeLeafArray(rpForest) { exports.makeLeafArray = makeLeafArray; +/***/ }), +/* 8 */ +/***/ (function(module, exports, __webpack_require__) { + +"use strict"; + + +var mlMatrix = __webpack_require__(9); + +/** + * Calculate current error + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} parameters - Array of current parameter values + * @param {function} parameterizedFunction - The parameters and returns a function with the independent variable as a parameter + * @return {number} + */ +function errorCalculation( + data, + parameters, + parameterizedFunction +) { + var error = 0; + const func = parameterizedFunction(parameters); + + for (var i = 0; i < data.x.length; i++) { + error += Math.abs(data.y[i] - func(data.x[i])); + } + + return error; +} + +/** + * Difference of the matrix function over the parameters + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} evaluatedData - Array of previous evaluated function values + * @param {Array} params - Array of previous parameter values + * @param {number} gradientDifference - Adjustment for decrease the damping parameter + * @param {function} paramFunction - The parameters and returns a function with the independent variable as a parameter + * @return {Matrix} + */ +function gradientFunction( + data, + evaluatedData, + params, + gradientDifference, + paramFunction +) { + const n = params.length; + const m = data.x.length; + + var ans = new Array(n); + + for (var param = 0; param < n; param++) { + ans[param] = new Array(m); + var auxParams = params.concat(); + auxParams[param] += gradientDifference; + var funcParam = paramFunction(auxParams); + + for (var point = 0; point < m; point++) { + ans[param][point] = evaluatedData[point] - funcParam(data.x[point]); + } + } + return new mlMatrix.Matrix(ans); +} + +/** + * Matrix function over the samples + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} evaluatedData - Array of previous evaluated function values + * @return {Matrix} + */ +function matrixFunction(data, evaluatedData) { + const m = data.x.length; + + var ans = new Array(m); + + for (var point = 0; point < m; point++) { + ans[point] = data.y[point] - evaluatedData[point]; + } + + return new mlMatrix.Matrix([ans]); +} + +/** + * Iteration for Levenberg-Marquardt + * @ignore + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {Array} params - Array of previous parameter values + * @param {number} damping - Levenberg-Marquardt parameter + * @param {number} gradientDifference - Adjustment for decrease the damping parameter + * @param {function} parameterizedFunction - The parameters and returns a function with the independent variable as a parameter + * @return {Array} + */ +function step( + data, + params, + damping, + gradientDifference, + parameterizedFunction +) { + var identity = mlMatrix.Matrix.eye(params.length).mul( + damping * gradientDifference * gradientDifference + ); + + var l = data.x.length; + var evaluatedData = new Array(l); + const func = parameterizedFunction(params); + for (var i = 0; i < l; i++) { + evaluatedData[i] = func(data.x[i]); + } + var gradientFunc = gradientFunction( + data, + evaluatedData, + params, + gradientDifference, + parameterizedFunction + ); + var matrixFunc = matrixFunction(data, evaluatedData).transposeView(); + var inverseMatrix = mlMatrix.inverse( + identity.add(gradientFunc.mmul(gradientFunc.transposeView())) + ); + params = new mlMatrix.Matrix([params]); + params = params.sub( + inverseMatrix + .mmul(gradientFunc) + .mmul(matrixFunc) + .mul(gradientDifference) + .transposeView() + ); + + return params.to1DArray(); +} + +/** + * Curve fitting algorithm + * @param {{x:Array, y:Array}} data - Array of points to fit in the format [x1, x2, ... ], [y1, y2, ... ] + * @param {function} parameterizedFunction - The parameters and returns a function with the independent variable as a parameter + * @param {object} [options] - Options object + * @param {number} [options.damping] - Levenberg-Marquardt parameter + * @param {number} [options.gradientDifference = 10e-2] - Adjustment for decrease the damping parameter + * @param {Array} [options.initialValues] - Array of initial parameter values + * @param {number} [options.maxIterations = 100] - Maximum of allowed iterations + * @param {number} [options.errorTolerance = 10e-3] - Minimum uncertainty allowed for each point + * @return {{parameterValues: Array, parameterError: number, iterations: number}} + */ +function levenbergMarquardt( + data, + parameterizedFunction, + options = {} +) { + let { + maxIterations = 100, + gradientDifference = 10e-2, + damping = 0, + errorTolerance = 10e-3, + initialValues + } = options; + + if (damping <= 0) { + throw new Error('The damping option must be a positive number'); + } else if (!data.x || !data.y) { + throw new Error('The data parameter must have x and y elements'); + } else if ( + !Array.isArray(data.x) || + data.x.length < 2 || + !Array.isArray(data.y) || + data.y.length < 2 + ) { + throw new Error( + 'The data parameter elements must be an array with more than 2 points' + ); + } else { + let dataLen = data.x.length; + if (dataLen !== data.y.length) { + throw new Error('The data parameter elements must have the same size'); + } + } + + var parameters = + initialValues || new Array(parameterizedFunction.length).fill(1); + + if (!Array.isArray(parameters)) { + throw new Error('initialValues must be an array'); + } + + var error = errorCalculation(data, parameters, parameterizedFunction); + + var converged = error <= errorTolerance; + + for ( + var iteration = 0; + iteration < maxIterations && !converged; + iteration++ + ) { + parameters = step( + data, + parameters, + damping, + gradientDifference, + parameterizedFunction + ); + error = errorCalculation(data, parameters, parameterizedFunction); + converged = error <= errorTolerance; + } + + return { + parameterValues: parameters, + parameterError: error, + iterations: iteration + }; +} + +module.exports = levenbergMarquardt; + + +/***/ }), +/* 9 */ +/***/ (function(module, __webpack_exports__, __webpack_require__) { + +"use strict"; +__webpack_require__.r(__webpack_exports__); + +// EXTERNAL MODULE: ./node_modules/is-any-array/src/index.js +var src = __webpack_require__(0); +var src_default = /*#__PURE__*/__webpack_require__.n(src); + +// CONCATENATED MODULE: ./node_modules/ml-array-max/lib-es6/index.js + + +/** + * Computes the maximum of the given values + * @param {Array} input + * @return {number} + */ + +function lib_es6_max(input) { + if (!src_default()(input)) { + throw new TypeError('input must be an array'); + } + + if (input.length === 0) { + throw new TypeError('input must not be empty'); + } + + var max = input[0]; + + for (var i = 1; i < input.length; i++) { + if (input[i] > max) max = input[i]; + } + + return max; +} + +/* harmony default export */ var lib_es6 = (lib_es6_max); + +// CONCATENATED MODULE: ./node_modules/ml-array-min/lib-es6/index.js + + +/** + * Computes the minimum of the given values + * @param {Array} input + * @return {number} + */ + +function lib_es6_min(input) { + if (!src_default()(input)) { + throw new TypeError('input must be an array'); + } + + if (input.length === 0) { + throw new TypeError('input must not be empty'); + } + + var min = input[0]; + + for (var i = 1; i < input.length; i++) { + if (input[i] < min) min = input[i]; + } + + return min; +} + +/* harmony default export */ var ml_array_min_lib_es6 = (lib_es6_min); + +// CONCATENATED MODULE: ./node_modules/ml-array-rescale/lib-es6/index.js + + + + +function rescale(input) { + var options = arguments.length > 1 && arguments[1] !== undefined ? arguments[1] : {}; + + if (!src_default()(input)) { + throw new TypeError('input must be an array'); + } else if (input.length === 0) { + throw new TypeError('input must not be empty'); + } + + var output; + + if (options.output !== undefined) { + if (!src_default()(options.output)) { + throw new TypeError('output option must be an array if specified'); + } + + output = options.output; + } else { + output = new Array(input.length); + } + + var currentMin = ml_array_min_lib_es6(input); + var currentMax = lib_es6(input); + + if (currentMin === currentMax) { + throw new RangeError('minimum and maximum input values are equal. Cannot rescale a constant array'); + } + + var _options$min = options.min, + minValue = _options$min === void 0 ? options.autoMinMax ? currentMin : 0 : _options$min, + _options$max = options.max, + maxValue = _options$max === void 0 ? options.autoMinMax ? currentMax : 1 : _options$max; + + if (minValue >= maxValue) { + throw new RangeError('min option must be smaller than max option'); + } + + var factor = (maxValue - minValue) / (currentMax - currentMin); + + for (var i = 0; i < input.length; i++) { + output[i] = (input[i] - currentMin) * factor + minValue; + } + + return output; +} + +/* harmony default export */ var ml_array_rescale_lib_es6 = (rescale); + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/dc/lu.js + + +/** + * @class LuDecomposition + * @link https://github.com/lutzroeder/Mapack/blob/master/Source/LuDecomposition.cs + * @param {Matrix} matrix + */ +class lu_LuDecomposition { + constructor(matrix) { + matrix = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(matrix); + + var lu = matrix.clone(); + var rows = lu.rows; + var columns = lu.columns; + var pivotVector = new Array(rows); + var pivotSign = 1; + var i, j, k, p, s, t, v; + var 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.get(i, j); + } + + for (i = 0; i < rows; i++) { + kmax = Math.min(i, j); + s = 0; + for (k = 0; k < kmax; k++) { + s += lu.get(i, k) * LUcolj[k]; + } + LUcolj[i] -= s; + lu.set(i, j, LUcolj[i]); + } + + 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.get(p, k); + lu.set(p, k, lu.get(j, k)); + lu.set(j, k, t); + } + + v = pivotVector[p]; + pivotVector[p] = pivotVector[j]; + pivotVector[j] = v; + + pivotSign = -pivotSign; + } + + if (j < rows && lu.get(j, j) !== 0) { + for (i = j + 1; i < rows; i++) { + lu.set(i, j, lu.get(i, j) / lu.get(j, j)); + } + } + } + + this.LU = lu; + this.pivotVector = pivotVector; + this.pivotSign = pivotSign; + } + + /** + * + * @return {boolean} + */ + isSingular() { + var data = this.LU; + var col = data.columns; + for (var j = 0; j < col; j++) { + if (data[j][j] === 0) { + return true; + } + } + return false; + } + + /** + * + * @param {Matrix} value + * @return {Matrix} + */ + solve(value) { + value = matrix_Matrix.checkMatrix(value); + + var lu = this.LU; + var 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; + } + + /** + * + * @return {number} + */ + get determinant() { + var data = this.LU; + if (!data.isSquare()) { + throw new Error('Matrix must be square'); + } + var determinant = this.pivotSign; + var col = data.columns; + for (var j = 0; j < col; j++) { + determinant *= data[j][j]; + } + return determinant; + } + + /** + * + * @return {Matrix} + */ + get lowerTriangularMatrix() { + var data = this.LU; + var rows = data.rows; + var columns = data.columns; + var 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; + } + + /** + * + * @return {Matrix} + */ + get upperTriangularMatrix() { + var data = this.LU; + var rows = data.rows; + var columns = data.columns; + var 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; + } + + /** + * + * @return {Array} + */ + get pivotPermutationVector() { + return this.pivotVector.slice(); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/dc/util.js +function hypotenuse(a, b) { + var r = 0; + 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; +} + +function getFilled2DArray(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; +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/dc/svd.js + + + + +/** + * @class SingularValueDecomposition + * @see https://github.com/accord-net/framework/blob/development/Sources/Accord.Math/Decompositions/SingularValueDecomposition.cs + * @param {Matrix} value + * @param {object} [options] + * @param {boolean} [options.computeLeftSingularVectors=true] + * @param {boolean} [options.computeRightSingularVectors=true] + * @param {boolean} [options.autoTranspose=false] + */ +class svd_SingularValueDecomposition { + constructor(value, options = {}) { + value = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(value); + + var m = value.rows; + var n = value.columns; + + const { + computeLeftSingularVectors = true, + computeRightSingularVectors = true, + autoTranspose = false + } = options; + + var wantu = Boolean(computeLeftSingularVectors); + var wantv = Boolean(computeRightSingularVectors); + + 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 nu = Math.min(m, n); + var ni = Math.min(m + 1, n); + var s = new Array(ni); + var U = getFilled2DArray(m, nu, 0); + var V = getFilled2DArray(n, n, 0); + + var e = new Array(n); + var work = new Array(m); + + var si = new Array(ni); + for (let i = 0; i < ni; i++) si[i] = i; + + var nct = Math.min(m - 1, n); + var nrt = Math.max(0, Math.min(n - 2, m)); + var mrc = Math.max(nct, nrt); + + for (let k = 0; k < mrc; k++) { + if (k < nct) { + s[k] = 0; + for (let 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 (let i = k; i < m; i++) { + a[i][k] /= s[k]; + } + a[k][k] += 1; + } + s[k] = -s[k]; + } + + for (let j = k + 1; j < n; j++) { + if (k < nct && s[k] !== 0) { + let t = 0; + for (let i = k; i < m; i++) { + t += a[i][k] * a[i][j]; + } + t = -t / a[k][k]; + for (let i = k; i < m; i++) { + a[i][j] += t * a[i][k]; + } + } + e[j] = a[k][j]; + } + + if (wantu && k < nct) { + for (let i = k; i < m; i++) { + U[i][k] = a[i][k]; + } + } + + if (k < nrt) { + e[k] = 0; + for (let 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 (let 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 (let i = k + 1; i < m; i++) { + work[i] = 0; + } + for (let i = k + 1; i < m; i++) { + for (let j = k + 1; j < n; j++) { + work[i] += e[j] * a[i][j]; + } + } + for (let j = k + 1; j < n; j++) { + let t = -e[j] / e[k + 1]; + for (let i = k + 1; i < m; i++) { + a[i][j] += t * work[i]; + } + } + } + if (wantv) { + for (let i = k + 1; i < n; i++) { + V[i][k] = e[i]; + } + } + } + } + + let 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 (let j = nct; j < nu; j++) { + for (let i = 0; i < m; i++) { + U[i][j] = 0; + } + U[j][j] = 1; + } + for (let k = nct - 1; k >= 0; k--) { + if (s[k] !== 0) { + for (let j = k + 1; j < nu; j++) { + let t = 0; + for (let i = k; i < m; i++) { + t += U[i][k] * U[i][j]; + } + t = -t / U[k][k]; + for (let i = k; i < m; i++) { + U[i][j] += t * U[i][k]; + } + } + for (let i = k; i < m; i++) { + U[i][k] = -U[i][k]; + } + U[k][k] = 1 + U[k][k]; + for (let i = 0; i < k - 1; i++) { + U[i][k] = 0; + } + } else { + for (let i = 0; i < m; i++) { + U[i][k] = 0; + } + U[k][k] = 1; + } + } + } + + if (wantv) { + for (let k = n - 1; k >= 0; k--) { + if (k < nrt && e[k] !== 0) { + for (let j = k + 1; j < n; j++) { + let t = 0; + for (let i = k + 1; i < n; i++) { + t += V[i][k] * V[i][j]; + } + t = -t / V[k + 1][k]; + for (let i = k + 1; i < n; i++) { + V[i][j] += t * V[i][k]; + } + } + } + for (let i = 0; i < n; i++) { + V[i][k] = 0; + } + V[k][k] = 1; + } + } + + var pp = p - 1; + var iter = 0; + var eps = Number.EPSILON; + while (p > 0) { + let k, kase; + for (k = p - 2; k >= -1; k--) { + if (k === -1) { + break; + } + const alpha = + Number.MIN_VALUE + eps * Math.abs(s[k] + Math.abs(s[k + 1])); + if (Math.abs(e[k]) <= alpha || Number.isNaN(e[k])) { + e[k] = 0; + break; + } + } + if (k === p - 2) { + kase = 4; + } else { + let ks; + for (ks = p - 1; ks >= k; ks--) { + if (ks === k) { + break; + } + let 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: { + let f = e[p - 2]; + e[p - 2] = 0; + for (let j = p - 2; j >= k; j--) { + let t = hypotenuse(s[j], f); + let cs = s[j] / t; + let sn = f / t; + s[j] = t; + if (j !== k) { + f = -sn * e[j - 1]; + e[j - 1] = cs * e[j - 1]; + } + if (wantv) { + for (let 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: { + let f = e[k - 1]; + e[k - 1] = 0; + for (let j = k; j < p; j++) { + let t = hypotenuse(s[j], f); + let cs = s[j] / t; + let sn = f / t; + s[j] = t; + f = -sn * e[j]; + e[j] = cs * e[j]; + if (wantu) { + for (let 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: { + const scale = 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]) + ); + const sp = s[p - 1] / scale; + const spm1 = s[p - 2] / scale; + const epm1 = e[p - 2] / scale; + const sk = s[k] / scale; + const ek = e[k] / scale; + const b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2; + const c = sp * epm1 * (sp * epm1); + let shift = 0; + if (b !== 0 || c !== 0) { + if (b < 0) { + shift = 0 - Math.sqrt(b * b + c); + } else { + shift = Math.sqrt(b * b + c); + } + shift = c / (b + shift); + } + let f = (sk + sp) * (sk - sp) + shift; + let g = sk * ek; + for (let j = k; j < p - 1; j++) { + let t = hypotenuse(f, g); + if (t === 0) t = Number.MIN_VALUE; + let cs = f / t; + let 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 (let 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); + if (t === 0) t = Number.MIN_VALUE; + 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 (let 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 (let i = 0; i <= pp; i++) { + V[i][k] = -V[i][k]; + } + } + } + while (k < pp) { + if (s[k] >= s[k + 1]) { + break; + } + let t = s[k]; + s[k] = s[k + 1]; + s[k + 1] = t; + if (wantv && k < n - 1) { + for (let 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 (let 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; + } + + /** + * Solve a problem of least square (Ax=b) by using the SVD. Useful when A is singular. When A is not singular, it would be better to use qr.solve(value). + * Example : We search to approximate x, with A matrix shape m*n, x vector size n, b vector size m (m > n). We will use : + * var svd = SingularValueDecomposition(A); + * var x = svd.solve(b); + * @param {Matrix} value - Matrix 1D which is the vector b (in the equation Ax = b) + * @return {Matrix} - The vector x + */ + solve(value) { + var Y = value; + var e = this.threshold; + var scols = this.s.length; + var Ls = matrix_Matrix.zeros(scols, scols); + + for (let 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); + var vrows = V.rows; + var urows = U.length; + var VLU = matrix_Matrix.zeros(vrows, urows); + + for (let i = 0; i < vrows; i++) { + for (let j = 0; j < urows; j++) { + let sum = 0; + for (let k = 0; k < scols; k++) { + sum += VL[i][k] * U[j][k]; + } + VLU[i][j] = sum; + } + } + + return VLU.mmul(Y); + } + + /** + * + * @param {Array} value + * @return {Matrix} + */ + solveForDiagonal(value) { + return this.solve(matrix_Matrix.diag(value)); + } + + /** + * Get the inverse of the matrix. We compute the inverse of a matrix using SVD when this matrix is singular or ill-conditioned. Example : + * var svd = SingularValueDecomposition(A); + * var inverseA = svd.inverse(); + * @return {Matrix} - The approximation of the inverse of the matrix + */ + inverse() { + var V = this.V; + var e = this.threshold; + var vrows = V.length; + var vcols = V[0].length; + var X = new matrix_Matrix(vrows, this.s.length); + + for (let i = 0; i < vrows; i++) { + for (let 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; + var ucols = U[0].length; + var Y = new matrix_Matrix(vrows, urows); + + for (let i = 0; i < vrows; i++) { + for (let j = 0; j < urows; j++) { + let sum = 0; + for (let k = 0; k < ucols; k++) { + sum += X[i][k] * U[j][k]; + } + Y[i][j] = sum; + } + } + + return Y; + } + + /** + * + * @return {number} + */ + get condition() { + return this.s[0] / this.s[Math.min(this.m, this.n) - 1]; + } + + /** + * + * @return {number} + */ + get norm2() { + return this.s[0]; + } + + /** + * + * @return {number} + */ + get rank() { + var tol = Math.max(this.m, this.n) * this.s[0] * Number.EPSILON; + var r = 0; + var s = this.s; + for (var i = 0, ii = s.length; i < ii; i++) { + if (s[i] > tol) { + r++; + } + } + return r; + } + + /** + * + * @return {Array} + */ + get diagonal() { + return this.s; + } + + /** + * + * @return {number} + */ + get threshold() { + return Number.EPSILON / 2 * Math.max(this.m, this.n) * this.s[0]; + } + + /** + * + * @return {Matrix} + */ + get leftSingularVectors() { + if (!matrix_Matrix.isMatrix(this.U)) { + this.U = new matrix_Matrix(this.U); + } + return this.U; + } + + /** + * + * @return {Matrix} + */ + get rightSingularVectors() { + if (!matrix_Matrix.isMatrix(this.V)) { + this.V = new matrix_Matrix(this.V); + } + return this.V; + } + + /** + * + * @return {Matrix} + */ + get diagonalMatrix() { + return matrix_Matrix.diag(this.s); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/util.js + + +/** + * @private + * Check that a row index is not out of bounds + * @param {Matrix} matrix + * @param {number} index + * @param {boolean} [outer] + */ +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] + */ +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} + */ +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} + */ +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; +} + +function checkIndices(matrix, rowIndices, columnIndices) { + return { + row: checkRowIndices(matrix, rowIndices), + column: checkColumnIndices(matrix, columnIndices) + }; +} + +function checkRowIndices(matrix, rowIndices) { + if (typeof rowIndices !== 'object') { + throw new TypeError('unexpected type for row indices'); + } + + var rowOut = rowIndices.some((r) => { + return r < 0 || r >= matrix.rows; + }); + + if (rowOut) { + throw new RangeError('row indices are out of range'); + } + + if (!Array.isArray(rowIndices)) rowIndices = Array.from(rowIndices); + + return rowIndices; +} + +function checkColumnIndices(matrix, columnIndices) { + if (typeof columnIndices !== 'object') { + throw new TypeError('unexpected type for column indices'); + } + + var columnOut = columnIndices.some((c) => { + return c < 0 || c >= matrix.columns; + }); + + if (columnOut) { + throw new RangeError('column indices are out of range'); + } + if (!Array.isArray(columnIndices)) columnIndices = Array.from(columnIndices); + + return columnIndices; +} + +function checkRange(matrix, startRow, endRow, startColumn, endColumn) { + if (arguments.length !== 5) { + throw new RangeError('expected 4 arguments'); + } + checkNumber('startRow', startRow); + checkNumber('endRow', endRow); + checkNumber('startColumn', startColumn); + checkNumber('endColumn', endColumn); + 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'); + } +} + +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; +} + +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; +} + +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; +} + +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; +} + +function checkNumber(name, value) { + if (typeof value !== 'number') { + throw new TypeError(`${name} must be a number`); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/base.js + + + +class base_BaseView extends AbstractMatrix() { + constructor(matrix, rows, columns) { + super(); + this.matrix = matrix; + this.rows = rows; + this.columns = columns; + } + + static get [Symbol.species]() { + return matrix_Matrix; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/transpose.js + + +class transpose_MatrixTransposeView extends base_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); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/row.js + + +class row_MatrixRowView extends base_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); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/sub.js + + + + +class sub_MatrixSubView extends base_BaseView { + constructor(matrix, startRow, endRow, startColumn, endColumn) { + 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 + ); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/selection.js + + + + +class selection_MatrixSelectionView extends base_BaseView { + constructor(matrix, rowIndices, columnIndices) { + var indices = 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] + ); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/rowSelection.js + + + + +class rowSelection_MatrixRowSelectionView extends base_BaseView { + constructor(matrix, rowIndices) { + rowIndices = checkRowIndices(matrix, rowIndices); + super(matrix, rowIndices.length, matrix.columns); + this.rowIndices = rowIndices; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(this.rowIndices[rowIndex], columnIndex, value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(this.rowIndices[rowIndex], columnIndex); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/columnSelection.js + + + + +class columnSelection_MatrixColumnSelectionView extends base_BaseView { + constructor(matrix, columnIndices) { + columnIndices = checkColumnIndices(matrix, columnIndices); + super(matrix, matrix.rows, columnIndices.length); + this.columnIndices = columnIndices; + } + + set(rowIndex, columnIndex, value) { + this.matrix.set(rowIndex, this.columnIndices[columnIndex], value); + return this; + } + + get(rowIndex, columnIndex) { + return this.matrix.get(rowIndex, this.columnIndices[columnIndex]); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/column.js + + +class column_MatrixColumnView extends base_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); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/flipRow.js + + +class flipRow_MatrixFlipRowView extends base_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); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/views/flipColumn.js + + +class flipColumn_MatrixFlipColumnView extends base_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); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/abstractMatrix.js + + + + + + + + + + + + + + + +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) { + 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) { + checkRowIndex(this, index); + array = 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) { + checkRowIndex(this, row1); + 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) { + 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) { + checkColumnIndex(this, index); + array = 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) { + checkColumnIndex(this, column1); + 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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 = 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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) { + 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 sumByRow(this); + case 'column': + return sumByColumn(this); + default: + return 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; + } + + /** + * Returns the norm of a matrix. + * @param {string} type - "frobenius" (default) or "max" return resp. the Frobenius norm and the max norm. + * @return {number} + */ + norm(type = 'frobenius') { + var result = 0; + if (type === 'max') { + return this.max(); + } else if (type === 'frobenius') { + for (var i = 0; i < this.rows; i++) { + for (var j = 0; j < this.columns; j++) { + result = result + this.get(i, j) * this.get(i, j); + } + } + return Math.sqrt(result); + } else { + throw new RangeError(`unknown norm type: ${type}`); + } + } + + /** + * 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, 10); + var halfCols = parseInt(a.columns / 2, 10); + // 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 = ml_array_rescale_lib_es6(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 = ml_array_rescale_lib_es6(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) { + 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; + 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} rowIndices - The row indices to select. Order matters and an index can be more than once. + * @param {Array} 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 = 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 transpose_MatrixTransposeView(this); + } + + /** + * Returns a view of the row vector with the given index + * @param {number} row - row index of the vector + * @return {MatrixRowView} + */ + rowView(row) { + checkRowIndex(this, row); + return new row_MatrixRowView(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) { + checkColumnIndex(this, column); + return new column_MatrixColumnView(this, column); + } + + /** + * Returns a view of the matrix flipped in the row axis + * @return {MatrixFlipRowView} + */ + flipRowView() { + return new flipRow_MatrixFlipRowView(this); + } + + /** + * Returns a view of the matrix flipped in the column axis + * @return {MatrixFlipColumnView} + */ + flipColumnView() { + return new flipColumn_MatrixFlipColumnView(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 sub_MatrixSubView(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} rowIndices + * @param {Array} columnIndices + * @return {MatrixSelectionView} + */ + selectionView(rowIndices, columnIndices) { + return new selection_MatrixSelectionView(this, rowIndices, columnIndices); + } + + /** + * Returns a view of the row indices + * @example + * // resulting vector is [[1,2,3], [1,2,3]] + * var matrix = new Matrix([[1,2,3], [4,5,6]]).rowSelectionView([0, 0]) + * @param {Array} rowIndices + * @return {MatrixRowSelectionView} + */ + rowSelectionView(rowIndices) { + return new rowSelection_MatrixRowSelectionView(this, rowIndices); + } + + /** + * Returns a view of the column indices + * @example + * // resulting vector is [[2, 2], [5, 5]] + * var matrix = new Matrix([[1,2,3], [4,5,6]]).columnSelectionView([1, 1]) + * @param {Array} columnIndices + * @return {MatrixColumnSelectionView} + */ + columnSelectionView(columnIndices) { + return new columnSelection_MatrixColumnSelectionView(this, 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 lu_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 svd_SingularValueDecomposition(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())); + } + + /** + * 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; + } + } + + Matrix.prototype.klass = 'Matrix'; + + 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); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + 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); + if (this.rows !== matrix.rows || + this.columns !== matrix.columns) { + throw new RangeError('Matrices dimensions must be equal'); + } + 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; + var eval2 = eval; // eslint-disable-line no-eval + for (var operator of operators) { + var inplaceOp = eval2(fillTemplateFunction(inplaceOperator, { name: operator[1], op: operator[0] })); + var inplaceOpS = eval2(fillTemplateFunction(inplaceOperatorScalar, { name: `${operator[1]}S`, op: operator[0] })); + var inplaceOpM = eval2(fillTemplateFunction(inplaceOperatorMatrix, { name: `${operator[1]}M`, op: operator[0] })); + var staticOp = eval2(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 = eval2(fillTemplateFunction(inplaceMethod, { name: method[1], method: method[0] })); + var staticMeth = eval2(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 = eval2(fillTemplateFunction(inplaceMethodWithArgs, { + name: methodWithArg[2], + method: methodWithArg[0], + args: args + })); + var staticMethWithArgs = eval2(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 = eval2(fillTemplateFunction(inplaceMethodWithOneArg, tmplVar)); + var inplaceMethodS = eval2(fillTemplateFunction(inplaceMethodWithOneArgScalar, tmplVar)); + var inplaceMethodM = eval2(fillTemplateFunction(inplaceMethodWithOneArgMatrix, tmplVar)); + var staticMethod2 = eval2(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; +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/matrix.js + + + +class matrix_Matrix extends AbstractMatrix(Array) { + constructor(nRows, nColumns) { + var i; + if (arguments.length === 1 && typeof nRows === 'number') { + return new Array(nRows); + } + if (matrix_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]; + } + + /** + * Removes a row from the given index + * @param {number} index - Row index + * @return {Matrix} this + */ + removeRow(index) { + 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; + } + checkRowIndex(this, index, true); + array = 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) { + 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; + } + checkColumnIndex(this, index, true); + array = checkColumnVector(this, array); + for (var i = 0; i < this.rows; i++) { + this[i].splice(index, 0, array[i]); + } + this.columns += 1; + return this; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/wrap/WrapperMatrix1D.js + + + +class WrapperMatrix1D_WrapperMatrix1D extends AbstractMatrix() { + /** + * @class WrapperMatrix1D + * @param {Array} data + * @param {object} [options] + * @param {object} [options.rows = 1] + */ + constructor(data, options = {}) { + const { rows = 1 } = options; + + if (data.length % rows !== 0) { + throw new Error('the data length is not divisible by the number of rows'); + } + super(); + this.rows = rows; + this.columns = data.length / rows; + this.data = data; + } + + set(rowIndex, columnIndex, value) { + var index = this._calculateIndex(rowIndex, columnIndex); + this.data[index] = value; + return this; + } + + get(rowIndex, columnIndex) { + var index = this._calculateIndex(rowIndex, columnIndex); + return this.data[index]; + } + + _calculateIndex(row, column) { + return row * this.columns + column; + } + + static get [Symbol.species]() { + return matrix_Matrix; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/wrap/WrapperMatrix2D.js + + + +class WrapperMatrix2D_WrapperMatrix2D extends AbstractMatrix() { + /** + * @class WrapperMatrix2D + * @param {Array>} data + */ + constructor(data) { + super(); + this.data = data; + this.rows = data.length; + this.columns = data[0].length; + } + + set(rowIndex, columnIndex, value) { + this.data[rowIndex][columnIndex] = value; + return this; + } + + get(rowIndex, columnIndex) { + return this.data[rowIndex][columnIndex]; + } + + static get [Symbol.species]() { + return matrix_Matrix; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/wrap/wrap.js + + + +/** + * @param {Array>|Array} array + * @param {object} [options] + * @param {object} [options.rows = 1] + * @return {WrapperMatrix1D|WrapperMatrix2D} + */ +function wrap(array, options) { + if (Array.isArray(array)) { + if (array[0] && Array.isArray(array[0])) { + return new WrapperMatrix2D_WrapperMatrix2D(array); + } else { + return new WrapperMatrix1D_WrapperMatrix1D(array, options); + } + } else { + throw new Error('the argument is not an array'); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/dc/qr.js + + + + +/** + * @class QrDecomposition + * @link https://github.com/lutzroeder/Mapack/blob/master/Source/QrDecomposition.cs + * @param {Matrix} value + */ +class qr_QrDecomposition { + constructor(value) { + value = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(value); + + var qr = value.clone(); + var m = value.rows; + var n = value.columns; + var rdiag = new Array(n); + var i, j, k, s; + + for (k = 0; k < n; k++) { + var nrm = 0; + for (i = k; i < m; i++) { + nrm = hypotenuse(nrm, qr.get(i, k)); + } + if (nrm !== 0) { + if (qr.get(k, k) < 0) { + nrm = -nrm; + } + for (i = k; i < m; i++) { + qr.set(i, k, qr.get(i, k) / nrm); + } + qr.set(k, k, qr.get(k, k) + 1); + for (j = k + 1; j < n; j++) { + s = 0; + for (i = k; i < m; i++) { + s += qr.get(i, k) * qr.get(i, j); + } + s = -s / qr.get(k, k); + for (i = k; i < m; i++) { + qr.set(i, j, qr.get(i, j) + s * qr.get(i, k)); + } + } + } + rdiag[k] = -nrm; + } + + this.QR = qr; + this.Rdiag = rdiag; + } + + /** + * Solve a problem of least square (Ax=b) by using the QR decomposition. Useful when A is rectangular, but not working when A is singular. + * Example : We search to approximate x, with A matrix shape m*n, x vector size n, b vector size m (m > n). We will use : + * var qr = QrDecomposition(A); + * var x = qr.solve(b); + * @param {Matrix} value - Matrix 1D which is the vector b (in the equation Ax = b) + * @return {Matrix} - The vector x + */ + solve(value) { + value = matrix_Matrix.checkMatrix(value); + + var qr = this.QR; + var 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); + } + + /** + * + * @return {boolean} + */ + isFullRank() { + var columns = this.QR.columns; + for (var i = 0; i < columns; i++) { + if (this.Rdiag[i] === 0) { + return false; + } + } + return true; + } + + /** + * + * @return {Matrix} + */ + get upperTriangularMatrix() { + var qr = this.QR; + var n = qr.columns; + var X = new matrix_Matrix(n, n); + var 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; + } + + /** + * + * @return {Matrix} + */ + get orthogonalMatrix() { + var qr = this.QR; + var rows = qr.rows; + var columns = qr.columns; + var X = new matrix_Matrix(rows, columns); + var 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; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/decompositions.js + + + + + + +/** + * Computes the inverse of a Matrix + * @param {Matrix} matrix + * @param {boolean} [useSVD=false] + * @return {Matrix} + */ +function inverse(matrix, useSVD = false) { + matrix = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(matrix); + if (useSVD) { + return new svd_SingularValueDecomposition(matrix).inverse(); + } else { + return solve(matrix, matrix_Matrix.eye(matrix.rows)); + } +} + +/** + * + * @param {Matrix} leftHandSide + * @param {Matrix} rightHandSide + * @param {boolean} [useSVD = false] + * @return {Matrix} + */ +function solve(leftHandSide, rightHandSide, useSVD = false) { + leftHandSide = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(leftHandSide); + rightHandSide = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(rightHandSide); + if (useSVD) { + return new svd_SingularValueDecomposition(leftHandSide).solve(rightHandSide); + } else { + return leftHandSide.isSquare() + ? new lu_LuDecomposition(leftHandSide).solve(rightHandSide) + : new qr_QrDecomposition(leftHandSide).solve(rightHandSide); + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/linearDependencies.js + + + + + +// function used by rowsDependencies +function xrange(n, exception) { + var range = []; + for (var i = 0; i < n; i++) { + if (i !== exception) { + range.push(i); + } + } + return range; +} + +// function used by rowsDependencies +function dependenciesOneRow( + error, + matrix, + index, + thresholdValue = 10e-10, + thresholdError = 10e-10 +) { + if (error > thresholdError) { + return new Array(matrix.rows + 1).fill(0); + } else { + var returnArray = matrix.addRow(index, [0]); + for (var i = 0; i < returnArray.rows; i++) { + if (Math.abs(returnArray.get(i, 0)) < thresholdValue) { + returnArray.set(i, 0, 0); + } + } + return returnArray.to1DArray(); + } +} + +/** + * Creates a matrix which represents the dependencies between rows. + * If a row is a linear combination of others rows, the result will be a row with the coefficients of this combination. + * For example : for A = [[2, 0, 0, 1], [0, 1, 6, 0], [0, 3, 0, 1], [0, 0, 1, 0], [0, 1, 2, 0]], the result will be [[0, 0, 0, 0, 0], [0, 0, 0, 4, 1], [0, 0, 0, 0, 0], [0, 0.25, 0, 0, -0.25], [0, 1, 0, -4, 0]] + * @param {Matrix} matrix + * @param {Object} [options] includes thresholdValue and thresholdError. + * @param {number} [options.thresholdValue = 10e-10] If an absolute value is inferior to this threshold, it will equals zero. + * @param {number} [options.thresholdError = 10e-10] If the error is inferior to that threshold, the linear combination found is accepted and the row is dependent from other rows. + * @return {Matrix} the matrix which represents the dependencies between rows. + */ + +function linearDependencies(matrix, options = {}) { + const { thresholdValue = 10e-10, thresholdError = 10e-10 } = options; + + var n = matrix.rows; + var results = new matrix_Matrix(n, n); + + for (var i = 0; i < n; i++) { + var b = matrix_Matrix.columnVector(matrix.getRow(i)); + var Abis = matrix.subMatrixRow(xrange(n, i)).transposeView(); + var svd = new svd_SingularValueDecomposition(Abis); + var x = svd.solve(b); + var error = lib_es6( + matrix_Matrix.sub(b, Abis.mmul(x)) + .abs() + .to1DArray() + ); + results.setRow( + i, + dependenciesOneRow(error, x, i, thresholdValue, thresholdError) + ); + } + return results; +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/dc/evd.js + + + + +/** + * @class EigenvalueDecomposition + * @link https://github.com/lutzroeder/Mapack/blob/master/Source/EigenvalueDecomposition.cs + * @param {Matrix} matrix + * @param {object} [options] + * @param {boolean} [options.assumeSymmetric=false] + */ +class evd_EigenvalueDecomposition { + constructor(matrix, options = {}) { + const { assumeSymmetric = false } = options; + + matrix = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(matrix); + if (!matrix.isSquare()) { + throw new Error('Matrix is not a square matrix'); + } + + var n = matrix.columns; + var V = getFilled2DArray(n, n, 0); + var d = new Array(n); + var e = new Array(n); + var value = matrix; + var i, j; + + var isSymmetric = false; + if (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); + var 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; + } + + /** + * + * @return {Array} + */ + get realEigenvalues() { + return this.d; + } + + /** + * + * @return {Array} + */ + get imaginaryEigenvalues() { + return this.e; + } + + /** + * + * @return {Matrix} + */ + get eigenvectorMatrix() { + if (!matrix_Matrix.isMatrix(this.V)) { + this.V = new matrix_Matrix(this.V); + } + return this.V; + } + + /** + * + * @return {Matrix} + */ + get diagonalMatrix() { + var n = this.n; + var e = this.e; + var d = this.d; + var X = new matrix_Matrix(n, n); + var 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; + var tst1 = 0; + var eps = Number.EPSILON; + + 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; + var high = n - 1; + var f, g, h, i, j, m; + var 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; + var low = 0; + var high = nn - 1; + var eps = Number.EPSILON; + var exshift = 0; + var norm = 0; + var p = 0; + var q = 0; + var r = 0; + var s = 0; + var z = 0; + var iter = 0; + var i, j, k, l, m, t, w, x, y; + var ra, sa, vr, vi; + var 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]; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/dc/cholesky.js + + +/** + * @class CholeskyDecomposition + * @link https://github.com/lutzroeder/Mapack/blob/master/Source/CholeskyDecomposition.cs + * @param {Matrix} value + */ +class cholesky_CholeskyDecomposition { + constructor(value) { + value = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(value); + if (!value.isSymmetric()) { + throw new Error('Matrix is not symmetric'); + } + + var a = value; + var dimension = a.rows; + var l = new matrix_Matrix(dimension, dimension); + var positiveDefinite = true; + var 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.get(j, k) - s) / l[k][k]; + d = d + s * s; + } + + d = a.get(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; + } + + /** + * + * @param {Matrix} value + * @return {Matrix} + */ + solve(value) { + value = WrapperMatrix2D_WrapperMatrix2D.checkMatrix(value); + + var l = this.L; + var dimension = l.rows; + + if (value.rows !== dimension) { + throw new Error('Matrix dimensions do not match'); + } + + var count = value.columns; + var B = value.clone(); + var 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; + } + + /** + * + * @return {Matrix} + */ + get lowerTriangularMatrix() { + return this.L; + } +} + +// CONCATENATED MODULE: ./node_modules/ml-matrix/src/index.js +/* concated harmony reexport default */__webpack_require__.d(__webpack_exports__, "default", function() { return matrix_Matrix; }); +/* concated harmony reexport Matrix */__webpack_require__.d(__webpack_exports__, "Matrix", function() { return matrix_Matrix; }); +/* concated harmony reexport abstractMatrix */__webpack_require__.d(__webpack_exports__, "abstractMatrix", function() { return AbstractMatrix; }); +/* concated harmony reexport wrap */__webpack_require__.d(__webpack_exports__, "wrap", function() { return wrap; }); +/* concated harmony reexport WrapperMatrix2D */__webpack_require__.d(__webpack_exports__, "WrapperMatrix2D", function() { return WrapperMatrix2D_WrapperMatrix2D; }); +/* concated harmony reexport WrapperMatrix1D */__webpack_require__.d(__webpack_exports__, "WrapperMatrix1D", function() { return WrapperMatrix1D_WrapperMatrix1D; }); +/* concated harmony reexport solve */__webpack_require__.d(__webpack_exports__, "solve", function() { return solve; }); +/* concated harmony reexport inverse */__webpack_require__.d(__webpack_exports__, "inverse", function() { return inverse; }); +/* concated harmony reexport linearDependencies */__webpack_require__.d(__webpack_exports__, "linearDependencies", function() { return linearDependencies; }); +/* concated harmony reexport SingularValueDecomposition */__webpack_require__.d(__webpack_exports__, "SingularValueDecomposition", function() { return svd_SingularValueDecomposition; }); +/* concated harmony reexport SVD */__webpack_require__.d(__webpack_exports__, "SVD", function() { return svd_SingularValueDecomposition; }); +/* concated harmony reexport EigenvalueDecomposition */__webpack_require__.d(__webpack_exports__, "EigenvalueDecomposition", function() { return evd_EigenvalueDecomposition; }); +/* concated harmony reexport EVD */__webpack_require__.d(__webpack_exports__, "EVD", function() { return evd_EigenvalueDecomposition; }); +/* concated harmony reexport CholeskyDecomposition */__webpack_require__.d(__webpack_exports__, "CholeskyDecomposition", function() { return cholesky_CholeskyDecomposition; }); +/* concated harmony reexport CHO */__webpack_require__.d(__webpack_exports__, "CHO", function() { return cholesky_CholeskyDecomposition; }); +/* concated harmony reexport LuDecomposition */__webpack_require__.d(__webpack_exports__, "LuDecomposition", function() { return lu_LuDecomposition; }); +/* concated harmony reexport LU */__webpack_require__.d(__webpack_exports__, "LU", function() { return lu_LuDecomposition; }); +/* concated harmony reexport QrDecomposition */__webpack_require__.d(__webpack_exports__, "QrDecomposition", function() { return qr_QrDecomposition; }); +/* concated harmony reexport QR */__webpack_require__.d(__webpack_exports__, "QR", function() { return qr_QrDecomposition; }); + + + + + + + + + + + + + + + + /***/ }) /******/ ]); \ No newline at end of file diff --git a/lib/umap-js.min.js b/lib/umap-js.min.js index a03476a..6494325 100644 --- a/lib/umap-js.min.js +++ b/lib/umap-js.min.js @@ -1 +1 @@ -!function(t){var e={};function n(r){if(e[r])return e[r].exports;var i=e[r]={i:r,l:!1,exports:{}};return t[r].call(i.exports,i,i.exports,n),i.l=!0,i.exports}n.m=t,n.c=e,n.d=function(t,e,r){n.o(t,e)||Object.defineProperty(t,e,{enumerable:!0,get:r})},n.r=function(t){"undefined"!=typeof Symbol&&Symbol.toStringTag&&Object.defineProperty(t,Symbol.toStringTag,{value:"Module"}),Object.defineProperty(t,"__esModule",{value:!0})},n.t=function(t,e){if(1&e&&(t=n(t)),8&e)return t;if(4&e&&"object"==typeof t&&t&&t.__esModule)return t;var r=Object.create(null);if(n.r(r),Object.defineProperty(r,"default",{enumerable:!0,value:t}),2&e&&"string"!=typeof t)for(var i in t)n.d(r,i,function(e){return t[e]}.bind(null,i));return r},n.n=function(t){var e=t&&t.__esModule?function(){return t.default}:function(){return t};return n.d(e,"a",e),e},n.o=function(t,e){return Object.prototype.hasOwnProperty.call(t,e)},n.p="",n(n.s=1)}([function(t,e,n){"use strict";function r(t){return void 0===t&&(t=Math.random),t()}function i(t){for(var e=[],n=0;ne?t[n]:e;return e},e.max2d=function(t){for(var e=0,n=0;ne?t[n][r]:e;return e}},function(t,e,n){"use strict";Object.defineProperty(e,"__esModule",{value:!0});var r=n(2);window.UMAP=r.UMAP},function(t,e,n){"use strict";var r=this&&this.__awaiter||function(t,e,n,r){return new(n||(n=Promise))(function(i,o){function a(t){try{h(r.next(t))}catch(t){o(t)}}function s(t){try{h(r.throw(t))}catch(t){o(t)}}function h(t){t.done?i(t.value):new n(function(e){e(t.value)}).then(a,s)}h((r=r.apply(t,e||[])).next())})},i=this&&this.__generator||function(t,e){var n,r,i,o,a={label:0,sent:function(){if(1&i[0])throw i[1];return i[1]},trys:[],ops:[]};return o={next:s(0),throw:s(1),return:s(2)},"function"==typeof Symbol&&(o[Symbol.iterator]=function(){return this}),o;function s(o){return function(s){return function(o){if(n)throw new TypeError("Generator is already executing.");for(;a;)try{if(n=1,r&&(i=2&o[0]?r.return:o[0]?r.throw||((i=r.return)&&i.call(r),0):r.next)&&!(i=i.call(r,o[1])).done)return i;switch(r=0,i&&(o=[2&o[0],i.value]),o[0]){case 0:case 1:i=o;break;case 4:return a.label++,{value:o[1],done:!1};case 5:a.label++,r=o[1],o=[0];continue;case 7:o=a.ops.pop(),a.trys.pop();continue;default:if(!(i=(i=a.trys).length>0&&i[i.length-1])&&(6===o[0]||2===o[0])){a=0;continue}if(3===o[0]&&(!i||o[1]>i[0]&&o[1]0});if(d.length>=n){var v=Math.floor(n),m=n-v;v>0?(a[u]=d[v-1],m>1e-5&&(a[u]+=m*(d[v]-d[v-1]))):a[u]=m*d[0]}else d.length>0&&(a[u]=h.max(d));for(var g=0;g0?Math.exp(-w/f):1}if(Math.abs(y-o)<1e-5)break;y>o?f=(l+(c=f))/2:(l=f,c===1/0?f*=2:f=(l+c)/2)}if(s[u]=f,a[u]>0){var M=h.mean(p);s[u]<.001*M&&(s[u]=.001*M)}else{var S=h.mean(t.map(h.mean));s[u]<.001*S&&(s[u]=.001*S)}}return{sigmas:s,rhos:a}},t.prototype.computeMembershipStrengths=function(t,e,n,r){for(var i=t.length,o=t[0].length,a=h.zeros(i*o),s=h.zeros(i*o),u=h.zeros(i*o),l=0;l0&&(n[r]=e/i[r])}),n},t.prototype.initializeOptimization=function(){var t=this.embedding,e=this.embedding,n=this.optimizationState,r=n.head,i=n.tail,o=n.epochsPerSample,a=this.getNEpochs(),s=this.graph.nCols,h=t[0].length,u=t.length===e.length,l=o.map(function(t){return t/5}),c=l.slice(),f=o.slice();Object.assign(this.optimizationState,{isInitialized:!0,headEmbedding:t,tailEmbedding:e,head:r,tail:i,epochsPerSample:o,epochOfNextSample:f,epochOfNextNegativeSample:c,epochsPerNegativeSample:l,moveOther:u,initialAlpha:1,alpha:1,gamma:1,a:1.5769434603113077,b:.8950608779109733,dim:h,nEpochs:a,nVertices:s})},t.prototype.optimizeLayoutStep=function(t){for(var e=this.optimizationState,n=e.head,r=e.tail,i=e.headEmbedding,o=e.tailEmbedding,a=e.epochsPerSample,s=e.epochOfNextSample,u=e.epochOfNextNegativeSample,l=e.epochsPerNegativeSample,c=e.moveOther,d=e.initialAlpha,v=e.alpha,m=e.gamma,g=e.a,y=e.b,b=e.dim,w=e.nEpochs,M=e.nVertices,S=0;St)){var z=n[S],N=r[S],k=i[z],E=o[N],P=p(k,E),x=0;P>0&&(x=-2*g*y*Math.pow(P,y-1),x/=g*Math.pow(P,y)+1);for(var C=0;C0)A=2*m*y,A/=(.001+L)*(g*Math.pow(L,y)+1);else if(z===I)continue;for(C=0;C0&&(O=f(A*(k[C]-j[C]),4)),k[C]+=O*v}}u[S]+=R*l[S]}return e.alpha=d*(1-t/w),e.currentEpoch+=1,this.embedding=i,e.currentEpoch},t.prototype.optimizeLayout=function(t){var e=this;return void 0===t&&(t=function(){return!0}),this.optimizationState.isInitialized||this.initializeOptimization(),new Promise(function(n,o){var a=function(){return r(e,void 0,void 0,function(){var e,r,s,h,u,l;return i(this,function(i){try{if(e=this.optimizationState,r=e.nEpochs,s=e.currentEpoch,h=this.optimizeLayoutStep(s),u=!1===t(h),l=h===r,u||l)return[2,n(l)];a()}catch(t){o(t)}return[2]})})};a()})},t.prototype.getNEpochs=function(){var t=this.graph;if(this.nEpochs>0)return this.nEpochs;var e=t.nRows;return e<=2500?500:e<=5e3?400:e<=7500?300:200},t}();function l(t,e){for(var n=0,r=0;re?e:t<-e?-e:t}function p(t,e){for(var n=0,r=0;r=a[0])return 0;for(var h=0;h=p)break;if(f>=p){if(!(a[c]>n))break;l=c}else if(a[c]>=a[f]){if(!(ni){var s=function(t,e,n){var i=t[0].length,o=r.tauRandInt(e.length,n),a=r.tauRandInt(e.length,n);a=(a+=o===a?1:0)%e.length;for(var s=e[o],h=e[a],u=0,l=r.zeros(i),c=0;c0?(d[c]=0,f+=1):(d[c]=1,p+=1)}var g=r.zeros(f),y=r.zeros(p);for(var c in f=0,p=0,r.range(d.length))0===d[c]?(g[f]=e[c],f+=1):(y[p]=e[c],p+=1);return{indicesLeft:g,indicesRight:y,hyperplane:l,offset:u}}(e,n,a),h=s.indicesLeft,u=s.indicesRight,l=s.hyperplane,c=s.offset,f=t(e,h,i,o+1,a),p=t(e,u,i,o+1,a),d={leftChild:f,rightChild:p,isLeaf:!1,hyperplane:l,offset:c};return d}var d={indices:n,isLeaf:!0};return d}(t,o,e,n,i)}(t,a,n,o)}).map(function(t){return function(t,e){var n=function t(e){return e.isLeaf?1:1+t(e.leftChild)+t(e.rightChild)}(t),o=function t(e){return e.isLeaf?1:t(e.leftChild)+t(e.rightChild)}(t),a=r.range(n).map(function(){return r.zeros(t.hyperplane.length)}),s=r.zeros(n),h=r.range(n).map(function(){return[-1,-1]}),u=r.range(o).map(function(){return r.range(e).map(function(){return-1})});return function t(e,n,r,i,o,a,s){var h;if(e.isLeaf)return i[a][0]=-s,(h=o[s]).splice.apply(h,[0,e.indices.length].concat(e.indices)),{nodeNum:a,leafNum:s+=1};n[a]=e.hyperplane,r[a]=e.offset,i[a][0]=a+1;var u=a,l=t(e.leftChild,n,r,i,o,a+1,s);return a=l.nodeNum,s=l.leafNum,i[u][1]=a+1,{nodeNum:(l=t(e.rightChild,n,r,i,o,a+1,s)).nodeNum,leafNum:l.leafNum}}(t,a,s,h,u,0,0),new i(a,s,h,u)}(t,a)})},e.makeLeafArray=function(t){if(t.length>0){for(var e=[],n=0,r=t;n=t.length&&(t=void 0),{value:t&&t[e++],done:!t}}}};function i(t){return void 0===t&&(t=Math.random),t()}function o(t){for(var r=[],e=0;er?t[e]:r;return r},r.max2d=function(t){for(var r=0,e=0;er?t[e][n]:r;return r}},function(t,r,e){"use strict";Object.defineProperty(r,"__esModule",{value:!0});var n=e(3);window.UMAP=n.UMAP},function(t,r,e){"use strict";var n=this&&this.__awaiter||function(t,r,e,n){return new(e||(e=Promise))(function(i,o){function s(t){try{h(n.next(t))}catch(t){o(t)}}function a(t){try{h(n.throw(t))}catch(t){o(t)}}function h(t){t.done?i(t.value):new e(function(r){r(t.value)}).then(s,a)}h((n=n.apply(t,r||[])).next())})},i=this&&this.__generator||function(t,r){var e,n,i,o,s={label:0,sent:function(){if(1&i[0])throw i[1];return i[1]},trys:[],ops:[]};return o={next:a(0),throw:a(1),return:a(2)},"function"==typeof Symbol&&(o[Symbol.iterator]=function(){return this}),o;function a(o){return function(a){return function(o){if(e)throw new TypeError("Generator is already executing.");for(;s;)try{if(e=1,n&&(i=2&o[0]?n.return:o[0]?n.throw||((i=n.return)&&i.call(n),0):n.next)&&!(i=i.call(n,o[1])).done)return i;switch(n=0,i&&(o=[2&o[0],i.value]),o[0]){case 0:case 1:i=o;break;case 4:return s.label++,{value:o[1],done:!1};case 5:s.label++,n=o[1],o=[0];continue;case 7:o=s.ops.pop(),s.trys.pop();continue;default:if(!(i=(i=s.trys).length>0&&i[i.length-1])&&(6===o[0]||2===o[0])){s=0;continue}if(3===o[0]&&(!i||o[1]>i[0]&&o[1]0)&&!(n=o.next()).done;)s.push(n.value)}catch(t){i={error:t}}finally{try{n&&!n.done&&(e=o.return)&&e.call(o)}finally{if(i)throw i.error}}return s},s=this&&this.__spread||function(){for(var t=[],r=0;r0});if(g.length>=e){var v=Math.floor(e),p=e-v;v>0?(s[h]=g[v-1],p>1e-5&&(s[h]+=p*(g[v]-g[v-1]))):s[h]=p*g[0]}else g.length>0&&(s[h]=l.max(g));for(var w=0;w0?Math.exp(-b/c):1}if(Math.abs(d-o)<1e-5)break;d>o?c=(u+(f=c))/2:(u=c,f===1/0?c*=2:c=(u+f)/2)}if(a[h]=c,s[h]>0){var M=l.mean(m);a[h]<.001*M&&(a[h]=.001*M)}else{var x=l.mean(t.map(l.mean));a[h]<.001*x&&(a[h]=.001*x)}}return{sigmas:a,rhos:s}},t.prototype.computeMembershipStrengths=function(t,r,e,n){for(var i=t.length,o=t[0].length,s=l.zeros(i*o),a=l.zeros(i*o),h=l.zeros(i*o),u=0;u0&&(e[n]=r/i[n])}),e},t.prototype.initializeOptimization=function(){var t=this.embedding,r=this.embedding,e=this.optimizationState,n=e.head,i=e.tail,o=e.epochsPerSample,a=this.getNEpochs(),h=this.graph.nCols,u=w(this.spread,this.minDist),l=u.a,f=u.b,c=t[0].length,m=t.length===r.length,g=o.map(function(t){return t/5}),v=s(g),p=s(o);Object.assign(this.optimizationState,{isInitialized:!0,headEmbedding:t,tailEmbedding:r,head:n,tail:i,epochsPerSample:o,epochOfNextSample:p,epochOfNextNegativeSample:v,epochsPerNegativeSample:g,moveOther:m,initialAlpha:1,alpha:1,gamma:1,a:l,b:f,dim:c,nEpochs:a,nVertices:h})},t.prototype.optimizeLayoutStep=function(t){for(var r=this.optimizationState,e=r.head,n=r.tail,i=r.headEmbedding,o=r.tailEmbedding,s=r.epochsPerSample,a=r.epochOfNextSample,h=r.epochOfNextNegativeSample,u=r.epochsPerNegativeSample,f=r.moveOther,c=r.initialAlpha,m=r.alpha,g=r.gamma,w=r.a,d=r.b,y=r.dim,b=r.nEpochs,M=r.nVertices,x=0;xt)){var S=e[x],E=n[x],R=i[S],k=o[E],A=p(R,k),N=0;A>0&&(N=-2*w*d*Math.pow(A,d-1),N/=w*Math.pow(A,d)+1);for(var V=0;V0)_=2*g*d,_/=(.001+D)*(w*Math.pow(D,d)+1);else if(S===j)continue;for(V=0;V0&&(z=v(_*(R[V]-P[V]),4)),R[V]+=z*m}}h[x]+=C*u[x]}return r.alpha=c*(1-t/b),r.currentEpoch+=1,this.embedding=i,r.currentEpoch},t.prototype.optimizeLayout=function(t){var r=this;return void 0===t&&(t=function(){return!0}),this.optimizationState.isInitialized||this.initializeOptimization(),new Promise(function(e,o){var s=function(){return n(r,void 0,void 0,function(){var r,n,a,h,u,l;return i(this,function(i){try{if(r=this.optimizationState,n=r.nEpochs,a=r.currentEpoch,h=this.optimizeLayoutStep(a),u=!1===t(h),l=h===n,u||l)return[2,e(l)];s()}catch(t){o(t)}return[2]})})};s()})},t.prototype.getNEpochs=function(){var t=this.graph;if(this.nEpochs>0)return this.nEpochs;var r=t.nRows;return r<=2500?500:r<=5e3?400:r<=7500?300:200},t}();function m(t,r){for(var e=0,n=0;nr?r:t<-r?-r:t}function p(t,r){for(var e=0,n=0;n=r?Math.exp(-(e[i]-r)/t):n}),i=f({x:e,y:n},function(t){var r=o(t,2),e=r[0],n=r[1];return function(t){return 1/(1+e*Math.pow(t,2*n))}},{damping:1.5,initialValues:[.5,.5],gradientDifference:.1,maxIterations:100,errorTolerance:.01}).parameterValues,s=o(i,2);return{a:s[0],b:s[1]}}function d(t,r,e,n){return void 0===e&&(e=1),void 0===n&&(n=5),t.map(function(t,i,o){return-1===r[i]||-1===r[o]?t*Math.exp(-e):r[i]!==r[o]?t*Math.exp(-n):t})}function y(t){t=a.normalize(t,"max");var r=a.transpose(t),e=a.pairwiseMultiply(r,t);return t=a.add(t,a.subtract(r,e)),a.eliminateZeros(t)}r.findABParams=w,r.fastIntersection=d,r.resetLocalConnectivity=y},function(t,r,e){"use strict";var n,i=this&&this.__read||function(t,r){var e="function"==typeof Symbol&&t[Symbol.iterator];if(!e)return t;var n,i,o=e.call(t),s=[];try{for(;(void 0===r||r-- >0)&&!(n=o.next()).done;)s.push(n.value)}catch(t){i={error:t}}finally{try{n&&!n.done&&(e=o.return)&&e.call(o)}finally{if(i)throw i.error}}return s},o=this&&this.__spread||function(){for(var t=[],r=0;r=t.length&&(t=void 0),{value:t&&t[e++],done:!t}}}};Object.defineProperty(r,"__esModule",{value:!0});var a=e(1),h=function(){function t(t,r,e,n){this.entries=new Map,this.nRows=0,this.nCols=0,this.rows=o(t),this.cols=o(r),this.values=o(e);for(var i=0;ir?t[e]:r;return t.map(function(t){return t/r})},n.l1=function(t){for(var r=0,e=0;e=s[0])return 0;for(var h=0;h=m)break;if(c>=m){if(!(s[f]>e))break;l=f}else if(s[f]>=s[c]){if(!(e0)&&!(n=o.next()).done;)s.push(n.value)}catch(t){i={error:t}}finally{try{n&&!n.done&&(e=o.return)&&e.call(o)}finally{if(i)throw i.error}}return s},i=this&&this.__spread||function(){for(var t=[],r=0;r=t.length&&(t=void 0),{value:t&&t[e++],done:!t}}}};Object.defineProperty(r,"__esModule",{value:!0});var s=e(1),a=function(){return function(t,r,e,n){this.hyperplanes=t,this.offsets=r,this.children=e,this.indices=n}}();r.FlatTree=a,r.makeForest=function(t,r,e,n){var o=Math.max(10,r);return s.range(e).map(function(r,e){return function(t,r,e,n){void 0===r&&(r=30);var i=s.range(t.length);return function t(r,e,n,i,o){if(void 0===n&&(n=30),e.length>n){var a=function(t,r,e){var n=t[0].length,i=s.tauRandInt(r.length,e),o=s.tauRandInt(r.length,e);o=(o+=i===o?1:0)%r.length;for(var a=r[i],h=r[o],u=0,l=s.zeros(n),f=0;f0?(g[f]=0,c+=1):(g[f]=1,m+=1)}var w=s.zeros(c),d=s.zeros(m);for(var f in c=0,m=0,s.range(g.length))0===g[f]?(w[c]=r[f],c+=1):(d[m]=r[f],m+=1);return{indicesLeft:w,indicesRight:d,hyperplane:l,offset:u}}(r,e,o),h=a.indicesLeft,u=a.indicesRight,l=a.hyperplane,f=a.offset,c=t(r,h,n,i+1,o),m=t(r,u,n,i+1,o),g={leftChild:c,rightChild:m,isLeaf:!1,hyperplane:l,offset:f};return g}var g={indices:e,isLeaf:!0};return g}(t,i,r,e,n)}(t,o,e,n)}).map(function(t){return function(t,r){var e=function t(r){return r.isLeaf?1:1+t(r.leftChild)+t(r.rightChild)}(t),n=function t(r){return r.isLeaf?1:t(r.leftChild)+t(r.rightChild)}(t),o=s.range(e).map(function(){return s.zeros(t.hyperplane.length)}),h=s.zeros(e),u=s.range(e).map(function(){return[-1,-1]}),l=s.range(n).map(function(){return s.range(r).map(function(){return-1})});return function t(r,e,n,o,s,a,h){var u;if(r.isLeaf)return o[a][0]=-h,(u=s[h]).splice.apply(u,i([0,r.indices.length],r.indices)),{nodeNum:a,leafNum:h+=1};e[a]=r.hyperplane,n[a]=r.offset,o[a][0]=a+1;var l=a,f=t(r.leftChild,e,n,o,s,a+1,h);return a=f.nodeNum,h=f.leafNum,o[l][1]=a+1,{nodeNum:(f=t(r.rightChild,e,n,o,s,a+1,h)).nodeNum,leafNum:f.leafNum}}(t,o,h,u,l,0,0),new a(o,h,u,l)}(t,o)})},r.makeLeafArray=function(t){var r,e;if(t.length>0){var n=[];try{for(var s=o(t),a=s.next();!a.done;a=s.next()){var h=a.value;n.push.apply(n,i(h.indices))}}catch(t){r={error:t}}finally{try{a&&!a.done&&(e=s.return)&&e.call(s)}finally{if(r)throw r.error}}return n}return[[-1]]}},function(t,r,e){"use strict";var n=e(9);function i(t,r,e){var n=0;const i=e(r);for(var o=0;or&&(r=t[e]);return r};var s=function(t){if(!i()(t))throw new TypeError("input must be an array");if(0===t.length)throw new TypeError("input must not be empty");for(var r=t[0],e=1;e1&&void 0!==arguments[1]?arguments[1]:{};if(!i()(t))throw new TypeError("input must be an array");if(0===t.length)throw new TypeError("input must not be empty");if(void 0!==e.output){if(!i()(e.output))throw new TypeError("output option must be an array if specified");r=e.output}else r=new Array(t.length);var n=s(t),a=o(t);if(n===a)throw new RangeError("minimum and maximum input values are equal. Cannot rescale a constant array");var h=e.min,u=void 0===h?e.autoMinMax?n:0:h,l=e.max,f=void 0===l?e.autoMinMax?a:1:l;if(u>=f)throw new RangeError("min option must be smaller than max option");for(var c=(f-u)/(a-n),m=0;mMath.abs(h[i])&&(i=r);if(i!==e){for(n=0;n=0;i--){for(n=0;no?t[i][o]:i===o?1:0;return n}get upperTriangularMatrix(){for(var t=this.LU,r=t.rows,e=t.columns,n=new I(r,e),i=0;iMath.abs(r)?(e=r/t,Math.abs(t)*Math.sqrt(1+e*e)):0!==r?(e=t/r,Math.abs(r)*Math.sqrt(1+e*e)):0}function l(t,r,e){for(var n=new Array(t),i=0;i=0;t--)if(0!==p[t]){for(let r=t+1;r=0;t--){if(t0;){let t,r;for(t=R-2;t>=-1&&-1!==t;t--){const r=Number.MIN_VALUE+A*Math.abs(p[t]+Math.abs(p[t+1]));if(Math.abs(y[t])<=r||Number.isNaN(y[t])){y[t]=0;break}}if(t===R-2)r=4;else{let e;for(e=R-1;e>=t&&e!==t;e--){let r=(e!==R?Math.abs(y[e]):0)+(e!==t+1?Math.abs(y[e-1]):0);if(Math.abs(p[e])<=A*r){p[e]=0;break}}e===t?r=3:e===R-1?r=1:(r=2,t=e)}switch(t++,r){case 1:{let r=y[R-2];y[R-2]=0;for(let e=R-2;e>=t;e--){let i=u(p[e],r),o=p[e]/i,s=r/i;if(p[e]=i,e!==t&&(r=-s*y[e-1],y[e-1]=o*y[e-1]),f)for(let t=0;t=p[t+1]);){let r=p[t];if(p[t]=p[t+1],p[t+1]=r,f&&tr?i[o][e]=t[o][e]/this.s[e]:i[o][e]=0;var o=this.U,s=o.length,a=o[0].length,h=new I(e,s);for(let t=0;tt&&r++;return r}get diagonal(){return this.s}get threshold(){return Number.EPSILON/2*Math.max(this.m,this.n)*this.s[0]}get leftSingularVectors(){return I.isMatrix(this.U)||(this.U=new I(this.U)),this.U}get rightSingularVectors(){return I.isMatrix(this.V)||(this.V=new I(this.V)),this.V}get diagonalMatrix(){return I.diag(this.s)}}function c(t,r,e){var n=e?t.rows:t.rows-1;if(r<0||r>n)throw new RangeError("Row index out of range")}function m(t,r,e){var n=e?t.columns:t.columns-1;if(r<0||r>n)throw new RangeError("Column index out of range")}function g(t,r){if(r.to1DArray&&(r=r.to1DArray()),r.length!==t.columns)throw new RangeError("vector size must be the same as the number of columns");return r}function v(t,r){if(r.to1DArray&&(r=r.to1DArray()),r.length!==t.rows)throw new RangeError("vector size must be the same as the number of rows");return r}function p(t,r,e){return{row:w(t,r),column:d(t,e)}}function w(t,r){if("object"!=typeof r)throw new TypeError("unexpected type for row indices");if(r.some(r=>r<0||r>=t.rows))throw new RangeError("row indices are out of range");return Array.isArray(r)||(r=Array.from(r)),r}function d(t,r){if("object"!=typeof r)throw new TypeError("unexpected type for column indices");if(r.some(r=>r<0||r>=t.columns))throw new RangeError("column indices are out of range");return Array.isArray(r)||(r=Array.from(r)),r}function y(t,r,e,n,i){if(5!==arguments.length)throw new RangeError("expected 4 arguments");if(b("startRow",r),b("endRow",e),b("startColumn",n),b("endColumn",i),r>e||n>i||r<0||r>=t.rows||e<0||e>=t.rows||n<0||n>=t.columns||i<0||i>=t.columns)throw new RangeError("Submatrix indices are out of range")}function b(t,r){if("number"!=typeof r)throw new TypeError(`${t} must be a number`)}class M extends(C()){constructor(t,r,e){super(),this.matrix=t,this.rows=r,this.columns=e}static get[Symbol.species](){return I}}class x extends M{constructor(t){super(t,t.columns,t.rows)}set(t,r,e){return this.matrix.set(r,t,e),this}get(t,r){return this.matrix.get(r,t)}}class S extends M{constructor(t,r){super(t,1,t.columns),this.row=r}set(t,r,e){return this.matrix.set(this.row,r,e),this}get(t,r){return this.matrix.get(this.row,r)}}class E extends M{constructor(t,r,e,n,i){y(t,r,e,n,i),super(t,e-r+1,i-n+1),this.startRow=r,this.startColumn=n}set(t,r,e){return this.matrix.set(this.startRow+t,this.startColumn+r,e),this}get(t,r){return this.matrix.get(this.startRow+t,this.startColumn+r)}}class R extends M{constructor(t,r,e){var n=p(t,r,e);super(t,n.row.length,n.column.length),this.rowIndices=n.row,this.columnIndices=n.column}set(t,r,e){return this.matrix.set(this.rowIndices[t],this.columnIndices[r],e),this}get(t,r){return this.matrix.get(this.rowIndices[t],this.columnIndices[r])}}class k extends M{constructor(t,r){super(t,(r=w(t,r)).length,t.columns),this.rowIndices=r}set(t,r,e){return this.matrix.set(this.rowIndices[t],r,e),this}get(t,r){return this.matrix.get(this.rowIndices[t],r)}}class A extends M{constructor(t,r){r=d(t,r),super(t,t.rows,r.length),this.columnIndices=r}set(t,r,e){return this.matrix.set(t,this.columnIndices[r],e),this}get(t,r){return this.matrix.get(t,this.columnIndices[r])}}class N extends M{constructor(t,r){super(t,t.rows,1),this.column=r}set(t,r,e){return this.matrix.set(t,this.column,e),this}get(t){return this.matrix.get(t,this.column)}}class V extends M{constructor(t){super(t,t.rows,t.columns)}set(t,r,e){return this.matrix.set(this.rows-t-1,r,e),this}get(t,r){return this.matrix.get(this.rows-t-1,r)}}class z extends M{constructor(t){super(t,t.rows,t.columns)}set(t,r,e){return this.matrix.set(t,this.columns-r-1,e),this}get(t,r){return this.matrix.get(t,this.columns-r-1)}}function C(t){void 0===t&&(t=Object);class r extends t{static get[Symbol.species](){return this}static from1DArray(t,r,e){if(t*r!==e.length)throw new RangeError("Data length does not match given dimensions");for(var n=new this(t,r),i=0;it&&(t=this.get(r,e));return t}maxIndex(){for(var t=this.get(0,0),r=[0,0],e=0;et&&(t=this.get(e,n),r[0]=e,r[1]=n);return r}min(){for(var t=this.get(0,0),r=0;rr&&(r=this.get(t,e));return r}maxRowIndex(t){c(this,t);for(var r=this.get(t,0),e=[t,0],n=1;nr&&(r=this.get(t,n),e[1]=n);return e}minRow(t){c(this,t);for(var r=this.get(t,0),e=1;er&&(r=this.get(e,t));return r}maxColumnIndex(t){m(this,t);for(var r=this.get(0,t),e=[0,t],n=1;nr&&(r=this.get(n,t),e[0]=n);return e}minColumn(t){m(this,t);for(var r=this.get(0,t),e=1;e=(r=void 0===r?1:r))throw new RangeError("min should be strictly smaller than max");for(var e=this.constructor.empty(this.rows,this.columns),n=0;n=(r=void 0===r?1:r))throw new RangeError("min should be strictly smaller than max");for(var e=this.constructor.empty(this.rows,this.columns),n=0;ne||r<0||r>=this.columns||e<0||e>=this.columns)throw new RangeError("Argument out of range");for(var n=new this.constructor[Symbol.species](t.length,e-r+1),i=0;i=this.rows)throw new RangeError(`Row index out of range: ${t[i]}`);n.set(i,o-r,this.get(t[i],o))}return n}subMatrixColumn(t,r,e){if(void 0===r&&(r=0),void 0===e&&(e=this.rows-1),r>e||r<0||r>=this.rows||e<0||e>=this.rows)throw new RangeError("Argument out of range");for(var n=new this.constructor[Symbol.species](e-r+1,t.length),i=0;i=this.columns)throw new RangeError(`Column index out of range: ${t[i]}`);n.set(o-r,i,this.get(o,t[i]))}return n}setSubMatrix(t,r,e){y(this,r,r+(t=this.constructor.checkMatrix(t)).rows-1,e,e+t.columns-1);for(var n=0;nt?i[o]=1/i[o]:i[o]=0;return i=this.constructor[Symbol.species].diag(i),n.mmul(i.mmul(e.transposeView()))}clone(){for(var t=new this.constructor[Symbol.species](this.rows,this.columns),r=0;r>","signPropagatingRightShift"],[">>>","rightShift","zeroFillRightShift"]]){var u=o(F("\n(function %name%(value) {\n if (typeof value === 'number') return this.%name%S(value);\n return this.%name%M(value);\n})\n",{name:s[1],op:s[0]})),l=o(F("\n(function %name%S(value) {\n for (var i = 0; i < this.rows; i++) {\n for (var j = 0; j < this.columns; j++) {\n this.set(i, j, this.get(i, j) %op% value);\n }\n }\n return this;\n})\n",{name:`${s[1]}S`,op:s[0]})),w=o(F("\n(function %name%M(matrix) {\n matrix = this.constructor.checkMatrix(matrix);\n if (this.rows !== matrix.rows ||\n this.columns !== matrix.columns) {\n throw new RangeError('Matrices dimensions must be equal');\n }\n for (var i = 0; i < this.rows; i++) {\n for (var j = 0; j < this.columns; j++) {\n this.set(i, j, this.get(i, j) %op% matrix.get(i, j));\n }\n }\n return this;\n})\n",{name:`${s[1]}M`,op:s[0]})),d=o(F("\n(function %name%(matrix, value) {\n var newMatrix = new this[Symbol.species](matrix);\n return newMatrix.%name%(value);\n})\n",{name:s[1]}));for(n=1;n0){if(super(t),!(Number.isInteger(r)&&r>0))throw new TypeError("nColumns must be a positive integer");for(e=0;e=0;o--){for(i=0;i=0;e--){for(t=0;ti)return new Array(r.rows+1).fill(0);for(var o=r.addRow(e,[0]),s=0;s0;a--){for(f=0,s=0,u=0;u0&&(o=-o),r[a]=f*o,s-=i*o,e[a-1]=i-o,h=0;hl){0;do{for(1,i=e[l],c=(e[l+1]-i)/(2*r[l]),m=u(c,1),c<0&&(m=-m),e[l]=r[l]/(c+m),e[l+1]=r[l]*(c+m),g=e[l+1],o=i-e[l],s=l+2;s=l;s--)for(w=p,p=v,b=y,i=v*r[s],o=v*c,m=u(c,r[s]),r[s+1]=y*m,y=r[s]/m,c=(v=c/m)*e[s]-y*i,e[s+1]=o+y*(v*i+y*e[s]),h=0;hS*x)}e[l]=e[l]+M,r[l]=0}for(s=0;s=u;a--)e[a]=r[a][u-1]/l,s+=e[a]*e[a];for(o=Math.sqrt(s),e[u]>0&&(o=-o),s-=e[u]*o,e[u]=e[u]-o,h=u;h=u;a--)i+=e[a]*r[a][h];for(i/=s,a=u;a<=f;a++)r[a][h]-=i*e[a]}for(a=0;a<=f;a++){for(i=0,h=f;h>=u;h--)i+=e[h]*r[a][h];for(i/=s,h=u;h<=f;h++)r[a][h]-=i*e[h]}e[u]=l*e[u],r[u][u-1]=l*o}}for(a=0;a=1;u--)if(0!==r[u][u-1]){for(a=u+1;a<=f;a++)e[a]=r[a][u-1];for(h=u;h<=f;h++){for(o=0,a=u;a<=f;a++)o+=e[a]*n[a][h];for(o=o/e[u]/r[u][u-1],a=u;a<=f;a++)n[a][h]+=o*e[a]}}}(o,c,m,s),function(t,r,e,n,i){var o,s,a,h,u,l,f,c,m,g,v,p,w,d,y,b=t-1,M=t-1,x=Number.EPSILON,S=0,E=0,R=0,k=0,A=0,N=0,V=0,z=0;for(o=0;oM)&&(e[o]=i[o][o],r[o]=0),s=Math.max(o-1,0);s=0;){for(h=b;h>0&&(0===(N=Math.abs(i[h-1][h-1])+Math.abs(i[h][h]))&&(N=E),!(Math.abs(i[h][h-1])=0){for(V=R>=0?R+V:R-V,e[b-1]=c+V,e[b]=e[b-1],0!==V&&(e[b]=c-f/V),r[b-1]=0,r[b]=0,c=i[b][b-1],N=Math.abs(c)+Math.abs(V),R=c/N,k=V/N,A=Math.sqrt(R*R+k*k),R/=A,k/=A,s=b-1;s0){for(N=Math.sqrt(N),m=h&&(V=i[u][u],R=((A=c-V)*(N=m-V)-f)/i[u+1][u]+i[u][u+1],k=i[u+1][u+1]-V-A-N,A=i[u+2][u+1],N=Math.abs(R)+Math.abs(k)+Math.abs(A),R/=N,k/=N,A/=N,u!==h)&&!(Math.abs(i[u][u-1])*(Math.abs(k)+Math.abs(A))u+2&&(i[o][o-3]=0);for(a=u;a<=b-1&&(d=a!==b-1,a!==u&&(R=i[a][a-1],k=i[a+1][a-1],A=d?i[a+2][a-1]:0,0!==(c=Math.abs(R)+Math.abs(k)+Math.abs(A))&&(R/=c,k/=c,A/=c)),0!==c);a++)if(N=Math.sqrt(R*R+k*k+A*A),R<0&&(N=-N),0!==N){for(a!==u?i[a][a-1]=-N*c:h!==u&&(i[a][a-1]=-i[a][a-1]),c=(R+=N)/N,m=k/N,V=A/N,k/=R,A/=R,s=a;s=0;b--)if(R=e[b],0===(k=r[b]))for(h=b,i[b][b]=1,o=b-1;o>=0;o--){for(f=i[o][o]-R,A=0,s=h;s<=b;s++)A+=i[o][s]*i[s][b];if(r[o]<0)V=f,N=A;else if(h=o,0===r[o]?i[o][b]=0!==f?-A/f:-A/(x*E):(c=i[o][o+1],m=i[o+1][o],k=(e[o]-R)*(e[o]-R)+r[o]*r[o],l=(c*N-V*A)/k,i[o][b]=l,i[o+1][b]=Math.abs(c)>Math.abs(V)?(-A-f*l)/c:(-N-m*l)/V),l=Math.abs(i[o][b]),x*l*l>1)for(s=o;s<=b;s++)i[s][b]=i[s][b]/l}else if(k<0)for(h=b-1,Math.abs(i[b][b-1])>Math.abs(i[b-1][b])?(i[b-1][b-1]=k/i[b][b-1],i[b-1][b]=-(i[b][b]-R)/i[b][b-1]):(y=F(0,-i[b-1][b],i[b-1][b-1]-R,k),i[b-1][b-1]=y[0],i[b-1][b]=y[1]),i[b][b-1]=0,i[b][b]=1,o=b-2;o>=0;o--){for(g=0,v=0,s=h;s<=b;s++)g+=i[o][s]*i[s][b-1],v+=i[o][s]*i[s][b];if(f=i[o][o]-R,r[o]<0)V=f,A=g,N=v;else if(h=o,0===r[o]?(y=F(-g,-v,f,k),i[o][b-1]=y[0],i[o][b]=y[1]):(c=i[o][o+1],m=i[o+1][o],p=(e[o]-R)*(e[o]-R)+r[o]*r[o]-k*k,w=2*(e[o]-R)*k,0===p&&0===w&&(p=x*E*(Math.abs(f)+Math.abs(k)+Math.abs(c)+Math.abs(m)+Math.abs(V))),y=F(c*A-V*g+k*v,c*N-V*v-k*g,p,w),i[o][b-1]=y[0],i[o][b]=y[1],Math.abs(c)>Math.abs(V)+Math.abs(k)?(i[o+1][b-1]=(-g-f*i[o][b-1]+k*i[o][b])/c,i[o+1][b]=(-v-f*i[o][b]-k*i[o][b-1])/c):(y=F(-A-m*i[o][b-1],-N-m*i[o][b],V,k),i[o+1][b-1]=y[0],i[o+1][b]=y[1])),l=Math.max(Math.abs(i[o][b-1]),Math.abs(i[o][b])),x*l*l>1)for(s=o;s<=b;s++)i[s][b-1]=i[s][b-1]/l,i[s][b]=i[s][b]/l}for(o=0;oM)for(s=o;s=0;s--)for(o=0;o<=M;o++){for(V=0,a=0;a<=Math.min(s,M);a++)V+=n[o][a]*i[a][s];n[o][s]=V}}(o,h,a,s,c)}this.n=o,this.e=h,this.d=a,this.V=s}get realEigenvalues(){return this.d}get imaginaryEigenvalues(){return this.e}get eigenvectorMatrix(){return I.isMatrix(this.V)||(this.V=new I(this.V)),this.V}get diagonalMatrix(){var t,r,e=this.n,n=this.e,i=this.d,o=new I(e,e);for(t=0;t0?o[t][t+1]=n[t]:n[t]<0&&(o[t][t-1]=n[t])}return o}}function F(t,r,e,n){var i,o;return Math.abs(e)>Math.abs(n)?[(t+(i=n/e)*r)/(o=e+i*n),(r-i*t)/o]:[((i=e/n)*t+r)/(o=n+i*e),(i*r-t)/o]}class W{constructor(t){if(!(t=P.checkMatrix(t)).isSymmetric())throw new Error("Matrix is not symmetric");var r,e,n,i=t,o=i.rows,s=new I(o,o),a=!0;for(e=0;e0,s[e][e]=Math.sqrt(Math.max(u,0)),n=e+1;n=0;o--)for(i=0;i