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 }