simple-squiggle

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

eigs.js (6984B)

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