Line data Source code
1 1 : import * as utils from "../utils.ts";
2 : /*
3 : Copyright (c) 2016 Chi Feng
4 :
5 : Permission is hereby granted, free of charge, to any person obtaining a copy
6 : of this software and associated documentation files (the "Software"), to deal
7 : in the Software without restriction, including without limitation the rights
8 : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9 : copies of the Software, and to permit persons to whom the Software is
10 : furnished to do so, subject to the following conditions:
11 :
12 : The above copyright notice and this permission notice shall be included in
13 : all copies or substantial portions of the Software.
14 :
15 : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18 : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20 : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21 : THE SOFTWARE.
22 : */
23 :
24 : /**
25 : * INFO: port of the _Bayes for physicists_ linear algebra library to Deno
26 : */
27 1 : export let precision = 1e-6;
28 :
29 :
30 : export interface LUResult {
31 : L: Matrix;
32 : U: Matrix;
33 : P: Matrix;
34 : }
35 :
36 : export interface SVDResult {
37 : U: Matrix;
38 : S: Vector;
39 : V: Matrix;
40 : }
41 :
42 1 : export class Vector {
43 : dim: number;
44 : data: Float64Array;
45 0 : constructor(array: Array<number> | Float64Array) {
46 0 : this.dim = array.length;
47 0 : if (array instanceof Array) {
48 0 : this.data = new Float64Array(this.dim);
49 0 : for (let i = 0; i < this.dim; i++) {
50 0 : this.data[i] = array[i];
51 0 : }
52 0 : } else {
53 0 : this.data = array;
54 0 : }
55 0 : }
56 0 : map(f: Function): Vector {
57 0 : this.data = this.data.map((x) => f.apply(x));
58 0 : return this;
59 0 : }
60 : /**
61 : * Create diagonal matrix from a vector
62 : * @return {Float64Array} a diagonal matrix
63 : */
64 :
65 0 : asDiagonal(): Matrix {
66 0 : var D = Zeros(this.data.length, this.data.length);
67 0 : for (var i = 0; i < this.data.length; ++i) {
68 0 : D.data[i * this.data.length + i] = this.data[i];
69 0 : }
70 0 : return D;
71 0 : }
72 0 : sum(): number {
73 0 : var sum = 0.0;
74 0 : for (let i = 0; i < this.dim; i++) {
75 0 : sum += this.data[i];
76 0 : }
77 0 : return sum;
78 0 : }
79 0 : static zeros(n: number): Vector {
80 0 : let matrix = new Float64Array(n);
81 0 : return new Vector(matrix);
82 0 : }
83 0 : add(other: Vector | number): Vector {
84 0 : if (other instanceof Vector) {
85 0 : for (let i = 0; i < this.dim; i++) {
86 0 : this.data[i] += other.data[i];
87 0 : }
88 0 : } else {
89 0 : for (let i = 0; i < this.dim; i++) {
90 0 : this.data[i] += other;
91 0 : }
92 0 : }
93 0 : return this;
94 0 : }
95 0 : copy(): Vector {
96 0 : return new Vector(this.data);
97 0 : }
98 1 : }
99 :
100 : /**
101 : * Creates a "matrix" from an existing array-like object with optional dimensions
102 : * @param {ArrayLike} array existing array, can be nested or flat (row-major)
103 : * @param {int} rows number of rows
104 : * @param {int} cols number of columns
105 : * @return {Float64Array}
106 : */
107 :
108 1 : export class Matrix {
109 : data: Float64Array;
110 : rows: number;
111 : cols: number;
112 : dim: number;
113 :
114 1 : constructor(
115 1 : array: Array<any> | Float64Array,
116 1 : rows: number,
117 1 : cols: number,
118 : ) {
119 0 : if (array instanceof Array) {
120 0 : if (Array.isArray(array[0])) { // flatten nested arrays
121 0 : rows = array.length;
122 0 : cols = array[0].length;
123 0 : this.data = new Float64Array(rows * cols);
124 0 : this.rows = rows;
125 0 : this.cols = cols;
126 0 : for (var i = 0; i < rows; ++i) {
127 0 : for (var j = 0; j < cols; ++j) {
128 0 : this.data[i * cols + j] = array[i][j];
129 0 : }
130 0 : }
131 0 : } else {
132 0 : this.data = new Float64Array(array);
133 0 : this.rows = rows || array.length;
134 0 : this.cols = cols || 1;
135 0 : }
136 0 : } else {
137 2 : this.data = array;
138 0 : this.rows = rows || array.length;
139 0 : this.cols = cols || 1;
140 2 : }
141 2 : this.dim = this.data.length;
142 1 : }
143 : /**
144 : * String representation
145 : * @param {int} precision (optional)
146 : * @return {string}
147 : */
148 0 : toString(precision: number): string {
149 0 : precision = precision || 4;
150 0 : var str = "";
151 0 : for (let i = 0; i < this.rows; ++i) {
152 0 : str += (i == 0) ? "[[ " : " [ ";
153 0 : str += this.data[i * this.cols + 0].toPrecision(precision);
154 0 : for (var j = 1; j < this.cols; ++j) {
155 0 : str += ", " + this.data[i * this.cols + j].toPrecision(precision);
156 0 : }
157 0 : str += (i == this.rows - 1) ? " ]]" : " ],\n";
158 0 : }
159 0 : return str;
160 0 : }
161 :
162 0 : asNestedArray(): Array<Array<number>> {
163 0 : let result: Array<Array<number>> = [];
164 0 : for (let i = 0; i < this.rows; i++) {
165 0 : result.push(Array.from(this.row(i).data));
166 0 : }
167 0 : return result;
168 0 : }
169 :
170 : /**
171 : * Returns a copy of a "matrix"
172 : * @return {Float64Array}
173 : */
174 0 : copy(): Matrix {
175 0 : var copy = new Float64Array(this.data);
176 :
177 0 : let rows = this.rows || this.data.length;
178 0 : let cols = this.cols || 1;
179 0 : return new Matrix(copy, rows, cols);
180 0 : }
181 :
182 : /**
183 : * (in place) Each element is replaced by a function applied to the element index
184 : * @param {function} f takes ([i, [j]]) as arguments
185 : * @return {Float64Array}
186 : */
187 0 : rebuild(f: Function): Matrix {
188 0 : if (this.cols == 1) {
189 0 : for (var i = 0; i < this.rows; ++i) {
190 0 : this.data[i] = f(i, i);
191 0 : }
192 0 : } else {
193 0 : for (var i = 0; i < this.rows; ++i) {
194 0 : for (var j = 0; j < this.cols; ++j) {
195 0 : this.data[i * this.cols + j] = f(i, j);
196 0 : }
197 0 : }
198 0 : }
199 0 : return this;
200 0 : }
201 :
202 : /**
203 : * Matrix transpose (copy)
204 : * @return {Float64Array}
205 : */
206 0 : transpose(): Matrix {
207 0 : var m = this.rows, n = this.cols;
208 0 : var transposed = Zeros(n, m);
209 0 : for (var i = 0; i < m; ++i) {
210 0 : for (var j = 0; j < n; ++j) {
211 0 : transposed.data[j * m + i] = this.data[i * n + j];
212 0 : }
213 0 : }
214 0 : return transposed;
215 0 : }
216 : /**
217 : * cwise map function onto matrix copy
218 : * @param {function} f arguments (A[i], i)
219 : * @return {Float64Array}
220 : */
221 0 : map(f: Function): Matrix {
222 0 : var A = Zeros(this.rows, this.cols);
223 0 : for (var i = 0; i < this.data.length; ++i) {
224 0 : A.data[i] = f(this.data[i], i);
225 0 : }
226 0 : return A;
227 0 : }
228 :
229 : /**
230 : * Add two matrices and return sum
231 : * @param {Float64Array} other
232 : * @return {Float64Array}
233 : */
234 0 : add(other: Matrix): Matrix {
235 0 : if (
236 0 : this.cols != other.cols || this.rows != other.rows
237 0 : ) {
238 0 : throw "matrix dimension mismatch";
239 0 : }
240 0 : var sum = Zeros(this.rows, this.cols);
241 0 : for (var i = 0; i < this.rows; ++i) {
242 0 : for (var j = 0; j < this.cols; ++j) {
243 0 : sum.data[i * this.cols + j] = this.data[i * this.cols + j] +
244 0 : other.data[i * this.cols + j];
245 0 : }
246 0 : }
247 0 : return sum;
248 0 : }
249 : /**
250 : * Increment matrix (in place)
251 : * @param {Float64Array} other
252 : * @return {Float64Array}
253 : */
254 0 : increment(other: Matrix): Matrix {
255 0 : if (
256 0 : this.cols != other.cols || this.rows != other.rows
257 0 : ) {
258 0 : throw "matrix dimension mismatch";
259 0 : }
260 0 : for (var i = 0; i < this.rows; ++i) {
261 0 : for (var j = 0; j < this.cols; ++j) {
262 0 : this.data[i * this.cols + j] += other.data[i * this.cols + j];
263 0 : }
264 0 : }
265 0 : return this;
266 0 : }
267 :
268 : /**
269 : * Decrement matrix (in place)
270 : * @param {Float64Array} other
271 : * @return {Float64Array}
272 : */
273 0 : decrement(other: Matrix): Matrix {
274 0 : if (
275 0 : this.cols != other.cols || this.rows != other.rows
276 0 : ) {
277 0 : throw "matrix dimension mismatch";
278 0 : }
279 0 : for (var i = 0; i < this.rows; ++i) {
280 0 : for (var j = 0; j < this.cols; ++j) {
281 0 : this.data[i * this.cols + j] -= other.data[i * this.cols + j];
282 0 : }
283 0 : }
284 0 : return this;
285 0 : }
286 :
287 : /**
288 : * Subtract two matrices and return difference
289 : * @param {Float64Array} other
290 : * @return {Float64Array}
291 : */
292 0 : subtract(other: Matrix): Matrix {
293 0 : if (
294 0 : this.cols != other.cols || this.rows != other.rows
295 0 : ) {
296 0 : throw "matrix dimension mismatch";
297 0 : }
298 0 : var difference = Zeros(this.rows, this.cols);
299 0 : for (var i = 0; i < this.rows; ++i) {
300 0 : for (var j = 0; j < this.cols; ++j) {
301 0 : difference.data[i * this.cols + j] = this.data[i * this.cols + j] -
302 0 : other.data[i * this.cols + j];
303 0 : }
304 0 : }
305 0 : return difference;
306 0 : }
307 :
308 : /**
309 : * Compute squared euclidian distance
310 : * @param {Float64Array} other
311 : * @return {float} square euclidian distance
312 : */
313 0 : dist2(other: Matrix): number {
314 0 : var d2 = 0;
315 0 : for (var i = 0; i < this.data.length; ++i) {
316 0 : d2 += Math.pow(this.data[i] - other.data[i], 2);
317 0 : }
318 0 : return d2;
319 0 : }
320 :
321 : /**
322 : * Compute euclidian distance
323 : * @param {Float64Array} other
324 : * @return {float} euclidian distance
325 : */
326 0 : dist(other: Matrix): number {
327 0 : return Math.sqrt(this.dist2(other));
328 0 : }
329 :
330 : /**
331 : * Multiply by scalar and return copy
332 : * @param {Float64Array} scalar
333 : * @return {Float64Array}
334 : */
335 0 : scale(scalar: number): Matrix {
336 0 : var scaled = Zeros(this.rows, this.cols);
337 0 : for (var i = 0; i < this.rows; ++i) {
338 0 : for (var j = 0; j < this.cols; ++j) {
339 0 : scaled.data[i * this.cols + j] = scalar *
340 0 : this.data[i * this.cols + j];
341 0 : }
342 0 : }
343 0 : return scaled;
344 0 : }
345 : /**
346 : * cwise negate matrix copy
347 : * @return {Float64Array}
348 : */
349 0 : negate(): Matrix {
350 0 : var A = Zeros(this.rows, this.cols);
351 0 : for (var i = 0; i < this.data.length; ++i) {
352 0 : A.data[i] = -this.data[i];
353 0 : }
354 0 : return A;
355 0 : }
356 :
357 : /**
358 : * Trace of a "matrix," i.e. sum along diagonal
359 : * @return {float}
360 : */
361 0 : trace(): number {
362 0 : let diagonal = this.diagonal();
363 0 : var trace = 0;
364 0 : for (let i = 0; i < diagonal.dim; ++i) {
365 0 : trace += diagonal.data[i];
366 0 : }
367 0 : return trace;
368 0 : }
369 :
370 : /**
371 : * Element-wise 2-norm (Frobenius-norm for matrices)
372 : * @return {float}
373 : */
374 0 : norm(): number {
375 0 : var norm = 0;
376 0 : for (var i = 0; i < this.data.length; ++i) {
377 0 : norm += this.data[i] * this.data[i];
378 0 : }
379 0 : return Math.sqrt(norm);
380 0 : }
381 :
382 : /**
383 : * Element-wise squared norm
384 : * @return {float}
385 : */
386 0 : norm2(): number {
387 0 : var norm = 0;
388 0 : for (var i = 0; i < this.data.length; ++i) {
389 0 : norm += this.data[i] * this.data[i];
390 0 : }
391 0 : return norm;
392 0 : }
393 :
394 : /**
395 : * Sum of all elements
396 : * @return {float}
397 : */
398 0 : sum(): number {
399 0 : var sum = 0;
400 0 : for (var i = 0; i < this.data.length; ++i) {
401 0 : sum += this.data[i];
402 0 : }
403 0 : return sum;
404 0 : }
405 :
406 : /**
407 : * Get diagonal as a column vector
408 : * @return {Float64Array}
409 : */
410 0 : diagonal(): Vector {
411 0 : var diagonal = Vector.zeros(Math.min(this.rows, this.cols));
412 0 : for (var i = 0; i < diagonal.dim; ++i) {
413 0 : diagonal.data[i] = this.data[i * this.cols + i];
414 0 : }
415 0 : return diagonal;
416 0 : }
417 :
418 : /**
419 : * Get row i as a row vector
420 : * @param {int} i row index
421 : * @return {Float64Array}
422 : */
423 0 : row(i: number): Vector {
424 0 : var row = Vector.zeros(this.cols);
425 0 : for (var j = 0; j < this.cols; ++j) {
426 0 : row.data[j] = this.data[i * this.cols + j];
427 0 : }
428 0 : return row;
429 0 : }
430 :
431 : /**
432 : * Get column j as as column vector
433 : * @param {int} j column index
434 : * @return {Float64Array}
435 : */
436 0 : col(j: number): Vector {
437 0 : var col = Vector.zeros(this.cols);
438 0 : for (var i = 0; i < this.rows; ++i) {
439 0 : col.data[i] = this.data[i * this.cols + j];
440 0 : }
441 0 : return col;
442 0 : }
443 :
444 0 : setRow(i: number, row: Vector | Array<number> | Float64Array): Matrix {
445 0 : if (row instanceof Vector) {
446 0 : row = row.data;
447 0 : }
448 0 : for (var j = 0; j < this.cols; ++i) {
449 0 : this.data[i * this.cols + j] = row[i];
450 0 : }
451 0 : return this;
452 0 : }
453 :
454 0 : setCol(j: number, col: Vector | Array<number> | Float64Array): Matrix {
455 0 : if (col instanceof Vector) {
456 0 : col = col.data;
457 0 : }
458 0 : for (var i = 0; i < this.rows; ++i) {
459 0 : this.data[i * this.cols + j] = col[i];
460 0 : }
461 0 : return this;
462 0 : }
463 :
464 : /**
465 : * Swap rows i and k
466 : * @param {int} i row index
467 : * @param {int} k row index
468 : * @return {Float64Array} (for chaining)
469 : */
470 0 : swap_rows(i: number, k: number): Matrix {
471 0 : for (var j = 0; j < this.cols; ++j) {
472 0 : var tmp = this.data[i * this.cols + j];
473 0 : this.data[i * this.cols + j] = this.data[k * this.cols + j];
474 0 : this.data[k * this.cols + j] = tmp;
475 0 : }
476 0 : return this;
477 0 : }
478 :
479 : /**
480 : * Computes determinant using upper triangulation
481 : * @return {float} NaN if uninvertible
482 : */
483 0 : det(): number {
484 0 : if (this.rows != this.cols) throw "det() requires square matrix";
485 0 : if (this.rows == 2 && this.cols == 2) {
486 0 : return this.data[0] * this.data[3] - this.data[1] * this.data[2];
487 0 : }
488 : // upper triangularize, then return product of diagonal
489 0 : var U = this.copy();
490 : // TODO: sure about this n?
491 0 : let n = this.rows;
492 0 : for (var i = 0; i < n; ++i) {
493 0 : var max = 0;
494 0 : for (var row = i; row < n; ++row) {
495 0 : if (Math.abs(U.data[row * n + i]) > Math.abs(U.data[max * n + i])) {
496 0 : max = row;
497 0 : }
498 0 : }
499 0 : if (max > 0) {
500 0 : U.swap_rows(i, max);
501 0 : }
502 0 : if (U.data[i * n + i] == 0) return NaN;
503 0 : for (var row = i + 1; row < n; ++row) {
504 0 : var r = U.data[row * n + i] / U.data[i * n + i];
505 0 : if (r == 0) continue;
506 0 : for (var col = i; col < n; ++col);
507 0 : U.data[row * n + col] -= U.data[i * n + col] * r;
508 0 : }
509 0 : }
510 0 : var det = 1;
511 0 : for (var i = 0; i < n; ++i) {
512 0 : det *= U.data[i * n + i];
513 0 : }
514 0 : return det;
515 0 : }
516 :
517 : /**
518 : * Generalized dot product (sum of element-wise multiplication)
519 : * @param {Float64Array} other another "matrix" of same size
520 : * @return {float}
521 : */
522 0 : dot(other: Matrix | Vector): number {
523 0 : var prod = 0;
524 0 : for (var i = 0; i < this.data.length; ++i) {
525 0 : prod += this.data[i] * other.data[i];
526 0 : }
527 0 : return prod;
528 0 : }
529 :
530 : /**
531 : * Matrix multiplication (naive implementation)
532 : * @param {Float64Array} other
533 : * @return {Float64Array}
534 : */
535 0 : multiply(other: Matrix | Vector): Matrix | Vector {
536 0 : let A = this;
537 0 : let B = other;
538 0 : var n = A.rows;
539 0 : var l = A.cols;
540 0 : if (B instanceof Matrix) {
541 0 : if (A.cols != B.rows) throw "multiply() dimension mismatch";
542 0 : let m = B.cols;
543 0 : var C = Zeros(n, m);
544 0 : for (var i = 0; i < n; ++i) {
545 0 : for (var j = 0; j < m; ++j) {
546 0 : var cij = 0;
547 0 : for (var k = 0; k < l; ++k) {
548 0 : cij += A.data[i * l + k] * B.data[k * m + j];
549 0 : }
550 0 : C.data[i * m + j] = cij;
551 0 : }
552 0 : }
553 0 : return C;
554 0 : } else {
555 : // matrix-vector product
556 0 : let V = Vector.zeros(n);
557 0 : for (var i = 0; i < n; ++i) {
558 0 : for (var j = 0; j < l; ++j) {
559 0 : V.data[i] += A.data[i * l + j] * B.data[j];
560 0 : }
561 0 : }
562 0 : return V;
563 0 : }
564 :
565 : // TODO: Add to Vector
566 : // // vector-matrix product
567 : // if (n == 1) {
568 : // for (var j = 0; j < m; ++j) {
569 : // for (var k = 0; k < l; ++k) {
570 : // C.matrix[j] += A.matrix[k] * B.matrix[k * m + j];
571 : // }
572 : // }
573 : // return C;
574 : // }
575 0 : }
576 :
577 : /**
578 : * Computes PA = LU decomposition
579 : * @return {LUResult} {L, U, P}
580 : */
581 0 : lu(): LUResult {
582 0 : if (this.rows != this.cols) throw "lu() requires square matrix";
583 0 : var n = this.rows;
584 0 : var L = Zeros(n, n);
585 0 : var U = Zeros(n, n);
586 0 : var P = Eye(n);
587 0 : for (var j = 0; j < n; ++j) {
588 0 : var max = j;
589 0 : for (var i = j; i < n; ++i) {
590 0 : if (
591 0 : Math.abs(this.data[i * n + j]) > Math.abs(this.data[max * n + j])
592 0 : ) {
593 0 : max = i;
594 0 : }
595 0 : }
596 0 : if (j != max) {
597 0 : P.swap_rows(j, max);
598 0 : }
599 0 : }
600 0 : var PA = P.multiply(this);
601 0 : for (var j = 0; j < n; ++j) {
602 0 : L.data[j * n + j] = 1;
603 0 : for (var i = 0; i < j + 1; ++i) {
604 0 : var s = 0;
605 0 : for (var k = 0; k < i; ++k) {
606 0 : s += U.data[k * n + j] * L.data[i * n + k];
607 0 : }
608 0 : U.data[i * n + j] = PA.data[i * n + j] - s;
609 0 : }
610 0 : for (var i = j; i < n; ++i) {
611 0 : var s = 0;
612 0 : for (var k = 0; k < i; ++k) {
613 0 : s += U.data[k * n + j] * L.data[i * n + k];
614 0 : }
615 0 : L.data[i * n + j] = (PA.data[i * n + j] - s) / U.data[j * n + j];
616 0 : }
617 0 : }
618 0 : return { L: L, U: U, P: P };
619 0 : }
620 :
621 : /**
622 : * Cholesky A = LL^T decomposition (in-place)
623 : * @return {[type]} [description]
624 : */
625 0 : chol_inplace(): Matrix {
626 0 : if (this.rows != this.cols) throw "chol_inplace() requires square matrix";
627 0 : let m = this.rows;
628 0 : let n = this.cols;
629 0 : var s = 0.0;
630 0 : for (let i = 0; i < n; ++i) {
631 0 : for (let j = 0; j < (i + 1); ++j) {
632 0 : s = 0.0;
633 0 : for (let k = 0; k < j; ++k) {
634 0 : s += this.data[i * n + k] * this.data[j * n + k];
635 0 : }
636 0 : if (i != j) this.data[j * n + i] = 0;
637 0 : if (
638 0 : i == j && this.data[i * n + i] - s < 0
639 0 : ) {
640 0 : throw "chol_inplace() matrix not positive definite";
641 0 : }
642 0 : this.data[i * n + j] = (i == j)
643 0 : ? Math.sqrt(this.data[i * n + i] - s)
644 0 : : ((this.data[i * n + j] - s) / this.data[j * n + j]);
645 0 : }
646 0 : }
647 0 : return this;
648 0 : }
649 :
650 : /**
651 : * Cholesky A = LL^T decomposition (returns copy)
652 : * @return {Float64Array}
653 : */
654 0 : chol(): Matrix {
655 0 : return this.copy().chol_inplace();
656 0 : }
657 :
658 : /**
659 : * Solves Lx = b using foward substitution, updates b
660 : * @param {Float64Array} b rhs
661 : * @return {Float64Array}
662 : */
663 0 : fsolve_inplace(b: Vector): Vector {
664 0 : var L = this;
665 0 : var m = L.rows, n = L.cols;
666 0 : for (var i = 0; i < n; ++i) {
667 0 : var s = 0.0;
668 0 : for (var j = 0; j < i; ++j) {
669 0 : s += L.data[i * n + j] * b.data[j];
670 0 : }
671 0 : b.data[i] = (b.data[i] - s) / L.data[i * n + i];
672 0 : }
673 0 : return b;
674 0 : }
675 :
676 : /**
677 : * Solves Lx = b using foward substitution
678 : * @param {Float64Array} b rhs
679 : * @return {Float64Array}
680 : */
681 0 : fsolve(b: Vector): Vector {
682 0 : return this.fsolve_inplace(b.copy());
683 0 : }
684 :
685 : // /**
686 : // * Solves Ux = b using backward substitution, updates b
687 : // * @param {Float64Array} b rhs
688 : // * @param {object} options {transpose: false}
689 : // * @return {Float64Array}
690 : // */
691 : // Float64Array.prototype.bsolve_inplace = function(b, options) {
692 : // var U = this;
693 : // var m = U.rows, n = U.cols;
694 : // options = options || {};
695 : // var transpose = options.hasOwnProperty('transpose') ? options.transpose : false;
696 : // for (var i = n - 1; i >= 0; --i) {
697 : // var s = 0.0;
698 : // for (var j = i + 1; j < n; ++j)
699 : // s += (transpose ? U[j * n + i] : U[i * n + j]) * b[j];
700 : // b[i] = (b[i] - s) / U[i * n + i];
701 : // }
702 : // return b;
703 : // };
704 :
705 : // /**
706 : // * Solves Ux = b using backward substitution
707 : // * @param {Float64Array} b rhs
708 : // * @param {object} options {transpose: false}
709 : // * @return {Float64Array}
710 : // */
711 : // Float64Array.prototype.bsolve = function(b, options) {
712 : // return this.bsolve_inplace(b.copy(), options);
713 : // };
714 :
715 : // /**
716 : // * Solve Ax = b using PA = LU decomposition
717 : // * @param {Float64Array} b rhs
718 : // * @return {Float64Array} x
719 : // */
720 : // Float64Array.prototype.lu_solve = function(b) {
721 : // var res = this.lu(), P = res.P, L = res.L, U = res.U;
722 : // return U.bsolve(L.fsolve(P.multiply(b)));
723 : // };
724 :
725 : /**
726 : * Computes the matrix inverse using PA = LU decomposition
727 : * @return {Float64Array} A^-1
728 : */
729 0 : lu_inverse(): Matrix {
730 0 : let res = this.lu();
731 0 : let P = res.P;
732 0 : let L = res.L;
733 0 : let U = res.U;
734 0 : var inverse = Zeros(this.rows, this.cols);
735 0 : var eye = Eye(this.rows);
736 0 : for (var j = 0; j < this.cols; ++j) {
737 0 : let mult = P.multiply(eye.col(j)) as Vector;
738 0 : inverse.setCol(j, U.bsolve(L.fsolve(mult)));
739 0 : }
740 0 : return inverse;
741 0 : }
742 :
743 : /**
744 : * Solves Ux = b using backward substitution, updates b
745 : * @param {Float64Array} b rhs
746 : * @param {object} options {transpose: false}
747 : * @return {Float64Array}
748 : */
749 0 : bsolve_inplace(b: Vector, options?: any): Vector {
750 0 : var U = this;
751 0 : var m = U.rows, n = U.cols;
752 0 : options = options || {};
753 0 : var transpose = options.hasOwnProperty("transpose")
754 0 : ? options.transpose
755 0 : : false;
756 0 : for (var i = n - 1; i >= 0; --i) {
757 0 : var s = 0.0;
758 0 : for (var j = i + 1; j < n; ++j) {
759 0 : s += (transpose ? U.data[j * n + i] : U.data[i * n + j]) * b.data[j];
760 0 : }
761 0 : b.data[i] = (b.data[i] - s) / U.data[i * n + i];
762 0 : }
763 0 : return b;
764 0 : }
765 :
766 : /**
767 : * Solves Ux = b using backward substitution
768 : * @param {Float64Array} b rhs
769 : * @param {object} options {transpose: false}
770 : * @return {Float64Array}
771 : */
772 0 : bsolve(b: Vector, options?: any): Vector {
773 0 : return this.bsolve_inplace(b.copy(), options);
774 0 : }
775 :
776 : // /**
777 : // * Solve Ax = b using A = LL^T decomposition
778 : // * @param {Float64Array} b rhs
779 : // * @return {Float64Array} x
780 : // */
781 : // Float64Array.prototype.llt_solve = function(b) {
782 : // var L = this.chol();
783 : // return L.bsolve(L.fsolve(b), {transpose: true});
784 : // };
785 :
786 : /**
787 : * Computes the matrix inverse using LL^T decomposition
788 : * @return {Float64Array} A^-1
789 : */
790 0 : llt_inverse(): Matrix {
791 0 : var L = this.chol();
792 0 : var inverse = Zeros(this.rows, this.cols);
793 0 : var eye = Eye(this.rows);
794 0 : for (var j = 0; j < this.cols; ++j) {
795 0 : inverse.setCol(j, L.bsolve(L.fsolve(eye.col(j)), { transpose: true }));
796 0 : }
797 0 : return inverse;
798 0 : }
799 :
800 : // /**
801 : // * Solve Ax = b using A = LL^T decomposition (in-place)
802 : // * @param {Float64Array} b rhs
803 : // * @return {Float64Array} x
804 : // */
805 : // Float64Array.prototype.llt_solve_inplace = function(b) {
806 : // var L = this.chol_inplace();
807 : // return L.bsolve_inplace(L.fsolve_inplace(b), {transpose: true});
808 : // };
809 :
810 : // /**
811 : // * Computes diagonal matrix D of eigenvalues and matrix V whose columns are the corresponding
812 : // * right eigenvectors so that AV = VD
813 : // * @param {object} options tolerance and maxIter
814 : // * @return {object} V:V D:D
815 : // */
816 : // Float64Array.prototype.jacobiRotation = function(options) {
817 :
818 : // if (this.cols != this.rows) throw 'matrix must be square';
819 :
820 : // if (arguments.length < 1)
821 : // options = {};
822 :
823 : // var maxIter = options.maxIter || 100;
824 : // var tolerance = options.tolerance || 1e-5;
825 :
826 : // var n = this.rows;
827 : // var D = this.copy();
828 : // var V = Float64Array.eye(n, n);
829 :
830 : // var iter, maxOffDiag, p, q;
831 : // for (iter = 0; iter < maxIter; ++iter) {
832 :
833 : // // find max off diagonal term at (p, q)
834 : // maxOffDiag = 0;
835 : // for (var i = 0; i < n - 1; ++i) {
836 : // for (var j = i + 1; j < n; ++j) {
837 : // if (Math.abs(D[i * n + j]) > maxOffDiag) {
838 : // maxOffDiag = Math.abs(D[i * n + j]);
839 : // p = i; q = j;
840 : // }
841 : // }
842 : // }
843 :
844 : // if (maxOffDiag < tolerance)
845 : // break;
846 :
847 : // // Rotates matrix D through theta in pq-plane to set D[p][q] = 0
848 : // // Rotation stored in matrix V whose columns are eigenvectors of D
849 : // // d = cot 2 * theta, t = tan theta, c = cos theta, s = sin theta
850 : // var d = (D[p * n + p] - D[q * n + q]) / (2.0 * D[p * n + q]);
851 : // var t = Math.sign(d) / (Math.abs(d) + Math.sqrt(d * d + 1));
852 : // var c = 1.0 / Math.sqrt(t * t + 1);
853 : // var s = t * c;
854 : // D[p * n + p] += t * D[p * n + q];
855 : // D[q * n + q] -= t * D[p * n + q];
856 : // D[p * n + q] = D[q * n + p] = 0.0;
857 : // for (var k = 0; k < n; k++) { // Transform D
858 : // if (k != p && k != q) {
859 : // var akp = c * D[k * n + p] + s * D[k * n + q];
860 : // var akq = -s * D[k * n + p] + c * D[k * n + q];
861 : // D[k * n + p] = akp;
862 : // D[p * n + k] = akp;
863 : // D[k * n + q] = akq;
864 : // D[q * n + k] = akq;
865 : // }
866 : // }
867 : // for (var k = 0; k < n; k++) { // Store V
868 : // var rkp = c * V[k * n + p] + s * V[k * n + q];
869 : // var rkq = -s * V[k * n + p] + c * V[k * n + q];
870 : // V[k * n + p] = rkp;
871 : // V[k * n + q] = rkq;
872 : // }
873 : // }
874 :
875 : // if (iter == maxIter) {
876 : // console.log('Hit maxIter: ', maxOffDiag, ' > ', tolerance);
877 : // }
878 :
879 : // return {V:V, D:D, eigenvalues: D.diagonal(), eigenvectors: V};
880 :
881 : // };
882 :
883 : // Float64Array.prototype.maxCoeff = function() {
884 : // var max = this[0];
885 : // for (var i = 0; i < this.length; ++i) {
886 : // if (this[i] > max)
887 : // max = this[i];
888 : // }
889 : // return max;
890 : // };
891 :
892 : // Float64Array.prototype.cwiseProduct = function(other) {
893 : // var A = this.copy();
894 : // for (var i = 0; i < this.length; ++i) {
895 : // A[i] = this[i] * other[i];
896 : // }
897 : // return A;
898 : // };
899 :
900 : // Float64Array.prototype.cwiseQuotient = function(other) {
901 : // var A = this.copy();
902 : // for (var i = 0; i < this.length; ++i) {
903 : // A[i] = this[i] / other[i];
904 : // }
905 : // return A;
906 : // };
907 :
908 : // Float64Array.prototype.cwiseInverse = function() {
909 : // var A = this.copy();
910 : // for (var i = 0; i < this.length; ++i) {
911 : // A[i] = 1.0 / this[i];
912 : // }
913 : // return A;
914 : // };
915 :
916 : // Float64Array.prototype.cwiseSqrt = function() {
917 : // var A = this.copy();
918 : // for (var i = 0; i < this.length; ++i) {
919 : // A[i] = Math.sqrt(this[i]);
920 : // }
921 : // return A;
922 : // };
923 :
924 : /*
925 : Shanti Rao sent me this routine by private email. I had to modify it
926 : slightly to work on Arrays instead of using a Matrix object.
927 : It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py
928 : */
929 :
930 0 : svd(): SVDResult {
931 0 : let A = this;
932 0 : var temp;
933 : //Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970)
934 0 : var prec = Math.pow(2, -52); // assumes double prec
935 0 : var tolerance = 1.e-64 / prec;
936 0 : var itmax = 50;
937 0 : var c = 0;
938 0 : var i = 0;
939 0 : var j = 0;
940 0 : var k = 0;
941 0 : var l = 0;
942 :
943 0 : var u = A.copy().asNestedArray();
944 0 : var m = u.length;
945 :
946 0 : var n = u[0].length;
947 :
948 0 : if (m < n) throw "Need more rows than columns";
949 :
950 0 : var e = new Array(n);
951 0 : var q: Array<number> = new Array(n);
952 0 : for (i = 0; i < n; i++) e[i] = q[i] = 0.0;
953 0 : var v: Array<Array<number>> = utils.rep([n, n], 0);
954 : // v.zero();
955 :
956 0 : let pythag = (a: number, b: number) => {
957 0 : a = Math.abs(a);
958 0 : b = Math.abs(b);
959 0 : if (a > b) {
960 0 : return a * Math.sqrt(1.0 + (b * b / a / a));
961 0 : } else if (b == 0.0) {
962 0 : return a;
963 0 : }
964 0 : return b * Math.sqrt(1.0 + (a * a / b / b));
965 0 : };
966 :
967 : //Householder's reduction to bidiagonal form
968 :
969 0 : var f = 0.0;
970 0 : var g = 0.0;
971 0 : var h = 0.0;
972 0 : var x = 0.0;
973 0 : var y = 0.0;
974 0 : var z = 0.0;
975 0 : var s = 0.0;
976 :
977 0 : for (i = 0; i < n; i++) {
978 0 : e[i] = g;
979 0 : s = 0.0;
980 0 : l = i + 1;
981 0 : for (j = i; j < m; j++) {
982 0 : s += (u[j][i] * u[j][i]);
983 0 : }
984 0 : if (s <= tolerance) {
985 0 : g = 0.0;
986 0 : } else {
987 0 : f = u[i][i];
988 0 : g = Math.sqrt(s);
989 0 : if (f >= 0.0) g = -g;
990 0 : h = f * g - s;
991 0 : u[i][i] = f - g;
992 0 : for (j = l; j < n; j++) {
993 0 : s = 0.0;
994 0 : for (k = i; k < m; k++) {
995 0 : s += u[k][i] * u[k][j];
996 0 : }
997 0 : f = s / h;
998 0 : for (k = i; k < m; k++) {
999 0 : u[k][j] += f * u[k][i];
1000 0 : }
1001 0 : }
1002 0 : }
1003 0 : q[i] = g;
1004 0 : s = 0.0;
1005 0 : for (j = l; j < n; j++) {
1006 0 : s = s + u[i][j] * u[i][j];
1007 0 : }
1008 0 : if (s <= tolerance) {
1009 0 : g = 0.0;
1010 0 : } else {
1011 0 : f = u[i][i + 1];
1012 0 : g = Math.sqrt(s);
1013 0 : if (f >= 0.0) g = -g;
1014 0 : h = f * g - s;
1015 0 : u[i][i + 1] = f - g;
1016 0 : for (j = l; j < n; j++) e[j] = u[i][j] / h;
1017 0 : for (j = l; j < m; j++) {
1018 0 : s = 0.0;
1019 0 : for (k = l; k < n; k++) {
1020 0 : s += (u[j][k] * u[i][k]);
1021 0 : }
1022 0 : for (k = l; k < n; k++) {
1023 0 : u[j][k] += s * e[k];
1024 0 : }
1025 0 : }
1026 0 : }
1027 0 : y = Math.abs(q[i]) + Math.abs(e[i]);
1028 0 : if (y > x) {
1029 0 : x = y;
1030 0 : }
1031 0 : }
1032 :
1033 : // accumulation of right hand gtransformations
1034 0 : for (i = n - 1; i != -1; i += -1) {
1035 0 : if (g != 0.0) {
1036 0 : h = g * u[i][i + 1];
1037 0 : for (j = l; j < n; j++) {
1038 0 : v[j][i] = u[i][j] / h;
1039 0 : }
1040 0 : for (j = l; j < n; j++) {
1041 0 : s = 0.0;
1042 0 : for (k = l; k < n; k++) {
1043 0 : s += u[i][k] * v[k][j];
1044 0 : }
1045 0 : for (k = l; k < n; k++) {
1046 0 : v[k][j] += (s * v[k][i]);
1047 0 : }
1048 0 : }
1049 0 : }
1050 0 : for (j = l; j < n; j++) {
1051 0 : v[i][j] = 0;
1052 0 : v[j][i] = 0;
1053 0 : }
1054 0 : v[i][i] = 1;
1055 0 : g = e[i];
1056 0 : l = i;
1057 0 : }
1058 :
1059 : // accumulation of left hand transformations
1060 0 : for (i = n - 1; i != -1; i += -1) {
1061 0 : l = i + 1;
1062 0 : g = q[i];
1063 0 : for (j = l; j < n; j++) {
1064 0 : u[i][j] = 0;
1065 0 : }
1066 0 : if (g != 0.0) {
1067 0 : h = u[i][i] * g;
1068 0 : for (j = l; j < n; j++) {
1069 0 : s = 0.0;
1070 0 : for (k = l; k < m; k++) s += u[k][i] * u[k][j];
1071 0 : f = s / h;
1072 0 : for (k = i; k < m; k++) u[k][j] += f * u[k][i];
1073 0 : }
1074 0 : for (j = i; j < m; j++) u[j][i] = u[j][i] / g;
1075 0 : } else {
1076 0 : for (j = i; j < m; j++) u[j][i] = 0;
1077 0 : }
1078 0 : u[i][i] += 1;
1079 0 : }
1080 :
1081 : // diagonalization of the bidiagonal form
1082 0 : prec = prec * x;
1083 0 : for (k = n - 1; k != -1; k += -1) {
1084 0 : for (var iteration = 0; iteration < itmax; iteration++) { // test f splitting
1085 0 : var test_convergence = false;
1086 0 : for (l = k; l != -1; l += -1) {
1087 0 : if (Math.abs(e[l]) <= prec) {
1088 0 : test_convergence = true;
1089 0 : break;
1090 0 : }
1091 0 : if (Math.abs(q[l - 1]) <= prec) {
1092 0 : break;
1093 0 : }
1094 0 : }
1095 0 : if (!test_convergence) { // cancellation of e[l] if l>0
1096 0 : c = 0.0;
1097 0 : s = 1.0;
1098 0 : var l1 = l - 1;
1099 0 : for (i = l; i < k + 1; i++) {
1100 0 : f = s * e[i];
1101 0 : e[i] = c * e[i];
1102 0 : if (Math.abs(f) <= prec) {
1103 0 : break;
1104 0 : }
1105 0 : g = q[i];
1106 0 : h = pythag(f, g);
1107 0 : q[i] = h;
1108 0 : c = g / h;
1109 0 : s = -f / h;
1110 0 : for (j = 0; j < m; j++) {
1111 0 : y = u[j][l1];
1112 0 : z = u[j][i];
1113 0 : u[j][l1] = y * c + (z * s);
1114 0 : u[j][i] = -y * s + (z * c);
1115 0 : }
1116 0 : }
1117 0 : }
1118 : // test f convergence
1119 0 : z = q[k];
1120 0 : if (l == k) { //convergence
1121 0 : if (z < 0.0) { //q[k] is made non-negative
1122 0 : q[k] = -z;
1123 0 : for (j = 0; j < n; j++) {
1124 0 : v[j][k] = -v[j][k];
1125 0 : }
1126 0 : }
1127 0 : break; //break out of iteration loop and move on to next k value
1128 0 : }
1129 0 : if (iteration >= itmax - 1) {
1130 0 : throw "Error: no convergence.";
1131 0 : }
1132 : // shift from bottom 2x2 minor
1133 0 : x = q[l];
1134 0 : y = q[k - 1];
1135 0 : g = e[k - 1];
1136 0 : h = e[k];
1137 0 : f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
1138 0 : g = pythag(f, 1.0);
1139 0 : if (f < 0.0) {
1140 0 : f = ((x - z) * (x + z) + h * (y / (f - g) - h)) / x;
1141 0 : } else {
1142 0 : f = ((x - z) * (x + z) + h * (y / (f + g) - h)) / x;
1143 0 : }
1144 : // next QR transformation
1145 0 : c = 1.0;
1146 0 : s = 1.0;
1147 0 : for (i = l + 1; i < k + 1; i++) {
1148 0 : g = e[i];
1149 0 : y = q[i];
1150 0 : h = s * g;
1151 0 : g = c * g;
1152 0 : z = pythag(f, h);
1153 0 : e[i - 1] = z;
1154 0 : c = f / z;
1155 0 : s = h / z;
1156 0 : f = x * c + g * s;
1157 0 : g = -x * s + g * c;
1158 0 : h = y * s;
1159 0 : y = y * c;
1160 0 : for (j = 0; j < n; j++) {
1161 0 : x = v[j][i - 1];
1162 0 : z = v[j][i];
1163 0 : v[j][i - 1] = x * c + z * s;
1164 0 : v[j][i] = -x * s + z * c;
1165 0 : }
1166 0 : z = pythag(f, h);
1167 0 : q[i - 1] = z;
1168 0 : c = f / z;
1169 0 : s = h / z;
1170 0 : f = c * g + s * y;
1171 0 : x = -s * g + c * y;
1172 0 : for (j = 0; j < m; j++) {
1173 0 : y = u[j][i - 1];
1174 0 : z = u[j][i];
1175 0 : u[j][i - 1] = y * c + z * s;
1176 0 : u[j][i] = -y * s + z * c;
1177 0 : }
1178 0 : }
1179 0 : e[l] = 0.0;
1180 0 : e[k] = f;
1181 0 : q[k] = x;
1182 0 : }
1183 0 : }
1184 :
1185 : //vt= transpose(v)
1186 : //return (u,q,vt)
1187 0 : for (i = 0; i < q.length; i++) {
1188 0 : if (q[i] < prec) q[i] = 0;
1189 0 : }
1190 :
1191 : //sort eigenvalues
1192 0 : for (i = 0; i < n; i++) {
1193 : //writeln(q)
1194 0 : for (j = i - 1; j >= 0; j--) {
1195 0 : if (q[j] < q[i]) {
1196 : // writeln(i,'-',j)
1197 0 : c = q[j];
1198 0 : q[j] = q[i];
1199 0 : q[i] = c;
1200 0 : for (k = 0; k < u.length; k++) {
1201 0 : temp = u[k][i];
1202 0 : u[k][i] = u[k][j];
1203 0 : u[k][j] = temp;
1204 0 : }
1205 0 : for (k = 0; k < v.length; k++) {
1206 0 : temp = v[k][i];
1207 0 : v[k][i] = v[k][j];
1208 0 : v[k][j] = temp;
1209 0 : }
1210 : // u.swapCols(i,j)
1211 : // v.swapCols(i,j)
1212 0 : i = j;
1213 0 : }
1214 0 : }
1215 0 : }
1216 :
1217 0 : let U = new Matrix(u, this.rows, this.cols);
1218 0 : console.log(`u: ${u}`);
1219 0 : console.log(`q: ${q}`);
1220 0 : console.log(`v: ${v}`);
1221 0 : let S = new Vector(q);
1222 0 : let V = new Matrix(v, this.rows, this.cols);
1223 :
1224 0 : return { U: U, S: S, V: V };
1225 0 : }
1226 1 : }
1227 :
1228 : /**
1229 : * Creates a "matrix" with all entries set to zero
1230 : * @param {int} rows number of rows
1231 : * @param {int} cols number of columns
1232 : * @return {Float64Array}
1233 : */
1234 1 : export function Zeros(rows: number, cols: number): Matrix {
1235 0 : cols = cols || 1;
1236 2 : var matrix = new Float64Array(rows * cols);
1237 2 : return new Matrix(matrix, rows, cols);
1238 1 : }
1239 :
1240 : /**
1241 : * Creates a column vector with linearly-spaced elements
1242 : * @param {float} min minimum value (inclusive)
1243 : * @param {float} max maximum value (inclusive)
1244 : * @param {int} n number of elements
1245 : * @return {Float64Array}
1246 : */
1247 0 : export function Linspace(min: number, max: number, n: number): Float64Array {
1248 0 : var matrix = new Float64Array(n);
1249 0 : var dx = (max - min) / (n - 1);
1250 0 : for (var i = 0; i < n; ++i) {
1251 0 : matrix[i] = i * dx + min;
1252 0 : }
1253 0 : return matrix;
1254 0 : }
1255 :
1256 : /**
1257 : * Creates a n x n identity matrix
1258 : * @param {int} n number of rows and columns
1259 : * @return {Float64Array}
1260 : *
1261 : * const n = 3;
1262 : * const diagonal = Eye(n);
1263 : */
1264 0 : export function Eye(n: number): Matrix {
1265 0 : var matrix = new Float64Array(n * n);
1266 0 : for (var i = 0; i < n; ++i) {
1267 0 : matrix[i * n + i] = 1;
1268 0 : }
1269 0 : return new Matrix(matrix, n, n);
1270 0 : } /**
1271 : * Creates a "matrix" filled with ones
1272 : * @param {int} rows number of rows
1273 : * @param {int} cols number of columns
1274 : * @return {Float64Array}
1275 : */
1276 :
1277 0 : export function Ones(rows: number, cols: number): Matrix {
1278 0 : cols = cols || 1;
1279 0 : var matrix = new Float64Array(rows * cols);
1280 0 : for (var i = 0; i < matrix.length; ++i) {
1281 0 : matrix[i] = 1;
1282 0 : }
1283 0 : return new Matrix(matrix, rows, cols);
1284 0 : } /**
1285 : * Creates a "matrix" filled with constant
1286 : * @param {float} const constant
1287 : * @param {int} rows number of rows
1288 : * @param {int} cols number of columns
1289 : * @return {Float64Array}
1290 : */
1291 :
1292 0 : export function Constant(constant: number, rows: number, cols: number): Matrix {
1293 0 : cols = cols || 1;
1294 0 : var matrix = new Float64Array(rows * cols);
1295 0 : for (var i = 0; i < matrix.length; ++i) {
1296 0 : matrix[i] = constant;
1297 0 : }
1298 0 : return new Matrix(matrix, rows, cols);
1299 0 : } /**
1300 : * Build a "matrix" where each element is a function applied to element index
1301 : * @param {function} f takes ([i, [j]]) as arguments
1302 : * @param {int} rows number of rows
1303 : * @param {int} cols number of cols
1304 : * @return {Float64Array}
1305 : */
1306 :
1307 0 : export function Build(f: Function, rows: number, cols: number) {
1308 0 : cols = cols || 1;
1309 0 : var matrix = Zeros(rows, cols);
1310 0 : for (var i = 0; i < rows; ++i) {
1311 0 : for (var j = 0; j < cols; ++j) {
1312 0 : matrix.data[i * cols + j] = f(i, j);
1313 0 : }
1314 0 : }
1315 0 : return new Matrix(matrix.data, rows, cols);
1316 0 : } /**
1317 : * Outer product (form matrix from vector tensor product)
1318 : * @param {Float64Array} u vector (column or row)
1319 : * @param {Float64Array} v vector (column or row)
1320 : * @return {Float64Array}
1321 : */
1322 :
1323 0 : export function Outer(u: Float64Array, v: Float64Array) {
1324 0 : var A = Zeros(u.length, v.length);
1325 0 : for (var i = 0; i < u.length; ++i) {
1326 0 : for (var j = 0; j < v.length; ++j) {
1327 0 : A.data[i * v.length + j] = u[i] * v[j];
1328 0 : }
1329 0 : }
1330 0 : return A;
1331 0 : }
|
