simple-squiggle

A restricted subset of Squiggle
Log | Files | Refs | README

eigs.js (8145B)


      1 "use strict";
      2 
      3 Object.defineProperty(exports, "__esModule", {
      4   value: true
      5 });
      6 exports.createEigs = void 0;
      7 
      8 var _factory = require("../../utils/factory.js");
      9 
     10 var _string = require("../../utils/string.js");
     11 
     12 var _complexEigs = require("./eigs/complexEigs.js");
     13 
     14 var _realSymetric = require("./eigs/realSymetric.js");
     15 
     16 var _is = require("../../utils/is.js");
     17 
     18 var name = 'eigs'; // The absolute state of math.js's dependency system:
     19 
     20 var dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolve', 'usolveAll', 'im', 're', 'smaller', 'matrixFromColumns', 'dot'];
     21 var createEigs = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
     22   var config = _ref.config,
     23       typed = _ref.typed,
     24       matrix = _ref.matrix,
     25       addScalar = _ref.addScalar,
     26       subtract = _ref.subtract,
     27       equal = _ref.equal,
     28       abs = _ref.abs,
     29       atan = _ref.atan,
     30       cos = _ref.cos,
     31       sin = _ref.sin,
     32       multiplyScalar = _ref.multiplyScalar,
     33       divideScalar = _ref.divideScalar,
     34       inv = _ref.inv,
     35       bignumber = _ref.bignumber,
     36       multiply = _ref.multiply,
     37       add = _ref.add,
     38       larger = _ref.larger,
     39       column = _ref.column,
     40       flatten = _ref.flatten,
     41       number = _ref.number,
     42       complex = _ref.complex,
     43       sqrt = _ref.sqrt,
     44       diag = _ref.diag,
     45       qr = _ref.qr,
     46       usolve = _ref.usolve,
     47       usolveAll = _ref.usolveAll,
     48       im = _ref.im,
     49       re = _ref.re,
     50       smaller = _ref.smaller,
     51       matrixFromColumns = _ref.matrixFromColumns,
     52       dot = _ref.dot;
     53   var doRealSymetric = (0, _realSymetric.createRealSymmetric)({
     54     config: config,
     55     addScalar: addScalar,
     56     subtract: subtract,
     57     column: column,
     58     flatten: flatten,
     59     equal: equal,
     60     abs: abs,
     61     atan: atan,
     62     cos: cos,
     63     sin: sin,
     64     multiplyScalar: multiplyScalar,
     65     inv: inv,
     66     bignumber: bignumber,
     67     complex: complex,
     68     multiply: multiply,
     69     add: add
     70   });
     71   var doComplexEigs = (0, _complexEigs.createComplexEigs)({
     72     config: config,
     73     addScalar: addScalar,
     74     subtract: subtract,
     75     multiply: multiply,
     76     multiplyScalar: multiplyScalar,
     77     flatten: flatten,
     78     divideScalar: divideScalar,
     79     sqrt: sqrt,
     80     abs: abs,
     81     bignumber: bignumber,
     82     diag: diag,
     83     qr: qr,
     84     inv: inv,
     85     usolve: usolve,
     86     usolveAll: usolveAll,
     87     equal: equal,
     88     complex: complex,
     89     larger: larger,
     90     smaller: smaller,
     91     matrixFromColumns: matrixFromColumns,
     92     dot: dot
     93   });
     94   /**
     95    * Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending.
     96    * An eigenvalue with multiplicity k will be listed k times. The eigenvectors are returned as columns of a matrix –
     97    * the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
     98    * If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information
     99    * in `err.values` and `err.vectors`.
    100    *
    101    * Syntax:
    102    *
    103    *     math.eigs(x, [prec])
    104    *
    105    * Examples:
    106    *
    107    *     const { eigs, multiply, column, transpose } = math
    108    *     const H = [[5, 2.3], [2.3, 1]]
    109    *     const ans = eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
    110    *     const E = ans.values
    111    *     const U = ans.vectors
    112    *     multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0))
    113    *     const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H
    114    *     E[0] == UTxHxU[0][0]  // returns true
    115    *
    116    * See also:
    117    *
    118    *     inv
    119    *
    120    * @param {Array | Matrix} x  Matrix to be diagonalized
    121    *
    122    * @param {number | BigNumber} [prec] Precision, default value: 1e-15
    123    * @return {{values: Array|Matrix, vectors: Array|Matrix}} Object containing an array of eigenvalues and a matrix with eigenvectors as columns.
    124    *
    125    */
    126 
    127   return typed('eigs', {
    128     Array: function Array(x) {
    129       var mat = matrix(x);
    130       return computeValuesAndVectors(mat);
    131     },
    132     'Array, number|BigNumber': function ArrayNumberBigNumber(x, prec) {
    133       var mat = matrix(x);
    134       return computeValuesAndVectors(mat, prec);
    135     },
    136     Matrix: function Matrix(mat) {
    137       var _computeValuesAndVect = computeValuesAndVectors(mat),
    138           values = _computeValuesAndVect.values,
    139           vectors = _computeValuesAndVect.vectors;
    140 
    141       return {
    142         values: matrix(values),
    143         vectors: matrix(vectors)
    144       };
    145     },
    146     'Matrix, number|BigNumber': function MatrixNumberBigNumber(mat, prec) {
    147       var _computeValuesAndVect2 = computeValuesAndVectors(mat, prec),
    148           values = _computeValuesAndVect2.values,
    149           vectors = _computeValuesAndVect2.vectors;
    150 
    151       return {
    152         values: matrix(values),
    153         vectors: matrix(vectors)
    154       };
    155     }
    156   });
    157 
    158   function computeValuesAndVectors(mat, prec) {
    159     if (prec === undefined) {
    160       prec = config.epsilon;
    161     }
    162 
    163     var size = mat.size();
    164 
    165     if (size.length !== 2 || size[0] !== size[1]) {
    166       throw new RangeError('Matrix must be square (size: ' + (0, _string.format)(size) + ')');
    167     }
    168 
    169     var arr = mat.toArray();
    170     var N = size[0];
    171 
    172     if (isReal(arr, N, prec)) {
    173       coerceReal(arr, N);
    174 
    175       if (isSymmetric(arr, N, prec)) {
    176         var _type = coerceTypes(mat, arr, N);
    177 
    178         return doRealSymetric(arr, N, prec, _type);
    179       }
    180     }
    181 
    182     var type = coerceTypes(mat, arr, N);
    183     return doComplexEigs(arr, N, prec, type);
    184   }
    185   /** @return {boolean} */
    186 
    187 
    188   function isSymmetric(arr, N, prec) {
    189     for (var i = 0; i < N; i++) {
    190       for (var j = i; j < N; j++) {
    191         // TODO proper comparison of bignum and frac
    192         if (larger(bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec)) {
    193           return false;
    194         }
    195       }
    196     }
    197 
    198     return true;
    199   }
    200   /** @return {boolean} */
    201 
    202 
    203   function isReal(arr, N, prec) {
    204     for (var i = 0; i < N; i++) {
    205       for (var j = 0; j < N; j++) {
    206         // TODO proper comparison of bignum and frac
    207         if (larger(bignumber(abs(im(arr[i][j]))), prec)) {
    208           return false;
    209         }
    210       }
    211     }
    212 
    213     return true;
    214   }
    215 
    216   function coerceReal(arr, N) {
    217     for (var i = 0; i < N; i++) {
    218       for (var j = 0; j < N; j++) {
    219         arr[i][j] = re(arr[i][j]);
    220       }
    221     }
    222   }
    223   /** @return {'number' | 'BigNumber' | 'Complex'} */
    224 
    225 
    226   function coerceTypes(mat, arr, N) {
    227     /** @type {string} */
    228     var type = mat.datatype();
    229 
    230     if (type === 'number' || type === 'BigNumber' || type === 'Complex') {
    231       return type;
    232     }
    233 
    234     var hasNumber = false;
    235     var hasBig = false;
    236     var hasComplex = false;
    237 
    238     for (var i = 0; i < N; i++) {
    239       for (var j = 0; j < N; j++) {
    240         var el = arr[i][j];
    241 
    242         if ((0, _is.isNumber)(el) || (0, _is.isFraction)(el)) {
    243           hasNumber = true;
    244         } else if ((0, _is.isBigNumber)(el)) {
    245           hasBig = true;
    246         } else if ((0, _is.isComplex)(el)) {
    247           hasComplex = true;
    248         } else {
    249           throw TypeError('Unsupported type in Matrix: ' + (0, _is.typeOf)(el));
    250         }
    251       }
    252     }
    253 
    254     if (hasBig && hasComplex) {
    255       console.warn('Complex BigNumbers not supported, this operation will lose precission.');
    256     }
    257 
    258     if (hasComplex) {
    259       for (var _i = 0; _i < N; _i++) {
    260         for (var _j = 0; _j < N; _j++) {
    261           arr[_i][_j] = complex(arr[_i][_j]);
    262         }
    263       }
    264 
    265       return 'Complex';
    266     }
    267 
    268     if (hasBig) {
    269       for (var _i2 = 0; _i2 < N; _i2++) {
    270         for (var _j2 = 0; _j2 < N; _j2++) {
    271           arr[_i2][_j2] = bignumber(arr[_i2][_j2]);
    272         }
    273       }
    274 
    275       return 'BigNumber';
    276     }
    277 
    278     if (hasNumber) {
    279       for (var _i3 = 0; _i3 < N; _i3++) {
    280         for (var _j3 = 0; _j3 < N; _j3++) {
    281           arr[_i3][_j3] = number(arr[_i3][_j3]);
    282         }
    283       }
    284 
    285       return 'number';
    286     } else {
    287       throw TypeError('Matrix contains unsupported types only.');
    288     }
    289   }
    290 });
    291 exports.createEigs = createEigs;