simple-squiggle

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

realSymetric.js (8202B)


      1 "use strict";
      2 
      3 Object.defineProperty(exports, "__esModule", {
      4   value: true
      5 });
      6 exports.createRealSymmetric = createRealSymmetric;
      7 
      8 var _object = require("../../../utils/object.js");
      9 
     10 function createRealSymmetric(_ref) {
     11   var config = _ref.config,
     12       addScalar = _ref.addScalar,
     13       subtract = _ref.subtract,
     14       abs = _ref.abs,
     15       atan = _ref.atan,
     16       cos = _ref.cos,
     17       sin = _ref.sin,
     18       multiplyScalar = _ref.multiplyScalar,
     19       inv = _ref.inv,
     20       bignumber = _ref.bignumber,
     21       multiply = _ref.multiply,
     22       add = _ref.add;
     23 
     24   /**
     25    * @param {number[] | BigNumber[]} arr
     26    * @param {number} N
     27    * @param {number} prec
     28    * @param {'number' | 'BigNumber'} type
     29    */
     30   function main(arr, N) {
     31     var prec = arguments.length > 2 && arguments[2] !== undefined ? arguments[2] : config.epsilon;
     32     var type = arguments.length > 3 ? arguments[3] : undefined;
     33 
     34     if (type === 'number') {
     35       return diag(arr, prec);
     36     }
     37 
     38     if (type === 'BigNumber') {
     39       return diagBig(arr, prec);
     40     }
     41 
     42     throw TypeError('Unsupported data type: ' + type);
     43   } // diagonalization implementation for number (efficient)
     44 
     45 
     46   function diag(x, precision) {
     47     var N = x.length;
     48     var e0 = Math.abs(precision / N);
     49     var psi;
     50     var Sij = new Array(N); // Sij is Identity Matrix
     51 
     52     for (var i = 0; i < N; i++) {
     53       Sij[i] = createArray(N, 0);
     54       Sij[i][i] = 1.0;
     55     } // initial error
     56 
     57 
     58     var Vab = getAij(x);
     59 
     60     while (Math.abs(Vab[1]) >= Math.abs(e0)) {
     61       var _i = Vab[0][0];
     62       var j = Vab[0][1];
     63       psi = getTheta(x[_i][_i], x[j][j], x[_i][j]);
     64       x = x1(x, psi, _i, j);
     65       Sij = Sij1(Sij, psi, _i, j);
     66       Vab = getAij(x);
     67     }
     68 
     69     var Ei = createArray(N, 0); // eigenvalues
     70 
     71     for (var _i2 = 0; _i2 < N; _i2++) {
     72       Ei[_i2] = x[_i2][_i2];
     73     }
     74 
     75     return sorting((0, _object.clone)(Ei), (0, _object.clone)(Sij));
     76   } // diagonalization implementation for bigNumber
     77 
     78 
     79   function diagBig(x, precision) {
     80     var N = x.length;
     81     var e0 = abs(precision / N);
     82     var psi;
     83     var Sij = new Array(N); // Sij is Identity Matrix
     84 
     85     for (var i = 0; i < N; i++) {
     86       Sij[i] = createArray(N, 0);
     87       Sij[i][i] = 1.0;
     88     } // initial error
     89 
     90 
     91     var Vab = getAijBig(x);
     92 
     93     while (abs(Vab[1]) >= abs(e0)) {
     94       var _i3 = Vab[0][0];
     95       var j = Vab[0][1];
     96       psi = getThetaBig(x[_i3][_i3], x[j][j], x[_i3][j]);
     97       x = x1Big(x, psi, _i3, j);
     98       Sij = Sij1Big(Sij, psi, _i3, j);
     99       Vab = getAijBig(x);
    100     }
    101 
    102     var Ei = createArray(N, 0); // eigenvalues
    103 
    104     for (var _i4 = 0; _i4 < N; _i4++) {
    105       Ei[_i4] = x[_i4][_i4];
    106     } // return [clone(Ei), clone(Sij)]
    107 
    108 
    109     return sorting((0, _object.clone)(Ei), (0, _object.clone)(Sij));
    110   } // get angle
    111 
    112 
    113   function getTheta(aii, ajj, aij) {
    114     var denom = ajj - aii;
    115 
    116     if (Math.abs(denom) <= config.epsilon) {
    117       return Math.PI / 4.0;
    118     } else {
    119       return 0.5 * Math.atan(2.0 * aij / (ajj - aii));
    120     }
    121   } // get angle
    122 
    123 
    124   function getThetaBig(aii, ajj, aij) {
    125     var denom = subtract(ajj, aii);
    126 
    127     if (abs(denom) <= config.epsilon) {
    128       return bignumber(-1).acos().div(4);
    129     } else {
    130       return multiplyScalar(0.5, atan(multiply(2.0, aij, inv(denom))));
    131     }
    132   } // update eigvec
    133 
    134 
    135   function Sij1(Sij, theta, i, j) {
    136     var N = Sij.length;
    137     var c = Math.cos(theta);
    138     var s = Math.sin(theta);
    139     var Ski = createArray(N, 0);
    140     var Skj = createArray(N, 0);
    141 
    142     for (var k = 0; k < N; k++) {
    143       Ski[k] = c * Sij[k][i] - s * Sij[k][j];
    144       Skj[k] = s * Sij[k][i] + c * Sij[k][j];
    145     }
    146 
    147     for (var _k = 0; _k < N; _k++) {
    148       Sij[_k][i] = Ski[_k];
    149       Sij[_k][j] = Skj[_k];
    150     }
    151 
    152     return Sij;
    153   } // update eigvec for overlap
    154 
    155 
    156   function Sij1Big(Sij, theta, i, j) {
    157     var N = Sij.length;
    158     var c = cos(theta);
    159     var s = sin(theta);
    160     var Ski = createArray(N, bignumber(0));
    161     var Skj = createArray(N, bignumber(0));
    162 
    163     for (var k = 0; k < N; k++) {
    164       Ski[k] = subtract(multiplyScalar(c, Sij[k][i]), multiplyScalar(s, Sij[k][j]));
    165       Skj[k] = addScalar(multiplyScalar(s, Sij[k][i]), multiplyScalar(c, Sij[k][j]));
    166     }
    167 
    168     for (var _k2 = 0; _k2 < N; _k2++) {
    169       Sij[_k2][i] = Ski[_k2];
    170       Sij[_k2][j] = Skj[_k2];
    171     }
    172 
    173     return Sij;
    174   } // update matrix
    175 
    176 
    177   function x1Big(Hij, theta, i, j) {
    178     var N = Hij.length;
    179     var c = bignumber(cos(theta));
    180     var s = bignumber(sin(theta));
    181     var c2 = multiplyScalar(c, c);
    182     var s2 = multiplyScalar(s, s);
    183     var Aki = createArray(N, bignumber(0));
    184     var Akj = createArray(N, bignumber(0)); // 2cs Hij
    185 
    186     var csHij = multiply(bignumber(2), c, s, Hij[i][j]); //  Aii
    187 
    188     var Aii = addScalar(subtract(multiplyScalar(c2, Hij[i][i]), csHij), multiplyScalar(s2, Hij[j][j]));
    189     var Ajj = add(multiplyScalar(s2, Hij[i][i]), csHij, multiplyScalar(c2, Hij[j][j])); // 0  to i
    190 
    191     for (var k = 0; k < N; k++) {
    192       Aki[k] = subtract(multiplyScalar(c, Hij[i][k]), multiplyScalar(s, Hij[j][k]));
    193       Akj[k] = addScalar(multiplyScalar(s, Hij[i][k]), multiplyScalar(c, Hij[j][k]));
    194     } // Modify Hij
    195 
    196 
    197     Hij[i][i] = Aii;
    198     Hij[j][j] = Ajj;
    199     Hij[i][j] = bignumber(0);
    200     Hij[j][i] = bignumber(0); // 0  to i
    201 
    202     for (var _k3 = 0; _k3 < N; _k3++) {
    203       if (_k3 !== i && _k3 !== j) {
    204         Hij[i][_k3] = Aki[_k3];
    205         Hij[_k3][i] = Aki[_k3];
    206         Hij[j][_k3] = Akj[_k3];
    207         Hij[_k3][j] = Akj[_k3];
    208       }
    209     }
    210 
    211     return Hij;
    212   } // update matrix
    213 
    214 
    215   function x1(Hij, theta, i, j) {
    216     var N = Hij.length;
    217     var c = Math.cos(theta);
    218     var s = Math.sin(theta);
    219     var c2 = c * c;
    220     var s2 = s * s;
    221     var Aki = createArray(N, 0);
    222     var Akj = createArray(N, 0); //  Aii
    223 
    224     var Aii = c2 * Hij[i][i] - 2 * c * s * Hij[i][j] + s2 * Hij[j][j];
    225     var Ajj = s2 * Hij[i][i] + 2 * c * s * Hij[i][j] + c2 * Hij[j][j]; // 0  to i
    226 
    227     for (var k = 0; k < N; k++) {
    228       Aki[k] = c * Hij[i][k] - s * Hij[j][k];
    229       Akj[k] = s * Hij[i][k] + c * Hij[j][k];
    230     } // Modify Hij
    231 
    232 
    233     Hij[i][i] = Aii;
    234     Hij[j][j] = Ajj;
    235     Hij[i][j] = 0;
    236     Hij[j][i] = 0; // 0  to i
    237 
    238     for (var _k4 = 0; _k4 < N; _k4++) {
    239       if (_k4 !== i && _k4 !== j) {
    240         Hij[i][_k4] = Aki[_k4];
    241         Hij[_k4][i] = Aki[_k4];
    242         Hij[j][_k4] = Akj[_k4];
    243         Hij[_k4][j] = Akj[_k4];
    244       }
    245     }
    246 
    247     return Hij;
    248   } // get max off-diagonal value from Upper Diagonal
    249 
    250 
    251   function getAij(Mij) {
    252     var N = Mij.length;
    253     var maxMij = 0;
    254     var maxIJ = [0, 1];
    255 
    256     for (var i = 0; i < N; i++) {
    257       for (var j = i + 1; j < N; j++) {
    258         if (Math.abs(maxMij) < Math.abs(Mij[i][j])) {
    259           maxMij = Math.abs(Mij[i][j]);
    260           maxIJ = [i, j];
    261         }
    262       }
    263     }
    264 
    265     return [maxIJ, maxMij];
    266   } // get max off-diagonal value from Upper Diagonal
    267 
    268 
    269   function getAijBig(Mij) {
    270     var N = Mij.length;
    271     var maxMij = 0;
    272     var maxIJ = [0, 1];
    273 
    274     for (var i = 0; i < N; i++) {
    275       for (var j = i + 1; j < N; j++) {
    276         if (abs(maxMij) < abs(Mij[i][j])) {
    277           maxMij = abs(Mij[i][j]);
    278           maxIJ = [i, j];
    279         }
    280       }
    281     }
    282 
    283     return [maxIJ, maxMij];
    284   } // sort results
    285 
    286 
    287   function sorting(E, S) {
    288     var N = E.length;
    289     var values = Array(N);
    290     var vectors = Array(N);
    291 
    292     for (var k = 0; k < N; k++) {
    293       vectors[k] = Array(N);
    294     }
    295 
    296     for (var i = 0; i < N; i++) {
    297       var minID = 0;
    298       var minE = E[0];
    299 
    300       for (var j = 0; j < E.length; j++) {
    301         if (abs(E[j]) < abs(minE)) {
    302           minID = j;
    303           minE = E[minID];
    304         }
    305       }
    306 
    307       values[i] = E.splice(minID, 1)[0];
    308 
    309       for (var _k5 = 0; _k5 < N; _k5++) {
    310         vectors[_k5][i] = S[_k5][minID];
    311 
    312         S[_k5].splice(minID, 1);
    313       }
    314     }
    315 
    316     return {
    317       values: values,
    318       vectors: vectors
    319     };
    320   }
    321   /**
    322    * Create an array of a certain size and fill all items with an initial value
    323    * @param {number} size
    324    * @param {number} value
    325    * @return {number[]}
    326    */
    327 
    328 
    329   function createArray(size, value) {
    330     // TODO: as soon as all browsers support Array.fill, use that instead (IE doesn't support it)
    331     var array = new Array(size);
    332 
    333     for (var i = 0; i < size; i++) {
    334       array[i] = value;
    335     }
    336 
    337     return array;
    338   }
    339 
    340   return main;
    341 }