vendor: 48,793 bytes, lines 1-1387
1"use strict"; 2 3var numeric = (typeof exports === "undefined")?(function numeric() {}):(exports); 4if(typeof global !== "undefined") { global.numeric = numeric; } 5 6numeric.version = "1.2.6"; 7 8// 1. Utility functions 9numeric.bench = function bench (f,interval) { 10 var t1,t2,n,i; 11 if(typeof interval === "undefined") { interval = 15; } 12 n = 0.5; 13 t1 = new Date(); 14 while(1) { 15 n*=2; 16 for(i=n;i>3;i-=4) { f(); f(); f(); f(); } 17 while(i>0) { f(); i--; } 18 t2 = new Date(); 19 if(t2-t1 > interval) break; 20 } 21 for(i=n;i>3;i-=4) { f(); f(); f(); f(); } 22 while(i>0) { f(); i--; } 23 t2 = new Date(); 24 return 1000*(3*n-1)/(t2-t1); 25} 26 27numeric._myIndexOf = (function _myIndexOf(w) { 28 var n = this.length,k; 29 for(k=0;k<n;++k) if(this[k]===w) return k; 30 return -1; 31}); 32numeric.myIndexOf = (Array.prototype.indexOf)?Array.prototype.indexOf:numeric._myIndexOf; 33 34numeric.Function = Function; 35numeric.precision = 4; 36numeric.largeArray = 50; 37 38numeric.prettyPrint = function prettyPrint(x) { 39 function fmtnum(x) { 40 if(x === 0) { return '0'; } 41 if(isNaN(x)) { return 'NaN'; } 42 if(x<0) { return '-'+fmtnum(-x); } 43 if(isFinite(x)) { 44 var scale = Math.floor(Math.log(x) / Math.log(10)); 45 var normalized = x / Math.pow(10,scale); 46 var basic = normalized.toPrecision(numeric.precision); 47 if(parseFloat(basic) === 10) { scale++; normalized = 1; basic = normalized.toPrecision(numeric.precision); } 48 return parseFloat(basic).toString()+'e'+scale.toString(); 49 } 50 return 'Infinity'; 51 } 52 var ret = []; 53 function foo(x) { 54 var k; 55 if(typeof x === "undefined") { ret.push(Array(numeric.precision+8).join(' ')); return false; } 56 if(typeof x === "string") { ret.push('"'+x+'"'); return false; } 57 if(typeof x === "boolean") { ret.push(x.toString()); return false; } 58 if(typeof x === "number") { 59 var a = fmtnum(x); 60 var b = x.toPrecision(numeric.precision); 61 var c = parseFloat(x.toString()).toString(); 62 var d = [a,b,c,parseFloat(b).toString(),parseFloat(c).toString()]; 63 for(k=1;k<d.length;k++) { if(d[k].length < a.length) a = d[k]; } 64 ret.push(Array(numeric.precision+8-a.length).join(' ')+a); 65 return false; 66 } 67 if(x === null) { ret.push("null"); return false; } 68 if(typeof x === "function") { 69 ret.push(x.toString()); 70 var flag = false; 71 for(k in x) { if(x.hasOwnProperty(k)) { 72 if(flag) ret.push(',\n'); 73 else ret.push('\n{'); 74 flag = true; 75 ret.push(k); 76 ret.push(': \n'); 77 foo(x[k]); 78 } } 79 if(flag) ret.push('}\n'); 80 return true; 81 } 82 if(x instanceof Array) { 83 if(x.length > numeric.largeArray) { ret.push('...Large Array...'); return true; } 84 var flag = false; 85 ret.push('['); 86 for(k=0;k<x.length;k++) { if(k>0) { ret.push(','); if(flag) ret.push('\n '); } flag = foo(x[k]); } 87 ret.push(']'); 88 return true; 89 } 90 ret.push('{'); 91 var flag = false; 92 for(k in x) { if(x.hasOwnProperty(k)) { if(flag) ret.push(',\n'); flag = true; ret.push(k); ret.push(': \n'); foo(x[k]); } } 93 ret.push('}'); 94 return true; 95 } 96 foo(x); 97 return ret.join(''); 98} 99 100numeric.parseDate = function parseDate(d) { 101 function foo(d) { 102 if(typeof d === 'string') { return Date.parse(d.replace(/-/g,'/')); } 103 if(!(d instanceof Array)) { throw new Error("parseDate: parameter must be arrays of strings"); } 104 var ret = [],k; 105 for(k=0;k<d.length;k++) { ret[k] = foo(d[k]); } 106 return ret; 107 } 108 return foo(d); 109} 110 111numeric.parseFloat = function parseFloat_(d) { 112 function foo(d) { 113 if(typeof d === 'string') { return parseFloat(d); } 114 if(!(d instanceof Array)) { throw new Error("parseFloat: parameter must be arrays of strings"); } 115 var ret = [],k; 116 for(k=0;k<d.length;k++) { ret[k] = foo(d[k]); } 117 return ret; 118 } 119 return foo(d); 120} 121 122numeric.parseCSV = function parseCSV(t) { 123 var foo = t.split('\n'); 124 var j,k; 125 var ret = []; 126 var pat = /(([^'",]*)|('[^']*')|("[^"]*")),/g; 127 var patnum = /^\s*(([+-]?[0-9]+(\.[0-9]*)?(e[+-]?[0-9]+)?)|([+-]?[0-9]*(\.[0-9]+)?(e[+-]?[0-9]+)?))\s*$/; 128 var stripper = function(n) { return n.substr(0,n.length-1); } 129 var count = 0; 130 for(k=0;k<foo.length;k++) { 131 var bar = (foo[k]+",").match(pat),baz; 132 if(bar.length>0) { 133 ret[count] = []; 134 for(j=0;j<bar.length;j++) { 135 baz = stripper(bar[j]); 136 if(patnum.test(baz)) { ret[count][j] = parseFloat(baz); } 137 else ret[count][j] = baz; 138 } 139 count++; 140 } 141 } 142 return ret; 143} 144 145numeric.toCSV = function toCSV(A) { 146 var s = numeric.dim(A); 147 var i,j,m,n,row,ret; 148 m = s[0]; 149 n = s[1]; 150 ret = []; 151 for(i=0;i<m;i++) { 152 row = []; 153 for(j=0;j<m;j++) { row[j] = A[i][j].toString(); } 154 ret[i] = row.join(', '); 155 } 156 return ret.join('\n')+'\n'; 157} 158 159numeric.getURL = function getURL(url) { 160 var client = new XMLHttpRequest(); 161 client.open("GET",url,false); 162 client.send(); 163 return client; 164} 165 166numeric.imageURL = function imageURL(img) { 167 function base64(A) { 168 var n = A.length, i,x,y,z,p,q,r,s; 169 var key = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/="; 170 var ret = ""; 171 for(i=0;i<n;i+=3) { 172 x = A[i]; 173 y = A[i+1]; 174 z = A[i+2]; 175 p = x >> 2; 176 q = ((x & 3) << 4) + (y >> 4); 177 r = ((y & 15) << 2) + (z >> 6); 178 s = z & 63; 179 if(i+1>=n) { r = s = 64; } 180 else if(i+2>=n) { s = 64; } 181 ret += key.charAt(p) + key.charAt(q) + key.charAt(r) + key.charAt(s); 182 } 183 return ret; 184 } 185 function crc32Array (a,from,to) { 186 if(typeof from === "undefined") { from = 0; } 187 if(typeof to === "undefined") { to = a.length; } 188 var table = [0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3, 189 0x0EDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91, 190 0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7, 191 0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5, 192 0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B, 193 0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59, 194 0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F, 195 0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D, 196 0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433, 197 0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01, 198 0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457, 199 0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65, 200 0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB, 201 0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9, 202 0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F, 203 0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD, 204 0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683, 205 0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1, 206 0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7, 207 0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5, 208 0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B, 209 0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79, 210 0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F, 211 0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D, 212 0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713, 213 0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21, 214 0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777, 215 0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45, 216 0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB, 217 0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9, 218 0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF, 219 0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D]; 220 221 var crc = -1, y = 0, n = a.length,i; 222 223 for (i = from; i < to; i++) { 224 y = (crc ^ a[i]) & 0xFF; 225 crc = (crc >>> 8) ^ table[y]; 226 } 227 228 return crc ^ (-1); 229 } 230 231 var h = img[0].length, w = img[0][0].length, s1, s2, next,k,length,a,b,i,j,adler32,crc32; 232 var stream = [ 233 137, 80, 78, 71, 13, 10, 26, 10, // 0: PNG signature 234 0,0,0,13, // 8: IHDR Chunk length 235 73, 72, 68, 82, // 12: "IHDR" 236 (w >> 24) & 255, (w >> 16) & 255, (w >> 8) & 255, w&255, // 16: Width 237 (h >> 24) & 255, (h >> 16) & 255, (h >> 8) & 255, h&255, // 20: Height 238 8, // 24: bit depth 239 2, // 25: RGB 240 0, // 26: deflate 241 0, // 27: no filter 242 0, // 28: no interlace 243 -1,-2,-3,-4, // 29: CRC 244 -5,-6,-7,-8, // 33: IDAT Chunk length 245 73, 68, 65, 84, // 37: "IDAT" 246 // RFC 1950 header starts here 247 8, // 41: RFC1950 CMF 248 29 // 42: RFC1950 FLG 249 ]; 250 crc32 = crc32Array(stream,12,29); 251 stream[29] = (crc32>>24)&255; 252 stream[30] = (crc32>>16)&255; 253 stream[31] = (crc32>>8)&255; 254 stream[32] = (crc32)&255; 255 s1 = 1; 256 s2 = 0; 257 for(i=0;i<h;i++) { 258 if(i<h-1) { stream.push(0); } 259 else { stream.push(1); } 260 a = (3*w+1+(i===0))&255; b = ((3*w+1+(i===0))>>8)&255; 261 stream.push(a); stream.push(b); 262 stream.push((~a)&255); stream.push((~b)&255); 263 if(i===0) stream.push(0); 264 for(j=0;j<w;j++) { 265 for(k=0;k<3;k++) { 266 a = img[k][i][j]; 267 if(a>255) a = 255; 268 else if(a<0) a=0; 269 else a = Math.round(a); 270 s1 = (s1 + a )%65521; 271 s2 = (s2 + s1)%65521; 272 stream.push(a); 273 } 274 } 275 stream.push(0); 276 } 277 adler32 = (s2<<16)+s1; 278 stream.push((adler32>>24)&255); 279 stream.push((adler32>>16)&255); 280 stream.push((adler32>>8)&255); 281 stream.push((adler32)&255); 282 length = stream.length - 41; 283 stream[33] = (length>>24)&255; 284 stream[34] = (length>>16)&255; 285 stream[35] = (length>>8)&255; 286 stream[36] = (length)&255; 287 crc32 = crc32Array(stream,37); 288 stream.push((crc32>>24)&255); 289 stream.push((crc32>>16)&255); 290 stream.push((crc32>>8)&255); 291 stream.push((crc32)&255); 292 stream.push(0); 293 stream.push(0); 294 stream.push(0); 295 stream.push(0); 296// a = stream.length; 297 stream.push(73); // I 298 stream.push(69); // E 299 stream.push(78); // N 300 stream.push(68); // D 301 stream.push(174); // CRC1 302 stream.push(66); // CRC2 303 stream.push(96); // CRC3 304 stream.push(130); // CRC4 305 return 'data:image/png;base64,'+base64(stream); 306} 307 308// 2. Linear algebra with Arrays. 309numeric._dim = function _dim(x) { 310 var ret = []; 311 while(typeof x === "object") { ret.push(x.length); x = x[0]; } 312 return ret; 313} 314 315numeric.dim = function dim(x) { 316 var y,z; 317 if(typeof x === "object") { 318 y = x[0]; 319 if(typeof y === "object") { 320 z = y[0]; 321 if(typeof z === "object") { 322 return numeric._dim(x); 323 } 324 return [x.length,y.length]; 325 } 326 return [x.length]; 327 } 328 return []; 329} 330 331numeric.mapreduce = function mapreduce(body,init) { 332 return Function('x','accum','_s','_k', 333 'if(typeof accum === "undefined") accum = '+init+';\n'+ 334 'if(typeof x === "number") { var xi = x; '+body+'; return accum; }\n'+ 335 'if(typeof _s === "undefined") _s = numeric.dim(x);\n'+ 336 'if(typeof _k === "undefined") _k = 0;\n'+ 337 'var _n = _s[_k];\n'+ 338 'var i,xi;\n'+ 339 'if(_k < _s.length-1) {\n'+ 340 ' for(i=_n-1;i>=0;i--) {\n'+ 341 ' accum = arguments.callee(x[i],accum,_s,_k+1);\n'+ 342 ' }'+ 343 ' return accum;\n'+ 344 '}\n'+ 345 'for(i=_n-1;i>=1;i-=2) { \n'+ 346 ' xi = x[i];\n'+ 347 ' '+body+';\n'+ 348 ' xi = x[i-1];\n'+ 349 ' '+body+';\n'+ 350 '}\n'+ 351 'if(i === 0) {\n'+ 352 ' xi = x[i];\n'+ 353 ' '+body+'\n'+ 354 '}\n'+ 355 'return accum;' 356 ); 357} 358numeric.mapreduce2 = function mapreduce2(body,setup) { 359 return Function('x', 360 'var n = x.length;\n'+ 361 'var i,xi;\n'+setup+';\n'+ 362 'for(i=n-1;i!==-1;--i) { \n'+ 363 ' xi = x[i];\n'+ 364 ' '+body+';\n'+ 365 '}\n'+ 366 'return accum;' 367 ); 368} 369 370 371numeric.same = function same(x,y) { 372 var i,n; 373 if(!(x instanceof Array) || !(y instanceof Array)) { return false; } 374 n = x.length; 375 if(n !== y.length) { return false; } 376 for(i=0;i<n;i++) { 377 if(x[i] === y[i]) { continue; } 378 if(typeof x[i] === "object") { if(!same(x[i],y[i])) return false; } 379 else { return false; } 380 } 381 return true; 382} 383 384numeric.rep = function rep(s,v,k) { 385 if(typeof k === "undefined") { k=0; } 386 var n = s[k], ret = Array(n), i; 387 if(k === s.length-1) { 388 for(i=n-2;i>=0;i-=2) { ret[i+1] = v; ret[i] = v; } 389 if(i===-1) { ret[0] = v; } 390 return ret; 391 } 392 for(i=n-1;i>=0;i--) { ret[i] = numeric.rep(s,v,k+1); } 393 return ret; 394} 395 396 397numeric.dotMMsmall = function dotMMsmall(x,y) { 398 var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0; 399 p = x.length; q = y.length; r = y[0].length; 400 ret = Array(p); 401 for(i=p-1;i>=0;i--) { 402 foo = Array(r); 403 bar = x[i]; 404 for(k=r-1;k>=0;k--) { 405 woo = bar[q-1]*y[q-1][k]; 406 for(j=q-2;j>=1;j-=2) { 407 i0 = j-1; 408 woo += bar[j]*y[j][k] + bar[i0]*y[i0][k]; 409 } 410 if(j===0) { woo += bar[0]*y[0][k]; } 411 foo[k] = woo; 412 } 413 ret[i] = foo; 414 } 415 return ret; 416} 417numeric._getCol = function _getCol(A,j,x) { 418 var n = A.length, i; 419 for(i=n-1;i>0;--i) { 420 x[i] = A[i][j]; 421 --i; 422 x[i] = A[i][j]; 423 } 424 if(i===0) x[0] = A[0][j]; 425} 426numeric.dotMMbig = function dotMMbig(x,y){ 427 var gc = numeric._getCol, p = y.length, v = Array(p); 428 var m = x.length, n = y[0].length, A = new Array(m), xj; 429 var VV = numeric.dotVV; 430 var i,j,k,z; 431 --p; 432 --m; 433 for(i=m;i!==-1;--i) A[i] = Array(n); 434 --n; 435 for(i=n;i!==-1;--i) { 436 gc(y,i,v); 437 for(j=m;j!==-1;--j) { 438 z=0; 439 xj = x[j]; 440 A[j][i] = VV(xj,v); 441 } 442 } 443 return A; 444} 445 446numeric.dotMV = function dotMV(x,y) { 447 var p = x.length, q = y.length,i; 448 var ret = Array(p), dotVV = numeric.dotVV; 449 for(i=p-1;i>=0;i--) { ret[i] = dotVV(x[i],y); } 450 return ret; 451} 452 453numeric.dotVM = function dotVM(x,y) { 454 var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0,s1,s2,s3,baz,accum; 455 p = x.length; q = y[0].length; 456 ret = Array(q); 457 for(k=q-1;k>=0;k--) { 458 woo = x[p-1]*y[p-1][k]; 459 for(j=p-2;j>=1;j-=2) { 460 i0 = j-1; 461 woo += x[j]*y[j][k] + x[i0]*y[i0][k]; 462 } 463 if(j===0) { woo += x[0]*y[0][k]; } 464 ret[k] = woo; 465 } 466 return ret; 467} 468 469numeric.dotVV = function dotVV(x,y) { 470 var i,n=x.length,i1,ret = x[n-1]*y[n-1]; 471 for(i=n-2;i>=1;i-=2) { 472 i1 = i-1; 473 ret += x[i]*y[i] + x[i1]*y[i1]; 474 } 475 if(i===0) { ret += x[0]*y[0]; } 476 return ret; 477} 478 479numeric.dot = function dot(x,y) { 480 var d = numeric.dim; 481 switch(d(x).length*1000+d(y).length) { 482 case 2002: 483 if(y.length < 10) return numeric.dotMMsmall(x,y); 484 else return numeric.dotMMbig(x,y); 485 case 2001: return numeric.dotMV(x,y); 486 case 1002: return numeric.dotVM(x,y); 487 case 1001: return numeric.dotVV(x,y); 488 case 1000: return numeric.mulVS(x,y); 489 case 1: return numeric.mulSV(x,y); 490 case 0: return x*y; 491 default: throw new Error('numeric.dot only works on vectors and matrices'); 492 } 493} 494 495numeric.diag = function diag(d) { 496 var i,i1,j,n = d.length, A = Array(n), Ai; 497 for(i=n-1;i>=0;i--) { 498 Ai = Array(n); 499 i1 = i+2; 500 for(j=n-1;j>=i1;j-=2) { 501 Ai[j] = 0; 502 Ai[j-1] = 0; 503 } 504 if(j>i) { Ai[j] = 0; } 505 Ai[i] = d[i]; 506 for(j=i-1;j>=1;j-=2) { 507 Ai[j] = 0; 508 Ai[j-1] = 0; 509 } 510 if(j===0) { Ai[0] = 0; } 511 A[i] = Ai; 512 } 513 return A; 514} 515numeric.getDiag = function(A) { 516 var n = Math.min(A.length,A[0].length),i,ret = Array(n); 517 for(i=n-1;i>=1;--i) { 518 ret[i] = A[i][i]; 519 --i; 520 ret[i] = A[i][i]; 521 } 522 if(i===0) { 523 ret[0] = A[0][0]; 524 } 525 return ret; 526} 527 528numeric.identity = function identity(n) { return numeric.diag(numeric.rep([n],1)); } 529numeric.pointwise = function pointwise(params,body,setup) { 530 if(typeof setup === "undefined") { setup = ""; } 531 var fun = []; 532 var k; 533 var avec = /\[i\]$/,p,thevec = ''; 534 var haveret = false; 535 for(k=0;k<params.length;k++) { 536 if(avec.test(params[k])) { 537 p = params[k].substring(0,params[k].length-3); 538 thevec = p; 539 } else { p = params[k]; } 540 if(p==='ret') haveret = true; 541 fun.push(p); 542 } 543 fun[params.length] = '_s'; 544 fun[params.length+1] = '_k'; 545 fun[params.length+2] = ( 546 'if(typeof _s === "undefined") _s = numeric.dim('+thevec+');\n'+ 547 'if(typeof _k === "undefined") _k = 0;\n'+ 548 'var _n = _s[_k];\n'+ 549 'var i'+(haveret?'':', ret = Array(_n)')+';\n'+ 550 'if(_k < _s.length-1) {\n'+ 551 ' for(i=_n-1;i>=0;i--) ret[i] = arguments.callee('+params.join(',')+',_s,_k+1);\n'+ 552 ' return ret;\n'+ 553 '}\n'+ 554 setup+'\n'+ 555 'for(i=_n-1;i!==-1;--i) {\n'+ 556 ' '+body+'\n'+ 557 '}\n'+ 558 'return ret;' 559 ); 560 return Function.apply(null,fun); 561} 562numeric.pointwise2 = function pointwise2(params,body,setup) { 563 if(typeof setup === "undefined") { setup = ""; } 564 var fun = []; 565 var k; 566 var avec = /\[i\]$/,p,thevec = ''; 567 var haveret = false; 568 for(k=0;k<params.length;k++) { 569 if(avec.test(params[k])) { 570 p = params[k].substring(0,params[k].length-3); 571 thevec = p; 572 } else { p = params[k]; } 573 if(p==='ret') haveret = true; 574 fun.push(p); 575 } 576 fun[params.length] = ( 577 'var _n = '+thevec+'.length;\n'+ 578 'var i'+(haveret?'':', ret = Array(_n)')+';\n'+ 579 setup+'\n'+ 580 'for(i=_n-1;i!==-1;--i) {\n'+ 581 body+'\n'+ 582 '}\n'+ 583 'return ret;' 584 ); 585 return Function.apply(null,fun); 586} 587numeric._biforeach = (function _biforeach(x,y,s,k,f) { 588 if(k === s.length-1) { f(x,y); return; } 589 var i,n=s[k]; 590 for(i=n-1;i>=0;i--) { _biforeach(typeof x==="object"?x[i]:x,typeof y==="object"?y[i]:y,s,k+1,f); } 591}); 592numeric._biforeach2 = (function _biforeach2(x,y,s,k,f) { 593 if(k === s.length-1) { return f(x,y); } 594 var i,n=s[k],ret = Array(n); 595 for(i=n-1;i>=0;--i) { ret[i] = _biforeach2(typeof x==="object"?x[i]:x,typeof y==="object"?y[i]:y,s,k+1,f); } 596 return ret; 597}); 598numeric._foreach = (function _foreach(x,s,k,f) { 599 if(k === s.length-1) { f(x); return; } 600 var i,n=s[k]; 601 for(i=n-1;i>=0;i--) { _foreach(x[i],s,k+1,f); } 602}); 603numeric._foreach2 = (function _foreach2(x,s,k,f) { 604 if(k === s.length-1) { return f(x); } 605 var i,n=s[k], ret = Array(n); 606 for(i=n-1;i>=0;i--) { ret[i] = _foreach2(x[i],s,k+1,f); } 607 return ret; 608}); 609 610/*numeric.anyV = numeric.mapreduce('if(xi) return true;','false'); 611numeric.allV = numeric.mapreduce('if(!xi) return false;','true'); 612numeric.any = function(x) { if(typeof x.length === "undefined") return x; return numeric.anyV(x); } 613numeric.all = function(x) { if(typeof x.length === "undefined") return x; return numeric.allV(x); }*/ 614 615numeric.ops2 = { 616 add: '+', 617 sub: '-', 618 mul: '*', 619 div: '/', 620 mod: '%', 621 and: '&&', 622 or: '||', 623 eq: '===', 624 neq: '!==', 625 lt: '<', 626 gt: '>', 627 leq: '<=', 628 geq: '>=', 629 band: '&', 630 bor: '|', 631 bxor: '^', 632 lshift: '<<', 633 rshift: '>>', 634 rrshift: '>>>' 635}; 636numeric.opseq = { 637 addeq: '+=', 638 subeq: '-=', 639 muleq: '*=', 640 diveq: '/=', 641 modeq: '%=', 642 lshifteq: '<<=', 643 rshifteq: '>>=', 644 rrshifteq: '>>>=', 645 bandeq: '&=', 646 boreq: '|=', 647 bxoreq: '^=' 648}; 649numeric.mathfuns = ['abs','acos','asin','atan','ceil','cos', 650 'exp','floor','log','round','sin','sqrt','tan', 651 'isNaN','isFinite']; 652numeric.mathfuns2 = ['atan2','pow','max','min']; 653numeric.ops1 = { 654 neg: '-', 655 not: '!', 656 bnot: '~', 657 clone: '' 658}; 659numeric.mapreducers = { 660 any: ['if(xi) return true;','var accum = false;'], 661 all: ['if(!xi) return false;','var accum = true;'], 662 sum: ['accum += xi;','var accum = 0;'], 663 prod: ['accum *= xi;','var accum = 1;'], 664 norm2Squared: ['accum += xi*xi;','var accum = 0;'], 665 norminf: ['accum = max(accum,abs(xi));','var accum = 0, max = Math.max, abs = Math.abs;'], 666 norm1: ['accum += abs(xi)','var accum = 0, abs = Math.abs;'], 667 sup: ['accum = max(accum,xi);','var accum = -Infinity, max = Math.max;'], 668 inf: ['accum = min(accum,xi);','var accum = Infinity, min = Math.min;'] 669}; 670 671(function () { 672 var i,o; 673 for(i=0;i<numeric.mathfuns2.length;++i) { 674 o = numeric.mathfuns2[i]; 675 numeric.ops2[o] = o; 676 } 677 for(i in numeric.ops2) { 678 if(numeric.ops2.hasOwnProperty(i)) { 679 o = numeric.ops2[i]; 680 var code, codeeq, setup = ''; 681 if(numeric.myIndexOf.call(numeric.mathfuns2,i)!==-1) { 682 setup = 'var '+o+' = Math.'+o+';\n'; 683 code = function(r,x,y) { return r+' = '+o+'('+x+','+y+')'; }; 684 codeeq = function(x,y) { return x+' = '+o+'('+x+','+y+')'; }; 685 } else { 686 code = function(r,x,y) { return r+' = '+x+' '+o+' '+y; }; 687 if(numeric.opseq.hasOwnProperty(i+'eq')) { 688 codeeq = function(x,y) { return x+' '+o+'= '+y; }; 689 } else { 690 codeeq = function(x,y) { return x+' = '+x+' '+o+' '+y; }; 691 } 692 } 693 numeric[i+'VV'] = numeric.pointwise2(['x[i]','y[i]'],code('ret[i]','x[i]','y[i]'),setup); 694 numeric[i+'SV'] = numeric.pointwise2(['x','y[i]'],code('ret[i]','x','y[i]'),setup); 695 numeric[i+'VS'] = numeric.pointwise2(['x[i]','y'],code('ret[i]','x[i]','y'),setup); 696 numeric[i] = Function( 697 'var n = arguments.length, i, x = arguments[0], y;\n'+ 698 'var VV = numeric.'+i+'VV, VS = numeric.'+i+'VS, SV = numeric.'+i+'SV;\n'+ 699 'var dim = numeric.dim;\n'+ 700 'for(i=1;i!==n;++i) { \n'+ 701 ' y = arguments[i];\n'+ 702 ' if(typeof x === "object") {\n'+ 703 ' if(typeof y === "object") x = numeric._biforeach2(x,y,dim(x),0,VV);\n'+ 704 ' else x = numeric._biforeach2(x,y,dim(x),0,VS);\n'+ 705 ' } else if(typeof y === "object") x = numeric._biforeach2(x,y,dim(y),0,SV);\n'+ 706 ' else '+codeeq('x','y')+'\n'+ 707 '}\nreturn x;\n'); 708 numeric[o] = numeric[i]; 709 numeric[i+'eqV'] = numeric.pointwise2(['ret[i]','x[i]'], codeeq('ret[i]','x[i]'),setup); 710 numeric[i+'eqS'] = numeric.pointwise2(['ret[i]','x'], codeeq('ret[i]','x'),setup); 711 numeric[i+'eq'] = Function( 712 'var n = arguments.length, i, x = arguments[0], y;\n'+ 713 'var V = numeric.'+i+'eqV, S = numeric.'+i+'eqS\n'+ 714 'var s = numeric.dim(x);\n'+ 715 'for(i=1;i!==n;++i) { \n'+ 716 ' y = arguments[i];\n'+ 717 ' if(typeof y === "object") numeric._biforeach(x,y,s,0,V);\n'+ 718 ' else numeric._biforeach(x,y,s,0,S);\n'+ 719 '}\nreturn x;\n'); 720 } 721 } 722 for(i=0;i<numeric.mathfuns2.length;++i) { 723 o = numeric.mathfuns2[i]; 724 delete numeric.ops2[o]; 725 } 726 for(i=0;i<numeric.mathfuns.length;++i) { 727 o = numeric.mathfuns[i]; 728 numeric.ops1[o] = o; 729 } 730 for(i in numeric.ops1) { 731 if(numeric.ops1.hasOwnProperty(i)) { 732 setup = ''; 733 o = numeric.ops1[i]; 734 if(numeric.myIndexOf.call(numeric.mathfuns,i)!==-1) { 735 if(Math.hasOwnProperty(o)) setup = 'var '+o+' = Math.'+o+';\n'; 736 } 737 numeric[i+'eqV'] = numeric.pointwise2(['ret[i]'],'ret[i] = '+o+'(ret[i]);',setup); 738 numeric[i+'eq'] = Function('x', 739 'if(typeof x !== "object") return '+o+'x\n'+ 740 'var i;\n'+ 741 'var V = numeric.'+i+'eqV;\n'+ 742 'var s = numeric.dim(x);\n'+ 743 'numeric._foreach(x,s,0,V);\n'+ 744 'return x;\n'); 745 numeric[i+'V'] = numeric.pointwise2(['x[i]'],'ret[i] = '+o+'(x[i]);',setup); 746 numeric[i] = Function('x', 747 'if(typeof x !== "object") return '+o+'(x)\n'+ 748 'var i;\n'+ 749 'var V = numeric.'+i+'V;\n'+ 750 'var s = numeric.dim(x);\n'+ 751 'return numeric._foreach2(x,s,0,V);\n'); 752 } 753 } 754 for(i=0;i<numeric.mathfuns.length;++i) { 755 o = numeric.mathfuns[i]; 756 delete numeric.ops1[o]; 757 } 758 for(i in numeric.mapreducers) { 759 if(numeric.mapreducers.hasOwnProperty(i)) { 760 o = numeric.mapreducers[i]; 761 numeric[i+'V'] = numeric.mapreduce2(o[0],o[1]); 762 numeric[i] = Function('x','s','k', 763 o[1]+ 764 'if(typeof x !== "object") {'+ 765 ' xi = x;\n'+ 766 o[0]+';\n'+ 767 ' return accum;\n'+ 768 '}'+ 769 'if(typeof s === "undefined") s = numeric.dim(x);\n'+ 770 'if(typeof k === "undefined") k = 0;\n'+ 771 'if(k === s.length-1) return numeric.'+i+'V(x);\n'+ 772 'var xi;\n'+ 773 'var n = x.length, i;\n'+ 774 'for(i=n-1;i!==-1;--i) {\n'+ 775 ' xi = arguments.callee(x[i]);\n'+ 776 o[0]+';\n'+ 777 '}\n'+ 778 'return accum;\n'); 779 } 780 } 781}()); 782 783numeric.truncVV = numeric.pointwise(['x[i]','y[i]'],'ret[i] = round(x[i]/y[i])*y[i];','var round = Math.round;'); 784numeric.truncVS = numeric.pointwise(['x[i]','y'],'ret[i] = round(x[i]/y)*y;','var round = Math.round;'); 785numeric.truncSV = numeric.pointwise(['x','y[i]'],'ret[i] = round(x/y[i])*y[i];','var round = Math.round;'); 786numeric.trunc = function trunc(x,y) { 787 if(typeof x === "object") { 788 if(typeof y === "object") return numeric.truncVV(x,y); 789 return numeric.truncVS(x,y); 790 } 791 if (typeof y === "object") return numeric.truncSV(x,y); 792 return Math.round(x/y)*y; 793} 794 795numeric.inv = function inv(x) { 796 var s = numeric.dim(x), abs = Math.abs, m = s[0], n = s[1]; 797 var A = numeric.clone(x), Ai, Aj; 798 var I = numeric.identity(m), Ii, Ij; 799 var i,j,k,x; 800 for(j=0;j<n;++j) { 801 var i0 = -1; 802 var v0 = -1; 803 for(i=j;i!==m;++i) { k = abs(A[i][j]); if(k>v0) { i0 = i; v0 = k; } } 804 Aj = A[i0]; A[i0] = A[j]; A[j] = Aj; 805 Ij = I[i0]; I[i0] = I[j]; I[j] = Ij; 806 x = Aj[j]; 807 for(k=j;k!==n;++k) Aj[k] /= x; 808 for(k=n-1;k!==-1;--k) Ij[k] /= x; 809 for(i=m-1;i!==-1;--i) { 810 if(i!==j) { 811 Ai = A[i]; 812 Ii = I[i]; 813 x = Ai[j]; 814 for(k=j+1;k!==n;++k) Ai[k] -= Aj[k]*x; 815 for(k=n-1;k>0;--k) { Ii[k] -= Ij[k]*x; --k; Ii[k] -= Ij[k]*x; } 816 if(k===0) Ii[0] -= Ij[0]*x; 817 } 818 } 819 } 820 return I; 821} 822 823numeric.det = function det(x) { 824 var s = numeric.dim(x); 825 if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: det() only works on square matrices'); } 826 var n = s[0], ret = 1,i,j,k,A = numeric.clone(x),Aj,Ai,alpha,temp,k1,k2,k3; 827 for(j=0;j<n-1;j++) { 828 k=j; 829 for(i=j+1;i<n;i++) { if(Math.abs(A[i][j]) > Math.abs(A[k][j])) { k = i; } } 830 if(k !== j) { 831 temp = A[k]; A[k] = A[j]; A[j] = temp; 832 ret *= -1; 833 } 834 Aj = A[j]; 835 for(i=j+1;i<n;i++) { 836 Ai = A[i]; 837 alpha = Ai[j]/Aj[j]; 838 for(k=j+1;k<n-1;k+=2) { 839 k1 = k+1; 840 Ai[k] -= Aj[k]*alpha; 841 Ai[k1] -= Aj[k1]*alpha; 842 } 843 if(k!==n) { Ai[k] -= Aj[k]*alpha; } 844 } 845 if(Aj[j] === 0) { return 0; } 846 ret *= Aj[j]; 847 } 848 return ret*A[j][j]; 849} 850 851numeric.transpose = function transpose(x) { 852 var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj; 853 for(j=0;j<n;j++) ret[j] = Array(m); 854 for(i=m-1;i>=1;i-=2) { 855 A1 = x[i]; 856 A0 = x[i-1]; 857 for(j=n-1;j>=1;--j) { 858 Bj = ret[j]; Bj[i] = A1[j]; Bj[i-1] = A0[j]; 859 --j; 860 Bj = ret[j]; Bj[i] = A1[j]; Bj[i-1] = A0[j]; 861 } 862 if(j===0) { 863 Bj = ret[0]; Bj[i] = A1[0]; Bj[i-1] = A0[0]; 864 } 865 } 866 if(i===0) { 867 A0 = x[0]; 868 for(j=n-1;j>=1;--j) { 869 ret[j][0] = A0[j]; 870 --j; 871 ret[j][0] = A0[j]; 872 } 873 if(j===0) { ret[0][0] = A0[0]; } 874 } 875 return ret; 876} 877numeric.negtranspose = function negtranspose(x) { 878 var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj; 879 for(j=0;j<n;j++) ret[j] = Array(m); 880 for(i=m-1;i>=1;i-=2) { 881 A1 = x[i]; 882 A0 = x[i-1]; 883 for(j=n-1;j>=1;--j) { 884 Bj = ret[j]; Bj[i] = -A1[j]; Bj[i-1] = -A0[j]; 885 --j; 886 Bj = ret[j]; Bj[i] = -A1[j]; Bj[i-1] = -A0[j]; 887 } 888 if(j===0) { 889 Bj = ret[0]; Bj[i] = -A1[0]; Bj[i-1] = -A0[0]; 890 } 891 } 892 if(i===0) { 893 A0 = x[0]; 894 for(j=n-1;j>=1;--j) { 895 ret[j][0] = -A0[j]; 896 --j; 897 ret[j][0] = -A0[j]; 898 } 899 if(j===0) { ret[0][0] = -A0[0]; } 900 } 901 return ret; 902} 903 904numeric._random = function _random(s,k) { 905 var i,n=s[k],ret=Array(n), rnd; 906 if(k === s.length-1) { 907 rnd = Math.random; 908 for(i=n-1;i>=1;i-=2) { 909 ret[i] = rnd(); 910 ret[i-1] = rnd(); 911 } 912 if(i===0) { ret[0] = rnd(); } 913 return ret; 914 } 915 for(i=n-1;i>=0;i--) ret[i] = _random(s,k+1); 916 return ret; 917} 918numeric.random = function random(s) { return numeric._random(s,0); } 919 920numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); } 921 922numeric.linspace = function linspace(a,b,n) { 923 if(typeof n === "undefined") n = Math.max(Math.round(b-a)+1,1); 924 if(n<2) { return n===1?[a]:[]; } 925 var i,ret = Array(n); 926 n--; 927 for(i=n;i>=0;i--) { ret[i] = (i*b+(n-i)*a)/n; } 928 return ret; 929} 930 931numeric.getBlock = function getBlock(x,from,to) { 932 var s = numeric.dim(x); 933 function foo(x,k) { 934 var i,a = from[k], n = to[k]-a, ret = Array(n); 935 if(k === s.length-1) { 936 for(i=n;i>=0;i--) { ret[i] = x[i+a]; } 937 return ret; 938 } 939 for(i=n;i>=0;i--) { ret[i] = foo(x[i+a],k+1); } 940 return ret; 941 } 942 return foo(x,0); 943} 944 945numeric.setBlock = function setBlock(x,from,to,B) { 946 var s = numeric.dim(x); 947 function foo(x,y,k) { 948 var i,a = from[k], n = to[k]-a; 949 if(k === s.length-1) { for(i=n;i>=0;i--) { x[i+a] = y[i]; } } 950 for(i=n;i>=0;i--) { foo(x[i+a],y[i],k+1); } 951 } 952 foo(x,B,0); 953 return x; 954} 955 956numeric.getRange = function getRange(A,I,J) { 957 var m = I.length, n = J.length; 958 var i,j; 959 var B = Array(m), Bi, AI; 960 for(i=m-1;i!==-1;--i) { 961 B[i] = Array(n); 962 Bi = B[i]; 963 AI = A[I[i]]; 964 for(j=n-1;j!==-1;--j) Bi[j] = AI[J[j]]; 965 } 966 return B; 967} 968 969numeric.blockMatrix = function blockMatrix(X) { 970 var s = numeric.dim(X); 971 if(s.length<4) return numeric.blockMatrix([X]); 972 var m=s[0],n=s[1],M,N,i,j,Xij; 973 M = 0; N = 0; 974 for(i=0;i<m;++i) M+=X[i][0].length; 975 for(j=0;j<n;++j) N+=X[0][j][0].length; 976 var Z = Array(M); 977 for(i=0;i<M;++i) Z[i] = Array(N); 978 var I=0,J,ZI,k,l,Xijk; 979 for(i=0;i<m;++i) { 980 J=N; 981 for(j=n-1;j!==-1;--j) { 982 Xij = X[i][j]; 983 J -= Xij[0].length; 984 for(k=Xij.length-1;k!==-1;--k) { 985 Xijk = Xij[k]; 986 ZI = Z[I+k]; 987 for(l = Xijk.length-1;l!==-1;--l) ZI[J+l] = Xijk[l]; 988 } 989 } 990 I += X[i][0].length; 991 } 992 return Z; 993} 994 995numeric.tensor = function tensor(x,y) { 996 if(typeof x === "number" || typeof y === "number") return numeric.mul(x,y); 997 var s1 = numeric.dim(x), s2 = numeric.dim(y); 998 if(s1.length !== 1 || s2.length !== 1) { 999 throw new Error('numeric: tensor product is only defined for vectors'); 1000 } 1001 var m = s1[0], n = s2[0], A = Array(m), Ai, i,j,xi; 1002 for(i=m-1;i>=0;i--) { 1003 Ai = Array(n); 1004 xi = x[i]; 1005 for(j=n-1;j>=3;--j) { 1006 Ai[j] = xi * y[j]; 1007 --j; 1008 Ai[j] = xi * y[j]; 1009 --j; 1010 Ai[j] = xi * y[j]; 1011 --j; 1012 Ai[j] = xi * y[j]; 1013 } 1014 while(j>=0) { Ai[j] = xi * y[j]; --j; } 1015 A[i] = Ai; 1016 } 1017 return A; 1018} 1019 1020// 3. The Tensor type T 1021numeric.T = function T(x,y) { this.x = x; this.y = y; } 1022numeric.t = function t(x,y) { return new numeric.T(x,y); } 1023 1024numeric.Tbinop = function Tbinop(rr,rc,cr,cc,setup) { 1025 var io = numeric.indexOf; 1026 if(typeof setup !== "string") { 1027 var k; 1028 setup = ''; 1029 for(k in numeric) { 1030 if(numeric.hasOwnProperty(k) && (rr.indexOf(k)>=0 || rc.indexOf(k)>=0 || cr.indexOf(k)>=0 || cc.indexOf(k)>=0) && k.length>1) { 1031 setup += 'var '+k+' = numeric.'+k+';\n'; 1032 } 1033 } 1034 } 1035 return Function(['y'], 1036 'var x = this;\n'+ 1037 'if(!(y instanceof numeric.T)) { y = new numeric.T(y); }\n'+ 1038 setup+'\n'+ 1039 'if(x.y) {'+ 1040 ' if(y.y) {'+ 1041 ' return new numeric.T('+cc+');\n'+ 1042 ' }\n'+ 1043 ' return new numeric.T('+cr+');\n'+ 1044 '}\n'+ 1045 'if(y.y) {\n'+ 1046 ' return new numeric.T('+rc+');\n'+ 1047 '}\n'+ 1048 'return new numeric.T('+rr+');\n' 1049 ); 1050} 1051 1052numeric.T.prototype.add = numeric.Tbinop( 1053 'add(x.x,y.x)', 1054 'add(x.x,y.x),y.y', 1055 'add(x.x,y.x),x.y', 1056 'add(x.x,y.x),add(x.y,y.y)'); 1057numeric.T.prototype.sub = numeric.Tbinop( 1058 'sub(x.x,y.x)', 1059 'sub(x.x,y.x),neg(y.y)', 1060 'sub(x.x,y.x),x.y', 1061 'sub(x.x,y.x),sub(x.y,y.y)'); 1062numeric.T.prototype.mul = numeric.Tbinop( 1063 'mul(x.x,y.x)', 1064 'mul(x.x,y.x),mul(x.x,y.y)', 1065 'mul(x.x,y.x),mul(x.y,y.x)', 1066 'sub(mul(x.x,y.x),mul(x.y,y.y)),add(mul(x.x,y.y),mul(x.y,y.x))'); 1067 1068numeric.T.prototype.reciprocal = function reciprocal() { 1069 var mul = numeric.mul, div = numeric.div; 1070 if(this.y) { 1071 var d = numeric.add(mul(this.x,this.x),mul(this.y,this.y)); 1072 return new numeric.T(div(this.x,d),div(numeric.neg(this.y),d)); 1073 } 1074 return new T(div(1,this.x)); 1075} 1076numeric.T.prototype.div = function div(y) { 1077 if(!(y instanceof numeric.T)) y = new numeric.T(y); 1078 if(y.y) { return this.mul(y.reciprocal()); } 1079 var div = numeric.div; 1080 if(this.y) { return new numeric.T(div(this.x,y.x),div(this.y,y.x)); } 1081 return new numeric.T(div(this.x,y.x)); 1082} 1083numeric.T.prototype.dot = numeric.Tbinop( 1084 'dot(x.x,y.x)', 1085 'dot(x.x,y.x),dot(x.x,y.y)', 1086 'dot(x.x,y.x),dot(x.y,y.x)', 1087 'sub(dot(x.x,y.x),dot(x.y,y.y)),add(dot(x.x,y.y),dot(x.y,y.x))' 1088 ); 1089numeric.T.prototype.transpose = function transpose() { 1090 var t = numeric.transpose, x = this.x, y = this.y; 1091 if(y) { return new numeric.T(t(x),t(y)); } 1092 return new numeric.T(t(x)); 1093} 1094numeric.T.prototype.transjugate = function transjugate() { 1095 var t = numeric.transpose, x = this.x, y = this.y; 1096 if(y) { return new numeric.T(t(x),numeric.negtranspose(y)); } 1097 return new numeric.T(t(x)); 1098} 1099numeric.Tunop = function Tunop(r,c,s) { 1100 if(typeof s !== "string") { s = ''; } 1101 return Function( 1102 'var x = this;\n'+ 1103 s+'\n'+ 1104 'if(x.y) {'+ 1105 ' '+c+';\n'+ 1106 '}\n'+ 1107 r+';\n' 1108 ); 1109} 1110 1111numeric.T.prototype.exp = numeric.Tunop( 1112 'return new numeric.T(ex)', 1113 'return new numeric.T(mul(cos(x.y),ex),mul(sin(x.y),ex))', 1114 'var ex = numeric.exp(x.x), cos = numeric.cos, sin = numeric.sin, mul = numeric.mul;'); 1115numeric.T.prototype.conj = numeric.Tunop( 1116 'return new numeric.T(x.x);', 1117 'return new numeric.T(x.x,numeric.neg(x.y));'); 1118numeric.T.prototype.neg = numeric.Tunop( 1119 'return new numeric.T(neg(x.x));', 1120 'return new numeric.T(neg(x.x),neg(x.y));', 1121 'var neg = numeric.neg;'); 1122numeric.T.prototype.sin = numeric.Tunop( 1123 'return new numeric.T(numeric.sin(x.x))', 1124 'return x.exp().sub(x.neg().exp()).div(new numeric.T(0,2));'); 1125numeric.T.prototype.cos = numeric.Tunop( 1126 'return new numeric.T(numeric.cos(x.x))', 1127 'return x.exp().add(x.neg().exp()).div(2);'); 1128numeric.T.prototype.abs = numeric.Tunop( 1129 'return new numeric.T(numeric.abs(x.x));', 1130 'return new numeric.T(numeric.sqrt(numeric.add(mul(x.x,x.x),mul(x.y,x.y))));', 1131 'var mul = numeric.mul;'); 1132numeric.T.prototype.log = numeric.Tunop( 1133 'return new numeric.T(numeric.log(x.x));', 1134 'var theta = new numeric.T(numeric.atan2(x.y,x.x)), r = x.abs();\n'+ 1135 'return new numeric.T(numeric.log(r.x),theta.x);'); 1136numeric.T.prototype.norm2 = numeric.Tunop( 1137 'return numeric.norm2(x.x);', 1138 'var f = numeric.norm2Squared;\n'+ 1139 'return Math.sqrt(f(x.x)+f(x.y));'); 1140numeric.T.prototype.inv = function inv() { 1141 var A = this; 1142 if(typeof A.y === "undefined") { return new numeric.T(numeric.inv(A.x)); } 1143 var n = A.x.length, i, j, k; 1144 var Rx = numeric.identity(n),Ry = numeric.rep([n,n],0); 1145 var Ax = numeric.clone(A.x), Ay = numeric.clone(A.y); 1146 var Aix, Aiy, Ajx, Ajy, Rix, Riy, Rjx, Rjy; 1147 var i,j,k,d,d1,ax,ay,bx,by,temp; 1148 for(i=0;i<n;i++) { 1149 ax = Ax[i][i]; ay = Ay[i][i]; 1150 d = ax*ax+ay*ay; 1151 k = i; 1152 for(j=i+1;j<n;j++) { 1153 ax = Ax[j][i]; ay = Ay[j][i]; 1154 d1 = ax*ax+ay*ay; 1155 if(d1 > d) { k=j; d = d1; } 1156 } 1157 if(k!==i) { 1158 temp = Ax[i]; Ax[i] = Ax[k]; Ax[k] = temp; 1159 temp = Ay[i]; Ay[i] = Ay[k]; Ay[k] = temp; 1160 temp = Rx[i]; Rx[i] = Rx[k]; Rx[k] = temp; 1161 temp = Ry[i]; Ry[i] = Ry[k]; Ry[k] = temp; 1162 } 1163 Aix = Ax[i]; Aiy = Ay[i]; 1164 Rix = Rx[i]; Riy = Ry[i]; 1165 ax = Aix[i]; ay = Aiy[i]; 1166 for(j=i+1;j<n;j++) { 1167 bx = Aix[j]; by = Aiy[j]; 1168 Aix[j] = (bx*ax+by*ay)/d; 1169 Aiy[j] = (by*ax-bx*ay)/d; 1170 } 1171 for(j=0;j<n;j++) { 1172 bx = Rix[j]; by = Riy[j]; 1173 Rix[j] = (bx*ax+by*ay)/d; 1174 Riy[j] = (by*ax-bx*ay)/d; 1175 } 1176 for(j=i+1;j<n;j++) { 1177 Ajx = Ax[j]; Ajy = Ay[j]; 1178 Rjx = Rx[j]; Rjy = Ry[j]; 1179 ax = Ajx[i]; ay = Ajy[i]; 1180 for(k=i+1;k<n;k++) { 1181 bx = Aix[k]; by = Aiy[k]; 1182 Ajx[k] -= bx*ax-by*ay; 1183 Ajy[k] -= by*ax+bx*ay; 1184 } 1185 for(k=0;k<n;k++) { 1186 bx = Rix[k]; by = Riy[k]; 1187 Rjx[k] -= bx*ax-by*ay; 1188 Rjy[k] -= by*ax+bx*ay; 1189 } 1190 } 1191 } 1192 for(i=n-1;i>0;i--) { 1193 Rix = Rx[i]; Riy = Ry[i]; 1194 for(j=i-1;j>=0;j--) { 1195 Rjx = Rx[j]; Rjy = Ry[j]; 1196 ax = Ax[j][i]; ay = Ay[j][i]; 1197 for(k=n-1;k>=0;k--) { 1198 bx = Rix[k]; by = Riy[k]; 1199 Rjx[k] -= ax*bx - ay*by; 1200 Rjy[k] -= ax*by + ay*bx; 1201 } 1202 } 1203 } 1204 return new numeric.T(Rx,Ry); 1205} 1206numeric.T.prototype.get = function get(i) { 1207 var x = this.x, y = this.y, k = 0, ik, n = i.length; 1208 if(y) { 1209 while(k<n) { 1210 ik = i[k]; 1211 x = x[ik]; 1212 y = y[ik]; 1213 k++; 1214 } 1215 return new numeric.T(x,y); 1216 } 1217 while(k<n) { 1218 ik = i[k]; 1219 x = x[ik]; 1220 k++; 1221 } 1222 return new numeric.T(x); 1223} 1224numeric.T.prototype.set = function set(i,v) { 1225 var x = this.x, y = this.y, k = 0, ik, n = i.length, vx = v.x, vy = v.y; 1226 if(n===0) { 1227 if(vy) { this.y = vy; } 1228 else if(y) { this.y = undefined; } 1229 this.x = x; 1230 return this; 1231 } 1232 if(vy) { 1233 if(y) { /* ok */ } 1234 else { 1235 y = numeric.rep(numeric.dim(x),0); 1236 this.y = y; 1237 } 1238 while(k<n-1) { 1239 ik = i[k]; 1240 x = x[ik]; 1241 y = y[ik]; 1242 k++; 1243 } 1244 ik = i[k]; 1245 x[ik] = vx; 1246 y[ik] = vy; 1247 return this; 1248 } 1249 if(y) { 1250 while(k<n-1) { 1251 ik = i[k]; 1252 x = x[ik]; 1253 y = y[ik]; 1254 k++; 1255 } 1256 ik = i[k]; 1257 x[ik] = vx; 1258 if(vx instanceof Array) y[ik] = numeric.rep(numeric.dim(vx),0); 1259 else y[ik] = 0; 1260 return this; 1261 } 1262 while(k<n-1) { 1263 ik = i[k]; 1264 x = x[ik]; 1265 k++; 1266 } 1267 ik = i[k]; 1268 x[ik] = vx; 1269 return this; 1270} 1271numeric.T.prototype.getRows = function getRows(i0,i1) { 1272 var n = i1-i0+1, j; 1273 var rx = Array(n), ry, x = this.x, y = this.y; 1274 for(j=i0;j<=i1;j++) { rx[j-i0] = x[j]; } 1275 if(y) { 1276 ry = Array(n); 1277 for(j=i0;j<=i1;j++) { ry[j-i0] = y[j]; } 1278 return new numeric.T(rx,ry); 1279 } 1280 return new numeric.T(rx); 1281} 1282numeric.T.prototype.setRows = function setRows(i0,i1,A) { 1283 var j; 1284 var rx = this.x, ry = this.y, x = A.x, y = A.y; 1285 for(j=i0;j<=i1;j++) { rx[j] = x[j-i0]; } 1286 if(y) { 1287 if(!ry) { ry = numeric.rep(numeric.dim(rx),0); this.y = ry; } 1288 for(j=i0;j<=i1;j++) { ry[j] = y[j-i0]; } 1289 } else if(ry) { 1290 for(j=i0;j<=i1;j++) { ry[j] = numeric.rep([x[j-i0].length],0); } 1291 } 1292 return this; 1293} 1294numeric.T.prototype.getRow = function getRow(k) { 1295 var x = this.x, y = this.y; 1296 if(y) { return new numeric.T(x[k],y[k]); } 1297 return new numeric.T(x[k]); 1298} 1299numeric.T.prototype.setRow = function setRow(i,v) { 1300 var rx = this.x, ry = this.y, x = v.x, y = v.y; 1301 rx[i] = x; 1302 if(y) { 1303 if(!ry) { ry = numeric.rep(numeric.dim(rx),0); this.y = ry; } 1304 ry[i] = y; 1305 } else if(ry) { 1306 ry = numeric.rep([x.length],0); 1307 } 1308 return this; 1309} 1310 1311numeric.T.prototype.getBlock = function getBlock(from,to) { 1312 var x = this.x, y = this.y, b = numeric.getBlock; 1313 if(y) { return new numeric.T(b(x,from,to),b(y,from,to)); } 1314 return new numeric.T(b(x,from,to)); 1315} 1316numeric.T.prototype.setBlock = function setBlock(from,to,A) { 1317 if(!(A instanceof numeric.T)) A = new numeric.T(A); 1318 var x = this.x, y = this.y, b = numeric.setBlock, Ax = A.x, Ay = A.y; 1319 if(Ay) { 1320 if(!y) { this.y = numeric.rep(numeric.dim(this),0); y = this.y; } 1321 b(x,from,to,Ax); 1322 b(y,from,to,Ay); 1323 return this; 1324 } 1325 b(x,from,to,Ax); 1326 if(y) b(y,from,to,numeric.rep(numeric.dim(Ax),0)); 1327} 1328numeric.T.rep = function rep(s,v) { 1329 var T = numeric.T; 1330 if(!(v instanceof T)) v = new T(v); 1331 var x = v.x, y = v.y, r = numeric.rep; 1332 if(y) return new T(r(s,x),r(s,y)); 1333 return new T(r(s,x)); 1334} 1335numeric.T.diag = function diag(d) { 1336 if(!(d instanceof numeric.T)) d = new numeric.T(d); 1337 var x = d.x, y = d.y, diag = numeric.diag; 1338 if(y) return new numeric.T(diag(x),diag(y)); 1339 return new numeric.T(diag(x)); 1340} 1341numeric.T.eig = function eig() { 1342 if(this.y) { throw new Error('eig: not implemented for complex matrices.'); } 1343 return numeric.eig(this.x); 1344} 1345numeric.T.identity = function identity(n) { return new numeric.T(numeric.identity(n)); } 1346numeric.T.prototype.getDiag = function getDiag() { 1347 var n = numeric; 1348 var x = this.x, y = this.y; 1349 if(y) { return new n.T(n.getDiag(x),n.getDiag(y)); } 1350 return new n.T(n.getDiag(x)); 1351} 1352 1353// 4. Eigenvalues of real matrices 1354 1355numeric.house = function house(x) { 1356 var v = numeric.clone(x); 1357 var s = x[0] >= 0 ? 1 : -1; 1358 var alpha = s*numeric.norm2(x); 1359 v[0] += alpha; 1360 var foo = numeric.norm2(v); 1361 if(foo === 0) { /* this should not happen */ throw new Error('eig: internal error'); } 1362 return numeric.div(v,foo); 1363} 1364 1365numeric.toUpperHessenberg = function toUpperHessenberg(me) { 1366 var s = numeric.dim(me); 1367 if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: toUpperHessenberg() only works on square matrices'); } 1368 var m = s[0], i,j,k,x,v,A = numeric.clone(me),B,C,Ai,Ci,Q = numeric.identity(m),Qi; 1369 for(j=0;j<m-2;j++) { 1370 x = Array(m-j-1); 1371 for(i=j+1;i<m;i++) { x[i-j-1] = A[i][j]; } 1372 if(numeric.norm2(x)>0) { 1373 v = numeric.house(x); 1374 B = numeric.getBlock(A,[j+1,j],[m-1,m-1]); 1375 C = numeric.tensor(v,numeric.dot(v,B)); 1376 for(i=j+1;i<m;i++) { Ai = A[i]; Ci = C[i-j-1]; for(k=j;k<m;k++) Ai[k] -= 2*Ci[k-j]; } 1377 B = numeric.getBlock(A,[0,j+1],[m-1,m-1]); 1378 C = numeric.tensor(numeric.dot(B,v),v); 1379 for(i=0;i<m;i++) { Ai = A[i]; Ci = C[i]; for(k=j+1;k<m;k++) Ai[k] -= 2*Ci[k-j-1]; } 1380 B = Array(m-j-1); 1381 for(i=j+1;i<m;i++) B[i-j-1] = Q[i]; 1382 C = numeric.tensor(v,numeric.dot(v,B)); 1383 for(i=j+1;i<m;i++) { Qi = Q[i]; Ci = C[i-j-1]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; } 1384 } 1385 } 1386 return {H:A, Q:Q}; 1387}
vendor: 4,298 bytes, lines 1387-1476
1387 1388 1389numeric.epsilon = 2.220446049250313e-16; 1390 1391numeric.QRFrancis = function(H,maxiter) { 1392 if(typeof maxiter === "undefined") { maxiter = 10000; } 1393 H = numeric.clone(H); 1394 var H0 = numeric.clone(H); 1395 var s = numeric.dim(H),m=s[0],x,v,a,b,c,d,det,tr, Hloc, Q = numeric.identity(m), Qi, Hi, B, C, Ci,i,j,k,iter; 1396 if(m<3) { return {Q:Q, B:[ [0,m-1] ]}; } 1397 var epsilon = numeric.epsilon; 1398 for(iter=0;iter<maxiter;iter++) { 1399 for(j=0;j<m-1;j++) { 1400 if(Math.abs(H[j+1][j]) < epsilon*(Math.abs(H[j][j])+Math.abs(H[j+1][j+1]))) { 1401 var QH1 = numeric.QRFrancis(numeric.getBlock(H,[0,0],[j,j]),maxiter); 1402 var QH2 = numeric.QRFrancis(numeric.getBlock(H,[j+1,j+1],[m-1,m-1]),maxiter); 1403 B = Array(j+1); 1404 for(i=0;i<=j;i++) { B[i] = Q[i]; } 1405 C = numeric.dot(QH1.Q,B); 1406 for(i=0;i<=j;i++) { Q[i] = C[i]; } 1407 B = Array(m-j-1); 1408 for(i=j+1;i<m;i++) { B[i-j-1] = Q[i]; } 1409 C = numeric.dot(QH2.Q,B); 1410 for(i=j+1;i<m;i++) { Q[i] = C[i-j-1]; } 1411 return {Q:Q,B:QH1.B.concat(numeric.add(QH2.B,j+1))}; 1412 } 1413 } 1414 a = H[m-2][m-2]; b = H[m-2][m-1]; 1415 c = H[m-1][m-2]; d = H[m-1][m-1]; 1416 tr = a+d; 1417 det = (a*d-b*c); 1418 Hloc = numeric.getBlock(H, [0,0], [2,2]); 1419 if(tr*tr>=4*det) { 1420 var s1,s2; 1421 s1 = 0.5*(tr+Math.sqrt(tr*tr-4*det)); 1422 s2 = 0.5*(tr-Math.sqrt(tr*tr-4*det)); 1423 Hloc = numeric.add(numeric.sub(numeric.dot(Hloc,Hloc), 1424 numeric.mul(Hloc,s1+s2)), 1425 numeric.diag(numeric.rep([3],s1*s2))); 1426 } else { 1427 Hloc = numeric.add(numeric.sub(numeric.dot(Hloc,Hloc), 1428 numeric.mul(Hloc,tr)), 1429 numeric.diag(numeric.rep([3],det))); 1430 } 1431 x = [Hloc[0][0],Hloc[1][0],Hloc[2][0]]; 1432 v = numeric.house(x); 1433 B = [H[0],H[1],H[2]]; 1434 C = numeric.tensor(v,numeric.dot(v,B)); 1435 for(i=0;i<3;i++) { Hi = H[i]; Ci = C[i]; for(k=0;k<m;k++) Hi[k] -= 2*Ci[k]; } 1436 B = numeric.getBlock(H, [0,0],[m-1,2]); 1437 C = numeric.tensor(numeric.dot(B,v),v); 1438 for(i=0;i<m;i++) { Hi = H[i]; Ci = C[i]; for(k=0;k<3;k++) Hi[k] -= 2*Ci[k]; } 1439 B = [Q[0],Q[1],Q[2]]; 1440 C = numeric.tensor(v,numeric.dot(v,B)); 1441 for(i=0;i<3;i++) { Qi = Q[i]; Ci = C[i]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; } 1442 var J; 1443 for(j=0;j<m-2;j++) { 1444 for(k=j;k<=j+1;k++) { 1445 if(Math.abs(H[k+1][k]) < epsilon*(Math.abs(H[k][k])+Math.abs(H[k+1][k+1]))) { 1446 var QH1 = numeric.QRFrancis(numeric.getBlock(H,[0,0],[k,k]),maxiter); 1447 var QH2 = numeric.QRFrancis(numeric.getBlock(H,[k+1,k+1],[m-1,m-1]),maxiter); 1448 B = Array(k+1); 1449 for(i=0;i<=k;i++) { B[i] = Q[i]; } 1450 C = numeric.dot(QH1.Q,B); 1451 for(i=0;i<=k;i++) { Q[i] = C[i]; } 1452 B = Array(m-k-1); 1453 for(i=k+1;i<m;i++) { B[i-k-1] = Q[i]; } 1454 C = numeric.dot(QH2.Q,B); 1455 for(i=k+1;i<m;i++) { Q[i] = C[i-k-1]; } 1456 return {Q:Q,B:QH1.B.concat(numeric.add(QH2.B,k+1))}; 1457 } 1458 } 1459 J = Math.min(m-1,j+3); 1460 x = Array(J-j); 1461 for(i=j+1;i<=J;i++) { x[i-j-1] = H[i][j]; } 1462 v = numeric.house(x); 1463 B = numeric.getBlock(H, [j+1,j],[J,m-1]); 1464 C = numeric.tensor(v,numeric.dot(v,B)); 1465 for(i=j+1;i<=J;i++) { Hi = H[i]; Ci = C[i-j-1]; for(k=j;k<m;k++) Hi[k] -= 2*Ci[k-j]; } 1466 B = numeric.getBlock(H, [0,j+1],[m-1,J]); 1467 C = numeric.tensor(numeric.dot(B,v),v); 1468 for(i=0;i<m;i++) { Hi = H[i]; Ci = C[i]; for(k=j+1;k<=J;k++) Hi[k] -= 2*Ci[k-j-1]; } 1469 B = Array(J-j); 1470 for(i=j+1;i<=J;i++) B[i-j-1] = Q[i]; 1471 C = numeric.tensor(v,numeric.dot(v,B)); 1472 for(i=j+1;i<=J;i++) { Qi = Q[i]; Ci = C[i-j-1]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; } 1473 } 1474 } 1475 throw new Error('numeric: eigenvalue iteration does not converge -- increase maxiter?'); 1476}
vendor: 56,498 bytes, lines 1476-3263
1476 1477 1478numeric.eig = function eig(A,maxiter) { 1479 var QH = numeric.toUpperHessenberg(A); 1480 var QB = numeric.QRFrancis(QH.H,maxiter); 1481 var T = numeric.T; 1482 var n = A.length,i,k,flag = false,B = QB.B,H = numeric.dot(QB.Q,numeric.dot(QH.H,numeric.transpose(QB.Q))); 1483 var Q = new T(numeric.dot(QB.Q,QH.Q)),Q0; 1484 var m = B.length,j; 1485 var a,b,c,d,p1,p2,disc,x,y,p,q,n1,n2; 1486 var sqrt = Math.sqrt; 1487 for(k=0;k<m;k++) { 1488 i = B[k][0]; 1489 if(i === B[k][1]) { 1490 // nothing 1491 } else { 1492 j = i+1; 1493 a = H[i][i]; 1494 b = H[i][j]; 1495 c = H[j][i]; 1496 d = H[j][j]; 1497 if(b === 0 && c === 0) continue; 1498 p1 = -a-d; 1499 p2 = a*d-b*c; 1500 disc = p1*p1-4*p2; 1501 if(disc>=0) { 1502 if(p1<0) x = -0.5*(p1-sqrt(disc)); 1503 else x = -0.5*(p1+sqrt(disc)); 1504 n1 = (a-x)*(a-x)+b*b; 1505 n2 = c*c+(d-x)*(d-x); 1506 if(n1>n2) { 1507 n1 = sqrt(n1); 1508 p = (a-x)/n1; 1509 q = b/n1; 1510 } else { 1511 n2 = sqrt(n2); 1512 p = c/n2; 1513 q = (d-x)/n2; 1514 } 1515 Q0 = new T([[q,-p],[p,q]]); 1516 Q.setRows(i,j,Q0.dot(Q.getRows(i,j))); 1517 } else { 1518 x = -0.5*p1; 1519 y = 0.5*sqrt(-disc); 1520 n1 = (a-x)*(a-x)+b*b; 1521 n2 = c*c+(d-x)*(d-x); 1522 if(n1>n2) { 1523 n1 = sqrt(n1+y*y); 1524 p = (a-x)/n1; 1525 q = b/n1; 1526 x = 0; 1527 y /= n1; 1528 } else { 1529 n2 = sqrt(n2+y*y); 1530 p = c/n2; 1531 q = (d-x)/n2; 1532 x = y/n2; 1533 y = 0; 1534 } 1535 Q0 = new T([[q,-p],[p,q]],[[x,y],[y,-x]]); 1536 Q.setRows(i,j,Q0.dot(Q.getRows(i,j))); 1537 } 1538 } 1539 } 1540 var R = Q.dot(A).dot(Q.transjugate()), n = A.length, E = numeric.T.identity(n); 1541 for(j=0;j<n;j++) { 1542 if(j>0) { 1543 for(k=j-1;k>=0;k--) { 1544 var Rk = R.get([k,k]), Rj = R.get([j,j]); 1545 if(numeric.neq(Rk.x,Rj.x) || numeric.neq(Rk.y,Rj.y)) { 1546 x = R.getRow(k).getBlock([k],[j-1]); 1547 y = E.getRow(j).getBlock([k],[j-1]); 1548 E.set([j,k],(R.get([k,j]).neg().sub(x.dot(y))).div(Rk.sub(Rj))); 1549 } else { 1550 E.setRow(j,E.getRow(k)); 1551 continue; 1552 } 1553 } 1554 } 1555 } 1556 for(j=0;j<n;j++) { 1557 x = E.getRow(j); 1558 E.setRow(j,x.div(x.norm2())); 1559 } 1560 E = E.transpose(); 1561 E = Q.transjugate().dot(E); 1562 return { lambda:R.getDiag(), E:E }; 1563}; 1564 1565// 5. Compressed Column Storage matrices 1566numeric.ccsSparse = function ccsSparse(A) { 1567 var m = A.length,n,foo, i,j, counts = []; 1568 for(i=m-1;i!==-1;--i) { 1569 foo = A[i]; 1570 for(j in foo) { 1571 j = parseInt(j); 1572 while(j>=counts.length) counts[counts.length] = 0; 1573 if(foo[j]!==0) counts[j]++; 1574 } 1575 } 1576 var n = counts.length; 1577 var Ai = Array(n+1); 1578 Ai[0] = 0; 1579 for(i=0;i<n;++i) Ai[i+1] = Ai[i] + counts[i]; 1580 var Aj = Array(Ai[n]), Av = Array(Ai[n]); 1581 for(i=m-1;i!==-1;--i) { 1582 foo = A[i]; 1583 for(j in foo) { 1584 if(foo[j]!==0) { 1585 counts[j]--; 1586 Aj[Ai[j]+counts[j]] = i; 1587 Av[Ai[j]+counts[j]] = foo[j]; 1588 } 1589 } 1590 } 1591 return [Ai,Aj,Av]; 1592} 1593numeric.ccsFull = function ccsFull(A) { 1594 var Ai = A[0], Aj = A[1], Av = A[2], s = numeric.ccsDim(A), m = s[0], n = s[1], i,j,j0,j1,k; 1595 var B = numeric.rep([m,n],0); 1596 for(i=0;i<n;i++) { 1597 j0 = Ai[i]; 1598 j1 = Ai[i+1]; 1599 for(j=j0;j<j1;++j) { B[Aj[j]][i] = Av[j]; } 1600 } 1601 return B; 1602} 1603numeric.ccsTSolve = function ccsTSolve(A,b,x,bj,xj) { 1604 var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, max = Math.max,n=0; 1605 if(typeof bj === "undefined") x = numeric.rep([m],0); 1606 if(typeof bj === "undefined") bj = numeric.linspace(0,x.length-1); 1607 if(typeof xj === "undefined") xj = []; 1608 function dfs(j) { 1609 var k; 1610 if(x[j] !== 0) return; 1611 x[j] = 1; 1612 for(k=Ai[j];k<Ai[j+1];++k) dfs(Aj[k]); 1613 xj[n] = j; 1614 ++n; 1615 } 1616 var i,j,j0,j1,k,l,l0,l1,a; 1617 for(i=bj.length-1;i!==-1;--i) { dfs(bj[i]); } 1618 xj.length = n; 1619 for(i=xj.length-1;i!==-1;--i) { x[xj[i]] = 0; } 1620 for(i=bj.length-1;i!==-1;--i) { j = bj[i]; x[j] = b[j]; } 1621 for(i=xj.length-1;i!==-1;--i) { 1622 j = xj[i]; 1623 j0 = Ai[j]; 1624 j1 = max(Ai[j+1],j0); 1625 for(k=j0;k!==j1;++k) { if(Aj[k] === j) { x[j] /= Av[k]; break; } } 1626 a = x[j]; 1627 for(k=j0;k!==j1;++k) { 1628 l = Aj[k]; 1629 if(l !== j) x[l] -= a*Av[k]; 1630 } 1631 } 1632 return x; 1633} 1634numeric.ccsDFS = function ccsDFS(n) { 1635 this.k = Array(n); 1636 this.k1 = Array(n); 1637 this.j = Array(n); 1638} 1639numeric.ccsDFS.prototype.dfs = function dfs(J,Ai,Aj,x,xj,Pinv) { 1640 var m = 0,foo,n=xj.length; 1641 var k = this.k, k1 = this.k1, j = this.j,km,k11; 1642 if(x[J]!==0) return; 1643 x[J] = 1; 1644 j[0] = J; 1645 k[0] = km = Ai[J]; 1646 k1[0] = k11 = Ai[J+1]; 1647 while(1) { 1648 if(km >= k11) { 1649 xj[n] = j[m]; 1650 if(m===0) return; 1651 ++n; 1652 --m; 1653 km = k[m]; 1654 k11 = k1[m]; 1655 } else { 1656 foo = Pinv[Aj[km]]; 1657 if(x[foo] === 0) { 1658 x[foo] = 1; 1659 k[m] = km; 1660 ++m; 1661 j[m] = foo; 1662 km = Ai[foo]; 1663 k1[m] = k11 = Ai[foo+1]; 1664 } else ++km; 1665 } 1666 } 1667} 1668numeric.ccsLPSolve = function ccsLPSolve(A,B,x,xj,I,Pinv,dfs) { 1669 var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, n=0; 1670 var Bi = B[0], Bj = B[1], Bv = B[2]; 1671 1672 var i,i0,i1,j,J,j0,j1,k,l,l0,l1,a; 1673 i0 = Bi[I]; 1674 i1 = Bi[I+1]; 1675 xj.length = 0; 1676 for(i=i0;i<i1;++i) { dfs.dfs(Pinv[Bj[i]],Ai,Aj,x,xj,Pinv); } 1677 for(i=xj.length-1;i!==-1;--i) { x[xj[i]] = 0; } 1678 for(i=i0;i!==i1;++i) { j = Pinv[Bj[i]]; x[j] = Bv[i]; } 1679 for(i=xj.length-1;i!==-1;--i) { 1680 j = xj[i]; 1681 j0 = Ai[j]; 1682 j1 = Ai[j+1]; 1683 for(k=j0;k<j1;++k) { if(Pinv[Aj[k]] === j) { x[j] /= Av[k]; break; } } 1684 a = x[j]; 1685 for(k=j0;k<j1;++k) { 1686 l = Pinv[Aj[k]]; 1687 if(l !== j) x[l] -= a*Av[k]; 1688 } 1689 } 1690 return x; 1691} 1692numeric.ccsLUP1 = function ccsLUP1(A,threshold) { 1693 var m = A[0].length-1; 1694 var L = [numeric.rep([m+1],0),[],[]], U = [numeric.rep([m+1], 0),[],[]]; 1695 var Li = L[0], Lj = L[1], Lv = L[2], Ui = U[0], Uj = U[1], Uv = U[2]; 1696 var x = numeric.rep([m],0), xj = numeric.rep([m],0); 1697 var i,j,k,j0,j1,a,e,c,d,K; 1698 var sol = numeric.ccsLPSolve, max = Math.max, abs = Math.abs; 1699 var P = numeric.linspace(0,m-1),Pinv = numeric.linspace(0,m-1); 1700 var dfs = new numeric.ccsDFS(m); 1701 if(typeof threshold === "undefined") { threshold = 1; } 1702 for(i=0;i<m;++i) { 1703 sol(L,A,x,xj,i,Pinv,dfs); 1704 a = -1; 1705 e = -1; 1706 for(j=xj.length-1;j!==-1;--j) { 1707 k = xj[j]; 1708 if(k <= i) continue; 1709 c = abs(x[k]); 1710 if(c > a) { e = k; a = c; } 1711 } 1712 if(abs(x[i])<threshold*a) { 1713 j = P[i]; 1714 a = P[e]; 1715 P[i] = a; Pinv[a] = i; 1716 P[e] = j; Pinv[j] = e; 1717 a = x[i]; x[i] = x[e]; x[e] = a; 1718 } 1719 a = Li[i]; 1720 e = Ui[i]; 1721 d = x[i]; 1722 Lj[a] = P[i]; 1723 Lv[a] = 1; 1724 ++a; 1725 for(j=xj.length-1;j!==-1;--j) { 1726 k = xj[j]; 1727 c = x[k]; 1728 xj[j] = 0; 1729 x[k] = 0; 1730 if(k<=i) { Uj[e] = k; Uv[e] = c; ++e; } 1731 else { Lj[a] = P[k]; Lv[a] = c/d; ++a; } 1732 } 1733 Li[i+1] = a; 1734 Ui[i+1] = e; 1735 } 1736 for(j=Lj.length-1;j!==-1;--j) { Lj[j] = Pinv[Lj[j]]; } 1737 return {L:L, U:U, P:P, Pinv:Pinv}; 1738} 1739numeric.ccsDFS0 = function ccsDFS0(n) { 1740 this.k = Array(n); 1741 this.k1 = Array(n); 1742 this.j = Array(n); 1743} 1744numeric.ccsDFS0.prototype.dfs = function dfs(J,Ai,Aj,x,xj,Pinv,P) { 1745 var m = 0,foo,n=xj.length; 1746 var k = this.k, k1 = this.k1, j = this.j,km,k11; 1747 if(x[J]!==0) return; 1748 x[J] = 1; 1749 j[0] = J; 1750 k[0] = km = Ai[Pinv[J]]; 1751 k1[0] = k11 = Ai[Pinv[J]+1]; 1752 while(1) { 1753 if(isNaN(km)) throw new Error("Ow!"); 1754 if(km >= k11) { 1755 xj[n] = Pinv[j[m]]; 1756 if(m===0) return; 1757 ++n; 1758 --m; 1759 km = k[m]; 1760 k11 = k1[m]; 1761 } else { 1762 foo = Aj[km]; 1763 if(x[foo] === 0) { 1764 x[foo] = 1; 1765 k[m] = km; 1766 ++m; 1767 j[m] = foo; 1768 foo = Pinv[foo]; 1769 km = Ai[foo]; 1770 k1[m] = k11 = Ai[foo+1]; 1771 } else ++km; 1772 } 1773 } 1774} 1775numeric.ccsLPSolve0 = function ccsLPSolve0(A,B,y,xj,I,Pinv,P,dfs) { 1776 var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, n=0; 1777 var Bi = B[0], Bj = B[1], Bv = B[2]; 1778 1779 var i,i0,i1,j,J,j0,j1,k,l,l0,l1,a; 1780 i0 = Bi[I]; 1781 i1 = Bi[I+1]; 1782 xj.length = 0; 1783 for(i=i0;i<i1;++i) { dfs.dfs(Bj[i],Ai,Aj,y,xj,Pinv,P); } 1784 for(i=xj.length-1;i!==-1;--i) { j = xj[i]; y[P[j]] = 0; } 1785 for(i=i0;i!==i1;++i) { j = Bj[i]; y[j] = Bv[i]; } 1786 for(i=xj.length-1;i!==-1;--i) { 1787 j = xj[i]; 1788 l = P[j]; 1789 j0 = Ai[j]; 1790 j1 = Ai[j+1]; 1791 for(k=j0;k<j1;++k) { if(Aj[k] === l) { y[l] /= Av[k]; break; } } 1792 a = y[l]; 1793 for(k=j0;k<j1;++k) y[Aj[k]] -= a*Av[k]; 1794 y[l] = a; 1795 } 1796} 1797numeric.ccsLUP0 = function ccsLUP0(A,threshold) { 1798 var m = A[0].length-1; 1799 var L = [numeric.rep([m+1],0),[],[]], U = [numeric.rep([m+1], 0),[],[]]; 1800 var Li = L[0], Lj = L[1], Lv = L[2], Ui = U[0], Uj = U[1], Uv = U[2]; 1801 var y = numeric.rep([m],0), xj = numeric.rep([m],0); 1802 var i,j,k,j0,j1,a,e,c,d,K; 1803 var sol = numeric.ccsLPSolve0, max = Math.max, abs = Math.abs; 1804 var P = numeric.linspace(0,m-1),Pinv = numeric.linspace(0,m-1); 1805 var dfs = new numeric.ccsDFS0(m); 1806 if(typeof threshold === "undefined") { threshold = 1; } 1807 for(i=0;i<m;++i) { 1808 sol(L,A,y,xj,i,Pinv,P,dfs); 1809 a = -1; 1810 e = -1; 1811 for(j=xj.length-1;j!==-1;--j) { 1812 k = xj[j]; 1813 if(k <= i) continue; 1814 c = abs(y[P[k]]); 1815 if(c > a) { e = k; a = c; } 1816 } 1817 if(abs(y[P[i]])<threshold*a) { 1818 j = P[i]; 1819 a = P[e]; 1820 P[i] = a; Pinv[a] = i; 1821 P[e] = j; Pinv[j] = e; 1822 } 1823 a = Li[i]; 1824 e = Ui[i]; 1825 d = y[P[i]]; 1826 Lj[a] = P[i]; 1827 Lv[a] = 1; 1828 ++a; 1829 for(j=xj.length-1;j!==-1;--j) { 1830 k = xj[j]; 1831 c = y[P[k]]; 1832 xj[j] = 0; 1833 y[P[k]] = 0; 1834 if(k<=i) { Uj[e] = k; Uv[e] = c; ++e; } 1835 else { Lj[a] = P[k]; Lv[a] = c/d; ++a; } 1836 } 1837 Li[i+1] = a; 1838 Ui[i+1] = e; 1839 } 1840 for(j=Lj.length-1;j!==-1;--j) { Lj[j] = Pinv[Lj[j]]; } 1841 return {L:L, U:U, P:P, Pinv:Pinv}; 1842} 1843numeric.ccsLUP = numeric.ccsLUP0; 1844 1845numeric.ccsDim = function ccsDim(A) { return [numeric.sup(A[1])+1,A[0].length-1]; } 1846numeric.ccsGetBlock = function ccsGetBlock(A,i,j) { 1847 var s = numeric.ccsDim(A),m=s[0],n=s[1]; 1848 if(typeof i === "undefined") { i = numeric.linspace(0,m-1); } 1849 else if(typeof i === "number") { i = [i]; } 1850 if(typeof j === "undefined") { j = numeric.linspace(0,n-1); } 1851 else if(typeof j === "number") { j = [j]; } 1852 var p,p0,p1,P = i.length,q,Q = j.length,r,jq,ip; 1853 var Bi = numeric.rep([n],0), Bj=[], Bv=[], B = [Bi,Bj,Bv]; 1854 var Ai = A[0], Aj = A[1], Av = A[2]; 1855 var x = numeric.rep([m],0),count=0,flags = numeric.rep([m],0); 1856 for(q=0;q<Q;++q) { 1857 jq = j[q]; 1858 var q0 = Ai[jq]; 1859 var q1 = Ai[jq+1]; 1860 for(p=q0;p<q1;++p) { 1861 r = Aj[p]; 1862 flags[r] = 1; 1863 x[r] = Av[p]; 1864 } 1865 for(p=0;p<P;++p) { 1866 ip = i[p]; 1867 if(flags[ip]) { 1868 Bj[count] = p; 1869 Bv[count] = x[i[p]]; 1870 ++count; 1871 } 1872 } 1873 for(p=q0;p<q1;++p) { 1874 r = Aj[p]; 1875 flags[r] = 0; 1876 } 1877 Bi[q+1] = count; 1878 } 1879 return B; 1880} 1881 1882numeric.ccsDot = function ccsDot(A,B) { 1883 var Ai = A[0], Aj = A[1], Av = A[2]; 1884 var Bi = B[0], Bj = B[1], Bv = B[2]; 1885 var sA = numeric.ccsDim(A), sB = numeric.ccsDim(B); 1886 var m = sA[0], n = sA[1], o = sB[1]; 1887 var x = numeric.rep([m],0), flags = numeric.rep([m],0), xj = Array(m); 1888 var Ci = numeric.rep([o],0), Cj = [], Cv = [], C = [Ci,Cj,Cv]; 1889 var i,j,k,j0,j1,i0,i1,l,p,a,b; 1890 for(k=0;k!==o;++k) { 1891 j0 = Bi[k]; 1892 j1 = Bi[k+1]; 1893 p = 0; 1894 for(j=j0;j<j1;++j) { 1895 a = Bj[j]; 1896 b = Bv[j]; 1897 i0 = Ai[a]; 1898 i1 = Ai[a+1]; 1899 for(i=i0;i<i1;++i) { 1900 l = Aj[i]; 1901 if(flags[l]===0) { 1902 xj[p] = l; 1903 flags[l] = 1; 1904 p = p+1; 1905 } 1906 x[l] = x[l] + Av[i]*b; 1907 } 1908 } 1909 j0 = Ci[k]; 1910 j1 = j0+p; 1911 Ci[k+1] = j1; 1912 for(j=p-1;j!==-1;--j) { 1913 b = j0+j; 1914 i = xj[j]; 1915 Cj[b] = i; 1916 Cv[b] = x[i]; 1917 flags[i] = 0; 1918 x[i] = 0; 1919 } 1920 Ci[k+1] = Ci[k]+p; 1921 } 1922 return C; 1923} 1924 1925numeric.ccsLUPSolve = function ccsLUPSolve(LUP,B) { 1926 var L = LUP.L, U = LUP.U, P = LUP.P; 1927 var Bi = B[0]; 1928 var flag = false; 1929 if(typeof Bi !== "object") { B = [[0,B.length],numeric.linspace(0,B.length-1),B]; Bi = B[0]; flag = true; } 1930 var Bj = B[1], Bv = B[2]; 1931 var n = L[0].length-1, m = Bi.length-1; 1932 var x = numeric.rep([n],0), xj = Array(n); 1933 var b = numeric.rep([n],0), bj = Array(n); 1934 var Xi = numeric.rep([m+1],0), Xj = [], Xv = []; 1935 var sol = numeric.ccsTSolve; 1936 var i,j,j0,j1,k,J,N=0; 1937 for(i=0;i<m;++i) { 1938 k = 0; 1939 j0 = Bi[i]; 1940 j1 = Bi[i+1]; 1941 for(j=j0;j<j1;++j) { 1942 J = LUP.Pinv[Bj[j]]; 1943 bj[k] = J; 1944 b[J] = Bv[j]; 1945 ++k; 1946 } 1947 bj.length = k; 1948 sol(L,b,x,bj,xj); 1949 for(j=bj.length-1;j!==-1;--j) b[bj[j]] = 0; 1950 sol(U,x,b,xj,bj); 1951 if(flag) return b; 1952 for(j=xj.length-1;j!==-1;--j) x[xj[j]] = 0; 1953 for(j=bj.length-1;j!==-1;--j) { 1954 J = bj[j]; 1955 Xj[N] = J; 1956 Xv[N] = b[J]; 1957 b[J] = 0; 1958 ++N; 1959 } 1960 Xi[i+1] = N; 1961 } 1962 return [Xi,Xj,Xv]; 1963} 1964 1965numeric.ccsbinop = function ccsbinop(body,setup) { 1966 if(typeof setup === "undefined") setup=''; 1967 return Function('X','Y', 1968 'var Xi = X[0], Xj = X[1], Xv = X[2];\n'+ 1969 'var Yi = Y[0], Yj = Y[1], Yv = Y[2];\n'+ 1970 'var n = Xi.length-1,m = Math.max(numeric.sup(Xj),numeric.sup(Yj))+1;\n'+ 1971 'var Zi = numeric.rep([n+1],0), Zj = [], Zv = [];\n'+ 1972 'var x = numeric.rep([m],0),y = numeric.rep([m],0);\n'+ 1973 'var xk,yk,zk;\n'+ 1974 'var i,j,j0,j1,k,p=0;\n'+ 1975 setup+ 1976 'for(i=0;i<n;++i) {\n'+ 1977 ' j0 = Xi[i]; j1 = Xi[i+1];\n'+ 1978 ' for(j=j0;j!==j1;++j) {\n'+ 1979 ' k = Xj[j];\n'+ 1980 ' x[k] = 1;\n'+ 1981 ' Zj[p] = k;\n'+ 1982 ' ++p;\n'+ 1983 ' }\n'+ 1984 ' j0 = Yi[i]; j1 = Yi[i+1];\n'+ 1985 ' for(j=j0;j!==j1;++j) {\n'+ 1986 ' k = Yj[j];\n'+ 1987 ' y[k] = Yv[j];\n'+ 1988 ' if(x[k] === 0) {\n'+ 1989 ' Zj[p] = k;\n'+ 1990 ' ++p;\n'+ 1991 ' }\n'+ 1992 ' }\n'+ 1993 ' Zi[i+1] = p;\n'+ 1994 ' j0 = Xi[i]; j1 = Xi[i+1];\n'+ 1995 ' for(j=j0;j!==j1;++j) x[Xj[j]] = Xv[j];\n'+ 1996 ' j0 = Zi[i]; j1 = Zi[i+1];\n'+ 1997 ' for(j=j0;j!==j1;++j) {\n'+ 1998 ' k = Zj[j];\n'+ 1999 ' xk = x[k];\n'+ 2000 ' yk = y[k];\n'+ 2001 body+'\n'+ 2002 ' Zv[j] = zk;\n'+ 2003 ' }\n'+ 2004 ' j0 = Xi[i]; j1 = Xi[i+1];\n'+ 2005 ' for(j=j0;j!==j1;++j) x[Xj[j]] = 0;\n'+ 2006 ' j0 = Yi[i]; j1 = Yi[i+1];\n'+ 2007 ' for(j=j0;j!==j1;++j) y[Yj[j]] = 0;\n'+ 2008 '}\n'+ 2009 'return [Zi,Zj,Zv];' 2010 ); 2011}; 2012 2013(function() { 2014 var k,A,B,C; 2015 for(k in numeric.ops2) { 2016 if(isFinite(eval('1'+numeric.ops2[k]+'0'))) A = '[Y[0],Y[1],numeric.'+k+'(X,Y[2])]'; 2017 else A = 'NaN'; 2018 if(isFinite(eval('0'+numeric.ops2[k]+'1'))) B = '[X[0],X[1],numeric.'+k+'(X[2],Y)]'; 2019 else B = 'NaN'; 2020 if(isFinite(eval('1'+numeric.ops2[k]+'0')) && isFinite(eval('0'+numeric.ops2[k]+'1'))) C = 'numeric.ccs'+k+'MM(X,Y)'; 2021 else C = 'NaN'; 2022 numeric['ccs'+k+'MM'] = numeric.ccsbinop('zk = xk '+numeric.ops2[k]+'yk;'); 2023 numeric['ccs'+k] = Function('X','Y', 2024 'if(typeof X === "number") return '+A+';\n'+ 2025 'if(typeof Y === "number") return '+B+';\n'+ 2026 'return '+C+';\n' 2027 ); 2028 } 2029}()); 2030 2031numeric.ccsScatter = function ccsScatter(A) { 2032 var Ai = A[0], Aj = A[1], Av = A[2]; 2033 var n = numeric.sup(Aj)+1,m=Ai.length; 2034 var Ri = numeric.rep([n],0),Rj=Array(m), Rv = Array(m); 2035 var counts = numeric.rep([n],0),i; 2036 for(i=0;i<m;++i) counts[Aj[i]]++; 2037 for(i=0;i<n;++i) Ri[i+1] = Ri[i] + counts[i]; 2038 var ptr = Ri.slice(0),k,Aii; 2039 for(i=0;i<m;++i) { 2040 Aii = Aj[i]; 2041 k = ptr[Aii]; 2042 Rj[k] = Ai[i]; 2043 Rv[k] = Av[i]; 2044 ptr[Aii]=ptr[Aii]+1; 2045 } 2046 return [Ri,Rj,Rv]; 2047} 2048 2049numeric.ccsGather = function ccsGather(A) { 2050 var Ai = A[0], Aj = A[1], Av = A[2]; 2051 var n = Ai.length-1,m = Aj.length; 2052 var Ri = Array(m), Rj = Array(m), Rv = Array(m); 2053 var i,j,j0,j1,p; 2054 p=0; 2055 for(i=0;i<n;++i) { 2056 j0 = Ai[i]; 2057 j1 = Ai[i+1]; 2058 for(j=j0;j!==j1;++j) { 2059 Rj[p] = i; 2060 Ri[p] = Aj[j]; 2061 Rv[p] = Av[j]; 2062 ++p; 2063 } 2064 } 2065 return [Ri,Rj,Rv]; 2066} 2067 2068// The following sparse linear algebra routines are deprecated. 2069 2070numeric.sdim = function dim(A,ret,k) { 2071 if(typeof ret === "undefined") { ret = []; } 2072 if(typeof A !== "object") return ret; 2073 if(typeof k === "undefined") { k=0; } 2074 if(!(k in ret)) { ret[k] = 0; } 2075 if(A.length > ret[k]) ret[k] = A.length; 2076 var i; 2077 for(i in A) { 2078 if(A.hasOwnProperty(i)) dim(A[i],ret,k+1); 2079 } 2080 return ret; 2081}; 2082 2083numeric.sclone = function clone(A,k,n) { 2084 if(typeof k === "undefined") { k=0; } 2085 if(typeof n === "undefined") { n = numeric.sdim(A).length; } 2086 var i,ret = Array(A.length); 2087 if(k === n-1) { 2088 for(i in A) { if(A.hasOwnProperty(i)) ret[i] = A[i]; } 2089 return ret; 2090 } 2091 for(i in A) { 2092 if(A.hasOwnProperty(i)) ret[i] = clone(A[i],k+1,n); 2093 } 2094 return ret; 2095} 2096 2097numeric.sdiag = function diag(d) { 2098 var n = d.length,i,ret = Array(n),i1,i2,i3; 2099 for(i=n-1;i>=1;i-=2) { 2100 i1 = i-1; 2101 ret[i] = []; ret[i][i] = d[i]; 2102 ret[i1] = []; ret[i1][i1] = d[i1]; 2103 } 2104 if(i===0) { ret[0] = []; ret[0][0] = d[i]; } 2105 return ret; 2106} 2107 2108numeric.sidentity = function identity(n) { return numeric.sdiag(numeric.rep([n],1)); } 2109 2110numeric.stranspose = function transpose(A) { 2111 var ret = [], n = A.length, i,j,Ai; 2112 for(i in A) { 2113 if(!(A.hasOwnProperty(i))) continue; 2114 Ai = A[i]; 2115 for(j in Ai) { 2116 if(!(Ai.hasOwnProperty(j))) continue; 2117 if(typeof ret[j] !== "object") { ret[j] = []; } 2118 ret[j][i] = Ai[j]; 2119 } 2120 } 2121 return ret; 2122} 2123 2124numeric.sLUP = function LUP(A,tol) { 2125 throw new Error("The function numeric.sLUP had a bug in it and has been removed. Please use the new numeric.ccsLUP function instead."); 2126}; 2127 2128numeric.sdotMM = function dotMM(A,B) { 2129 var p = A.length, q = B.length, BT = numeric.stranspose(B), r = BT.length, Ai, BTk; 2130 var i,j,k,accum; 2131 var ret = Array(p),reti; 2132 for(i=p-1;i>=0;i--) { 2133 reti = []; 2134 Ai = A[i]; 2135 for(k=r-1;k>=0;k--) { 2136 accum = 0; 2137 BTk = BT[k]; 2138 for(j in Ai) { 2139 if(!(Ai.hasOwnProperty(j))) continue; 2140 if(j in BTk) { accum += Ai[j]*BTk[j]; } 2141 } 2142 if(accum) reti[k] = accum; 2143 } 2144 ret[i] = reti; 2145 } 2146 return ret; 2147} 2148 2149numeric.sdotMV = function dotMV(A,x) { 2150 var p = A.length, Ai, i,j; 2151 var ret = Array(p), accum; 2152 for(i=p-1;i>=0;i--) { 2153 Ai = A[i]; 2154 accum = 0; 2155 for(j in Ai) { 2156 if(!(Ai.hasOwnProperty(j))) continue; 2157 if(x[j]) accum += Ai[j]*x[j]; 2158 } 2159 if(accum) ret[i] = accum; 2160 } 2161 return ret; 2162} 2163 2164numeric.sdotVM = function dotMV(x,A) { 2165 var i,j,Ai,alpha; 2166 var ret = [], accum; 2167 for(i in x) { 2168 if(!x.hasOwnProperty(i)) continue; 2169 Ai = A[i]; 2170 alpha = x[i]; 2171 for(j in Ai) { 2172 if(!Ai.hasOwnProperty(j)) continue; 2173 if(!ret[j]) { ret[j] = 0; } 2174 ret[j] += alpha*Ai[j]; 2175 } 2176 } 2177 return ret; 2178} 2179 2180numeric.sdotVV = function dotVV(x,y) { 2181 var i,ret=0; 2182 for(i in x) { if(x[i] && y[i]) ret+= x[i]*y[i]; } 2183 return ret; 2184} 2185 2186numeric.sdot = function dot(A,B) { 2187 var m = numeric.sdim(A).length, n = numeric.sdim(B).length; 2188 var k = m*1000+n; 2189 switch(k) { 2190 case 0: return A*B; 2191 case 1001: return numeric.sdotVV(A,B); 2192 case 2001: return numeric.sdotMV(A,B); 2193 case 1002: return numeric.sdotVM(A,B); 2194 case 2002: return numeric.sdotMM(A,B); 2195 default: throw new Error('numeric.sdot not implemented for tensors of order '+m+' and '+n); 2196 } 2197} 2198 2199numeric.sscatter = function scatter(V) { 2200 var n = V[0].length, Vij, i, j, m = V.length, A = [], Aj; 2201 for(i=n-1;i>=0;--i) { 2202 if(!V[m-1][i]) continue; 2203 Aj = A; 2204 for(j=0;j<m-2;j++) { 2205 Vij = V[j][i]; 2206 if(!Aj[Vij]) Aj[Vij] = []; 2207 Aj = Aj[Vij]; 2208 } 2209 Aj[V[j][i]] = V[j+1][i]; 2210 } 2211 return A; 2212} 2213 2214numeric.sgather = function gather(A,ret,k) { 2215 if(typeof ret === "undefined") ret = []; 2216 if(typeof k === "undefined") k = []; 2217 var n,i,Ai; 2218 n = k.length; 2219 for(i in A) { 2220 if(A.hasOwnProperty(i)) { 2221 k[n] = parseInt(i); 2222 Ai = A[i]; 2223 if(typeof Ai === "number") { 2224 if(Ai) { 2225 if(ret.length === 0) { 2226 for(i=n+1;i>=0;--i) ret[i] = []; 2227 } 2228 for(i=n;i>=0;--i) ret[i].push(k[i]); 2229 ret[n+1].push(Ai); 2230 } 2231 } else gather(Ai,ret,k); 2232 } 2233 } 2234 if(k.length>n) k.pop(); 2235 return ret; 2236} 2237 2238// 6. Coordinate matrices 2239numeric.cLU = function LU(A) { 2240 var I = A[0], J = A[1], V = A[2]; 2241 var p = I.length, m=0, i,j,k,a,b,c; 2242 for(i=0;i<p;i++) if(I[i]>m) m=I[i]; 2243 m++; 2244 var L = Array(m), U = Array(m), left = numeric.rep([m],Infinity), right = numeric.rep([m],-Infinity); 2245 var Ui, Uj,alpha; 2246 for(k=0;k<p;k++) { 2247 i = I[k]; 2248 j = J[k]; 2249 if(j<left[i]) left[i] = j; 2250 if(j>right[i]) right[i] = j; 2251 } 2252 for(i=0;i<m-1;i++) { if(right[i] > right[i+1]) right[i+1] = right[i]; } 2253 for(i=m-1;i>=1;i--) { if(left[i]<left[i-1]) left[i-1] = left[i]; } 2254 var countL = 0, countU = 0; 2255 for(i=0;i<m;i++) { 2256 U[i] = numeric.rep([right[i]-left[i]+1],0); 2257 L[i] = numeric.rep([i-left[i]],0); 2258 countL += i-left[i]+1; 2259 countU += right[i]-i+1; 2260 } 2261 for(k=0;k<p;k++) { i = I[k]; U[i][J[k]-left[i]] = V[k]; } 2262 for(i=0;i<m-1;i++) { 2263 a = i-left[i]; 2264 Ui = U[i]; 2265 for(j=i+1;left[j]<=i && j<m;j++) { 2266 b = i-left[j]; 2267 c = right[i]-i; 2268 Uj = U[j]; 2269 alpha = Uj[b]/Ui[a]; 2270 if(alpha) { 2271 for(k=1;k<=c;k++) { Uj[k+b] -= alpha*Ui[k+a]; } 2272 L[j][i-left[j]] = alpha; 2273 } 2274 } 2275 } 2276 var Ui = [], Uj = [], Uv = [], Li = [], Lj = [], Lv = []; 2277 var p,q,foo; 2278 p=0; q=0; 2279 for(i=0;i<m;i++) { 2280 a = left[i]; 2281 b = right[i]; 2282 foo = U[i]; 2283 for(j=i;j<=b;j++) { 2284 if(foo[j-a]) { 2285 Ui[p] = i; 2286 Uj[p] = j; 2287 Uv[p] = foo[j-a]; 2288 p++; 2289 } 2290 } 2291 foo = L[i]; 2292 for(j=a;j<i;j++) { 2293 if(foo[j-a]) { 2294 Li[q] = i; 2295 Lj[q] = j; 2296 Lv[q] = foo[j-a]; 2297 q++; 2298 } 2299 } 2300 Li[q] = i; 2301 Lj[q] = i; 2302 Lv[q] = 1; 2303 q++; 2304 } 2305 return {U:[Ui,Uj,Uv], L:[Li,Lj,Lv]}; 2306}; 2307 2308numeric.cLUsolve = function LUsolve(lu,b) { 2309 var L = lu.L, U = lu.U, ret = numeric.clone(b); 2310 var Li = L[0], Lj = L[1], Lv = L[2]; 2311 var Ui = U[0], Uj = U[1], Uv = U[2]; 2312 var p = Ui.length, q = Li.length; 2313 var m = ret.length,i,j,k; 2314 k = 0; 2315 for(i=0;i<m;i++) { 2316 while(Lj[k] < i) { 2317 ret[i] -= Lv[k]*ret[Lj[k]]; 2318 k++; 2319 } 2320 k++; 2321 } 2322 k = p-1; 2323 for(i=m-1;i>=0;i--) { 2324 while(Uj[k] > i) { 2325 ret[i] -= Uv[k]*ret[Uj[k]]; 2326 k--; 2327 } 2328 ret[i] /= Uv[k]; 2329 k--; 2330 } 2331 return ret; 2332}; 2333 2334numeric.cgrid = function grid(n,shape) { 2335 if(typeof n === "number") n = [n,n]; 2336 var ret = numeric.rep(n,-1); 2337 var i,j,count; 2338 if(typeof shape !== "function") { 2339 switch(shape) { 2340 case 'L': 2341 shape = function(i,j) { return (i>=n[0]/2 || j<n[1]/2); } 2342 break; 2343 default: 2344 shape = function(i,j) { return true; }; 2345 break; 2346 } 2347 } 2348 count=0; 2349 for(i=1;i<n[0]-1;i++) for(j=1;j<n[1]-1;j++) 2350 if(shape(i,j)) { 2351 ret[i][j] = count; 2352 count++; 2353 } 2354 return ret; 2355} 2356 2357numeric.cdelsq = function delsq(g) { 2358 var dir = [[-1,0],[0,-1],[0,1],[1,0]]; 2359 var s = numeric.dim(g), m = s[0], n = s[1], i,j,k,p,q; 2360 var Li = [], Lj = [], Lv = []; 2361 for(i=1;i<m-1;i++) for(j=1;j<n-1;j++) { 2362 if(g[i][j]<0) continue; 2363 for(k=0;k<4;k++) { 2364 p = i+dir[k][0]; 2365 q = j+dir[k][1]; 2366 if(g[p][q]<0) continue; 2367 Li.push(g[i][j]); 2368 Lj.push(g[p][q]); 2369 Lv.push(-1); 2370 } 2371 Li.push(g[i][j]); 2372 Lj.push(g[i][j]); 2373 Lv.push(4); 2374 } 2375 return [Li,Lj,Lv]; 2376} 2377 2378numeric.cdotMV = function dotMV(A,x) { 2379 var ret, Ai = A[0], Aj = A[1], Av = A[2],k,p=Ai.length,N; 2380 N=0; 2381 for(k=0;k<p;k++) { if(Ai[k]>N) N = Ai[k]; } 2382 N++; 2383 ret = numeric.rep([N],0); 2384 for(k=0;k<p;k++) { ret[Ai[k]]+=Av[k]*x[Aj[k]]; } 2385 return ret; 2386} 2387 2388// 7. Splines 2389 2390numeric.Spline = function Spline(x,yl,yr,kl,kr) { this.x = x; this.yl = yl; this.yr = yr; this.kl = kl; this.kr = kr; } 2391numeric.Spline.prototype._at = function _at(x1,p) { 2392 var x = this.x; 2393 var yl = this.yl; 2394 var yr = this.yr; 2395 var kl = this.kl; 2396 var kr = this.kr; 2397 var x1,a,b,t; 2398 var add = numeric.add, sub = numeric.sub, mul = numeric.mul; 2399 a = sub(mul(kl[p],x[p+1]-x[p]),sub(yr[p+1],yl[p])); 2400 b = add(mul(kr[p+1],x[p]-x[p+1]),sub(yr[p+1],yl[p])); 2401 t = (x1-x[p])/(x[p+1]-x[p]); 2402 var s = t*(1-t); 2403 return add(add(add(mul(1-t,yl[p]),mul(t,yr[p+1])),mul(a,s*(1-t))),mul(b,s*t)); 2404} 2405numeric.Spline.prototype.at = function at(x0) { 2406 if(typeof x0 === "number") { 2407 var x = this.x; 2408 var n = x.length; 2409 var p,q,mid,floor = Math.floor,a,b,t; 2410 p = 0; 2411 q = n-1; 2412 while(q-p>1) { 2413 mid = floor((p+q)/2); 2414 if(x[mid] <= x0) p = mid; 2415 else q = mid; 2416 } 2417 return this._at(x0,p); 2418 } 2419 var n = x0.length, i, ret = Array(n); 2420 for(i=n-1;i!==-1;--i) ret[i] = this.at(x0[i]); 2421 return ret; 2422} 2423numeric.Spline.prototype.diff = function diff() { 2424 var x = this.x; 2425 var yl = this.yl; 2426 var yr = this.yr; 2427 var kl = this.kl; 2428 var kr = this.kr; 2429 var n = yl.length; 2430 var i,dx,dy; 2431 var zl = kl, zr = kr, pl = Array(n), pr = Array(n); 2432 var add = numeric.add, mul = numeric.mul, div = numeric.div, sub = numeric.sub; 2433 for(i=n-1;i!==-1;--i) { 2434 dx = x[i+1]-x[i]; 2435 dy = sub(yr[i+1],yl[i]); 2436 pl[i] = div(add(mul(dy, 6),mul(kl[i],-4*dx),mul(kr[i+1],-2*dx)),dx*dx); 2437 pr[i+1] = div(add(mul(dy,-6),mul(kl[i], 2*dx),mul(kr[i+1], 4*dx)),dx*dx); 2438 } 2439 return new numeric.Spline(x,zl,zr,pl,pr); 2440} 2441numeric.Spline.prototype.roots = function roots() { 2442 function sqr(x) { return x*x; } 2443 function heval(y0,y1,k0,k1,x) { 2444 var A = k0*2-(y1-y0); 2445 var B = -k1*2+(y1-y0); 2446 var t = (x+1)*0.5; 2447 var s = t*(1-t); 2448 return (1-t)*y0+t*y1+A*s*(1-t)+B*s*t; 2449 } 2450 var ret = []; 2451 var x = this.x, yl = this.yl, yr = this.yr, kl = this.kl, kr = this.kr; 2452 if(typeof yl[0] === "number") { 2453 yl = [yl]; 2454 yr = [yr]; 2455 kl = [kl]; 2456 kr = [kr]; 2457 } 2458 var m = yl.length,n=x.length-1,i,j,k,y,s,t; 2459 var ai,bi,ci,di, ret = Array(m),ri,k0,k1,y0,y1,A,B,D,dx,cx,stops,z0,z1,zm,t0,t1,tm; 2460 var sqrt = Math.sqrt; 2461 for(i=0;i!==m;++i) { 2462 ai = yl[i]; 2463 bi = yr[i]; 2464 ci = kl[i]; 2465 di = kr[i]; 2466 ri = []; 2467 for(j=0;j!==n;j++) { 2468 if(j>0 && bi[j]*ai[j]<0) ri.push(x[j]); 2469 dx = (x[j+1]-x[j]); 2470 cx = x[j]; 2471 y0 = ai[j]; 2472 y1 = bi[j+1]; 2473 k0 = ci[j]/dx; 2474 k1 = di[j+1]/dx; 2475 D = sqr(k0-k1+3*(y0-y1)) + 12*k1*y0; 2476 A = k1+3*y0+2*k0-3*y1; 2477 B = 3*(k1+k0+2*(y0-y1)); 2478 if(D<=0) { 2479 z0 = A/B; 2480 if(z0>x[j] && z0<x[j+1]) stops = [x[j],z0,x[j+1]]; 2481 else stops = [x[j],x[j+1]]; 2482 } else { 2483 z0 = (A-sqrt(D))/B; 2484 z1 = (A+sqrt(D))/B; 2485 stops = [x[j]]; 2486 if(z0>x[j] && z0<x[j+1]) stops.push(z0); 2487 if(z1>x[j] && z1<x[j+1]) stops.push(z1); 2488 stops.push(x[j+1]); 2489 } 2490 t0 = stops[0]; 2491 z0 = this._at(t0,j); 2492 for(k=0;k<stops.length-1;k++) { 2493 t1 = stops[k+1]; 2494 z1 = this._at(t1,j); 2495 if(z0 === 0) { 2496 ri.push(t0); 2497 t0 = t1; 2498 z0 = z1; 2499 continue; 2500 } 2501 if(z1 === 0 || z0*z1>0) { 2502 t0 = t1; 2503 z0 = z1; 2504 continue; 2505 } 2506 var side = 0; 2507 while(1) { 2508 tm = (z0*t1-z1*t0)/(z0-z1); 2509 if(tm <= t0 || tm >= t1) { break; } 2510 zm = this._at(tm,j); 2511 if(zm*z1>0) { 2512 t1 = tm; 2513 z1 = zm; 2514 if(side === -1) z0*=0.5; 2515 side = -1; 2516 } else if(zm*z0>0) { 2517 t0 = tm; 2518 z0 = zm; 2519 if(side === 1) z1*=0.5; 2520 side = 1; 2521 } else break; 2522 } 2523 ri.push(tm); 2524 t0 = stops[k+1]; 2525 z0 = this._at(t0, j); 2526 } 2527 if(z1 === 0) ri.push(t1); 2528 } 2529 ret[i] = ri; 2530 } 2531 if(typeof this.yl[0] === "number") return ret[0]; 2532 return ret; 2533} 2534numeric.spline = function spline(x,y,k1,kn) { 2535 var n = x.length, b = [], dx = [], dy = []; 2536 var i; 2537 var sub = numeric.sub,mul = numeric.mul,add = numeric.add; 2538 for(i=n-2;i>=0;i--) { dx[i] = x[i+1]-x[i]; dy[i] = sub(y[i+1],y[i]); } 2539 if(typeof k1 === "string" || typeof kn === "string") { 2540 k1 = kn = "periodic"; 2541 } 2542 // Build sparse tridiagonal system 2543 var T = [[],[],[]]; 2544 switch(typeof k1) { 2545 case "undefined": 2546 b[0] = mul(3/(dx[0]*dx[0]),dy[0]); 2547 T[0].push(0,0); 2548 T[1].push(0,1); 2549 T[2].push(2/dx[0],1/dx[0]); 2550 break; 2551 case "string": 2552 b[0] = add(mul(3/(dx[n-2]*dx[n-2]),dy[n-2]),mul(3/(dx[0]*dx[0]),dy[0])); 2553 T[0].push(0,0,0); 2554 T[1].push(n-2,0,1); 2555 T[2].push(1/dx[n-2],2/dx[n-2]+2/dx[0],1/dx[0]); 2556 break; 2557 default: 2558 b[0] = k1; 2559 T[0].push(0); 2560 T[1].push(0); 2561 T[2].push(1); 2562 break; 2563 } 2564 for(i=1;i<n-1;i++) { 2565 b[i] = add(mul(3/(dx[i-1]*dx[i-1]),dy[i-1]),mul(3/(dx[i]*dx[i]),dy[i])); 2566 T[0].push(i,i,i); 2567 T[1].push(i-1,i,i+1); 2568 T[2].push(1/dx[i-1],2/dx[i-1]+2/dx[i],1/dx[i]); 2569 } 2570 switch(typeof kn) { 2571 case "undefined": 2572 b[n-1] = mul(3/(dx[n-2]*dx[n-2]),dy[n-2]); 2573 T[0].push(n-1,n-1); 2574 T[1].push(n-2,n-1); 2575 T[2].push(1/dx[n-2],2/dx[n-2]); 2576 break; 2577 case "string": 2578 T[1][T[1].length-1] = 0; 2579 break; 2580 default: 2581 b[n-1] = kn; 2582 T[0].push(n-1); 2583 T[1].push(n-1); 2584 T[2].push(1); 2585 break; 2586 } 2587 if(typeof b[0] !== "number") b = numeric.transpose(b); 2588 else b = [b]; 2589 var k = Array(b.length); 2590 if(typeof k1 === "string") { 2591 for(i=k.length-1;i!==-1;--i) { 2592 k[i] = numeric.ccsLUPSolve(numeric.ccsLUP(numeric.ccsScatter(T)),b[i]); 2593 k[i][n-1] = k[i][0]; 2594 } 2595 } else { 2596 for(i=k.length-1;i!==-1;--i) { 2597 k[i] = numeric.cLUsolve(numeric.cLU(T),b[i]); 2598 } 2599 } 2600 if(typeof y[0] === "number") k = k[0]; 2601 else k = numeric.transpose(k); 2602 return new numeric.Spline(x,y,y,k,k); 2603} 2604 2605// 8. FFT 2606numeric.fftpow2 = function fftpow2(x,y) { 2607 var n = x.length; 2608 if(n === 1) return; 2609 var cos = Math.cos, sin = Math.sin, i,j; 2610 var xe = Array(n/2), ye = Array(n/2), xo = Array(n/2), yo = Array(n/2); 2611 j = n/2; 2612 for(i=n-1;i!==-1;--i) { 2613 --j; 2614 xo[j] = x[i]; 2615 yo[j] = y[i]; 2616 --i; 2617 xe[j] = x[i]; 2618 ye[j] = y[i]; 2619 } 2620 fftpow2(xe,ye); 2621 fftpow2(xo,yo); 2622 j = n/2; 2623 var t,k = (-6.2831853071795864769252867665590057683943387987502116419/n),ci,si; 2624 for(i=n-1;i!==-1;--i) { 2625 --j; 2626 if(j === -1) j = n/2-1; 2627 t = k*i; 2628 ci = cos(t); 2629 si = sin(t); 2630 x[i] = xe[j] + ci*xo[j] - si*yo[j]; 2631 y[i] = ye[j] + ci*yo[j] + si*xo[j]; 2632 } 2633} 2634numeric._ifftpow2 = function _ifftpow2(x,y) { 2635 var n = x.length; 2636 if(n === 1) return; 2637 var cos = Math.cos, sin = Math.sin, i,j; 2638 var xe = Array(n/2), ye = Array(n/2), xo = Array(n/2), yo = Array(n/2); 2639 j = n/2; 2640 for(i=n-1;i!==-1;--i) { 2641 --j; 2642 xo[j] = x[i]; 2643 yo[j] = y[i]; 2644 --i; 2645 xe[j] = x[i]; 2646 ye[j] = y[i]; 2647 } 2648 _ifftpow2(xe,ye); 2649 _ifftpow2(xo,yo); 2650 j = n/2; 2651 var t,k = (6.2831853071795864769252867665590057683943387987502116419/n),ci,si; 2652 for(i=n-1;i!==-1;--i) { 2653 --j; 2654 if(j === -1) j = n/2-1; 2655 t = k*i; 2656 ci = cos(t); 2657 si = sin(t); 2658 x[i] = xe[j] + ci*xo[j] - si*yo[j]; 2659 y[i] = ye[j] + ci*yo[j] + si*xo[j]; 2660 } 2661} 2662numeric.ifftpow2 = function ifftpow2(x,y) { 2663 numeric._ifftpow2(x,y); 2664 numeric.diveq(x,x.length); 2665 numeric.diveq(y,y.length); 2666} 2667numeric.convpow2 = function convpow2(ax,ay,bx,by) { 2668 numeric.fftpow2(ax,ay); 2669 numeric.fftpow2(bx,by); 2670 var i,n = ax.length,axi,bxi,ayi,byi; 2671 for(i=n-1;i!==-1;--i) { 2672 axi = ax[i]; ayi = ay[i]; bxi = bx[i]; byi = by[i]; 2673 ax[i] = axi*bxi-ayi*byi; 2674 ay[i] = axi*byi+ayi*bxi; 2675 } 2676 numeric.ifftpow2(ax,ay); 2677} 2678numeric.T.prototype.fft = function fft() { 2679 var x = this.x, y = this.y; 2680 var n = x.length, log = Math.log, log2 = log(2), 2681 p = Math.ceil(log(2*n-1)/log2), m = Math.pow(2,p); 2682 var cx = numeric.rep([m],0), cy = numeric.rep([m],0), cos = Math.cos, sin = Math.sin; 2683 var k, c = (-3.141592653589793238462643383279502884197169399375105820/n),t; 2684 var a = numeric.rep([m],0), b = numeric.rep([m],0),nhalf = Math.floor(n/2); 2685 for(k=0;k<n;k++) a[k] = x[k]; 2686 if(typeof y !== "undefined") for(k=0;k<n;k++) b[k] = y[k]; 2687 cx[0] = 1; 2688 for(k=1;k<=m/2;k++) { 2689 t = c*k*k; 2690 cx[k] = cos(t); 2691 cy[k] = sin(t); 2692 cx[m-k] = cos(t); 2693 cy[m-k] = sin(t) 2694 } 2695 var X = new numeric.T(a,b), Y = new numeric.T(cx,cy); 2696 X = X.mul(Y); 2697 numeric.convpow2(X.x,X.y,numeric.clone(Y.x),numeric.neg(Y.y)); 2698 X = X.mul(Y); 2699 X.x.length = n; 2700 X.y.length = n; 2701 return X; 2702} 2703numeric.T.prototype.ifft = function ifft() { 2704 var x = this.x, y = this.y; 2705 var n = x.length, log = Math.log, log2 = log(2), 2706 p = Math.ceil(log(2*n-1)/log2), m = Math.pow(2,p); 2707 var cx = numeric.rep([m],0), cy = numeric.rep([m],0), cos = Math.cos, sin = Math.sin; 2708 var k, c = (3.141592653589793238462643383279502884197169399375105820/n),t; 2709 var a = numeric.rep([m],0), b = numeric.rep([m],0),nhalf = Math.floor(n/2); 2710 for(k=0;k<n;k++) a[k] = x[k]; 2711 if(typeof y !== "undefined") for(k=0;k<n;k++) b[k] = y[k]; 2712 cx[0] = 1; 2713 for(k=1;k<=m/2;k++) { 2714 t = c*k*k; 2715 cx[k] = cos(t); 2716 cy[k] = sin(t); 2717 cx[m-k] = cos(t); 2718 cy[m-k] = sin(t) 2719 } 2720 var X = new numeric.T(a,b), Y = new numeric.T(cx,cy); 2721 X = X.mul(Y); 2722 numeric.convpow2(X.x,X.y,numeric.clone(Y.x),numeric.neg(Y.y)); 2723 X = X.mul(Y); 2724 X.x.length = n; 2725 X.y.length = n; 2726 return X.div(n); 2727} 2728 2729//9. Unconstrained optimization 2730numeric.gradient = function gradient(f,x) { 2731 var n = x.length; 2732 var f0 = f(x); 2733 if(isNaN(f0)) throw new Error('gradient: f(x) is a NaN!'); 2734 var max = Math.max; 2735 var i,x0 = numeric.clone(x),f1,f2, J = Array(n); 2736 var div = numeric.div, sub = numeric.sub,errest,roundoff,max = Math.max,eps = 1e-3,abs = Math.abs, min = Math.min; 2737 var t0,t1,t2,it=0,d1,d2,N; 2738 for(i=0;i<n;i++) { 2739 var h = max(1e-6*f0,1e-8); 2740 while(1) { 2741 ++it; 2742 if(it>20) { throw new Error("Numerical gradient fails"); } 2743 x0[i] = x[i]+h; 2744 f1 = f(x0); 2745 x0[i] = x[i]-h; 2746 f2 = f(x0); 2747 x0[i] = x[i]; 2748 if(isNaN(f1) || isNaN(f2)) { h/=16; continue; } 2749 J[i] = (f1-f2)/(2*h); 2750 t0 = x[i]-h; 2751 t1 = x[i]; 2752 t2 = x[i]+h; 2753 d1 = (f1-f0)/h; 2754 d2 = (f0-f2)/h; 2755 N = max(abs(J[i]),abs(f0),abs(f1),abs(f2),abs(t0),abs(t1),abs(t2),1e-8); 2756 errest = min(max(abs(d1-J[i]),abs(d2-J[i]),abs(d1-d2))/N,h/N); 2757 if(errest>eps) { h/=16; } 2758 else break; 2759 } 2760 } 2761 return J; 2762} 2763 2764numeric.uncmin = function uncmin(f,x0,tol,gradient,maxit,callback,options) { 2765 var grad = numeric.gradient; 2766 if(typeof options === "undefined") { options = {}; } 2767 if(typeof tol === "undefined") { tol = 1e-8; } 2768 if(typeof gradient === "undefined") { gradient = function(x) { return grad(f,x); }; } 2769 if(typeof maxit === "undefined") maxit = 1000; 2770 x0 = numeric.clone(x0); 2771 var n = x0.length; 2772 var f0 = f(x0),f1,df0; 2773 if(isNaN(f0)) throw new Error('uncmin: f(x0) is a NaN!'); 2774 var max = Math.max, norm2 = numeric.norm2; 2775 tol = max(tol,numeric.epsilon); 2776 var step,g0,g1,H1 = options.Hinv || numeric.identity(n); 2777 var dot = numeric.dot, inv = numeric.inv, sub = numeric.sub, add = numeric.add, ten = numeric.tensor, div = numeric.div, mul = numeric.mul; 2778 var all = numeric.all, isfinite = numeric.isFinite, neg = numeric.neg; 2779 var it=0,i,s,x1,y,Hy,Hs,ys,i0,t,nstep,t1,t2; 2780 var msg = ""; 2781 g0 = gradient(x0); 2782 while(it<maxit) { 2783 if(typeof callback === "function") { if(callback(it,x0,f0,g0,H1)) { msg = "Callback returned true"; break; } } 2784 if(!all(isfinite(g0))) { msg = "Gradient has Infinity or NaN"; break; } 2785 step = neg(dot(H1,g0)); 2786 if(!all(isfinite(step))) { msg = "Search direction has Infinity or NaN"; break; } 2787 nstep = norm2(step); 2788 if(nstep < tol) { msg="Newton step smaller than tol"; break; } 2789 t = 1; 2790 df0 = dot(g0,step); 2791 // line search 2792 x1 = x0; 2793 while(it < maxit) { 2794 if(t*nstep < tol) { break; } 2795 s = mul(step,t); 2796 x1 = add(x0,s); 2797 f1 = f(x1); 2798 if(f1-f0 >= 0.1*t*df0 || isNaN(f1)) { 2799 t *= 0.5; 2800 ++it; 2801 continue; 2802 } 2803 break; 2804 } 2805 if(t*nstep < tol) { msg = "Line search step size smaller than tol"; break; } 2806 if(it === maxit) { msg = "maxit reached during line search"; break; } 2807 g1 = gradient(x1); 2808 y = sub(g1,g0); 2809 ys = dot(y,s); 2810 Hy = dot(H1,y); 2811 H1 = sub(add(H1, 2812 mul( 2813 (ys+dot(y,Hy))/(ys*ys), 2814 ten(s,s) )), 2815 div(add(ten(Hy,s),ten(s,Hy)),ys)); 2816 x0 = x1; 2817 f0 = f1; 2818 g0 = g1; 2819 ++it; 2820 } 2821 return {solution: x0, f: f0, gradient: g0, invHessian: H1, iterations:it, message: msg}; 2822} 2823 2824// 10. Ode solver (Dormand-Prince) 2825numeric.Dopri = function Dopri(x,y,f,ymid,iterations,msg,events) { 2826 this.x = x; 2827 this.y = y; 2828 this.f = f; 2829 this.ymid = ymid; 2830 this.iterations = iterations; 2831 this.events = events; 2832 this.message = msg; 2833} 2834numeric.Dopri.prototype._at = function _at(xi,j) { 2835 function sqr(x) { return x*x; } 2836 var sol = this; 2837 var xs = sol.x; 2838 var ys = sol.y; 2839 var k1 = sol.f; 2840 var ymid = sol.ymid; 2841 var n = xs.length; 2842 var x0,x1,xh,y0,y1,yh,xi; 2843 var floor = Math.floor,h; 2844 var c = 0.5; 2845 var add = numeric.add, mul = numeric.mul,sub = numeric.sub, p,q,w; 2846 x0 = xs[j]; 2847 x1 = xs[j+1]; 2848 y0 = ys[j]; 2849 y1 = ys[j+1]; 2850 h = x1-x0; 2851 xh = x0+c*h; 2852 yh = ymid[j]; 2853 p = sub(k1[j ],mul(y0,1/(x0-xh)+2/(x0-x1))); 2854 q = sub(k1[j+1],mul(y1,1/(x1-xh)+2/(x1-x0))); 2855 w = [sqr(xi - x1) * (xi - xh) / sqr(x0 - x1) / (x0 - xh), 2856 sqr(xi - x0) * sqr(xi - x1) / sqr(x0 - xh) / sqr(x1 - xh), 2857 sqr(xi - x0) * (xi - xh) / sqr(x1 - x0) / (x1 - xh), 2858 (xi - x0) * sqr(xi - x1) * (xi - xh) / sqr(x0-x1) / (x0 - xh), 2859 (xi - x1) * sqr(xi - x0) * (xi - xh) / sqr(x0-x1) / (x1 - xh)]; 2860 return add(add(add(add(mul(y0,w[0]), 2861 mul(yh,w[1])), 2862 mul(y1,w[2])), 2863 mul( p,w[3])), 2864 mul( q,w[4])); 2865} 2866numeric.Dopri.prototype.at = function at(x) { 2867 var i,j,k,floor = Math.floor; 2868 if(typeof x !== "number") { 2869 var n = x.length, ret = Array(n); 2870 for(i=n-1;i!==-1;--i) { 2871 ret[i] = this.at(x[i]); 2872 } 2873 return ret; 2874 } 2875 var x0 = this.x; 2876 i = 0; j = x0.length-1; 2877 while(j-i>1) { 2878 k = floor(0.5*(i+j)); 2879 if(x0[k] <= x) i = k; 2880 else j = k; 2881 } 2882 return this._at(x,i); 2883} 2884 2885numeric.dopri = function dopri(x0,x1,y0,f,tol,maxit,event) { 2886 if(typeof tol === "undefined") { tol = 1e-6; } 2887 if(typeof maxit === "undefined") { maxit = 1000; } 2888 var xs = [x0], ys = [y0], k1 = [f(x0,y0)], k2,k3,k4,k5,k6,k7, ymid = []; 2889 var A2 = 1/5; 2890 var A3 = [3/40,9/40]; 2891 var A4 = [44/45,-56/15,32/9]; 2892 var A5 = [19372/6561,-25360/2187,64448/6561,-212/729]; 2893 var A6 = [9017/3168,-355/33,46732/5247,49/176,-5103/18656]; 2894 var b = [35/384,0,500/1113,125/192,-2187/6784,11/84]; 2895 var bm = [0.5*6025192743/30085553152, 2896 0, 2897 0.5*51252292925/65400821598, 2898 0.5*-2691868925/45128329728, 2899 0.5*187940372067/1594534317056, 2900 0.5*-1776094331/19743644256, 2901 0.5*11237099/235043384]; 2902 var c = [1/5,3/10,4/5,8/9,1,1]; 2903 var e = [-71/57600,0,71/16695,-71/1920,17253/339200,-22/525,1/40]; 2904 var i = 0,er,j; 2905 var h = (x1-x0)/10; 2906 var it = 0; 2907 var add = numeric.add, mul = numeric.mul, y1,erinf; 2908 var max = Math.max, min = Math.min, abs = Math.abs, norminf = numeric.norminf,pow = Math.pow; 2909 var any = numeric.any, lt = numeric.lt, and = numeric.and, sub = numeric.sub; 2910 var e0, e1, ev; 2911 var ret = new numeric.Dopri(xs,ys,k1,ymid,-1,""); 2912 if(typeof event === "function") e0 = event(x0,y0); 2913 while(x0<x1 && it<maxit) { 2914 ++it; 2915 if(x0+h>x1) h = x1-x0; 2916 k2 = f(x0+c[0]*h, add(y0,mul( A2*h,k1[i]))); 2917 k3 = f(x0+c[1]*h, add(add(y0,mul(A3[0]*h,k1[i])),mul(A3[1]*h,k2))); 2918 k4 = f(x0+c[2]*h, add(add(add(y0,mul(A4[0]*h,k1[i])),mul(A4[1]*h,k2)),mul(A4[2]*h,k3))); 2919 k5 = f(x0+c[3]*h, add(add(add(add(y0,mul(A5[0]*h,k1[i])),mul(A5[1]*h,k2)),mul(A5[2]*h,k3)),mul(A5[3]*h,k4))); 2920 k6 = f(x0+c[4]*h,add(add(add(add(add(y0,mul(A6[0]*h,k1[i])),mul(A6[1]*h,k2)),mul(A6[2]*h,k3)),mul(A6[3]*h,k4)),mul(A6[4]*h,k5))); 2921 y1 = add(add(add(add(add(y0,mul(k1[i],h*b[0])),mul(k3,h*b[2])),mul(k4,h*b[3])),mul(k5,h*b[4])),mul(k6,h*b[5])); 2922 k7 = f(x0+h,y1); 2923 er = add(add(add(add(add(mul(k1[i],h*e[0]),mul(k3,h*e[2])),mul(k4,h*e[3])),mul(k5,h*e[4])),mul(k6,h*e[5])),mul(k7,h*e[6])); 2924 if(typeof er === "number") erinf = abs(er); 2925 else erinf = norminf(er); 2926 if(erinf > tol) { // reject 2927 h = 0.2*h*pow(tol/erinf,0.25); 2928 if(x0+h === x0) { 2929 ret.msg = "Step size became too small"; 2930 break; 2931 } 2932 continue; 2933 } 2934 ymid[i] = add(add(add(add(add(add(y0, 2935 mul(k1[i],h*bm[0])), 2936 mul(k3 ,h*bm[2])), 2937 mul(k4 ,h*bm[3])), 2938 mul(k5 ,h*bm[4])), 2939 mul(k6 ,h*bm[5])), 2940 mul(k7 ,h*bm[6])); 2941 ++i; 2942 xs[i] = x0+h; 2943 ys[i] = y1; 2944 k1[i] = k7; 2945 if(typeof event === "function") { 2946 var yi,xl = x0,xr = x0+0.5*h,xi; 2947 e1 = event(xr,ymid[i-1]); 2948 ev = and(lt(e0,0),lt(0,e1)); 2949 if(!any(ev)) { xl = xr; xr = x0+h; e0 = e1; e1 = event(xr,y1); ev = and(lt(e0,0),lt(0,e1)); } 2950 if(any(ev)) { 2951 var xc, yc, en,ei; 2952 var side=0, sl = 1.0, sr = 1.0; 2953 while(1) { 2954 if(typeof e0 === "number") xi = (sr*e1*xl-sl*e0*xr)/(sr*e1-sl*e0); 2955 else { 2956 xi = xr; 2957 for(j=e0.length-1;j!==-1;--j) { 2958 if(e0[j]<0 && e1[j]>0) xi = min(xi,(sr*e1[j]*xl-sl*e0[j]*xr)/(sr*e1[j]-sl*e0[j])); 2959 } 2960 } 2961 if(xi <= xl || xi >= xr) break; 2962 yi = ret._at(xi, i-1); 2963 ei = event(xi,yi); 2964 en = and(lt(e0,0),lt(0,ei)); 2965 if(any(en)) { 2966 xr = xi; 2967 e1 = ei; 2968 ev = en; 2969 sr = 1.0; 2970 if(side === -1) sl *= 0.5; 2971 else sl = 1.0; 2972 side = -1; 2973 } else { 2974 xl = xi; 2975 e0 = ei; 2976 sl = 1.0; 2977 if(side === 1) sr *= 0.5; 2978 else sr = 1.0; 2979 side = 1; 2980 } 2981 } 2982 y1 = ret._at(0.5*(x0+xi),i-1); 2983 ret.f[i] = f(xi,yi); 2984 ret.x[i] = xi; 2985 ret.y[i] = yi; 2986 ret.ymid[i-1] = y1; 2987 ret.events = ev; 2988 ret.iterations = it; 2989 return ret; 2990 } 2991 } 2992 x0 += h; 2993 y0 = y1; 2994 e0 = e1; 2995 h = min(0.8*h*pow(tol/erinf,0.25),4*h); 2996 } 2997 ret.iterations = it; 2998 return ret; 2999} 3000 3001// 11. Ax = b 3002numeric.LU = function(A, fast) { 3003 fast = fast || false; 3004 3005 var abs = Math.abs; 3006 var i, j, k, absAjk, Akk, Ak, Pk, Ai; 3007 var max; 3008 var n = A.length, n1 = n-1; 3009 var P = new Array(n); 3010 if(!fast) A = numeric.clone(A); 3011 3012 for (k = 0; k < n; ++k) { 3013 Pk = k; 3014 Ak = A[k]; 3015 max = abs(Ak[k]); 3016 for (j = k + 1; j < n; ++j) { 3017 absAjk = abs(A[j][k]); 3018 if (max < absAjk) { 3019 max = absAjk; 3020 Pk = j; 3021 } 3022 } 3023 P[k] = Pk; 3024 3025 if (Pk != k) { 3026 A[k] = A[Pk]; 3027 A[Pk] = Ak; 3028 Ak = A[k]; 3029 } 3030 3031 Akk = Ak[k]; 3032 3033 for (i = k + 1; i < n; ++i) { 3034 A[i][k] /= Akk; 3035 } 3036 3037 for (i = k + 1; i < n; ++i) { 3038 Ai = A[i]; 3039 for (j = k + 1; j < n1; ++j) { 3040 Ai[j] -= Ai[k] * Ak[j]; 3041 ++j; 3042 Ai[j] -= Ai[k] * Ak[j]; 3043 } 3044 if(j===n1) Ai[j] -= Ai[k] * Ak[j]; 3045 } 3046 } 3047 3048 return { 3049 LU: A, 3050 P: P 3051 }; 3052} 3053 3054numeric.LUsolve = function LUsolve(LUP, b) { 3055 var i, j; 3056 var LU = LUP.LU; 3057 var n = LU.length; 3058 var x = numeric.clone(b); 3059 var P = LUP.P; 3060 var Pi, LUi, LUii, tmp; 3061 3062 for (i=n-1;i!==-1;--i) x[i] = b[i]; 3063 for (i = 0; i < n; ++i) { 3064 Pi = P[i]; 3065 if (P[i] !== i) { 3066 tmp = x[i]; 3067 x[i] = x[Pi]; 3068 x[Pi] = tmp; 3069 } 3070 3071 LUi = LU[i]; 3072 for (j = 0; j < i; ++j) { 3073 x[i] -= x[j] * LUi[j]; 3074 } 3075 } 3076 3077 for (i = n - 1; i >= 0; --i) { 3078 LUi = LU[i]; 3079 for (j = i + 1; j < n; ++j) { 3080 x[i] -= x[j] * LUi[j]; 3081 } 3082 3083 x[i] /= LUi[i]; 3084 } 3085 3086 return x; 3087} 3088 3089numeric.solve = function solve(A,b,fast) { return numeric.LUsolve(numeric.LU(A,fast), b); } 3090 3091// 12. Linear programming 3092numeric.echelonize = function echelonize(A) { 3093 var s = numeric.dim(A), m = s[0], n = s[1]; 3094 var I = numeric.identity(m); 3095 var P = Array(m); 3096 var i,j,k,l,Ai,Ii,Z,a; 3097 var abs = Math.abs; 3098 var diveq = numeric.diveq; 3099 A = numeric.clone(A); 3100 for(i=0;i<m;++i) { 3101 k = 0; 3102 Ai = A[i]; 3103 Ii = I[i]; 3104 for(j=1;j<n;++j) if(abs(Ai[k])<abs(Ai[j])) k=j; 3105 P[i] = k; 3106 diveq(Ii,Ai[k]); 3107 diveq(Ai,Ai[k]); 3108 for(j=0;j<m;++j) if(j!==i) { 3109 Z = A[j]; a = Z[k]; 3110 for(l=n-1;l!==-1;--l) Z[l] -= Ai[l]*a; 3111 Z = I[j]; 3112 for(l=m-1;l!==-1;--l) Z[l] -= Ii[l]*a; 3113 } 3114 } 3115 return {I:I, A:A, P:P}; 3116} 3117 3118numeric.__solveLP = function __solveLP(c,A,b,tol,maxit,x,flag) { 3119 var sum = numeric.sum, log = numeric.log, mul = numeric.mul, sub = numeric.sub, dot = numeric.dot, div = numeric.div, add = numeric.add; 3120 var m = c.length, n = b.length,y; 3121 var unbounded = false, cb,i0=0; 3122 var alpha = 1.0; 3123 var f0,df0,AT = numeric.transpose(A), svd = numeric.svd,transpose = numeric.transpose,leq = numeric.leq, sqrt = Math.sqrt, abs = Math.abs; 3124 var muleq = numeric.muleq; 3125 var norm = numeric.norminf, any = numeric.any,min = Math.min; 3126 var all = numeric.all, gt = numeric.gt; 3127 var p = Array(m), A0 = Array(n),e=numeric.rep([n],1), H; 3128 var solve = numeric.solve, z = sub(b,dot(A,x)),count; 3129 var dotcc = dot(c,c); 3130 var g; 3131 for(count=i0;count<maxit;++count) { 3132 var i,j,d; 3133 for(i=n-1;i!==-1;--i) A0[i] = div(A[i],z[i]); 3134 var A1 = transpose(A0); 3135 for(i=m-1;i!==-1;--i) p[i] = (/*x[i]+*/sum(A1[i])); 3136 alpha = 0.25*abs(dotcc/dot(c,p)); 3137 var a1 = 100*sqrt(dotcc/dot(p,p)); 3138 if(!isFinite(alpha) || alpha>a1) alpha = a1; 3139 g = add(c,mul(alpha,p)); 3140 H = dot(A1,A0); 3141 for(i=m-1;i!==-1;--i) H[i][i] += 1; 3142 d = solve(H,div(g,alpha),true); 3143 var t0 = div(z,dot(A,d)); 3144 var t = 1.0; 3145 for(i=n-1;i!==-1;--i) if(t0[i]<0) t = min(t,-0.999*t0[i]); 3146 y = sub(x,mul(d,t)); 3147 z = sub(b,dot(A,y)); 3148 if(!all(gt(z,0))) return { solution: x, message: "", iterations: count }; 3149 x = y; 3150 if(alpha<tol) return { solution: y, message: "", iterations: count }; 3151 if(flag) { 3152 var s = dot(c,g), Ag = dot(A,g); 3153 unbounded = true; 3154 for(i=n-1;i!==-1;--i) if(s*Ag[i]<0) { unbounded = false; break; } 3155 } else { 3156 if(x[m-1]>=0) unbounded = false; 3157 else unbounded = true; 3158 } 3159 if(unbounded) return { solution: y, message: "Unbounded", iterations: count }; 3160 } 3161 return { solution: x, message: "maximum iteration count exceeded", iterations:count }; 3162} 3163 3164numeric._solveLP = function _solveLP(c,A,b,tol,maxit) { 3165 var m = c.length, n = b.length,y; 3166 var sum = numeric.sum, log = numeric.log, mul = numeric.mul, sub = numeric.sub, dot = numeric.dot, div = numeric.div, add = numeric.add; 3167 var c0 = numeric.rep([m],0).concat([1]); 3168 var J = numeric.rep([n,1],-1); 3169 var A0 = numeric.blockMatrix([[A , J ]]); 3170 var b0 = b; 3171 var y = numeric.rep([m],0).concat(Math.max(0,numeric.sup(numeric.neg(b)))+1); 3172 var x0 = numeric.__solveLP(c0,A0,b0,tol,maxit,y,false); 3173 var x = numeric.clone(x0.solution); 3174 x.length = m; 3175 var foo = numeric.inf(sub(b,dot(A,x))); 3176 if(foo<0) { return { solution: NaN, message: "Infeasible", iterations: x0.iterations }; } 3177 var ret = numeric.__solveLP(c, A, b, tol, maxit-x0.iterations, x, true); 3178 ret.iterations += x0.iterations; 3179 return ret; 3180}; 3181 3182numeric.solveLP = function solveLP(c,A,b,Aeq,beq,tol,maxit) { 3183 if(typeof maxit === "undefined") maxit = 1000; 3184 if(typeof tol === "undefined") tol = numeric.epsilon; 3185 if(typeof Aeq === "undefined") return numeric._solveLP(c,A,b,tol,maxit); 3186 var m = Aeq.length, n = Aeq[0].length, o = A.length; 3187 var B = numeric.echelonize(Aeq); 3188 var flags = numeric.rep([n],0); 3189 var P = B.P; 3190 var Q = []; 3191 var i; 3192 for(i=P.length-1;i!==-1;--i) flags[P[i]] = 1; 3193 for(i=n-1;i!==-1;--i) if(flags[i]===0) Q.push(i); 3194 var g = numeric.getRange; 3195 var I = numeric.linspace(0,m-1), J = numeric.linspace(0,o-1); 3196 var Aeq2 = g(Aeq,I,Q), A1 = g(A,J,P), A2 = g(A,J,Q), dot = numeric.dot, sub = numeric.sub; 3197 var A3 = dot(A1,B.I); 3198 var A4 = sub(A2,dot(A3,Aeq2)), b4 = sub(b,dot(A3,beq)); 3199 var c1 = Array(P.length), c2 = Array(Q.length); 3200 for(i=P.length-1;i!==-1;--i) c1[i] = c[P[i]]; 3201 for(i=Q.length-1;i!==-1;--i) c2[i] = c[Q[i]]; 3202 var c4 = sub(c2,dot(c1,dot(B.I,Aeq2))); 3203 var S = numeric._solveLP(c4,A4,b4,tol,maxit); 3204 var x2 = S.solution; 3205 if(x2!==x2) return S; 3206 var x1 = dot(B.I,sub(beq,dot(Aeq2,x2))); 3207 var x = Array(c.length); 3208 for(i=P.length-1;i!==-1;--i) x[P[i]] = x1[i]; 3209 for(i=Q.length-1;i!==-1;--i) x[Q[i]] = x2[i]; 3210 return { solution: x, message:S.message, iterations: S.iterations }; 3211} 3212 3213numeric.MPStoLP = function MPStoLP(MPS) { 3214 if(MPS instanceof String) { MPS.split('\n'); } 3215 var state = 0; 3216 var states = ['Initial state','NAME','ROWS','COLUMNS','RHS','BOUNDS','ENDATA']; 3217 var n = MPS.length; 3218 var i,j,z,N=0,rows = {}, sign = [], rl = 0, vars = {}, nv = 0; 3219 var name; 3220 var c = [], A = [], b = []; 3221 function err(e) { throw new Error('MPStoLP: '+e+'\nLine '+i+': '+MPS[i]+'\nCurrent state: '+states[state]+'\n'); } 3222 for(i=0;i<n;++i) { 3223 z = MPS[i]; 3224 var w0 = z.match(/\S*/g); 3225 var w = []; 3226 for(j=0;j<w0.length;++j) if(w0[j]!=="") w.push(w0[j]); 3227 if(w.length === 0) continue; 3228 for(j=0;j<states.length;++j) if(z.substr(0,states[j].length) === states[j]) break; 3229 if(j<states.length) { 3230 state = j; 3231 if(j===1) { name = w[1]; } 3232 if(j===6) return { name:name, c:c, A:numeric.transpose(A), b:b, rows:rows, vars:vars }; 3233 continue; 3234 } 3235 switch(state) { 3236 case 0: case 1: err('Unexpected line'); 3237 case 2: 3238 switch(w[0]) { 3239 case 'N': if(N===0) N = w[1]; else err('Two or more N rows'); break; 3240 case 'L': rows[w[1]] = rl; sign[rl] = 1; b[rl] = 0; ++rl; break; 3241 case 'G': rows[w[1]] = rl; sign[rl] = -1;b[rl] = 0; ++rl; break; 3242 case 'E': rows[w[1]] = rl; sign[rl] = 0;b[rl] = 0; ++rl; break; 3243 default: err('Parse error '+numeric.prettyPrint(w)); 3244 } 3245 break; 3246 case 3: 3247 if(!vars.hasOwnProperty(w[0])) { vars[w[0]] = nv; c[nv] = 0; A[nv] = numeric.rep([rl],0); ++nv; } 3248 var p = vars[w[0]]; 3249 for(j=1;j<w.length;j+=2) { 3250 if(w[j] === N) { c[p] = parseFloat(w[j+1]); continue; } 3251 var q = rows[w[j]]; 3252 A[p][q] = (sign[q]<0?-1:1)*parseFloat(w[j+1]); 3253 } 3254 break; 3255 case 4: 3256 for(j=1;j<w.length;j+=2) b[rows[w[j]]] = (sign[rows[w[j]]]<0?-1:1)*parseFloat(w[j+1]); 3257 break; 3258 case 5: /*FIXME*/ break; 3259 case 6: err('Internal error'); 3260 } 3261 } 3262 err('Reached end of file without ENDATA'); 3263}
vendor: 26,806 bytes, lines 3263-4129
3263 3264// seedrandom.js version 2.0. 3265// Author: David Bau 4/2/2011 3266// 3267// Defines a method Math.seedrandom() that, when called, substitutes 3268// an explicitly seeded RC4-based algorithm for Math.random(). Also 3269// supports automatic seeding from local or network sources of entropy. 3270// 3271// Usage: 3272// 3273// <script src=http://davidbau.com/encode/seedrandom-min.js></script> 3274// 3275// Math.seedrandom('yipee'); Sets Math.random to a function that is 3276// initialized using the given explicit seed. 3277// 3278// Math.seedrandom(); Sets Math.random to a function that is 3279// seeded using the current time, dom state, 3280// and other accumulated local entropy. 3281// The generated seed string is returned. 3282// 3283// Math.seedrandom('yowza', true); 3284// Seeds using the given explicit seed mixed 3285// together with accumulated entropy. 3286// 3287// <script src="http://bit.ly/srandom-512"></script> 3288// Seeds using physical random bits downloaded 3289// from random.org. 3290// 3291// <script src="https://jsonlib.appspot.com/urandom?callback=Math.seedrandom"> 3292// </script> Seeds using urandom bits from call.jsonlib.com, 3293// which is faster than random.org. 3294// 3295// Examples: 3296// 3297// Math.seedrandom("hello"); // Use "hello" as the seed. 3298// document.write(Math.random()); // Always 0.5463663768140734 3299// document.write(Math.random()); // Always 0.43973793770592234 3300// var rng1 = Math.random; // Remember the current prng. 3301// 3302// var autoseed = Math.seedrandom(); // New prng with an automatic seed. 3303// document.write(Math.random()); // Pretty much unpredictable. 3304// 3305// Math.random = rng1; // Continue "hello" prng sequence. 3306// document.write(Math.random()); // Always 0.554769432473455 3307// 3308// Math.seedrandom(autoseed); // Restart at the previous seed. 3309// document.write(Math.random()); // Repeat the 'unpredictable' value. 3310// 3311// Notes: 3312// 3313// Each time seedrandom('arg') is called, entropy from the passed seed 3314// is accumulated in a pool to help generate future seeds for the 3315// zero-argument form of Math.seedrandom, so entropy can be injected over 3316// time by calling seedrandom with explicit data repeatedly. 3317// 3318// On speed - This javascript implementation of Math.random() is about 3319// 3-10x slower than the built-in Math.random() because it is not native 3320// code, but this is typically fast enough anyway. Seeding is more expensive, 3321// especially if you use auto-seeding. Some details (timings on Chrome 4): 3322// 3323// Our Math.random() - avg less than 0.002 milliseconds per call 3324// seedrandom('explicit') - avg less than 0.5 milliseconds per call 3325// seedrandom('explicit', true) - avg less than 2 milliseconds per call 3326// seedrandom() - avg about 38 milliseconds per call 3327// 3328// LICENSE (BSD): 3329// 3330// Copyright 2010 David Bau, all rights reserved. 3331// 3332// Redistribution and use in source and binary forms, with or without 3333// modification, are permitted provided that the following conditions are met: 3334// 3335// 1. Redistributions of source code must retain the above copyright 3336// notice, this list of conditions and the following disclaimer. 3337// 3338// 2. Redistributions in binary form must reproduce the above copyright 3339// notice, this list of conditions and the following disclaimer in the 3340// documentation and/or other materials provided with the distribution. 3341// 3342// 3. Neither the name of this module nor the names of its contributors may 3343// be used to endorse or promote products derived from this software 3344// without specific prior written permission. 3345// 3346// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 3347// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 3348// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 3349// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 3350// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 3351// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 3352// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 3353// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 3354// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 3355// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 3356// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 3357// 3358/** 3359 * All code is in an anonymous closure to keep the global namespace clean. 3360 * 3361 * @param {number=} overflow 3362 * @param {number=} startdenom 3363 */ 3364 3365// Patched by Seb so that seedrandom.js does not pollute the Math object. 3366// My tests suggest that doing Math.trouble = 1 makes Math lookups about 5% 3367// slower. 3368numeric.seedrandom = { pow:Math.pow, random:Math.random }; 3369 3370(function (pool, math, width, chunks, significance, overflow, startdenom) { 3371 3372 3373// 3374// seedrandom() 3375// This is the seedrandom function described above. 3376// 3377math['seedrandom'] = function seedrandom(seed, use_entropy) { 3378 var key = []; 3379 var arc4; 3380 3381 // Flatten the seed string or build one from local entropy if needed. 3382 seed = mixkey(flatten( 3383 use_entropy ? [seed, pool] : 3384 arguments.length ? seed : 3385 [new Date().getTime(), pool, window], 3), key); 3386 3387 // Use the seed to initialize an ARC4 generator. 3388 arc4 = new ARC4(key); 3389 3390 // Mix the randomness into accumulated entropy. 3391 mixkey(arc4.S, pool); 3392 3393 // Override Math.random 3394 3395 // This function returns a random double in [0, 1) that contains 3396 // randomness in every bit of the mantissa of the IEEE 754 value. 3397 3398 math['random'] = function random() { // Closure to return a random double: 3399 var n = arc4.g(chunks); // Start with a numerator n < 2 ^ 48 3400 var d = startdenom; // and denominator d = 2 ^ 48. 3401 var x = 0; // and no 'extra last byte'. 3402 while (n < significance) { // Fill up all significant digits by 3403 n = (n + x) * width; // shifting numerator and 3404 d *= width; // denominator and generating a 3405 x = arc4.g(1); // new least-significant-byte. 3406 } 3407 while (n >= overflow) { // To avoid rounding up, before adding 3408 n /= 2; // last byte, shift everything 3409 d /= 2; // right using integer math until 3410 x >>>= 1; // we have exactly the desired bits. 3411 } 3412 return (n + x) / d; // Form the number within [0, 1). 3413 }; 3414 3415 // Return the seed that was used 3416 return seed; 3417}; 3418 3419// 3420// ARC4 3421// 3422// An ARC4 implementation. The constructor takes a key in the form of 3423// an array of at most (width) integers that should be 0 <= x < (width). 3424// 3425// The g(count) method returns a pseudorandom integer that concatenates 3426// the next (count) outputs from ARC4. Its return value is a number x 3427// that is in the range 0 <= x < (width ^ count). 3428// 3429/** @constructor */ 3430function ARC4(key) { 3431 var t, u, me = this, keylen = key.length; 3432 var i = 0, j = me.i = me.j = me.m = 0; 3433 me.S = []; 3434 me.c = []; 3435 3436 // The empty key [] is treated as [0]. 3437 if (!keylen) { key = [keylen++]; } 3438 3439 // Set up S using the standard key scheduling algorithm. 3440 while (i < width) { me.S[i] = i++; } 3441 for (i = 0; i < width; i++) { 3442 t = me.S[i]; 3443 j = lowbits(j + t + key[i % keylen]); 3444 u = me.S[j]; 3445 me.S[i] = u; 3446 me.S[j] = t; 3447 } 3448 3449 // The "g" method returns the next (count) outputs as one number. 3450 me.g = function getnext(count) { 3451 var s = me.S; 3452 var i = lowbits(me.i + 1); var t = s[i]; 3453 var j = lowbits(me.j + t); var u = s[j]; 3454 s[i] = u; 3455 s[j] = t; 3456 var r = s[lowbits(t + u)]; 3457 while (--count) { 3458 i = lowbits(i + 1); t = s[i]; 3459 j = lowbits(j + t); u = s[j]; 3460 s[i] = u; 3461 s[j] = t; 3462 r = r * width + s[lowbits(t + u)]; 3463 } 3464 me.i = i; 3465 me.j = j; 3466 return r; 3467 }; 3468 // For robust unpredictability discard an initial batch of values. 3469 // See http://www.rsa.com/rsalabs/node.asp?id=2009 3470 me.g(width); 3471} 3472 3473// 3474// flatten() 3475// Converts an object tree to nested arrays of strings. 3476// 3477/** @param {Object=} result 3478 * @param {string=} prop 3479 * @param {string=} typ */ 3480function flatten(obj, depth, result, prop, typ) { 3481 result = []; 3482 typ = typeof(obj); 3483 if (depth && typ == 'object') { 3484 for (prop in obj) { 3485 if (prop.indexOf('S') < 5) { // Avoid FF3 bug (local/sessionStorage) 3486 try { result.push(flatten(obj[prop], depth - 1)); } catch (e) {} 3487 } 3488 } 3489 } 3490 return (result.length ? result : obj + (typ != 'string' ? '\0' : '')); 3491} 3492 3493// 3494// mixkey() 3495// Mixes a string seed into a key that is an array of integers, and 3496// returns a shortened string seed that is equivalent to the result key. 3497// 3498/** @param {number=} smear 3499 * @param {number=} j */ 3500function mixkey(seed, key, smear, j) { 3501 seed += ''; // Ensure the seed is a string 3502 smear = 0; 3503 for (j = 0; j < seed.length; j++) { 3504 key[lowbits(j)] = 3505 lowbits((smear ^= key[lowbits(j)] * 19) + seed.charCodeAt(j)); 3506 } 3507 seed = ''; 3508 for (j in key) { seed += String.fromCharCode(key[j]); } 3509 return seed; 3510} 3511 3512// 3513// lowbits() 3514// A quick "n mod width" for width a power of 2. 3515// 3516function lowbits(n) { return n & (width - 1); } 3517 3518// 3519// The following constants are related to IEEE 754 limits. 3520// 3521startdenom = math.pow(width, chunks); 3522significance = math.pow(2, significance); 3523overflow = significance * 2; 3524 3525// 3526// When seedrandom.js is loaded, we immediately mix a few bits 3527// from the built-in RNG into the entropy pool. Because we do 3528// not want to intefere with determinstic PRNG state later, 3529// seedrandom will not call math.random on its own again after 3530// initialization. 3531// 3532mixkey(math.random(), pool); 3533 3534// End anonymous scope, and pass initial values. 3535}( 3536 [], // pool: entropy pool starts empty 3537 numeric.seedrandom, // math: package containing random, pow, and seedrandom 3538 256, // width: each RC4 output is 0 <= x < 256 3539 6, // chunks: at least six RC4 outputs for each double 3540 52 // significance: there are 52 significant digits in a double 3541 )); 3542/* This file is a slightly modified version of quadprog.js from Alberto Santini. 3543 * It has been slightly modified by Sébastien Loisel to make sure that it handles 3544 * 0-based Arrays instead of 1-based Arrays. 3545 * License is in resources/LICENSE.quadprog */ 3546(function(exports) { 3547 3548function base0to1(A) { 3549 if(typeof A !== "object") { return A; } 3550 var ret = [], i,n=A.length; 3551 for(i=0;i<n;i++) ret[i+1] = base0to1(A[i]); 3552 return ret; 3553} 3554function base1to0(A) { 3555 if(typeof A !== "object") { return A; } 3556 var ret = [], i,n=A.length; 3557 for(i=1;i<n;i++) ret[i-1] = base1to0(A[i]); 3558 return ret; 3559} 3560 3561function dpori(a, lda, n) { 3562 var i, j, k, kp1, t; 3563 3564 for (k = 1; k <= n; k = k + 1) { 3565 a[k][k] = 1 / a[k][k]; 3566 t = -a[k][k]; 3567 //~ dscal(k - 1, t, a[1][k], 1); 3568 for (i = 1; i < k; i = i + 1) { 3569 a[i][k] = t * a[i][k]; 3570 } 3571 3572 kp1 = k + 1; 3573 if (n < kp1) { 3574 break; 3575 } 3576 for (j = kp1; j <= n; j = j + 1) { 3577 t = a[k][j]; 3578 a[k][j] = 0; 3579 //~ daxpy(k, t, a[1][k], 1, a[1][j], 1); 3580 for (i = 1; i <= k; i = i + 1) { 3581 a[i][j] = a[i][j] + (t * a[i][k]); 3582 } 3583 } 3584 } 3585 3586} 3587 3588function dposl(a, lda, n, b) { 3589 var i, k, kb, t; 3590 3591 for (k = 1; k <= n; k = k + 1) { 3592 //~ t = ddot(k - 1, a[1][k], 1, b[1], 1); 3593 t = 0; 3594 for (i = 1; i < k; i = i + 1) { 3595 t = t + (a[i][k] * b[i]); 3596 } 3597 3598 b[k] = (b[k] - t) / a[k][k]; 3599 } 3600 3601 for (kb = 1; kb <= n; kb = kb + 1) { 3602 k = n + 1 - kb; 3603 b[k] = b[k] / a[k][k]; 3604 t = -b[k]; 3605 //~ daxpy(k - 1, t, a[1][k], 1, b[1], 1); 3606 for (i = 1; i < k; i = i + 1) { 3607 b[i] = b[i] + (t * a[i][k]); 3608 } 3609 } 3610} 3611 3612function dpofa(a, lda, n, info) { 3613 var i, j, jm1, k, t, s; 3614 3615 for (j = 1; j <= n; j = j + 1) { 3616 info[1] = j; 3617 s = 0; 3618 jm1 = j - 1; 3619 if (jm1 < 1) { 3620 s = a[j][j] - s; 3621 if (s <= 0) { 3622 break; 3623 } 3624 a[j][j] = Math.sqrt(s); 3625 } else { 3626 for (k = 1; k <= jm1; k = k + 1) { 3627 //~ t = a[k][j] - ddot(k - 1, a[1][k], 1, a[1][j], 1); 3628 t = a[k][j]; 3629 for (i = 1; i < k; i = i + 1) { 3630 t = t - (a[i][j] * a[i][k]); 3631 } 3632 t = t / a[k][k]; 3633 a[k][j] = t; 3634 s = s + t * t; 3635 } 3636 s = a[j][j] - s; 3637 if (s <= 0) { 3638 break; 3639 } 3640 a[j][j] = Math.sqrt(s); 3641 } 3642 info[1] = 0; 3643 } 3644} 3645 3646function qpgen2(dmat, dvec, fddmat, n, sol, crval, amat, 3647 bvec, fdamat, q, meq, iact, nact, iter, work, ierr) { 3648 3649 var i, j, l, l1, info, it1, iwzv, iwrv, iwrm, iwsv, iwuv, nvl, r, iwnbv, 3650 temp, sum, t1, tt, gc, gs, nu, 3651 t1inf, t2min, 3652 vsmall, tmpa, tmpb, 3653 go; 3654 3655 r = Math.min(n, q); 3656 l = 2 * n + (r * (r + 5)) / 2 + 2 * q + 1; 3657 3658 vsmall = 1.0e-60; 3659 do { 3660 vsmall = vsmall + vsmall; 3661 tmpa = 1 + 0.1 * vsmall; 3662 tmpb = 1 + 0.2 * vsmall; 3663 } while (tmpa <= 1 || tmpb <= 1); 3664 3665 for (i = 1; i <= n; i = i + 1) { 3666 work[i] = dvec[i]; 3667 } 3668 for (i = n + 1; i <= l; i = i + 1) { 3669 work[i] = 0; 3670 } 3671 for (i = 1; i <= q; i = i + 1) { 3672 iact[i] = 0; 3673 } 3674 3675 info = []; 3676 3677 if (ierr[1] === 0) { 3678 dpofa(dmat, fddmat, n, info); 3679 if (info[1] !== 0) { 3680 ierr[1] = 2; 3681 return; 3682 } 3683 dposl(dmat, fddmat, n, dvec); 3684 dpori(dmat, fddmat, n); 3685 } else { 3686 for (j = 1; j <= n; j = j + 1) { 3687 sol[j] = 0; 3688 for (i = 1; i <= j; i = i + 1) { 3689 sol[j] = sol[j] + dmat[i][j] * dvec[i]; 3690 } 3691 } 3692 for (j = 1; j <= n; j = j + 1) { 3693 dvec[j] = 0; 3694 for (i = j; i <= n; i = i + 1) { 3695 dvec[j] = dvec[j] + dmat[j][i] * sol[i]; 3696 } 3697 } 3698 } 3699 3700 crval[1] = 0; 3701 for (j = 1; j <= n; j = j + 1) { 3702 sol[j] = dvec[j]; 3703 crval[1] = crval[1] + work[j] * sol[j]; 3704 work[j] = 0; 3705 for (i = j + 1; i <= n; i = i + 1) { 3706 dmat[i][j] = 0; 3707 } 3708 } 3709 crval[1] = -crval[1] / 2; 3710 ierr[1] = 0; 3711 3712 iwzv = n; 3713 iwrv = iwzv + n; 3714 iwuv = iwrv + r; 3715 iwrm = iwuv + r + 1; 3716 iwsv = iwrm + (r * (r + 1)) / 2; 3717 iwnbv = iwsv + q; 3718 3719 for (i = 1; i <= q; i = i + 1) { 3720 sum = 0; 3721 for (j = 1; j <= n; j = j + 1) { 3722 sum = sum + amat[j][i] * amat[j][i]; 3723 } 3724 work[iwnbv + i] = Math.sqrt(sum); 3725 } 3726 nact = 0; 3727 iter[1] = 0; 3728 iter[2] = 0; 3729 3730 function fn_goto_50() { 3731 iter[1] = iter[1] + 1; 3732 3733 l = iwsv; 3734 for (i = 1; i <= q; i = i + 1) { 3735 l = l + 1; 3736 sum = -bvec[i]; 3737 for (j = 1; j <= n; j = j + 1) { 3738 sum = sum + amat[j][i] * sol[j]; 3739 } 3740 if (Math.abs(sum) < vsmall) { 3741 sum = 0; 3742 } 3743 if (i > meq) { 3744 work[l] = sum; 3745 } else { 3746 work[l] = -Math.abs(sum); 3747 if (sum > 0) { 3748 for (j = 1; j <= n; j = j + 1) { 3749 amat[j][i] = -amat[j][i]; 3750 } 3751 bvec[i] = -bvec[i]; 3752 } 3753 } 3754 } 3755 3756 for (i = 1; i <= nact; i = i + 1) { 3757 work[iwsv + iact[i]] = 0; 3758 } 3759 3760 nvl = 0; 3761 temp = 0; 3762 for (i = 1; i <= q; i = i + 1) { 3763 if (work[iwsv + i] < temp * work[iwnbv + i]) { 3764 nvl = i; 3765 temp = work[iwsv + i] / work[iwnbv + i]; 3766 } 3767 } 3768 if (nvl === 0) { 3769 return 999; 3770 } 3771 3772 return 0; 3773 } 3774 3775 function fn_goto_55() { 3776 for (i = 1; i <= n; i = i + 1) { 3777 sum = 0; 3778 for (j = 1; j <= n; j = j + 1) { 3779 sum = sum + dmat[j][i] * amat[j][nvl]; 3780 } 3781 work[i] = sum; 3782 } 3783 3784 l1 = iwzv; 3785 for (i = 1; i <= n; i = i + 1) { 3786 work[l1 + i] = 0; 3787 } 3788 for (j = nact + 1; j <= n; j = j + 1) { 3789 for (i = 1; i <= n; i = i + 1) { 3790 work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j]; 3791 } 3792 } 3793 3794 t1inf = true; 3795 for (i = nact; i >= 1; i = i - 1) { 3796 sum = work[i]; 3797 l = iwrm + (i * (i + 3)) / 2; 3798 l1 = l - i; 3799 for (j = i + 1; j <= nact; j = j + 1) { 3800 sum = sum - work[l] * work[iwrv + j]; 3801 l = l + j; 3802 } 3803 sum = sum / work[l1]; 3804 work[iwrv + i] = sum; 3805 if (iact[i] < meq) { 3806 // continue; 3807 break; 3808 } 3809 if (sum < 0) { 3810 // continue; 3811 break; 3812 } 3813 t1inf = false; 3814 it1 = i; 3815 } 3816 3817 if (!t1inf) { 3818 t1 = work[iwuv + it1] / work[iwrv + it1]; 3819 for (i = 1; i <= nact; i = i + 1) { 3820 if (iact[i] < meq) { 3821 // continue; 3822 break; 3823 } 3824 if (work[iwrv + i] < 0) { 3825 // continue; 3826 break; 3827 } 3828 temp = work[iwuv + i] / work[iwrv + i]; 3829 if (temp < t1) { 3830 t1 = temp; 3831 it1 = i; 3832 } 3833 } 3834 } 3835 3836 sum = 0; 3837 for (i = iwzv + 1; i <= iwzv + n; i = i + 1) { 3838 sum = sum + work[i] * work[i]; 3839 } 3840 if (Math.abs(sum) <= vsmall) { 3841 if (t1inf) { 3842 ierr[1] = 1; 3843 // GOTO 999 3844 return 999; 3845 } else { 3846 for (i = 1; i <= nact; i = i + 1) { 3847 work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i]; 3848 } 3849 work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1; 3850 // GOTO 700 3851 return 700; 3852 } 3853 } else { 3854 sum = 0; 3855 for (i = 1; i <= n; i = i + 1) { 3856 sum = sum + work[iwzv + i] * amat[i][nvl]; 3857 } 3858 tt = -work[iwsv + nvl] / sum; 3859 t2min = true; 3860 if (!t1inf) { 3861 if (t1 < tt) { 3862 tt = t1; 3863 t2min = false; 3864 } 3865 } 3866 3867 for (i = 1; i <= n; i = i + 1) { 3868 sol[i] = sol[i] + tt * work[iwzv + i]; 3869 if (Math.abs(sol[i]) < vsmall) { 3870 sol[i] = 0; 3871 } 3872 } 3873 3874 crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]); 3875 for (i = 1; i <= nact; i = i + 1) { 3876 work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i]; 3877 } 3878 work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt; 3879 3880 if (t2min) { 3881 nact = nact + 1; 3882 iact[nact] = nvl; 3883 3884 l = iwrm + ((nact - 1) * nact) / 2 + 1; 3885 for (i = 1; i <= nact - 1; i = i + 1) { 3886 work[l] = work[i]; 3887 l = l + 1; 3888 } 3889 3890 if (nact === n) { 3891 work[l] = work[n]; 3892 } else { 3893 for (i = n; i >= nact + 1; i = i - 1) { 3894 if (work[i] === 0) { 3895 // continue; 3896 break; 3897 } 3898 gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i])); 3899 gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i])); 3900 if (work[i - 1] >= 0) { 3901 temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); 3902 } else { 3903 temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); 3904 } 3905 gc = work[i - 1] / temp; 3906 gs = work[i] / temp; 3907 3908 if (gc === 1) { 3909 // continue; 3910 break; 3911 } 3912 if (gc === 0) { 3913 work[i - 1] = gs * temp; 3914 for (j = 1; j <= n; j = j + 1) { 3915 temp = dmat[j][i - 1]; 3916 dmat[j][i - 1] = dmat[j][i]; 3917 dmat[j][i] = temp; 3918 } 3919 } else { 3920 work[i - 1] = temp; 3921 nu = gs / (1 + gc); 3922 for (j = 1; j <= n; j = j + 1) { 3923 temp = gc * dmat[j][i - 1] + gs * dmat[j][i]; 3924 dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i]; 3925 dmat[j][i - 1] = temp; 3926 3927 } 3928 } 3929 } 3930 work[l] = work[nact]; 3931 } 3932 } else { 3933 sum = -bvec[nvl]; 3934 for (j = 1; j <= n; j = j + 1) { 3935 sum = sum + sol[j] * amat[j][nvl]; 3936 } 3937 if (nvl > meq) { 3938 work[iwsv + nvl] = sum; 3939 } else { 3940 work[iwsv + nvl] = -Math.abs(sum); 3941 if (sum > 0) { 3942 for (j = 1; j <= n; j = j + 1) { 3943 amat[j][nvl] = -amat[j][nvl]; 3944 } 3945 bvec[nvl] = -bvec[nvl]; 3946 } 3947 } 3948 // GOTO 700 3949 return 700; 3950 } 3951 } 3952 3953 return 0; 3954 } 3955 3956 function fn_goto_797() { 3957 l = iwrm + (it1 * (it1 + 1)) / 2 + 1; 3958 l1 = l + it1; 3959 if (work[l1] === 0) { 3960 // GOTO 798 3961 return 798; 3962 } 3963 gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1])); 3964 gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1])); 3965 if (work[l1 - 1] >= 0) { 3966 temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); 3967 } else { 3968 temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); 3969 } 3970 gc = work[l1 - 1] / temp; 3971 gs = work[l1] / temp; 3972 3973 if (gc === 1) { 3974 // GOTO 798 3975 return 798; 3976 } 3977 if (gc === 0) { 3978 for (i = it1 + 1; i <= nact; i = i + 1) { 3979 temp = work[l1 - 1]; 3980 work[l1 - 1] = work[l1]; 3981 work[l1] = temp; 3982 l1 = l1 + i; 3983 } 3984 for (i = 1; i <= n; i = i + 1) { 3985 temp = dmat[i][it1]; 3986 dmat[i][it1] = dmat[i][it1 + 1]; 3987 dmat[i][it1 + 1] = temp; 3988 } 3989 } else { 3990 nu = gs / (1 + gc); 3991 for (i = it1 + 1; i <= nact; i = i + 1) { 3992 temp = gc * work[l1 - 1] + gs * work[l1]; 3993 work[l1] = nu * (work[l1 - 1] + temp) - work[l1]; 3994 work[l1 - 1] = temp; 3995 l1 = l1 + i; 3996 } 3997 for (i = 1; i <= n; i = i + 1) { 3998 temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1]; 3999 dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1]; 4000 dmat[i][it1] = temp; 4001 } 4002 } 4003 4004 return 0; 4005 } 4006 4007 function fn_goto_798() { 4008 l1 = l - it1; 4009 for (i = 1; i <= it1; i = i + 1) { 4010 work[l1] = work[l]; 4011 l = l + 1; 4012 l1 = l1 + 1; 4013 } 4014 4015 work[iwuv + it1] = work[iwuv + it1 + 1]; 4016 iact[it1] = iact[it1 + 1]; 4017 it1 = it1 + 1; 4018 if (it1 < nact) { 4019 // GOTO 797 4020 return 797; 4021 } 4022 4023 return 0; 4024 } 4025 4026 function fn_goto_799() { 4027 work[iwuv + nact] = work[iwuv + nact + 1]; 4028 work[iwuv + nact + 1] = 0; 4029 iact[nact] = 0; 4030 nact = nact - 1; 4031 iter[2] = iter[2] + 1; 4032 4033 return 0; 4034 } 4035 4036 go = 0; 4037 while (true) { 4038 go = fn_goto_50(); 4039 if (go === 999) { 4040 return; 4041 } 4042 while (true) { 4043 go = fn_goto_55(); 4044 if (go === 0) { 4045 break; 4046 } 4047 if (go === 999) { 4048 return; 4049 } 4050 if (go === 700) { 4051 if (it1 === nact) { 4052 fn_goto_799(); 4053 } else { 4054 while (true) { 4055 fn_goto_797(); 4056 go = fn_goto_798(); 4057 if (go !== 797) { 4058 break; 4059 } 4060 } 4061 fn_goto_799(); 4062 } 4063 } 4064 } 4065 } 4066 4067} 4068 4069function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) { 4070 Dmat = base0to1(Dmat); 4071 dvec = base0to1(dvec); 4072 Amat = base0to1(Amat); 4073 var i, n, q, 4074 nact, r, 4075 crval = [], iact = [], sol = [], work = [], iter = [], 4076 message; 4077 4078 meq = meq || 0; 4079 factorized = factorized ? base0to1(factorized) : [undefined, 0]; 4080 bvec = bvec ? base0to1(bvec) : []; 4081 4082 // In Fortran the array index starts from 1 4083 n = Dmat.length - 1; 4084 q = Amat[1].length - 1; 4085 4086 if (!bvec) { 4087 for (i = 1; i <= q; i = i + 1) { 4088 bvec[i] = 0; 4089 } 4090 } 4091 for (i = 1; i <= q; i = i + 1) { 4092 iact[i] = 0; 4093 } 4094 nact = 0; 4095 r = Math.min(n, q); 4096 for (i = 1; i <= n; i = i + 1) { 4097 sol[i] = 0; 4098 } 4099 crval[1] = 0; 4100 for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) { 4101 work[i] = 0; 4102 } 4103 for (i = 1; i <= 2; i = i + 1) { 4104 iter[i] = 0; 4105 } 4106 4107 qpgen2(Dmat, dvec, n, n, sol, crval, Amat, 4108 bvec, n, q, meq, iact, nact, iter, work, factorized); 4109 4110 message = ""; 4111 if (factorized[1] === 1) { 4112 message = "constraints are inconsistent, no solution!"; 4113 } 4114 if (factorized[1] === 2) { 4115 message = "matrix D in quadratic function is not positive definite!"; 4116 } 4117 4118 return { 4119 solution: base1to0(sol), 4120 value: base1to0(crval), 4121 unconstrained_solution: base1to0(dvec), 4122 iterations: base1to0(iter), 4123 iact: base1to0(iact), 4124 message: message 4125 }; 4126} 4127exports.solveQP = solveQP; 4128}(numeric)); 4129/*
4129 4130Shanti Rao sent me this routine by private email. I had to modify it 4131slightly to work on Arrays instead of using a Matrix object. 4132It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py 4133*/ 4134 4135numeric.svd= function svd(A) { 4136 var temp; 4137//Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970) 4138 var prec= numeric.epsilon; //Math.pow(2,-52) // assumes double prec 4139 var tolerance= 1.e-64/prec; 4140 var itmax= 50; 4141 var c=0; 4142 var i=0; 4143 var j=0; 4144 var k=0; 4145 var l=0; 4146 4147 var u= numeric.clone(A); 4148 var m= u.length; 4149 4150 var n= u[0].length; 4151 4152 if (m < n) throw "Need more rows than columns" 4153 4154 var e = new Array(n); 4155 var q = new Array(n); 4156 for (i=0; i<n; i++) e[i] = q[i] = 0.0; 4157 var v = numeric.rep([n,n],0); 4158// v.zero(); 4159 4160 function pythag(a,b) 4161 { 4162 a = Math.abs(a) 4163 b = Math.abs(b) 4164 if (a > b) 4165 return a*Math.sqrt(1.0+(b*b/a/a)) 4166 else if (b == 0.0) 4167 return a 4168 return b*Math.sqrt(1.0+(a*a/b/b)) 4169 } 4170 4171 //Householder's reduction to bidiagonal form 4172 4173 var f= 0.0; 4174 var g= 0.0; 4175 var h= 0.0; 4176 var x= 0.0; 4177 var y= 0.0; 4178 var z= 0.0; 4179 var s= 0.0; 4180 4181 for (i=0; i < n; i++) 4182 { 4183 e[i]= g; 4184 s= 0.0; 4185 l= i+1; 4186 for (j=i; j < m; j++) 4187 s += (u[j][i]*u[j][i]); 4188 if (s <= tolerance) 4189 g= 0.0; 4190 else 4191 { 4192 f= u[i][i]; 4193 g= Math.sqrt(s); 4194 if (f >= 0.0) g= -g; 4195 h= f*g-s 4196 u[i][i]=f-g; 4197 for (j=l; j < n; j++) 4198 { 4199 s= 0.0 4200 for (k=i; k < m; k++) 4201 s += u[k][i]*u[k][j] 4202 f= s/h 4203 for (k=i; k < m; k++) 4204 u[k][j]+=f*u[k][i] 4205 } 4206 } 4207 q[i]= g 4208 s= 0.0 4209 for (j=l; j < n; j++) 4210 s= s + u[i][j]*u[i][j] 4211 if (s <= tolerance) 4212 g= 0.0 4213 else 4214 { 4215 f= u[i][i+1] 4216 g= Math.sqrt(s) 4217 if (f >= 0.0) g= -g 4218 h= f*g - s 4219 u[i][i+1] = f-g; 4220 for (j=l; j < n; j++) e[j]= u[i][j]/h 4221 for (j=l; j < m; j++) 4222 { 4223 s=0.0 4224 for (k=l; k < n; k++) 4225 s += (u[j][k]*u[i][k]) 4226 for (k=l; k < n; k++) 4227 u[j][k]+=s*e[k] 4228 } 4229 } 4230 y= Math.abs(q[i])+Math.abs(e[i]) 4231 if (y>x) 4232 x=y 4233 } 4234 4235 // accumulation of right hand gtransformations 4236 for (i=n-1; i != -1; i+= -1) 4237 { 4238 if (g != 0.0) 4239 { 4240 h= g*u[i][i+1] 4241 for (j=l; j < n; j++) 4242 v[j][i]=u[i][j]/h 4243 for (j=l; j < n; j++) 4244 { 4245 s=0.0 4246 for (k=l; k < n; k++) 4247 s += u[i][k]*v[k][j] 4248 for (k=l; k < n; k++) 4249 v[k][j]+=(s*v[k][i]) 4250 } 4251 } 4252 for (j=l; j < n; j++) 4253 { 4254 v[i][j] = 0; 4255 v[j][i] = 0; 4256 } 4257 v[i][i] = 1; 4258 g= e[i] 4259 l= i 4260 } 4261 4262 // accumulation of left hand transformations 4263 for (i=n-1; i != -1; i+= -1) 4264 { 4265 l= i+1 4266 g= q[i] 4267 for (j=l; j < n; j++) 4268 u[i][j] = 0; 4269 if (g != 0.0) 4270 { 4271 h= u[i][i]*g 4272 for (j=l; j < n; j++) 4273 { 4274 s=0.0 4275 for (k=l; k < m; k++) s += u[k][i]*u[k][j]; 4276 f= s/h 4277 for (k=i; k < m; k++) u[k][j]+=f*u[k][i]; 4278 } 4279 for (j=i; j < m; j++) u[j][i] = u[j][i]/g; 4280 } 4281 else 4282 for (j=i; j < m; j++) u[j][i] = 0; 4283 u[i][i] += 1; 4284 } 4285 4286 // diagonalization of the bidiagonal form 4287 prec= prec*x 4288 for (k=n-1; k != -1; k+= -1) 4289 { 4290 for (var iteration=0; iteration < itmax; iteration++) 4291 { // test f splitting 4292 var test_convergence = false 4293 for (l=k; l != -1; l+= -1) 4294 { 4295 if (Math.abs(e[l]) <= prec) 4296 { test_convergence= true 4297 break 4298 } 4299 if (Math.abs(q[l-1]) <= prec) 4300 break 4301 }
4302 if (!test_convergence) 4303 { // cancellation of e[l] if l>0 4304 c= 0.0 4305 s= 1.0 4306 var l1= l-1 4307 for (i =l; i<k+1; i++) 4308 { 4309 f= s*e[i] 4310 e[i]= c*e[i] 4311 if (Math.abs(f) <= prec) 4312 break 4313 g= q[i] 4314 h= pythag(f,g) 4315 q[i]= h 4316 c= g/h 4317 s= -f/h 4318 for (j=0; j < m; j++) 4319 { 4320 y= u[j][l1] 4321 z= u[j][i] 4322 u[j][l1] = y*c+(z*s) 4323 u[j][i] = -y*s+(z*c) 4324 } 4325 } 4326 } 4327 // test f convergence 4328 z= q[k] 4329 if (l== k) 4330 { //convergence 4331 if (z<0.0) 4332 { //q[k] is made non-negative 4333 q[k]= -z 4334 for (j=0; j < n; j++) 4335 v[j][k] = -v[j][k] 4336 } 4337 break //break out of iteration loop and move on to next k value 4338 } 4339 if (iteration >= itmax-1) 4340 throw 'Error: no convergence.' 4341 // shift from bottom 2x2 minor 4342 x= q[l] 4343 y= q[k-1] 4344 g= e[k-1] 4345 h= e[k] 4346 f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) 4347 g= pythag(f,1.0) 4348 if (f < 0.0) 4349 f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x 4350 else 4351 f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x 4352 // next QR transformation 4353 c= 1.0 4354 s= 1.0 4355 for (i=l+1; i< k+1; i++) 4356 { 4357 g= e[i] 4358 y= q[i] 4359 h= s*g 4360 g= c*g 4361 z= pythag(f,h) 4362 e[i-1]= z 4363 c= f/z 4364 s= h/z 4365 f= x*c+g*s 4366 g= -x*s+g*c 4367 h= y*s 4368 y= y*c 4369 for (j=0; j < n; j++) 4370 { 4371 x= v[j][i-1] 4372 z= v[j][i] 4373 v[j][i-1] = x*c+z*s 4374 v[j][i] = -x*s+z*c 4375 } 4376 z= pythag(f,h) 4377 q[i-1]= z 4378 c= f/z 4379 s= h/z 4380 f= c*g+s*y 4381 x= -s*g+c*y 4382 for (j=0; j < m; j++) 4383 { 4384 y= u[j][i-1] 4385 z= u[j][i] 4386 u[j][i-1] = y*c+z*s 4387 u[j][i] = -y*s+z*c 4388 } 4389 } 4390 e[l]= 0.0 4391 e[k]= f 4392 q[k]= x 4393 } 4394 } 4395 4396 //vt= transpose(v) 4397 //return (u,q,vt) 4398 for (i=0;i<q.length; i++) 4399 if (q[i] < prec) q[i] = 0 4400 4401 //sort eigenvalues 4402 for (i=0; i< n; i++) 4403 { 4404 //writeln(q) 4405 for (j=i-1; j >= 0; j--) 4406 { 4407 if (q[j] < q[i]) 4408 { 4409 // writeln(i,'-',j) 4410 c = q[j] 4411 q[j] = q[i] 4412 q[i] = c 4413 for(k=0;k<u.length;k++) { temp = u[k][i]; u[k][i] = u[k][j]; u[k][j] = temp; } 4414 for(k=0;k<v.length;k++) { temp = v[k][i]; v[k][i] = v[k][j]; v[k][j] = temp; } 4415// u.swapCols(i,j) 4416// v.swapCols(i,j) 4417 i = j 4418 } 4419 } 4420 } 4421 4422 return {U:u,S:q,V:v} 4423}; 4424
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.