123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184 |
- import { isMatrix } from '../../utils/is.js';
- import { arraySize } from '../../utils/array.js';
- import { factory } from '../../utils/factory.js';
- import { format } from '../../utils/string.js';
- var name = 'inv';
- var dependencies = ['typed', 'matrix', 'divideScalar', 'addScalar', 'multiply', 'unaryMinus', 'det', 'identity', 'abs'];
- export var createInv = /* #__PURE__ */factory(name, dependencies, _ref => {
- var {
- typed,
- matrix,
- divideScalar,
- addScalar,
- multiply,
- unaryMinus,
- det,
- identity,
- abs
- } = _ref;
- /**
- * Calculate the inverse of a square matrix.
- *
- * Syntax:
- *
- * math.inv(x)
- *
- * Examples:
- *
- * math.inv([[1, 2], [3, 4]]) // returns [[-2, 1], [1.5, -0.5]]
- * math.inv(4) // returns 0.25
- * 1 / 4 // returns 0.25
- *
- * See also:
- *
- * det, transpose
- *
- * @param {number | Complex | Array | Matrix} x Matrix to be inversed
- * @return {number | Complex | Array | Matrix} The inverse of `x`.
- */
- return typed(name, {
- 'Array | Matrix': function ArrayMatrix(x) {
- var size = isMatrix(x) ? x.size() : arraySize(x);
- switch (size.length) {
- case 1:
- // vector
- if (size[0] === 1) {
- if (isMatrix(x)) {
- return matrix([divideScalar(1, x.valueOf()[0])]);
- } else {
- return [divideScalar(1, x[0])];
- }
- } else {
- throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
- }
- case 2:
- // two dimensional array
- {
- var rows = size[0];
- var cols = size[1];
- if (rows === cols) {
- if (isMatrix(x)) {
- return matrix(_inv(x.valueOf(), rows, cols), x.storage());
- } else {
- // return an Array
- return _inv(x, rows, cols);
- }
- } else {
- throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
- }
- }
- default:
- // multi dimensional array
- throw new RangeError('Matrix must be two dimensional ' + '(size: ' + format(size) + ')');
- }
- },
- any: function any(x) {
- // scalar
- return divideScalar(1, x); // FIXME: create a BigNumber one when configured for bignumbers
- }
- });
- /**
- * Calculate the inverse of a square matrix
- * @param {Array[]} mat A square matrix
- * @param {number} rows Number of rows
- * @param {number} cols Number of columns, must equal rows
- * @return {Array[]} inv Inverse matrix
- * @private
- */
- function _inv(mat, rows, cols) {
- var r, s, f, value, temp;
- if (rows === 1) {
- // this is a 1 x 1 matrix
- value = mat[0][0];
- if (value === 0) {
- throw Error('Cannot calculate inverse, determinant is zero');
- }
- return [[divideScalar(1, value)]];
- } else if (rows === 2) {
- // this is a 2 x 2 matrix
- var d = det(mat);
- if (d === 0) {
- throw Error('Cannot calculate inverse, determinant is zero');
- }
- return [[divideScalar(mat[1][1], d), divideScalar(unaryMinus(mat[0][1]), d)], [divideScalar(unaryMinus(mat[1][0]), d), divideScalar(mat[0][0], d)]];
- } else {
- // this is a matrix of 3 x 3 or larger
- // calculate inverse using gauss-jordan elimination
- // https://en.wikipedia.org/wiki/Gaussian_elimination
- // http://mathworld.wolfram.com/MatrixInverse.html
- // http://math.uww.edu/~mcfarlat/inverse.htm
- // make a copy of the matrix (only the arrays, not of the elements)
- var A = mat.concat();
- for (r = 0; r < rows; r++) {
- A[r] = A[r].concat();
- }
- // create an identity matrix which in the end will contain the
- // matrix inverse
- var B = identity(rows).valueOf();
- // loop over all columns, and perform row reductions
- for (var c = 0; c < cols; c++) {
- // Pivoting: Swap row c with row r, where row r contains the largest element A[r][c]
- var ABig = abs(A[c][c]);
- var rBig = c;
- r = c + 1;
- while (r < rows) {
- if (abs(A[r][c]) > ABig) {
- ABig = abs(A[r][c]);
- rBig = r;
- }
- r++;
- }
- if (ABig === 0) {
- throw Error('Cannot calculate inverse, determinant is zero');
- }
- r = rBig;
- if (r !== c) {
- temp = A[c];
- A[c] = A[r];
- A[r] = temp;
- temp = B[c];
- B[c] = B[r];
- B[r] = temp;
- }
- // eliminate non-zero values on the other rows at column c
- var Ac = A[c];
- var Bc = B[c];
- for (r = 0; r < rows; r++) {
- var Ar = A[r];
- var Br = B[r];
- if (r !== c) {
- // eliminate value at column c and row r
- if (Ar[c] !== 0) {
- f = divideScalar(unaryMinus(Ar[c]), Ac[c]);
- // add (f * row c) to row r to eliminate the value
- // at column c
- for (s = c; s < cols; s++) {
- Ar[s] = addScalar(Ar[s], multiply(f, Ac[s]));
- }
- for (s = 0; s < cols; s++) {
- Br[s] = addScalar(Br[s], multiply(f, Bc[s]));
- }
- }
- } else {
- // normalize value at Acc to 1,
- // divide each value on row r with the value at Acc
- f = Ac[c];
- for (s = c; s < cols; s++) {
- Ar[s] = divideScalar(Ar[s], f);
- }
- for (s = 0; s < cols; s++) {
- Br[s] = divideScalar(Br[s], f);
- }
- }
- }
- }
- return B;
- }
- }
- });
|