From 1957077f2a0735c27754814f5444fb360365213b Mon Sep 17 00:00:00 2001 From: Franklin Ta Date: Tue, 27 May 2014 05:37:07 -0700 Subject: [PATCH] Add numeric 1.2.6 Source: https://www.npmjs.org/package/numeric --- ajax/libs/numeric/1.2.6/numeric.js | 4424 ++++++++++++++++++++++++ ajax/libs/numeric/1.2.6/numeric.min.js | 3 + ajax/libs/numeric/package.json | 22 + 3 files changed, 4449 insertions(+) create mode 100644 ajax/libs/numeric/1.2.6/numeric.js create mode 100644 ajax/libs/numeric/1.2.6/numeric.min.js create mode 100644 ajax/libs/numeric/package.json diff --git a/ajax/libs/numeric/1.2.6/numeric.js b/ajax/libs/numeric/1.2.6/numeric.js new file mode 100644 index 000000000..44edb9b8d --- /dev/null +++ b/ajax/libs/numeric/1.2.6/numeric.js @@ -0,0 +1,4424 @@ +"use strict"; + +var numeric = (typeof exports === "undefined")?(function numeric() {}):(exports); +if(typeof global !== "undefined") { global.numeric = numeric; } + +numeric.version = "1.2.6"; + +// 1. Utility functions +numeric.bench = function bench (f,interval) { + var t1,t2,n,i; + if(typeof interval === "undefined") { interval = 15; } + n = 0.5; + t1 = new Date(); + while(1) { + n*=2; + for(i=n;i>3;i-=4) { f(); f(); f(); f(); } + while(i>0) { f(); i--; } + t2 = new Date(); + if(t2-t1 > interval) break; + } + for(i=n;i>3;i-=4) { f(); f(); f(); f(); } + while(i>0) { f(); i--; } + t2 = new Date(); + return 1000*(3*n-1)/(t2-t1); +} + +numeric._myIndexOf = (function _myIndexOf(w) { + var n = this.length,k; + for(k=0;k numeric.largeArray) { ret.push('...Large Array...'); return true; } + var flag = false; + ret.push('['); + for(k=0;k0) { ret.push(','); if(flag) ret.push('\n '); } flag = foo(x[k]); } + ret.push(']'); + return true; + } + ret.push('{'); + var flag = false; + for(k in x) { if(x.hasOwnProperty(k)) { if(flag) ret.push(',\n'); flag = true; ret.push(k); ret.push(': \n'); foo(x[k]); } } + ret.push('}'); + return true; + } + foo(x); + return ret.join(''); +} + +numeric.parseDate = function parseDate(d) { + function foo(d) { + if(typeof d === 'string') { return Date.parse(d.replace(/-/g,'/')); } + if(!(d instanceof Array)) { throw new Error("parseDate: parameter must be arrays of strings"); } + var ret = [],k; + for(k=0;k0) { + ret[count] = []; + for(j=0;j> 2; + q = ((x & 3) << 4) + (y >> 4); + r = ((y & 15) << 2) + (z >> 6); + s = z & 63; + if(i+1>=n) { r = s = 64; } + else if(i+2>=n) { s = 64; } + ret += key.charAt(p) + key.charAt(q) + key.charAt(r) + key.charAt(s); + } + return ret; + } + function crc32Array (a,from,to) { + if(typeof from === "undefined") { from = 0; } + if(typeof to === "undefined") { to = a.length; } + var table = [0x00000000, 0x77073096, 0xEE0E612C, 0x990951BA, 0x076DC419, 0x706AF48F, 0xE963A535, 0x9E6495A3, + 0x0EDB8832, 0x79DCB8A4, 0xE0D5E91E, 0x97D2D988, 0x09B64C2B, 0x7EB17CBD, 0xE7B82D07, 0x90BF1D91, + 0x1DB71064, 0x6AB020F2, 0xF3B97148, 0x84BE41DE, 0x1ADAD47D, 0x6DDDE4EB, 0xF4D4B551, 0x83D385C7, + 0x136C9856, 0x646BA8C0, 0xFD62F97A, 0x8A65C9EC, 0x14015C4F, 0x63066CD9, 0xFA0F3D63, 0x8D080DF5, + 0x3B6E20C8, 0x4C69105E, 0xD56041E4, 0xA2677172, 0x3C03E4D1, 0x4B04D447, 0xD20D85FD, 0xA50AB56B, + 0x35B5A8FA, 0x42B2986C, 0xDBBBC9D6, 0xACBCF940, 0x32D86CE3, 0x45DF5C75, 0xDCD60DCF, 0xABD13D59, + 0x26D930AC, 0x51DE003A, 0xC8D75180, 0xBFD06116, 0x21B4F4B5, 0x56B3C423, 0xCFBA9599, 0xB8BDA50F, + 0x2802B89E, 0x5F058808, 0xC60CD9B2, 0xB10BE924, 0x2F6F7C87, 0x58684C11, 0xC1611DAB, 0xB6662D3D, + 0x76DC4190, 0x01DB7106, 0x98D220BC, 0xEFD5102A, 0x71B18589, 0x06B6B51F, 0x9FBFE4A5, 0xE8B8D433, + 0x7807C9A2, 0x0F00F934, 0x9609A88E, 0xE10E9818, 0x7F6A0DBB, 0x086D3D2D, 0x91646C97, 0xE6635C01, + 0x6B6B51F4, 0x1C6C6162, 0x856530D8, 0xF262004E, 0x6C0695ED, 0x1B01A57B, 0x8208F4C1, 0xF50FC457, + 0x65B0D9C6, 0x12B7E950, 0x8BBEB8EA, 0xFCB9887C, 0x62DD1DDF, 0x15DA2D49, 0x8CD37CF3, 0xFBD44C65, + 0x4DB26158, 0x3AB551CE, 0xA3BC0074, 0xD4BB30E2, 0x4ADFA541, 0x3DD895D7, 0xA4D1C46D, 0xD3D6F4FB, + 0x4369E96A, 0x346ED9FC, 0xAD678846, 0xDA60B8D0, 0x44042D73, 0x33031DE5, 0xAA0A4C5F, 0xDD0D7CC9, + 0x5005713C, 0x270241AA, 0xBE0B1010, 0xC90C2086, 0x5768B525, 0x206F85B3, 0xB966D409, 0xCE61E49F, + 0x5EDEF90E, 0x29D9C998, 0xB0D09822, 0xC7D7A8B4, 0x59B33D17, 0x2EB40D81, 0xB7BD5C3B, 0xC0BA6CAD, + 0xEDB88320, 0x9ABFB3B6, 0x03B6E20C, 0x74B1D29A, 0xEAD54739, 0x9DD277AF, 0x04DB2615, 0x73DC1683, + 0xE3630B12, 0x94643B84, 0x0D6D6A3E, 0x7A6A5AA8, 0xE40ECF0B, 0x9309FF9D, 0x0A00AE27, 0x7D079EB1, + 0xF00F9344, 0x8708A3D2, 0x1E01F268, 0x6906C2FE, 0xF762575D, 0x806567CB, 0x196C3671, 0x6E6B06E7, + 0xFED41B76, 0x89D32BE0, 0x10DA7A5A, 0x67DD4ACC, 0xF9B9DF6F, 0x8EBEEFF9, 0x17B7BE43, 0x60B08ED5, + 0xD6D6A3E8, 0xA1D1937E, 0x38D8C2C4, 0x4FDFF252, 0xD1BB67F1, 0xA6BC5767, 0x3FB506DD, 0x48B2364B, + 0xD80D2BDA, 0xAF0A1B4C, 0x36034AF6, 0x41047A60, 0xDF60EFC3, 0xA867DF55, 0x316E8EEF, 0x4669BE79, + 0xCB61B38C, 0xBC66831A, 0x256FD2A0, 0x5268E236, 0xCC0C7795, 0xBB0B4703, 0x220216B9, 0x5505262F, + 0xC5BA3BBE, 0xB2BD0B28, 0x2BB45A92, 0x5CB36A04, 0xC2D7FFA7, 0xB5D0CF31, 0x2CD99E8B, 0x5BDEAE1D, + 0x9B64C2B0, 0xEC63F226, 0x756AA39C, 0x026D930A, 0x9C0906A9, 0xEB0E363F, 0x72076785, 0x05005713, + 0x95BF4A82, 0xE2B87A14, 0x7BB12BAE, 0x0CB61B38, 0x92D28E9B, 0xE5D5BE0D, 0x7CDCEFB7, 0x0BDBDF21, + 0x86D3D2D4, 0xF1D4E242, 0x68DDB3F8, 0x1FDA836E, 0x81BE16CD, 0xF6B9265B, 0x6FB077E1, 0x18B74777, + 0x88085AE6, 0xFF0F6A70, 0x66063BCA, 0x11010B5C, 0x8F659EFF, 0xF862AE69, 0x616BFFD3, 0x166CCF45, + 0xA00AE278, 0xD70DD2EE, 0x4E048354, 0x3903B3C2, 0xA7672661, 0xD06016F7, 0x4969474D, 0x3E6E77DB, + 0xAED16A4A, 0xD9D65ADC, 0x40DF0B66, 0x37D83BF0, 0xA9BCAE53, 0xDEBB9EC5, 0x47B2CF7F, 0x30B5FFE9, + 0xBDBDF21C, 0xCABAC28A, 0x53B39330, 0x24B4A3A6, 0xBAD03605, 0xCDD70693, 0x54DE5729, 0x23D967BF, + 0xB3667A2E, 0xC4614AB8, 0x5D681B02, 0x2A6F2B94, 0xB40BBE37, 0xC30C8EA1, 0x5A05DF1B, 0x2D02EF8D]; + + var crc = -1, y = 0, n = a.length,i; + + for (i = from; i < to; i++) { + y = (crc ^ a[i]) & 0xFF; + crc = (crc >>> 8) ^ table[y]; + } + + return crc ^ (-1); + } + + var h = img[0].length, w = img[0][0].length, s1, s2, next,k,length,a,b,i,j,adler32,crc32; + var stream = [ + 137, 80, 78, 71, 13, 10, 26, 10, // 0: PNG signature + 0,0,0,13, // 8: IHDR Chunk length + 73, 72, 68, 82, // 12: "IHDR" + (w >> 24) & 255, (w >> 16) & 255, (w >> 8) & 255, w&255, // 16: Width + (h >> 24) & 255, (h >> 16) & 255, (h >> 8) & 255, h&255, // 20: Height + 8, // 24: bit depth + 2, // 25: RGB + 0, // 26: deflate + 0, // 27: no filter + 0, // 28: no interlace + -1,-2,-3,-4, // 29: CRC + -5,-6,-7,-8, // 33: IDAT Chunk length + 73, 68, 65, 84, // 37: "IDAT" + // RFC 1950 header starts here + 8, // 41: RFC1950 CMF + 29 // 42: RFC1950 FLG + ]; + crc32 = crc32Array(stream,12,29); + stream[29] = (crc32>>24)&255; + stream[30] = (crc32>>16)&255; + stream[31] = (crc32>>8)&255; + stream[32] = (crc32)&255; + s1 = 1; + s2 = 0; + for(i=0;i>8)&255; + stream.push(a); stream.push(b); + stream.push((~a)&255); stream.push((~b)&255); + if(i===0) stream.push(0); + for(j=0;j255) a = 255; + else if(a<0) a=0; + else a = Math.round(a); + s1 = (s1 + a )%65521; + s2 = (s2 + s1)%65521; + stream.push(a); + } + } + stream.push(0); + } + adler32 = (s2<<16)+s1; + stream.push((adler32>>24)&255); + stream.push((adler32>>16)&255); + stream.push((adler32>>8)&255); + stream.push((adler32)&255); + length = stream.length - 41; + stream[33] = (length>>24)&255; + stream[34] = (length>>16)&255; + stream[35] = (length>>8)&255; + stream[36] = (length)&255; + crc32 = crc32Array(stream,37); + stream.push((crc32>>24)&255); + stream.push((crc32>>16)&255); + stream.push((crc32>>8)&255); + stream.push((crc32)&255); + stream.push(0); + stream.push(0); + stream.push(0); + stream.push(0); +// a = stream.length; + stream.push(73); // I + stream.push(69); // E + stream.push(78); // N + stream.push(68); // D + stream.push(174); // CRC1 + stream.push(66); // CRC2 + stream.push(96); // CRC3 + stream.push(130); // CRC4 + return 'data:image/png;base64,'+base64(stream); +} + +// 2. Linear algebra with Arrays. +numeric._dim = function _dim(x) { + var ret = []; + while(typeof x === "object") { ret.push(x.length); x = x[0]; } + return ret; +} + +numeric.dim = function dim(x) { + var y,z; + if(typeof x === "object") { + y = x[0]; + if(typeof y === "object") { + z = y[0]; + if(typeof z === "object") { + return numeric._dim(x); + } + return [x.length,y.length]; + } + return [x.length]; + } + return []; +} + +numeric.mapreduce = function mapreduce(body,init) { + return Function('x','accum','_s','_k', + 'if(typeof accum === "undefined") accum = '+init+';\n'+ + 'if(typeof x === "number") { var xi = x; '+body+'; return accum; }\n'+ + 'if(typeof _s === "undefined") _s = numeric.dim(x);\n'+ + 'if(typeof _k === "undefined") _k = 0;\n'+ + 'var _n = _s[_k];\n'+ + 'var i,xi;\n'+ + 'if(_k < _s.length-1) {\n'+ + ' for(i=_n-1;i>=0;i--) {\n'+ + ' accum = arguments.callee(x[i],accum,_s,_k+1);\n'+ + ' }'+ + ' return accum;\n'+ + '}\n'+ + 'for(i=_n-1;i>=1;i-=2) { \n'+ + ' xi = x[i];\n'+ + ' '+body+';\n'+ + ' xi = x[i-1];\n'+ + ' '+body+';\n'+ + '}\n'+ + 'if(i === 0) {\n'+ + ' xi = x[i];\n'+ + ' '+body+'\n'+ + '}\n'+ + 'return accum;' + ); +} +numeric.mapreduce2 = function mapreduce2(body,setup) { + return Function('x', + 'var n = x.length;\n'+ + 'var i,xi;\n'+setup+';\n'+ + 'for(i=n-1;i!==-1;--i) { \n'+ + ' xi = x[i];\n'+ + ' '+body+';\n'+ + '}\n'+ + 'return accum;' + ); +} + + +numeric.same = function same(x,y) { + var i,n; + if(!(x instanceof Array) || !(y instanceof Array)) { return false; } + n = x.length; + if(n !== y.length) { return false; } + for(i=0;i=0;i-=2) { ret[i+1] = v; ret[i] = v; } + if(i===-1) { ret[0] = v; } + return ret; + } + for(i=n-1;i>=0;i--) { ret[i] = numeric.rep(s,v,k+1); } + return ret; +} + + +numeric.dotMMsmall = function dotMMsmall(x,y) { + var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0; + p = x.length; q = y.length; r = y[0].length; + ret = Array(p); + for(i=p-1;i>=0;i--) { + foo = Array(r); + bar = x[i]; + for(k=r-1;k>=0;k--) { + woo = bar[q-1]*y[q-1][k]; + for(j=q-2;j>=1;j-=2) { + i0 = j-1; + woo += bar[j]*y[j][k] + bar[i0]*y[i0][k]; + } + if(j===0) { woo += bar[0]*y[0][k]; } + foo[k] = woo; + } + ret[i] = foo; + } + return ret; +} +numeric._getCol = function _getCol(A,j,x) { + var n = A.length, i; + for(i=n-1;i>0;--i) { + x[i] = A[i][j]; + --i; + x[i] = A[i][j]; + } + if(i===0) x[0] = A[0][j]; +} +numeric.dotMMbig = function dotMMbig(x,y){ + var gc = numeric._getCol, p = y.length, v = Array(p); + var m = x.length, n = y[0].length, A = new Array(m), xj; + var VV = numeric.dotVV; + var i,j,k,z; + --p; + --m; + for(i=m;i!==-1;--i) A[i] = Array(n); + --n; + for(i=n;i!==-1;--i) { + gc(y,i,v); + for(j=m;j!==-1;--j) { + z=0; + xj = x[j]; + A[j][i] = VV(xj,v); + } + } + return A; +} + +numeric.dotMV = function dotMV(x,y) { + var p = x.length, q = y.length,i; + var ret = Array(p), dotVV = numeric.dotVV; + for(i=p-1;i>=0;i--) { ret[i] = dotVV(x[i],y); } + return ret; +} + +numeric.dotVM = function dotVM(x,y) { + var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0,s1,s2,s3,baz,accum; + p = x.length; q = y[0].length; + ret = Array(q); + for(k=q-1;k>=0;k--) { + woo = x[p-1]*y[p-1][k]; + for(j=p-2;j>=1;j-=2) { + i0 = j-1; + woo += x[j]*y[j][k] + x[i0]*y[i0][k]; + } + if(j===0) { woo += x[0]*y[0][k]; } + ret[k] = woo; + } + return ret; +} + +numeric.dotVV = function dotVV(x,y) { + var i,n=x.length,i1,ret = x[n-1]*y[n-1]; + for(i=n-2;i>=1;i-=2) { + i1 = i-1; + ret += x[i]*y[i] + x[i1]*y[i1]; + } + if(i===0) { ret += x[0]*y[0]; } + return ret; +} + +numeric.dot = function dot(x,y) { + var d = numeric.dim; + switch(d(x).length*1000+d(y).length) { + case 2002: + if(y.length < 10) return numeric.dotMMsmall(x,y); + else return numeric.dotMMbig(x,y); + case 2001: return numeric.dotMV(x,y); + case 1002: return numeric.dotVM(x,y); + case 1001: return numeric.dotVV(x,y); + case 1000: return numeric.mulVS(x,y); + case 1: return numeric.mulSV(x,y); + case 0: return x*y; + default: throw new Error('numeric.dot only works on vectors and matrices'); + } +} + +numeric.diag = function diag(d) { + var i,i1,j,n = d.length, A = Array(n), Ai; + for(i=n-1;i>=0;i--) { + Ai = Array(n); + i1 = i+2; + for(j=n-1;j>=i1;j-=2) { + Ai[j] = 0; + Ai[j-1] = 0; + } + if(j>i) { Ai[j] = 0; } + Ai[i] = d[i]; + for(j=i-1;j>=1;j-=2) { + Ai[j] = 0; + Ai[j-1] = 0; + } + if(j===0) { Ai[0] = 0; } + A[i] = Ai; + } + return A; +} +numeric.getDiag = function(A) { + var n = Math.min(A.length,A[0].length),i,ret = Array(n); + for(i=n-1;i>=1;--i) { + ret[i] = A[i][i]; + --i; + ret[i] = A[i][i]; + } + if(i===0) { + ret[0] = A[0][0]; + } + return ret; +} + +numeric.identity = function identity(n) { return numeric.diag(numeric.rep([n],1)); } +numeric.pointwise = function pointwise(params,body,setup) { + if(typeof setup === "undefined") { setup = ""; } + var fun = []; + var k; + var avec = /\[i\]$/,p,thevec = ''; + var haveret = false; + for(k=0;k=0;i--) ret[i] = arguments.callee('+params.join(',')+',_s,_k+1);\n'+ + ' return ret;\n'+ + '}\n'+ + setup+'\n'+ + 'for(i=_n-1;i!==-1;--i) {\n'+ + ' '+body+'\n'+ + '}\n'+ + 'return ret;' + ); + return Function.apply(null,fun); +} +numeric.pointwise2 = function pointwise2(params,body,setup) { + if(typeof setup === "undefined") { setup = ""; } + var fun = []; + var k; + var avec = /\[i\]$/,p,thevec = ''; + var haveret = false; + for(k=0;k=0;i--) { _biforeach(typeof x==="object"?x[i]:x,typeof y==="object"?y[i]:y,s,k+1,f); } +}); +numeric._biforeach2 = (function _biforeach2(x,y,s,k,f) { + if(k === s.length-1) { return f(x,y); } + var i,n=s[k],ret = Array(n); + 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); } + return ret; +}); +numeric._foreach = (function _foreach(x,s,k,f) { + if(k === s.length-1) { f(x); return; } + var i,n=s[k]; + for(i=n-1;i>=0;i--) { _foreach(x[i],s,k+1,f); } +}); +numeric._foreach2 = (function _foreach2(x,s,k,f) { + if(k === s.length-1) { return f(x); } + var i,n=s[k], ret = Array(n); + for(i=n-1;i>=0;i--) { ret[i] = _foreach2(x[i],s,k+1,f); } + return ret; +}); + +/*numeric.anyV = numeric.mapreduce('if(xi) return true;','false'); +numeric.allV = numeric.mapreduce('if(!xi) return false;','true'); +numeric.any = function(x) { if(typeof x.length === "undefined") return x; return numeric.anyV(x); } +numeric.all = function(x) { if(typeof x.length === "undefined") return x; return numeric.allV(x); }*/ + +numeric.ops2 = { + add: '+', + sub: '-', + mul: '*', + div: '/', + mod: '%', + and: '&&', + or: '||', + eq: '===', + neq: '!==', + lt: '<', + gt: '>', + leq: '<=', + geq: '>=', + band: '&', + bor: '|', + bxor: '^', + lshift: '<<', + rshift: '>>', + rrshift: '>>>' +}; +numeric.opseq = { + addeq: '+=', + subeq: '-=', + muleq: '*=', + diveq: '/=', + modeq: '%=', + lshifteq: '<<=', + rshifteq: '>>=', + rrshifteq: '>>>=', + bandeq: '&=', + boreq: '|=', + bxoreq: '^=' +}; +numeric.mathfuns = ['abs','acos','asin','atan','ceil','cos', + 'exp','floor','log','round','sin','sqrt','tan', + 'isNaN','isFinite']; +numeric.mathfuns2 = ['atan2','pow','max','min']; +numeric.ops1 = { + neg: '-', + not: '!', + bnot: '~', + clone: '' +}; +numeric.mapreducers = { + any: ['if(xi) return true;','var accum = false;'], + all: ['if(!xi) return false;','var accum = true;'], + sum: ['accum += xi;','var accum = 0;'], + prod: ['accum *= xi;','var accum = 1;'], + norm2Squared: ['accum += xi*xi;','var accum = 0;'], + norminf: ['accum = max(accum,abs(xi));','var accum = 0, max = Math.max, abs = Math.abs;'], + norm1: ['accum += abs(xi)','var accum = 0, abs = Math.abs;'], + sup: ['accum = max(accum,xi);','var accum = -Infinity, max = Math.max;'], + inf: ['accum = min(accum,xi);','var accum = Infinity, min = Math.min;'] +}; + +(function () { + var i,o; + for(i=0;iv0) { i0 = i; v0 = k; } } + Aj = A[i0]; A[i0] = A[j]; A[j] = Aj; + Ij = I[i0]; I[i0] = I[j]; I[j] = Ij; + x = Aj[j]; + for(k=j;k!==n;++k) Aj[k] /= x; + for(k=n-1;k!==-1;--k) Ij[k] /= x; + for(i=m-1;i!==-1;--i) { + if(i!==j) { + Ai = A[i]; + Ii = I[i]; + x = Ai[j]; + for(k=j+1;k!==n;++k) Ai[k] -= Aj[k]*x; + for(k=n-1;k>0;--k) { Ii[k] -= Ij[k]*x; --k; Ii[k] -= Ij[k]*x; } + if(k===0) Ii[0] -= Ij[0]*x; + } + } + } + return I; +} + +numeric.det = function det(x) { + var s = numeric.dim(x); + if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: det() only works on square matrices'); } + var n = s[0], ret = 1,i,j,k,A = numeric.clone(x),Aj,Ai,alpha,temp,k1,k2,k3; + for(j=0;j Math.abs(A[k][j])) { k = i; } } + if(k !== j) { + temp = A[k]; A[k] = A[j]; A[j] = temp; + ret *= -1; + } + Aj = A[j]; + for(i=j+1;i=1;i-=2) { + A1 = x[i]; + A0 = x[i-1]; + for(j=n-1;j>=1;--j) { + Bj = ret[j]; Bj[i] = A1[j]; Bj[i-1] = A0[j]; + --j; + Bj = ret[j]; Bj[i] = A1[j]; Bj[i-1] = A0[j]; + } + if(j===0) { + Bj = ret[0]; Bj[i] = A1[0]; Bj[i-1] = A0[0]; + } + } + if(i===0) { + A0 = x[0]; + for(j=n-1;j>=1;--j) { + ret[j][0] = A0[j]; + --j; + ret[j][0] = A0[j]; + } + if(j===0) { ret[0][0] = A0[0]; } + } + return ret; +} +numeric.negtranspose = function negtranspose(x) { + var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj; + for(j=0;j=1;i-=2) { + A1 = x[i]; + A0 = x[i-1]; + for(j=n-1;j>=1;--j) { + Bj = ret[j]; Bj[i] = -A1[j]; Bj[i-1] = -A0[j]; + --j; + Bj = ret[j]; Bj[i] = -A1[j]; Bj[i-1] = -A0[j]; + } + if(j===0) { + Bj = ret[0]; Bj[i] = -A1[0]; Bj[i-1] = -A0[0]; + } + } + if(i===0) { + A0 = x[0]; + for(j=n-1;j>=1;--j) { + ret[j][0] = -A0[j]; + --j; + ret[j][0] = -A0[j]; + } + if(j===0) { ret[0][0] = -A0[0]; } + } + return ret; +} + +numeric._random = function _random(s,k) { + var i,n=s[k],ret=Array(n), rnd; + if(k === s.length-1) { + rnd = Math.random; + for(i=n-1;i>=1;i-=2) { + ret[i] = rnd(); + ret[i-1] = rnd(); + } + if(i===0) { ret[0] = rnd(); } + return ret; + } + for(i=n-1;i>=0;i--) ret[i] = _random(s,k+1); + return ret; +} +numeric.random = function random(s) { return numeric._random(s,0); } + +numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); } + +numeric.linspace = function linspace(a,b,n) { + if(typeof n === "undefined") n = Math.max(Math.round(b-a)+1,1); + if(n<2) { return n===1?[a]:[]; } + var i,ret = Array(n); + n--; + for(i=n;i>=0;i--) { ret[i] = (i*b+(n-i)*a)/n; } + return ret; +} + +numeric.getBlock = function getBlock(x,from,to) { + var s = numeric.dim(x); + function foo(x,k) { + var i,a = from[k], n = to[k]-a, ret = Array(n); + if(k === s.length-1) { + for(i=n;i>=0;i--) { ret[i] = x[i+a]; } + return ret; + } + for(i=n;i>=0;i--) { ret[i] = foo(x[i+a],k+1); } + return ret; + } + return foo(x,0); +} + +numeric.setBlock = function setBlock(x,from,to,B) { + var s = numeric.dim(x); + function foo(x,y,k) { + var i,a = from[k], n = to[k]-a; + if(k === s.length-1) { for(i=n;i>=0;i--) { x[i+a] = y[i]; } } + for(i=n;i>=0;i--) { foo(x[i+a],y[i],k+1); } + } + foo(x,B,0); + return x; +} + +numeric.getRange = function getRange(A,I,J) { + var m = I.length, n = J.length; + var i,j; + var B = Array(m), Bi, AI; + for(i=m-1;i!==-1;--i) { + B[i] = Array(n); + Bi = B[i]; + AI = A[I[i]]; + for(j=n-1;j!==-1;--j) Bi[j] = AI[J[j]]; + } + return B; +} + +numeric.blockMatrix = function blockMatrix(X) { + var s = numeric.dim(X); + if(s.length<4) return numeric.blockMatrix([X]); + var m=s[0],n=s[1],M,N,i,j,Xij; + M = 0; N = 0; + for(i=0;i=0;i--) { + Ai = Array(n); + xi = x[i]; + for(j=n-1;j>=3;--j) { + Ai[j] = xi * y[j]; + --j; + Ai[j] = xi * y[j]; + --j; + Ai[j] = xi * y[j]; + --j; + Ai[j] = xi * y[j]; + } + while(j>=0) { Ai[j] = xi * y[j]; --j; } + A[i] = Ai; + } + return A; +} + +// 3. The Tensor type T +numeric.T = function T(x,y) { this.x = x; this.y = y; } +numeric.t = function t(x,y) { return new numeric.T(x,y); } + +numeric.Tbinop = function Tbinop(rr,rc,cr,cc,setup) { + var io = numeric.indexOf; + if(typeof setup !== "string") { + var k; + setup = ''; + for(k in numeric) { + if(numeric.hasOwnProperty(k) && (rr.indexOf(k)>=0 || rc.indexOf(k)>=0 || cr.indexOf(k)>=0 || cc.indexOf(k)>=0) && k.length>1) { + setup += 'var '+k+' = numeric.'+k+';\n'; + } + } + } + return Function(['y'], + 'var x = this;\n'+ + 'if(!(y instanceof numeric.T)) { y = new numeric.T(y); }\n'+ + setup+'\n'+ + 'if(x.y) {'+ + ' if(y.y) {'+ + ' return new numeric.T('+cc+');\n'+ + ' }\n'+ + ' return new numeric.T('+cr+');\n'+ + '}\n'+ + 'if(y.y) {\n'+ + ' return new numeric.T('+rc+');\n'+ + '}\n'+ + 'return new numeric.T('+rr+');\n' + ); +} + +numeric.T.prototype.add = numeric.Tbinop( + 'add(x.x,y.x)', + 'add(x.x,y.x),y.y', + 'add(x.x,y.x),x.y', + 'add(x.x,y.x),add(x.y,y.y)'); +numeric.T.prototype.sub = numeric.Tbinop( + 'sub(x.x,y.x)', + 'sub(x.x,y.x),neg(y.y)', + 'sub(x.x,y.x),x.y', + 'sub(x.x,y.x),sub(x.y,y.y)'); +numeric.T.prototype.mul = numeric.Tbinop( + 'mul(x.x,y.x)', + 'mul(x.x,y.x),mul(x.x,y.y)', + 'mul(x.x,y.x),mul(x.y,y.x)', + 'sub(mul(x.x,y.x),mul(x.y,y.y)),add(mul(x.x,y.y),mul(x.y,y.x))'); + +numeric.T.prototype.reciprocal = function reciprocal() { + var mul = numeric.mul, div = numeric.div; + if(this.y) { + var d = numeric.add(mul(this.x,this.x),mul(this.y,this.y)); + return new numeric.T(div(this.x,d),div(numeric.neg(this.y),d)); + } + return new T(div(1,this.x)); +} +numeric.T.prototype.div = function div(y) { + if(!(y instanceof numeric.T)) y = new numeric.T(y); + if(y.y) { return this.mul(y.reciprocal()); } + var div = numeric.div; + if(this.y) { return new numeric.T(div(this.x,y.x),div(this.y,y.x)); } + return new numeric.T(div(this.x,y.x)); +} +numeric.T.prototype.dot = numeric.Tbinop( + 'dot(x.x,y.x)', + 'dot(x.x,y.x),dot(x.x,y.y)', + 'dot(x.x,y.x),dot(x.y,y.x)', + 'sub(dot(x.x,y.x),dot(x.y,y.y)),add(dot(x.x,y.y),dot(x.y,y.x))' + ); +numeric.T.prototype.transpose = function transpose() { + var t = numeric.transpose, x = this.x, y = this.y; + if(y) { return new numeric.T(t(x),t(y)); } + return new numeric.T(t(x)); +} +numeric.T.prototype.transjugate = function transjugate() { + var t = numeric.transpose, x = this.x, y = this.y; + if(y) { return new numeric.T(t(x),numeric.negtranspose(y)); } + return new numeric.T(t(x)); +} +numeric.Tunop = function Tunop(r,c,s) { + if(typeof s !== "string") { s = ''; } + return Function( + 'var x = this;\n'+ + s+'\n'+ + 'if(x.y) {'+ + ' '+c+';\n'+ + '}\n'+ + r+';\n' + ); +} + +numeric.T.prototype.exp = numeric.Tunop( + 'return new numeric.T(ex)', + 'return new numeric.T(mul(cos(x.y),ex),mul(sin(x.y),ex))', + 'var ex = numeric.exp(x.x), cos = numeric.cos, sin = numeric.sin, mul = numeric.mul;'); +numeric.T.prototype.conj = numeric.Tunop( + 'return new numeric.T(x.x);', + 'return new numeric.T(x.x,numeric.neg(x.y));'); +numeric.T.prototype.neg = numeric.Tunop( + 'return new numeric.T(neg(x.x));', + 'return new numeric.T(neg(x.x),neg(x.y));', + 'var neg = numeric.neg;'); +numeric.T.prototype.sin = numeric.Tunop( + 'return new numeric.T(numeric.sin(x.x))', + 'return x.exp().sub(x.neg().exp()).div(new numeric.T(0,2));'); +numeric.T.prototype.cos = numeric.Tunop( + 'return new numeric.T(numeric.cos(x.x))', + 'return x.exp().add(x.neg().exp()).div(2);'); +numeric.T.prototype.abs = numeric.Tunop( + 'return new numeric.T(numeric.abs(x.x));', + 'return new numeric.T(numeric.sqrt(numeric.add(mul(x.x,x.x),mul(x.y,x.y))));', + 'var mul = numeric.mul;'); +numeric.T.prototype.log = numeric.Tunop( + 'return new numeric.T(numeric.log(x.x));', + 'var theta = new numeric.T(numeric.atan2(x.y,x.x)), r = x.abs();\n'+ + 'return new numeric.T(numeric.log(r.x),theta.x);'); +numeric.T.prototype.norm2 = numeric.Tunop( + 'return numeric.norm2(x.x);', + 'var f = numeric.norm2Squared;\n'+ + 'return Math.sqrt(f(x.x)+f(x.y));'); +numeric.T.prototype.inv = function inv() { + var A = this; + if(typeof A.y === "undefined") { return new numeric.T(numeric.inv(A.x)); } + var n = A.x.length, i, j, k; + var Rx = numeric.identity(n),Ry = numeric.rep([n,n],0); + var Ax = numeric.clone(A.x), Ay = numeric.clone(A.y); + var Aix, Aiy, Ajx, Ajy, Rix, Riy, Rjx, Rjy; + var i,j,k,d,d1,ax,ay,bx,by,temp; + for(i=0;i d) { k=j; d = d1; } + } + if(k!==i) { + temp = Ax[i]; Ax[i] = Ax[k]; Ax[k] = temp; + temp = Ay[i]; Ay[i] = Ay[k]; Ay[k] = temp; + temp = Rx[i]; Rx[i] = Rx[k]; Rx[k] = temp; + temp = Ry[i]; Ry[i] = Ry[k]; Ry[k] = temp; + } + Aix = Ax[i]; Aiy = Ay[i]; + Rix = Rx[i]; Riy = Ry[i]; + ax = Aix[i]; ay = Aiy[i]; + for(j=i+1;j0;i--) { + Rix = Rx[i]; Riy = Ry[i]; + for(j=i-1;j>=0;j--) { + Rjx = Rx[j]; Rjy = Ry[j]; + ax = Ax[j][i]; ay = Ay[j][i]; + for(k=n-1;k>=0;k--) { + bx = Rix[k]; by = Riy[k]; + Rjx[k] -= ax*bx - ay*by; + Rjy[k] -= ax*by + ay*bx; + } + } + } + return new numeric.T(Rx,Ry); +} +numeric.T.prototype.get = function get(i) { + var x = this.x, y = this.y, k = 0, ik, n = i.length; + if(y) { + while(k= 0 ? 1 : -1; + var alpha = s*numeric.norm2(x); + v[0] += alpha; + var foo = numeric.norm2(v); + if(foo === 0) { /* this should not happen */ throw new Error('eig: internal error'); } + return numeric.div(v,foo); +} + +numeric.toUpperHessenberg = function toUpperHessenberg(me) { + var s = numeric.dim(me); + if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: toUpperHessenberg() only works on square matrices'); } + var m = s[0], i,j,k,x,v,A = numeric.clone(me),B,C,Ai,Ci,Q = numeric.identity(m),Qi; + for(j=0;j0) { + v = numeric.house(x); + B = numeric.getBlock(A,[j+1,j],[m-1,m-1]); + C = numeric.tensor(v,numeric.dot(v,B)); + for(i=j+1;i=4*det) { + var s1,s2; + s1 = 0.5*(tr+Math.sqrt(tr*tr-4*det)); + s2 = 0.5*(tr-Math.sqrt(tr*tr-4*det)); + Hloc = numeric.add(numeric.sub(numeric.dot(Hloc,Hloc), + numeric.mul(Hloc,s1+s2)), + numeric.diag(numeric.rep([3],s1*s2))); + } else { + Hloc = numeric.add(numeric.sub(numeric.dot(Hloc,Hloc), + numeric.mul(Hloc,tr)), + numeric.diag(numeric.rep([3],det))); + } + x = [Hloc[0][0],Hloc[1][0],Hloc[2][0]]; + v = numeric.house(x); + B = [H[0],H[1],H[2]]; + C = numeric.tensor(v,numeric.dot(v,B)); + for(i=0;i<3;i++) { Hi = H[i]; Ci = C[i]; for(k=0;k=0) { + if(p1<0) x = -0.5*(p1-sqrt(disc)); + else x = -0.5*(p1+sqrt(disc)); + n1 = (a-x)*(a-x)+b*b; + n2 = c*c+(d-x)*(d-x); + if(n1>n2) { + n1 = sqrt(n1); + p = (a-x)/n1; + q = b/n1; + } else { + n2 = sqrt(n2); + p = c/n2; + q = (d-x)/n2; + } + Q0 = new T([[q,-p],[p,q]]); + Q.setRows(i,j,Q0.dot(Q.getRows(i,j))); + } else { + x = -0.5*p1; + y = 0.5*sqrt(-disc); + n1 = (a-x)*(a-x)+b*b; + n2 = c*c+(d-x)*(d-x); + if(n1>n2) { + n1 = sqrt(n1+y*y); + p = (a-x)/n1; + q = b/n1; + x = 0; + y /= n1; + } else { + n2 = sqrt(n2+y*y); + p = c/n2; + q = (d-x)/n2; + x = y/n2; + y = 0; + } + Q0 = new T([[q,-p],[p,q]],[[x,y],[y,-x]]); + Q.setRows(i,j,Q0.dot(Q.getRows(i,j))); + } + } + } + var R = Q.dot(A).dot(Q.transjugate()), n = A.length, E = numeric.T.identity(n); + for(j=0;j0) { + for(k=j-1;k>=0;k--) { + var Rk = R.get([k,k]), Rj = R.get([j,j]); + if(numeric.neq(Rk.x,Rj.x) || numeric.neq(Rk.y,Rj.y)) { + x = R.getRow(k).getBlock([k],[j-1]); + y = E.getRow(j).getBlock([k],[j-1]); + E.set([j,k],(R.get([k,j]).neg().sub(x.dot(y))).div(Rk.sub(Rj))); + } else { + E.setRow(j,E.getRow(k)); + continue; + } + } + } + } + for(j=0;j=counts.length) counts[counts.length] = 0; + if(foo[j]!==0) counts[j]++; + } + } + var n = counts.length; + var Ai = Array(n+1); + Ai[0] = 0; + for(i=0;i= k11) { + xj[n] = j[m]; + if(m===0) return; + ++n; + --m; + km = k[m]; + k11 = k1[m]; + } else { + foo = Pinv[Aj[km]]; + if(x[foo] === 0) { + x[foo] = 1; + k[m] = km; + ++m; + j[m] = foo; + km = Ai[foo]; + k1[m] = k11 = Ai[foo+1]; + } else ++km; + } + } +} +numeric.ccsLPSolve = function ccsLPSolve(A,B,x,xj,I,Pinv,dfs) { + var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, n=0; + var Bi = B[0], Bj = B[1], Bv = B[2]; + + var i,i0,i1,j,J,j0,j1,k,l,l0,l1,a; + i0 = Bi[I]; + i1 = Bi[I+1]; + xj.length = 0; + for(i=i0;i a) { e = k; a = c; } + } + if(abs(x[i])= k11) { + xj[n] = Pinv[j[m]]; + if(m===0) return; + ++n; + --m; + km = k[m]; + k11 = k1[m]; + } else { + foo = Aj[km]; + if(x[foo] === 0) { + x[foo] = 1; + k[m] = km; + ++m; + j[m] = foo; + foo = Pinv[foo]; + km = Ai[foo]; + k1[m] = k11 = Ai[foo+1]; + } else ++km; + } + } +} +numeric.ccsLPSolve0 = function ccsLPSolve0(A,B,y,xj,I,Pinv,P,dfs) { + var Ai = A[0], Aj = A[1], Av = A[2],m = Ai.length-1, n=0; + var Bi = B[0], Bj = B[1], Bv = B[2]; + + var i,i0,i1,j,J,j0,j1,k,l,l0,l1,a; + i0 = Bi[I]; + i1 = Bi[I+1]; + xj.length = 0; + for(i=i0;i a) { e = k; a = c; } + } + if(abs(y[P[i]]) ret[k]) ret[k] = A.length; + var i; + for(i in A) { + if(A.hasOwnProperty(i)) dim(A[i],ret,k+1); + } + return ret; +}; + +numeric.sclone = function clone(A,k,n) { + if(typeof k === "undefined") { k=0; } + if(typeof n === "undefined") { n = numeric.sdim(A).length; } + var i,ret = Array(A.length); + if(k === n-1) { + for(i in A) { if(A.hasOwnProperty(i)) ret[i] = A[i]; } + return ret; + } + for(i in A) { + if(A.hasOwnProperty(i)) ret[i] = clone(A[i],k+1,n); + } + return ret; +} + +numeric.sdiag = function diag(d) { + var n = d.length,i,ret = Array(n),i1,i2,i3; + for(i=n-1;i>=1;i-=2) { + i1 = i-1; + ret[i] = []; ret[i][i] = d[i]; + ret[i1] = []; ret[i1][i1] = d[i1]; + } + if(i===0) { ret[0] = []; ret[0][0] = d[i]; } + return ret; +} + +numeric.sidentity = function identity(n) { return numeric.sdiag(numeric.rep([n],1)); } + +numeric.stranspose = function transpose(A) { + var ret = [], n = A.length, i,j,Ai; + for(i in A) { + if(!(A.hasOwnProperty(i))) continue; + Ai = A[i]; + for(j in Ai) { + if(!(Ai.hasOwnProperty(j))) continue; + if(typeof ret[j] !== "object") { ret[j] = []; } + ret[j][i] = Ai[j]; + } + } + return ret; +} + +numeric.sLUP = function LUP(A,tol) { + throw new Error("The function numeric.sLUP had a bug in it and has been removed. Please use the new numeric.ccsLUP function instead."); +}; + +numeric.sdotMM = function dotMM(A,B) { + var p = A.length, q = B.length, BT = numeric.stranspose(B), r = BT.length, Ai, BTk; + var i,j,k,accum; + var ret = Array(p),reti; + for(i=p-1;i>=0;i--) { + reti = []; + Ai = A[i]; + for(k=r-1;k>=0;k--) { + accum = 0; + BTk = BT[k]; + for(j in Ai) { + if(!(Ai.hasOwnProperty(j))) continue; + if(j in BTk) { accum += Ai[j]*BTk[j]; } + } + if(accum) reti[k] = accum; + } + ret[i] = reti; + } + return ret; +} + +numeric.sdotMV = function dotMV(A,x) { + var p = A.length, Ai, i,j; + var ret = Array(p), accum; + for(i=p-1;i>=0;i--) { + Ai = A[i]; + accum = 0; + for(j in Ai) { + if(!(Ai.hasOwnProperty(j))) continue; + if(x[j]) accum += Ai[j]*x[j]; + } + if(accum) ret[i] = accum; + } + return ret; +} + +numeric.sdotVM = function dotMV(x,A) { + var i,j,Ai,alpha; + var ret = [], accum; + for(i in x) { + if(!x.hasOwnProperty(i)) continue; + Ai = A[i]; + alpha = x[i]; + for(j in Ai) { + if(!Ai.hasOwnProperty(j)) continue; + if(!ret[j]) { ret[j] = 0; } + ret[j] += alpha*Ai[j]; + } + } + return ret; +} + +numeric.sdotVV = function dotVV(x,y) { + var i,ret=0; + for(i in x) { if(x[i] && y[i]) ret+= x[i]*y[i]; } + return ret; +} + +numeric.sdot = function dot(A,B) { + var m = numeric.sdim(A).length, n = numeric.sdim(B).length; + var k = m*1000+n; + switch(k) { + case 0: return A*B; + case 1001: return numeric.sdotVV(A,B); + case 2001: return numeric.sdotMV(A,B); + case 1002: return numeric.sdotVM(A,B); + case 2002: return numeric.sdotMM(A,B); + default: throw new Error('numeric.sdot not implemented for tensors of order '+m+' and '+n); + } +} + +numeric.sscatter = function scatter(V) { + var n = V[0].length, Vij, i, j, m = V.length, A = [], Aj; + for(i=n-1;i>=0;--i) { + if(!V[m-1][i]) continue; + Aj = A; + for(j=0;j=0;--i) ret[i] = []; + } + for(i=n;i>=0;--i) ret[i].push(k[i]); + ret[n+1].push(Ai); + } + } else gather(Ai,ret,k); + } + } + if(k.length>n) k.pop(); + return ret; +} + +// 6. Coordinate matrices +numeric.cLU = function LU(A) { + var I = A[0], J = A[1], V = A[2]; + var p = I.length, m=0, i,j,k,a,b,c; + for(i=0;im) m=I[i]; + m++; + var L = Array(m), U = Array(m), left = numeric.rep([m],Infinity), right = numeric.rep([m],-Infinity); + var Ui, Uj,alpha; + for(k=0;kright[i]) right[i] = j; + } + for(i=0;i right[i+1]) right[i+1] = right[i]; } + for(i=m-1;i>=1;i--) { if(left[i]=0;i--) { + while(Uj[k] > i) { + ret[i] -= Uv[k]*ret[Uj[k]]; + k--; + } + ret[i] /= Uv[k]; + k--; + } + return ret; +}; + +numeric.cgrid = function grid(n,shape) { + if(typeof n === "number") n = [n,n]; + var ret = numeric.rep(n,-1); + var i,j,count; + if(typeof shape !== "function") { + switch(shape) { + case 'L': + shape = function(i,j) { return (i>=n[0]/2 || jN) N = Ai[k]; } + N++; + ret = numeric.rep([N],0); + for(k=0;k1) { + mid = floor((p+q)/2); + if(x[mid] <= x0) p = mid; + else q = mid; + } + return this._at(x0,p); + } + var n = x0.length, i, ret = Array(n); + for(i=n-1;i!==-1;--i) ret[i] = this.at(x0[i]); + return ret; +} +numeric.Spline.prototype.diff = function diff() { + var x = this.x; + var yl = this.yl; + var yr = this.yr; + var kl = this.kl; + var kr = this.kr; + var n = yl.length; + var i,dx,dy; + var zl = kl, zr = kr, pl = Array(n), pr = Array(n); + var add = numeric.add, mul = numeric.mul, div = numeric.div, sub = numeric.sub; + for(i=n-1;i!==-1;--i) { + dx = x[i+1]-x[i]; + dy = sub(yr[i+1],yl[i]); + pl[i] = div(add(mul(dy, 6),mul(kl[i],-4*dx),mul(kr[i+1],-2*dx)),dx*dx); + pr[i+1] = div(add(mul(dy,-6),mul(kl[i], 2*dx),mul(kr[i+1], 4*dx)),dx*dx); + } + return new numeric.Spline(x,zl,zr,pl,pr); +} +numeric.Spline.prototype.roots = function roots() { + function sqr(x) { return x*x; } + function heval(y0,y1,k0,k1,x) { + var A = k0*2-(y1-y0); + var B = -k1*2+(y1-y0); + var t = (x+1)*0.5; + var s = t*(1-t); + return (1-t)*y0+t*y1+A*s*(1-t)+B*s*t; + } + var ret = []; + var x = this.x, yl = this.yl, yr = this.yr, kl = this.kl, kr = this.kr; + if(typeof yl[0] === "number") { + yl = [yl]; + yr = [yr]; + kl = [kl]; + kr = [kr]; + } + var m = yl.length,n=x.length-1,i,j,k,y,s,t; + 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; + var sqrt = Math.sqrt; + for(i=0;i!==m;++i) { + ai = yl[i]; + bi = yr[i]; + ci = kl[i]; + di = kr[i]; + ri = []; + for(j=0;j!==n;j++) { + if(j>0 && bi[j]*ai[j]<0) ri.push(x[j]); + dx = (x[j+1]-x[j]); + cx = x[j]; + y0 = ai[j]; + y1 = bi[j+1]; + k0 = ci[j]/dx; + k1 = di[j+1]/dx; + D = sqr(k0-k1+3*(y0-y1)) + 12*k1*y0; + A = k1+3*y0+2*k0-3*y1; + B = 3*(k1+k0+2*(y0-y1)); + if(D<=0) { + z0 = A/B; + if(z0>x[j] && z0x[j] && z0x[j] && z10) { + t0 = t1; + z0 = z1; + continue; + } + var side = 0; + while(1) { + tm = (z0*t1-z1*t0)/(z0-z1); + if(tm <= t0 || tm >= t1) { break; } + zm = this._at(tm,j); + if(zm*z1>0) { + t1 = tm; + z1 = zm; + if(side === -1) z0*=0.5; + side = -1; + } else if(zm*z0>0) { + t0 = tm; + z0 = zm; + if(side === 1) z1*=0.5; + side = 1; + } else break; + } + ri.push(tm); + t0 = stops[k+1]; + z0 = this._at(t0, j); + } + if(z1 === 0) ri.push(t1); + } + ret[i] = ri; + } + if(typeof this.yl[0] === "number") return ret[0]; + return ret; +} +numeric.spline = function spline(x,y,k1,kn) { + var n = x.length, b = [], dx = [], dy = []; + var i; + var sub = numeric.sub,mul = numeric.mul,add = numeric.add; + for(i=n-2;i>=0;i--) { dx[i] = x[i+1]-x[i]; dy[i] = sub(y[i+1],y[i]); } + if(typeof k1 === "string" || typeof kn === "string") { + k1 = kn = "periodic"; + } + // Build sparse tridiagonal system + var T = [[],[],[]]; + switch(typeof k1) { + case "undefined": + b[0] = mul(3/(dx[0]*dx[0]),dy[0]); + T[0].push(0,0); + T[1].push(0,1); + T[2].push(2/dx[0],1/dx[0]); + break; + case "string": + b[0] = add(mul(3/(dx[n-2]*dx[n-2]),dy[n-2]),mul(3/(dx[0]*dx[0]),dy[0])); + T[0].push(0,0,0); + T[1].push(n-2,0,1); + T[2].push(1/dx[n-2],2/dx[n-2]+2/dx[0],1/dx[0]); + break; + default: + b[0] = k1; + T[0].push(0); + T[1].push(0); + T[2].push(1); + break; + } + for(i=1;i20) { throw new Error("Numerical gradient fails"); } + x0[i] = x[i]+h; + f1 = f(x0); + x0[i] = x[i]-h; + f2 = f(x0); + x0[i] = x[i]; + if(isNaN(f1) || isNaN(f2)) { h/=16; continue; } + J[i] = (f1-f2)/(2*h); + t0 = x[i]-h; + t1 = x[i]; + t2 = x[i]+h; + d1 = (f1-f0)/h; + d2 = (f0-f2)/h; + N = max(abs(J[i]),abs(f0),abs(f1),abs(f2),abs(t0),abs(t1),abs(t2),1e-8); + errest = min(max(abs(d1-J[i]),abs(d2-J[i]),abs(d1-d2))/N,h/N); + if(errest>eps) { h/=16; } + else break; + } + } + return J; +} + +numeric.uncmin = function uncmin(f,x0,tol,gradient,maxit,callback,options) { + var grad = numeric.gradient; + if(typeof options === "undefined") { options = {}; } + if(typeof tol === "undefined") { tol = 1e-8; } + if(typeof gradient === "undefined") { gradient = function(x) { return grad(f,x); }; } + if(typeof maxit === "undefined") maxit = 1000; + x0 = numeric.clone(x0); + var n = x0.length; + var f0 = f(x0),f1,df0; + if(isNaN(f0)) throw new Error('uncmin: f(x0) is a NaN!'); + var max = Math.max, norm2 = numeric.norm2; + tol = max(tol,numeric.epsilon); + var step,g0,g1,H1 = options.Hinv || numeric.identity(n); + var dot = numeric.dot, inv = numeric.inv, sub = numeric.sub, add = numeric.add, ten = numeric.tensor, div = numeric.div, mul = numeric.mul; + var all = numeric.all, isfinite = numeric.isFinite, neg = numeric.neg; + var it=0,i,s,x1,y,Hy,Hs,ys,i0,t,nstep,t1,t2; + var msg = ""; + g0 = gradient(x0); + while(it= 0.1*t*df0 || isNaN(f1)) { + t *= 0.5; + ++it; + continue; + } + break; + } + if(t*nstep < tol) { msg = "Line search step size smaller than tol"; break; } + if(it === maxit) { msg = "maxit reached during line search"; break; } + g1 = gradient(x1); + y = sub(g1,g0); + ys = dot(y,s); + Hy = dot(H1,y); + H1 = sub(add(H1, + mul( + (ys+dot(y,Hy))/(ys*ys), + ten(s,s) )), + div(add(ten(Hy,s),ten(s,Hy)),ys)); + x0 = x1; + f0 = f1; + g0 = g1; + ++it; + } + return {solution: x0, f: f0, gradient: g0, invHessian: H1, iterations:it, message: msg}; +} + +// 10. Ode solver (Dormand-Prince) +numeric.Dopri = function Dopri(x,y,f,ymid,iterations,msg,events) { + this.x = x; + this.y = y; + this.f = f; + this.ymid = ymid; + this.iterations = iterations; + this.events = events; + this.message = msg; +} +numeric.Dopri.prototype._at = function _at(xi,j) { + function sqr(x) { return x*x; } + var sol = this; + var xs = sol.x; + var ys = sol.y; + var k1 = sol.f; + var ymid = sol.ymid; + var n = xs.length; + var x0,x1,xh,y0,y1,yh,xi; + var floor = Math.floor,h; + var c = 0.5; + var add = numeric.add, mul = numeric.mul,sub = numeric.sub, p,q,w; + x0 = xs[j]; + x1 = xs[j+1]; + y0 = ys[j]; + y1 = ys[j+1]; + h = x1-x0; + xh = x0+c*h; + yh = ymid[j]; + p = sub(k1[j ],mul(y0,1/(x0-xh)+2/(x0-x1))); + q = sub(k1[j+1],mul(y1,1/(x1-xh)+2/(x1-x0))); + w = [sqr(xi - x1) * (xi - xh) / sqr(x0 - x1) / (x0 - xh), + sqr(xi - x0) * sqr(xi - x1) / sqr(x0 - xh) / sqr(x1 - xh), + sqr(xi - x0) * (xi - xh) / sqr(x1 - x0) / (x1 - xh), + (xi - x0) * sqr(xi - x1) * (xi - xh) / sqr(x0-x1) / (x0 - xh), + (xi - x1) * sqr(xi - x0) * (xi - xh) / sqr(x0-x1) / (x1 - xh)]; + return add(add(add(add(mul(y0,w[0]), + mul(yh,w[1])), + mul(y1,w[2])), + mul( p,w[3])), + mul( q,w[4])); +} +numeric.Dopri.prototype.at = function at(x) { + var i,j,k,floor = Math.floor; + if(typeof x !== "number") { + var n = x.length, ret = Array(n); + for(i=n-1;i!==-1;--i) { + ret[i] = this.at(x[i]); + } + return ret; + } + var x0 = this.x; + i = 0; j = x0.length-1; + while(j-i>1) { + k = floor(0.5*(i+j)); + if(x0[k] <= x) i = k; + else j = k; + } + return this._at(x,i); +} + +numeric.dopri = function dopri(x0,x1,y0,f,tol,maxit,event) { + if(typeof tol === "undefined") { tol = 1e-6; } + if(typeof maxit === "undefined") { maxit = 1000; } + var xs = [x0], ys = [y0], k1 = [f(x0,y0)], k2,k3,k4,k5,k6,k7, ymid = []; + var A2 = 1/5; + var A3 = [3/40,9/40]; + var A4 = [44/45,-56/15,32/9]; + var A5 = [19372/6561,-25360/2187,64448/6561,-212/729]; + var A6 = [9017/3168,-355/33,46732/5247,49/176,-5103/18656]; + var b = [35/384,0,500/1113,125/192,-2187/6784,11/84]; + var bm = [0.5*6025192743/30085553152, + 0, + 0.5*51252292925/65400821598, + 0.5*-2691868925/45128329728, + 0.5*187940372067/1594534317056, + 0.5*-1776094331/19743644256, + 0.5*11237099/235043384]; + var c = [1/5,3/10,4/5,8/9,1,1]; + var e = [-71/57600,0,71/16695,-71/1920,17253/339200,-22/525,1/40]; + var i = 0,er,j; + var h = (x1-x0)/10; + var it = 0; + var add = numeric.add, mul = numeric.mul, y1,erinf; + var max = Math.max, min = Math.min, abs = Math.abs, norminf = numeric.norminf,pow = Math.pow; + var any = numeric.any, lt = numeric.lt, and = numeric.and, sub = numeric.sub; + var e0, e1, ev; + var ret = new numeric.Dopri(xs,ys,k1,ymid,-1,""); + if(typeof event === "function") e0 = event(x0,y0); + while(x0x1) h = x1-x0; + k2 = f(x0+c[0]*h, add(y0,mul( A2*h,k1[i]))); + k3 = f(x0+c[1]*h, add(add(y0,mul(A3[0]*h,k1[i])),mul(A3[1]*h,k2))); + 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))); + 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))); + 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))); + 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])); + k7 = f(x0+h,y1); + 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])); + if(typeof er === "number") erinf = abs(er); + else erinf = norminf(er); + if(erinf > tol) { // reject + h = 0.2*h*pow(tol/erinf,0.25); + if(x0+h === x0) { + ret.msg = "Step size became too small"; + break; + } + continue; + } + ymid[i] = add(add(add(add(add(add(y0, + mul(k1[i],h*bm[0])), + mul(k3 ,h*bm[2])), + mul(k4 ,h*bm[3])), + mul(k5 ,h*bm[4])), + mul(k6 ,h*bm[5])), + mul(k7 ,h*bm[6])); + ++i; + xs[i] = x0+h; + ys[i] = y1; + k1[i] = k7; + if(typeof event === "function") { + var yi,xl = x0,xr = x0+0.5*h,xi; + e1 = event(xr,ymid[i-1]); + ev = and(lt(e0,0),lt(0,e1)); + if(!any(ev)) { xl = xr; xr = x0+h; e0 = e1; e1 = event(xr,y1); ev = and(lt(e0,0),lt(0,e1)); } + if(any(ev)) { + var xc, yc, en,ei; + var side=0, sl = 1.0, sr = 1.0; + while(1) { + if(typeof e0 === "number") xi = (sr*e1*xl-sl*e0*xr)/(sr*e1-sl*e0); + else { + xi = xr; + for(j=e0.length-1;j!==-1;--j) { + if(e0[j]<0 && e1[j]>0) xi = min(xi,(sr*e1[j]*xl-sl*e0[j]*xr)/(sr*e1[j]-sl*e0[j])); + } + } + if(xi <= xl || xi >= xr) break; + yi = ret._at(xi, i-1); + ei = event(xi,yi); + en = and(lt(e0,0),lt(0,ei)); + if(any(en)) { + xr = xi; + e1 = ei; + ev = en; + sr = 1.0; + if(side === -1) sl *= 0.5; + else sl = 1.0; + side = -1; + } else { + xl = xi; + e0 = ei; + sl = 1.0; + if(side === 1) sr *= 0.5; + else sr = 1.0; + side = 1; + } + } + y1 = ret._at(0.5*(x0+xi),i-1); + ret.f[i] = f(xi,yi); + ret.x[i] = xi; + ret.y[i] = yi; + ret.ymid[i-1] = y1; + ret.events = ev; + ret.iterations = it; + return ret; + } + } + x0 += h; + y0 = y1; + e0 = e1; + h = min(0.8*h*pow(tol/erinf,0.25),4*h); + } + ret.iterations = it; + return ret; +} + +// 11. Ax = b +numeric.LU = function(A, fast) { + fast = fast || false; + + var abs = Math.abs; + var i, j, k, absAjk, Akk, Ak, Pk, Ai; + var max; + var n = A.length, n1 = n-1; + var P = new Array(n); + if(!fast) A = numeric.clone(A); + + for (k = 0; k < n; ++k) { + Pk = k; + Ak = A[k]; + max = abs(Ak[k]); + for (j = k + 1; j < n; ++j) { + absAjk = abs(A[j][k]); + if (max < absAjk) { + max = absAjk; + Pk = j; + } + } + P[k] = Pk; + + if (Pk != k) { + A[k] = A[Pk]; + A[Pk] = Ak; + Ak = A[k]; + } + + Akk = Ak[k]; + + for (i = k + 1; i < n; ++i) { + A[i][k] /= Akk; + } + + for (i = k + 1; i < n; ++i) { + Ai = A[i]; + for (j = k + 1; j < n1; ++j) { + Ai[j] -= Ai[k] * Ak[j]; + ++j; + Ai[j] -= Ai[k] * Ak[j]; + } + if(j===n1) Ai[j] -= Ai[k] * Ak[j]; + } + } + + return { + LU: A, + P: P + }; +} + +numeric.LUsolve = function LUsolve(LUP, b) { + var i, j; + var LU = LUP.LU; + var n = LU.length; + var x = numeric.clone(b); + var P = LUP.P; + var Pi, LUi, LUii, tmp; + + for (i=n-1;i!==-1;--i) x[i] = b[i]; + for (i = 0; i < n; ++i) { + Pi = P[i]; + if (P[i] !== i) { + tmp = x[i]; + x[i] = x[Pi]; + x[Pi] = tmp; + } + + LUi = LU[i]; + for (j = 0; j < i; ++j) { + x[i] -= x[j] * LUi[j]; + } + } + + for (i = n - 1; i >= 0; --i) { + LUi = LU[i]; + for (j = i + 1; j < n; ++j) { + x[i] -= x[j] * LUi[j]; + } + + x[i] /= LUi[i]; + } + + return x; +} + +numeric.solve = function solve(A,b,fast) { return numeric.LUsolve(numeric.LU(A,fast), b); } + +// 12. Linear programming +numeric.echelonize = function echelonize(A) { + var s = numeric.dim(A), m = s[0], n = s[1]; + var I = numeric.identity(m); + var P = Array(m); + var i,j,k,l,Ai,Ii,Z,a; + var abs = Math.abs; + var diveq = numeric.diveq; + A = numeric.clone(A); + for(i=0;ia1) alpha = a1; + g = add(c,mul(alpha,p)); + H = dot(A1,A0); + for(i=m-1;i!==-1;--i) H[i][i] += 1; + d = solve(H,div(g,alpha),true); + var t0 = div(z,dot(A,d)); + var t = 1.0; + for(i=n-1;i!==-1;--i) if(t0[i]<0) t = min(t,-0.999*t0[i]); + y = sub(x,mul(d,t)); + z = sub(b,dot(A,y)); + if(!all(gt(z,0))) return { solution: x, message: "", iterations: count }; + x = y; + if(alpha=0) unbounded = false; + else unbounded = true; + } + if(unbounded) return { solution: y, message: "Unbounded", iterations: count }; + } + return { solution: x, message: "maximum iteration count exceeded", iterations:count }; +} + +numeric._solveLP = function _solveLP(c,A,b,tol,maxit) { + var m = c.length, n = b.length,y; + var sum = numeric.sum, log = numeric.log, mul = numeric.mul, sub = numeric.sub, dot = numeric.dot, div = numeric.div, add = numeric.add; + var c0 = numeric.rep([m],0).concat([1]); + var J = numeric.rep([n,1],-1); + var A0 = numeric.blockMatrix([[A , J ]]); + var b0 = b; + var y = numeric.rep([m],0).concat(Math.max(0,numeric.sup(numeric.neg(b)))+1); + var x0 = numeric.__solveLP(c0,A0,b0,tol,maxit,y,false); + var x = numeric.clone(x0.solution); + x.length = m; + var foo = numeric.inf(sub(b,dot(A,x))); + if(foo<0) { return { solution: NaN, message: "Infeasible", iterations: x0.iterations }; } + var ret = numeric.__solveLP(c, A, b, tol, maxit-x0.iterations, x, true); + ret.iterations += x0.iterations; + return ret; +}; + +numeric.solveLP = function solveLP(c,A,b,Aeq,beq,tol,maxit) { + if(typeof maxit === "undefined") maxit = 1000; + if(typeof tol === "undefined") tol = numeric.epsilon; + if(typeof Aeq === "undefined") return numeric._solveLP(c,A,b,tol,maxit); + var m = Aeq.length, n = Aeq[0].length, o = A.length; + var B = numeric.echelonize(Aeq); + var flags = numeric.rep([n],0); + var P = B.P; + var Q = []; + var i; + for(i=P.length-1;i!==-1;--i) flags[P[i]] = 1; + for(i=n-1;i!==-1;--i) if(flags[i]===0) Q.push(i); + var g = numeric.getRange; + var I = numeric.linspace(0,m-1), J = numeric.linspace(0,o-1); + var Aeq2 = g(Aeq,I,Q), A1 = g(A,J,P), A2 = g(A,J,Q), dot = numeric.dot, sub = numeric.sub; + var A3 = dot(A1,B.I); + var A4 = sub(A2,dot(A3,Aeq2)), b4 = sub(b,dot(A3,beq)); + var c1 = Array(P.length), c2 = Array(Q.length); + for(i=P.length-1;i!==-1;--i) c1[i] = c[P[i]]; + for(i=Q.length-1;i!==-1;--i) c2[i] = c[Q[i]]; + var c4 = sub(c2,dot(c1,dot(B.I,Aeq2))); + var S = numeric._solveLP(c4,A4,b4,tol,maxit); + var x2 = S.solution; + if(x2!==x2) return S; + var x1 = dot(B.I,sub(beq,dot(Aeq2,x2))); + var x = Array(c.length); + for(i=P.length-1;i!==-1;--i) x[P[i]] = x1[i]; + for(i=Q.length-1;i!==-1;--i) x[Q[i]] = x2[i]; + return { solution: x, message:S.message, iterations: S.iterations }; +} + +numeric.MPStoLP = function MPStoLP(MPS) { + if(MPS instanceof String) { MPS.split('\n'); } + var state = 0; + var states = ['Initial state','NAME','ROWS','COLUMNS','RHS','BOUNDS','ENDATA']; + var n = MPS.length; + var i,j,z,N=0,rows = {}, sign = [], rl = 0, vars = {}, nv = 0; + var name; + var c = [], A = [], b = []; + function err(e) { throw new Error('MPStoLP: '+e+'\nLine '+i+': '+MPS[i]+'\nCurrent state: '+states[state]+'\n'); } + for(i=0;i +// +// Math.seedrandom('yipee'); Sets Math.random to a function that is +// initialized using the given explicit seed. +// +// Math.seedrandom(); Sets Math.random to a function that is +// seeded using the current time, dom state, +// and other accumulated local entropy. +// The generated seed string is returned. +// +// Math.seedrandom('yowza', true); +// Seeds using the given explicit seed mixed +// together with accumulated entropy. +// +// +// Seeds using physical random bits downloaded +// from random.org. +// +// Seeds using urandom bits from call.jsonlib.com, +// which is faster than random.org. +// +// Examples: +// +// Math.seedrandom("hello"); // Use "hello" as the seed. +// document.write(Math.random()); // Always 0.5463663768140734 +// document.write(Math.random()); // Always 0.43973793770592234 +// var rng1 = Math.random; // Remember the current prng. +// +// var autoseed = Math.seedrandom(); // New prng with an automatic seed. +// document.write(Math.random()); // Pretty much unpredictable. +// +// Math.random = rng1; // Continue "hello" prng sequence. +// document.write(Math.random()); // Always 0.554769432473455 +// +// Math.seedrandom(autoseed); // Restart at the previous seed. +// document.write(Math.random()); // Repeat the 'unpredictable' value. +// +// Notes: +// +// Each time seedrandom('arg') is called, entropy from the passed seed +// is accumulated in a pool to help generate future seeds for the +// zero-argument form of Math.seedrandom, so entropy can be injected over +// time by calling seedrandom with explicit data repeatedly. +// +// On speed - This javascript implementation of Math.random() is about +// 3-10x slower than the built-in Math.random() because it is not native +// code, but this is typically fast enough anyway. Seeding is more expensive, +// especially if you use auto-seeding. Some details (timings on Chrome 4): +// +// Our Math.random() - avg less than 0.002 milliseconds per call +// seedrandom('explicit') - avg less than 0.5 milliseconds per call +// seedrandom('explicit', true) - avg less than 2 milliseconds per call +// seedrandom() - avg about 38 milliseconds per call +// +// LICENSE (BSD): +// +// Copyright 2010 David Bau, all rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of this module nor the names of its contributors may +// be used to endorse or promote products derived from this software +// without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +/** + * All code is in an anonymous closure to keep the global namespace clean. + * + * @param {number=} overflow + * @param {number=} startdenom + */ + +// Patched by Seb so that seedrandom.js does not pollute the Math object. +// My tests suggest that doing Math.trouble = 1 makes Math lookups about 5% +// slower. +numeric.seedrandom = { pow:Math.pow, random:Math.random }; + +(function (pool, math, width, chunks, significance, overflow, startdenom) { + + +// +// seedrandom() +// This is the seedrandom function described above. +// +math['seedrandom'] = function seedrandom(seed, use_entropy) { + var key = []; + var arc4; + + // Flatten the seed string or build one from local entropy if needed. + seed = mixkey(flatten( + use_entropy ? [seed, pool] : + arguments.length ? seed : + [new Date().getTime(), pool, window], 3), key); + + // Use the seed to initialize an ARC4 generator. + arc4 = new ARC4(key); + + // Mix the randomness into accumulated entropy. + mixkey(arc4.S, pool); + + // Override Math.random + + // This function returns a random double in [0, 1) that contains + // randomness in every bit of the mantissa of the IEEE 754 value. + + math['random'] = function random() { // Closure to return a random double: + var n = arc4.g(chunks); // Start with a numerator n < 2 ^ 48 + var d = startdenom; // and denominator d = 2 ^ 48. + var x = 0; // and no 'extra last byte'. + while (n < significance) { // Fill up all significant digits by + n = (n + x) * width; // shifting numerator and + d *= width; // denominator and generating a + x = arc4.g(1); // new least-significant-byte. + } + while (n >= overflow) { // To avoid rounding up, before adding + n /= 2; // last byte, shift everything + d /= 2; // right using integer math until + x >>>= 1; // we have exactly the desired bits. + } + return (n + x) / d; // Form the number within [0, 1). + }; + + // Return the seed that was used + return seed; +}; + +// +// ARC4 +// +// An ARC4 implementation. The constructor takes a key in the form of +// an array of at most (width) integers that should be 0 <= x < (width). +// +// The g(count) method returns a pseudorandom integer that concatenates +// the next (count) outputs from ARC4. Its return value is a number x +// that is in the range 0 <= x < (width ^ count). +// +/** @constructor */ +function ARC4(key) { + var t, u, me = this, keylen = key.length; + var i = 0, j = me.i = me.j = me.m = 0; + me.S = []; + me.c = []; + + // The empty key [] is treated as [0]. + if (!keylen) { key = [keylen++]; } + + // Set up S using the standard key scheduling algorithm. + while (i < width) { me.S[i] = i++; } + for (i = 0; i < width; i++) { + t = me.S[i]; + j = lowbits(j + t + key[i % keylen]); + u = me.S[j]; + me.S[i] = u; + me.S[j] = t; + } + + // The "g" method returns the next (count) outputs as one number. + me.g = function getnext(count) { + var s = me.S; + var i = lowbits(me.i + 1); var t = s[i]; + var j = lowbits(me.j + t); var u = s[j]; + s[i] = u; + s[j] = t; + var r = s[lowbits(t + u)]; + while (--count) { + i = lowbits(i + 1); t = s[i]; + j = lowbits(j + t); u = s[j]; + s[i] = u; + s[j] = t; + r = r * width + s[lowbits(t + u)]; + } + me.i = i; + me.j = j; + return r; + }; + // For robust unpredictability discard an initial batch of values. + // See http://www.rsa.com/rsalabs/node.asp?id=2009 + me.g(width); +} + +// +// flatten() +// Converts an object tree to nested arrays of strings. +// +/** @param {Object=} result + * @param {string=} prop + * @param {string=} typ */ +function flatten(obj, depth, result, prop, typ) { + result = []; + typ = typeof(obj); + if (depth && typ == 'object') { + for (prop in obj) { + if (prop.indexOf('S') < 5) { // Avoid FF3 bug (local/sessionStorage) + try { result.push(flatten(obj[prop], depth - 1)); } catch (e) {} + } + } + } + return (result.length ? result : obj + (typ != 'string' ? '\0' : '')); +} + +// +// mixkey() +// Mixes a string seed into a key that is an array of integers, and +// returns a shortened string seed that is equivalent to the result key. +// +/** @param {number=} smear + * @param {number=} j */ +function mixkey(seed, key, smear, j) { + seed += ''; // Ensure the seed is a string + smear = 0; + for (j = 0; j < seed.length; j++) { + key[lowbits(j)] = + lowbits((smear ^= key[lowbits(j)] * 19) + seed.charCodeAt(j)); + } + seed = ''; + for (j in key) { seed += String.fromCharCode(key[j]); } + return seed; +} + +// +// lowbits() +// A quick "n mod width" for width a power of 2. +// +function lowbits(n) { return n & (width - 1); } + +// +// The following constants are related to IEEE 754 limits. +// +startdenom = math.pow(width, chunks); +significance = math.pow(2, significance); +overflow = significance * 2; + +// +// When seedrandom.js is loaded, we immediately mix a few bits +// from the built-in RNG into the entropy pool. Because we do +// not want to intefere with determinstic PRNG state later, +// seedrandom will not call math.random on its own again after +// initialization. +// +mixkey(math.random(), pool); + +// End anonymous scope, and pass initial values. +}( + [], // pool: entropy pool starts empty + numeric.seedrandom, // math: package containing random, pow, and seedrandom + 256, // width: each RC4 output is 0 <= x < 256 + 6, // chunks: at least six RC4 outputs for each double + 52 // significance: there are 52 significant digits in a double + )); +/* This file is a slightly modified version of quadprog.js from Alberto Santini. + * It has been slightly modified by Sébastien Loisel to make sure that it handles + * 0-based Arrays instead of 1-based Arrays. + * License is in resources/LICENSE.quadprog */ +(function(exports) { + +function base0to1(A) { + if(typeof A !== "object") { return A; } + var ret = [], i,n=A.length; + for(i=0;i meq) { + work[l] = sum; + } else { + work[l] = -Math.abs(sum); + if (sum > 0) { + for (j = 1; j <= n; j = j + 1) { + amat[j][i] = -amat[j][i]; + } + bvec[i] = -bvec[i]; + } + } + } + + for (i = 1; i <= nact; i = i + 1) { + work[iwsv + iact[i]] = 0; + } + + nvl = 0; + temp = 0; + for (i = 1; i <= q; i = i + 1) { + if (work[iwsv + i] < temp * work[iwnbv + i]) { + nvl = i; + temp = work[iwsv + i] / work[iwnbv + i]; + } + } + if (nvl === 0) { + return 999; + } + + return 0; + } + + function fn_goto_55() { + for (i = 1; i <= n; i = i + 1) { + sum = 0; + for (j = 1; j <= n; j = j + 1) { + sum = sum + dmat[j][i] * amat[j][nvl]; + } + work[i] = sum; + } + + l1 = iwzv; + for (i = 1; i <= n; i = i + 1) { + work[l1 + i] = 0; + } + for (j = nact + 1; j <= n; j = j + 1) { + for (i = 1; i <= n; i = i + 1) { + work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j]; + } + } + + t1inf = true; + for (i = nact; i >= 1; i = i - 1) { + sum = work[i]; + l = iwrm + (i * (i + 3)) / 2; + l1 = l - i; + for (j = i + 1; j <= nact; j = j + 1) { + sum = sum - work[l] * work[iwrv + j]; + l = l + j; + } + sum = sum / work[l1]; + work[iwrv + i] = sum; + if (iact[i] < meq) { + // continue; + break; + } + if (sum < 0) { + // continue; + break; + } + t1inf = false; + it1 = i; + } + + if (!t1inf) { + t1 = work[iwuv + it1] / work[iwrv + it1]; + for (i = 1; i <= nact; i = i + 1) { + if (iact[i] < meq) { + // continue; + break; + } + if (work[iwrv + i] < 0) { + // continue; + break; + } + temp = work[iwuv + i] / work[iwrv + i]; + if (temp < t1) { + t1 = temp; + it1 = i; + } + } + } + + sum = 0; + for (i = iwzv + 1; i <= iwzv + n; i = i + 1) { + sum = sum + work[i] * work[i]; + } + if (Math.abs(sum) <= vsmall) { + if (t1inf) { + ierr[1] = 1; + // GOTO 999 + return 999; + } else { + for (i = 1; i <= nact; i = i + 1) { + work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i]; + } + work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1; + // GOTO 700 + return 700; + } + } else { + sum = 0; + for (i = 1; i <= n; i = i + 1) { + sum = sum + work[iwzv + i] * amat[i][nvl]; + } + tt = -work[iwsv + nvl] / sum; + t2min = true; + if (!t1inf) { + if (t1 < tt) { + tt = t1; + t2min = false; + } + } + + for (i = 1; i <= n; i = i + 1) { + sol[i] = sol[i] + tt * work[iwzv + i]; + if (Math.abs(sol[i]) < vsmall) { + sol[i] = 0; + } + } + + crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]); + for (i = 1; i <= nact; i = i + 1) { + work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i]; + } + work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt; + + if (t2min) { + nact = nact + 1; + iact[nact] = nvl; + + l = iwrm + ((nact - 1) * nact) / 2 + 1; + for (i = 1; i <= nact - 1; i = i + 1) { + work[l] = work[i]; + l = l + 1; + } + + if (nact === n) { + work[l] = work[n]; + } else { + for (i = n; i >= nact + 1; i = i - 1) { + if (work[i] === 0) { + // continue; + break; + } + gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i])); + gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i])); + if (work[i - 1] >= 0) { + temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); + } else { + temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); + } + gc = work[i - 1] / temp; + gs = work[i] / temp; + + if (gc === 1) { + // continue; + break; + } + if (gc === 0) { + work[i - 1] = gs * temp; + for (j = 1; j <= n; j = j + 1) { + temp = dmat[j][i - 1]; + dmat[j][i - 1] = dmat[j][i]; + dmat[j][i] = temp; + } + } else { + work[i - 1] = temp; + nu = gs / (1 + gc); + for (j = 1; j <= n; j = j + 1) { + temp = gc * dmat[j][i - 1] + gs * dmat[j][i]; + dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i]; + dmat[j][i - 1] = temp; + + } + } + } + work[l] = work[nact]; + } + } else { + sum = -bvec[nvl]; + for (j = 1; j <= n; j = j + 1) { + sum = sum + sol[j] * amat[j][nvl]; + } + if (nvl > meq) { + work[iwsv + nvl] = sum; + } else { + work[iwsv + nvl] = -Math.abs(sum); + if (sum > 0) { + for (j = 1; j <= n; j = j + 1) { + amat[j][nvl] = -amat[j][nvl]; + } + bvec[nvl] = -bvec[nvl]; + } + } + // GOTO 700 + return 700; + } + } + + return 0; + } + + function fn_goto_797() { + l = iwrm + (it1 * (it1 + 1)) / 2 + 1; + l1 = l + it1; + if (work[l1] === 0) { + // GOTO 798 + return 798; + } + gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1])); + gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1])); + if (work[l1 - 1] >= 0) { + temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); + } else { + temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc))); + } + gc = work[l1 - 1] / temp; + gs = work[l1] / temp; + + if (gc === 1) { + // GOTO 798 + return 798; + } + if (gc === 0) { + for (i = it1 + 1; i <= nact; i = i + 1) { + temp = work[l1 - 1]; + work[l1 - 1] = work[l1]; + work[l1] = temp; + l1 = l1 + i; + } + for (i = 1; i <= n; i = i + 1) { + temp = dmat[i][it1]; + dmat[i][it1] = dmat[i][it1 + 1]; + dmat[i][it1 + 1] = temp; + } + } else { + nu = gs / (1 + gc); + for (i = it1 + 1; i <= nact; i = i + 1) { + temp = gc * work[l1 - 1] + gs * work[l1]; + work[l1] = nu * (work[l1 - 1] + temp) - work[l1]; + work[l1 - 1] = temp; + l1 = l1 + i; + } + for (i = 1; i <= n; i = i + 1) { + temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1]; + dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1]; + dmat[i][it1] = temp; + } + } + + return 0; + } + + function fn_goto_798() { + l1 = l - it1; + for (i = 1; i <= it1; i = i + 1) { + work[l1] = work[l]; + l = l + 1; + l1 = l1 + 1; + } + + work[iwuv + it1] = work[iwuv + it1 + 1]; + iact[it1] = iact[it1 + 1]; + it1 = it1 + 1; + if (it1 < nact) { + // GOTO 797 + return 797; + } + + return 0; + } + + function fn_goto_799() { + work[iwuv + nact] = work[iwuv + nact + 1]; + work[iwuv + nact + 1] = 0; + iact[nact] = 0; + nact = nact - 1; + iter[2] = iter[2] + 1; + + return 0; + } + + go = 0; + while (true) { + go = fn_goto_50(); + if (go === 999) { + return; + } + while (true) { + go = fn_goto_55(); + if (go === 0) { + break; + } + if (go === 999) { + return; + } + if (go === 700) { + if (it1 === nact) { + fn_goto_799(); + } else { + while (true) { + fn_goto_797(); + go = fn_goto_798(); + if (go !== 797) { + break; + } + } + fn_goto_799(); + } + } + } + } + +} + +function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) { + Dmat = base0to1(Dmat); + dvec = base0to1(dvec); + Amat = base0to1(Amat); + var i, n, q, + nact, r, + crval = [], iact = [], sol = [], work = [], iter = [], + message; + + meq = meq || 0; + factorized = factorized ? base0to1(factorized) : [undefined, 0]; + bvec = bvec ? base0to1(bvec) : []; + + // In Fortran the array index starts from 1 + n = Dmat.length - 1; + q = Amat[1].length - 1; + + if (!bvec) { + for (i = 1; i <= q; i = i + 1) { + bvec[i] = 0; + } + } + for (i = 1; i <= q; i = i + 1) { + iact[i] = 0; + } + nact = 0; + r = Math.min(n, q); + for (i = 1; i <= n; i = i + 1) { + sol[i] = 0; + } + crval[1] = 0; + for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) { + work[i] = 0; + } + for (i = 1; i <= 2; i = i + 1) { + iter[i] = 0; + } + + qpgen2(Dmat, dvec, n, n, sol, crval, Amat, + bvec, n, q, meq, iact, nact, iter, work, factorized); + + message = ""; + if (factorized[1] === 1) { + message = "constraints are inconsistent, no solution!"; + } + if (factorized[1] === 2) { + message = "matrix D in quadratic function is not positive definite!"; + } + + return { + solution: base1to0(sol), + value: base1to0(crval), + unconstrained_solution: base1to0(dvec), + iterations: base1to0(iter), + iact: base1to0(iact), + message: message + }; +} +exports.solveQP = solveQP; +}(numeric)); +/* +Shanti Rao sent me this routine by private email. I had to modify it +slightly to work on Arrays instead of using a Matrix object. +It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py +*/ + +numeric.svd= function svd(A) { + var temp; +//Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970) + var prec= numeric.epsilon; //Math.pow(2,-52) // assumes double prec + var tolerance= 1.e-64/prec; + var itmax= 50; + var c=0; + var i=0; + var j=0; + var k=0; + var l=0; + + var u= numeric.clone(A); + var m= u.length; + + var n= u[0].length; + + if (m < n) throw "Need more rows than columns" + + var e = new Array(n); + var q = new Array(n); + for (i=0; i b) + return a*Math.sqrt(1.0+(b*b/a/a)) + else if (b == 0.0) + return a + return b*Math.sqrt(1.0+(a*a/b/b)) + } + + //Householder's reduction to bidiagonal form + + var f= 0.0; + var g= 0.0; + var h= 0.0; + var x= 0.0; + var y= 0.0; + var z= 0.0; + var s= 0.0; + + for (i=0; i < n; i++) + { + e[i]= g; + s= 0.0; + l= i+1; + for (j=i; j < m; j++) + s += (u[j][i]*u[j][i]); + if (s <= tolerance) + g= 0.0; + else + { + f= u[i][i]; + g= Math.sqrt(s); + if (f >= 0.0) g= -g; + h= f*g-s + u[i][i]=f-g; + for (j=l; j < n; j++) + { + s= 0.0 + for (k=i; k < m; k++) + s += u[k][i]*u[k][j] + f= s/h + for (k=i; k < m; k++) + u[k][j]+=f*u[k][i] + } + } + q[i]= g + s= 0.0 + for (j=l; j < n; j++) + s= s + u[i][j]*u[i][j] + if (s <= tolerance) + g= 0.0 + else + { + f= u[i][i+1] + g= Math.sqrt(s) + if (f >= 0.0) g= -g + h= f*g - s + u[i][i+1] = f-g; + for (j=l; j < n; j++) e[j]= u[i][j]/h + for (j=l; j < m; j++) + { + s=0.0 + for (k=l; k < n; k++) + s += (u[j][k]*u[i][k]) + for (k=l; k < n; k++) + u[j][k]+=s*e[k] + } + } + y= Math.abs(q[i])+Math.abs(e[i]) + if (y>x) + x=y + } + + // accumulation of right hand gtransformations + for (i=n-1; i != -1; i+= -1) + { + if (g != 0.0) + { + h= g*u[i][i+1] + for (j=l; j < n; j++) + v[j][i]=u[i][j]/h + for (j=l; j < n; j++) + { + s=0.0 + for (k=l; k < n; k++) + s += u[i][k]*v[k][j] + for (k=l; k < n; k++) + v[k][j]+=(s*v[k][i]) + } + } + for (j=l; j < n; j++) + { + v[i][j] = 0; + v[j][i] = 0; + } + v[i][i] = 1; + g= e[i] + l= i + } + + // accumulation of left hand transformations + for (i=n-1; i != -1; i+= -1) + { + l= i+1 + g= q[i] + for (j=l; j < n; j++) + u[i][j] = 0; + if (g != 0.0) + { + h= u[i][i]*g + for (j=l; j < n; j++) + { + s=0.0 + for (k=l; k < m; k++) s += u[k][i]*u[k][j]; + f= s/h + for (k=i; k < m; k++) u[k][j]+=f*u[k][i]; + } + for (j=i; j < m; j++) u[j][i] = u[j][i]/g; + } + else + for (j=i; j < m; j++) u[j][i] = 0; + u[i][i] += 1; + } + + // diagonalization of the bidiagonal form + prec= prec*x + for (k=n-1; k != -1; k+= -1) + { + for (var iteration=0; iteration < itmax; iteration++) + { // test f splitting + var test_convergence = false + for (l=k; l != -1; l+= -1) + { + if (Math.abs(e[l]) <= prec) + { test_convergence= true + break + } + if (Math.abs(q[l-1]) <= prec) + break + } + if (!test_convergence) + { // cancellation of e[l] if l>0 + c= 0.0 + s= 1.0 + var l1= l-1 + for (i =l; i= itmax-1) + throw 'Error: no convergence.' + // shift from bottom 2x2 minor + x= q[l] + y= q[k-1] + g= e[k-1] + h= e[k] + f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) + g= pythag(f,1.0) + if (f < 0.0) + f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x + else + f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x + // next QR transformation + c= 1.0 + s= 1.0 + for (i=l+1; i< k+1; i++) + { + g= e[i] + y= q[i] + h= s*g + g= c*g + z= pythag(f,h) + e[i-1]= z + c= f/z + s= h/z + f= x*c+g*s + g= -x*s+g*c + h= y*s + y= y*c + for (j=0; j < n; j++) + { + x= v[j][i-1] + z= v[j][i] + v[j][i-1] = x*c+z*s + v[j][i] = -x*s+z*c + } + z= pythag(f,h) + q[i-1]= z + c= f/z + s= h/z + f= c*g+s*y + x= -s*g+c*y + for (j=0; j < m; j++) + { + y= u[j][i-1] + z= u[j][i] + u[j][i-1] = y*c+z*s + u[j][i] = -y*s+z*c + } + } + e[l]= 0.0 + e[k]= f + q[k]= x + } + } + + //vt= transpose(v) + //return (u,q,vt) + for (i=0;i= 0; j--) + { + if (q[j] < q[i]) + { + // writeln(i,'-',j) + c = q[j] + q[j] = q[i] + q[i] = c + for(k=0;k3;o-=4)t(),t(),t(),t();while(o>0)t(),o--;i=new Date;if(i-r>n)break}for(o=s;o>3;o-=4)t(),t(),t(),t();while(o>0)t(),o--;return i=new Date,1e3*(3*s-1)/(i-r)},numeric._myIndexOf=function(t){var n=this.length,r;for(r=0;rnumeric.largeArray)return r.push("...Large Array..."),!0;var f=!1;r.push("[");for(t=0;t0&&(r.push(","),f&&r.push("\n ")),f=i(e[t]);return r.push("]"),!0}r.push("{");var f=!1;for(t in e)e.hasOwnProperty(t)&&(f&&r.push(",\n"),f=!0,r.push(t),r.push(": \n"),i(e[t]));return r.push("}"),!0}var r=[];return i(t),r.join("")},numeric.parseDate=function(t){function n(e){if(typeof e=="string")return Date.parse(e.replace(/-/g,"/"));if(e instanceof Array){var t=[],r;for(r=0;r0){s[f]=[];for(r=0;r>2,u=((r&3)<<4)+(i>>4),a=((i&15)<<2)+(s>>6),f=s&63,n+1>=t?a=f=64:n+2>=t&&(f=64),c+=l.charAt(o)+l.charAt(u)+l.charAt(a)+l.charAt(f);return c}function r(e,t,n){typeof t=="undefined"&&(t=0),typeof n=="undefined"&&(n=e.length);var r=[0,1996959894,3993919788,2567524794,124634137,1886057615,3915621685,2657392035,249268274,2044508324,3772115230,2547177864,162941995,2125561021,3887607047,2428444049,498536548,1789927666,4089016648,2227061214,450548861,1843258603,4107580753,2211677639,325883990,1684777152,4251122042,2321926636,335633487,1661365465,4195302755,2366115317,997073096,1281953886,3579855332,2724688242,1006888145,1258607687,3524101629,2768942443,901097722,1119000684,3686517206,2898065728,853044451,1172266101,3705015759,2882616665,651767980,1373503546,3369554304,3218104598,565507253,1454621731,3485111705,3099436303,671266974,1594198024,3322730930,2970347812,795835527,1483230225,3244367275,3060149565,1994146192,31158534,2563907772,4023717930,1907459465,112637215,2680153253,3904427059,2013776290,251722036,2517215374,3775830040,2137656763,141376813,2439277719,3865271297,1802195444,476864866,2238001368,4066508878,1812370925,453092731,2181625025,4111451223,1706088902,314042704,2344532202,4240017532,1658658271,366619977,2362670323,4224994405,1303535960,984961486,2747007092,3569037538,1256170817,1037604311,2765210733,3554079995,1131014506,879679996,2909243462,3663771856,1141124467,855842277,2852801631,3708648649,1342533948,654459306,3188396048,3373015174,1466479909,544179635,3110523913,3462522015,1591671054,702138776,2966460450,3352799412,1504918807,783551873,3082640443,3233442989,3988292384,2596254646,62317068,1957810842,3939845945,2647816111,81470997,1943803523,3814918930,2489596804,225274430,2053790376,3826175755,2466906013,167816743,2097651377,4027552580,2265490386,503444072,1762050814,4150417245,2154129355,426522225,1852507879,4275313526,2312317920,282753626,1742555852,4189708143,2394877945,397917763,1622183637,3604390888,2714866558,953729732,1340076626,3518719985,2797360999,1068828381,1219638859,3624741850,2936675148,906185462,1090812512,3747672003,2825379669,829329135,1181335161,3412177804,3160834842,628085408,1382605366,3423369109,3138078467,570562233,1426400815,3317316542,2998733608,733239954,1555261956,3268935591,3050360625,752459403,1541320221,2607071920,3965973030,1969922972,40735498,2617837225,3943577151,1913087877,83908371,2512341634,3803740692,2075208622,213261112,2463272603,3855990285,2094854071,198958881,2262029012,4057260610,1759359992,534414190,2176718541,4139329115,1873836001,414664567,2282248934,4279200368,1711684554,285281116,2405801727,4167216745,1634467795,376229701,2685067896,3608007406,1308918612,956543938,2808555105,3495958263,1231636301,1047427035,2932959818,3654703836,1088359270,936918e3,2847714899,3736837829,1202900863,817233897,3183342108,3401237130,1404277552,615818150,3134207493,3453421203,1423857449,601450431,3009837614,3294710456,1567103746,711928724,3020668471,3272380065,1510334235,755167117],i=-1,s=0,o=e.length,u;for(u=t;u>>8^r[s];return i^-1}var i=t[0].length,s=t[0][0].length,o,u,a,f,l,c,h,p,d,v,m,g=[137,80,78,71,13,10,26,10,0,0,0,13,73,72,68,82,s>>24&255,s>>16&255,s>>8&255,s&255,i>>24&255,i>>16&255,i>>8&255,i&255,8,2,0,0,0,-1,-2,-3,-4,-5,-6,-7,-8,73,68,65,84,8,29];m=r(g,12,29),g[29]=m>>24&255,g[30]=m>>16&255,g[31]=m>>8&255,g[32]=m&255,o=1,u=0;for(p=0;p>8&255,g.push(c),g.push(h),g.push(~c&255),g.push(~h&255),p===0&&g.push(0);for(d=0;d255?c=255:c<0?c=0:c=Math.round(c),o=(o+c)%65521,u=(u+o)%65521,g.push(c);g.push(0)}return v=(u<<16)+o,g.push(v>>24&255),g.push(v>>16&255),g.push(v>>8&255),g.push(v&255),l=g.length-41,g[33]=l>>24&255,g[34]=l>>16&255,g[35]=l>>8&255,g[36]=l&255,m=r(g,37),g.push(m>>24&255),g.push(m>>16&255),g.push(m>>8&255),g.push(m&255),g.push(0),g.push(0),g.push(0),g.push(0),g.push(73),g.push(69),g.push(78),g.push(68),g.push(174),g.push(66),g.push(96),g.push(130),"data:image/png;base64,"+n(g)},numeric._dim=function(t){var n=[];while(typeof t=="object")n.push(t.length),t=t[0];return n},numeric.dim=function(t){var n,r;if(typeof t=="object")return n=t[0],typeof n=="object"?(r=n[0],typeof r=="object"?numeric._dim(t):[t.length,n.length]):[t.length];return[]},numeric.mapreduce=function(t,n){return Function("x","accum","_s","_k",'if(typeof accum === "undefined") accum = '+n+";\n"+'if(typeof x === "number") { var xi = x; '+t+"; return accum; }\n"+'if(typeof _s === "undefined") _s = numeric.dim(x);\n'+'if(typeof _k === "undefined") _k = 0;\n'+"var _n = _s[_k];\n"+"var i,xi;\n"+"if(_k < _s.length-1) {\n"+" for(i=_n-1;i>=0;i--) {\n"+" accum = arguments.callee(x[i],accum,_s,_k+1);\n"+" }"+" return accum;\n"+"}\n"+"for(i=_n-1;i>=1;i-=2) { \n"+" xi = x[i];\n"+" "+t+";\n"+" xi = x[i-1];\n"+" "+t+";\n"+"}\n"+"if(i === 0) {\n"+" xi = x[i];\n"+" "+t+"\n"+"}\n"+"return accum;")},numeric.mapreduce2=function(t,n){return Function("x","var n = x.length;\nvar i,xi;\n"+n+";\n"+"for(i=n-1;i!==-1;--i) { \n"+" xi = x[i];\n"+" "+t+";\n"+"}\n"+"return accum;")},numeric.same=function same(e,t){var n,r;if(e instanceof Array&&t instanceof Array){r=e.length;if(r!==t.length)return!1;for(n=0;n=0;o-=2)s[o+1]=n,s[o]=n;return o===-1&&(s[0]=n),s}for(o=i-1;o>=0;o--)s[o]=numeric.rep(t,n,r+1);return s},numeric.dotMMsmall=function(t,n){var r,i,s,o,u,a,f,l,c,h,p,d,v,m;o=t.length,u=n.length,a=n[0].length,f=Array(o);for(r=o-1;r>=0;r--){l=Array(a),c=t[r];for(s=a-1;s>=0;s--){h=c[u-1]*n[u-1][s];for(i=u-2;i>=1;i-=2)p=i-1,h+=c[i]*n[i][s]+c[p]*n[p][s];i===0&&(h+=c[0]*n[0][s]),l[s]=h}f[r]=l}return f},numeric._getCol=function(t,n,r){var i=t.length,s;for(s=i-1;s>0;--s)r[s]=t[s][n],--s,r[s]=t[s][n];s===0&&(r[0]=t[0][n])},numeric.dotMMbig=function(t,n){var r=numeric._getCol,i=n.length,s=Array(i),o=t.length,u=n[0].length,a=new Array(o),f,l=numeric.dotVV,c,h,p,d;--i,--o;for(c=o;c!==-1;--c)a[c]=Array(u);--u;for(c=u;c!==-1;--c){r(n,c,s);for(h=o;h!==-1;--h)d=0,f=t[h],a[h][c]=l(f,s)}return a},numeric.dotMV=function(t,n){var r=t.length,i=n.length,s,o=Array(r),u=numeric.dotVV;for(s=r-1;s>=0;s--)o[s]=u(t[s],n);return o},numeric.dotVM=function(t,n){var r,i,s,o,u,a,f,l,c,h,p,d,v,m,g,y,b,w,E;o=t.length,u=n[0].length,f=Array(u);for(s=u-1;s>=0;s--){h=t[o-1]*n[o-1][s];for(i=o-2;i>=1;i-=2)p=i-1,h+=t[i]*n[i][s]+t[p]*n[p][s];i===0&&(h+=t[0]*n[0][s]),f[s]=h}return f},numeric.dotVV=function(t,n){var r,i=t.length,s,o=t[i-1]*n[i-1];for(r=i-2;r>=1;r-=2)s=r-1,o+=t[r]*n[r]+t[s]*n[s];return r===0&&(o+=t[0]*n[0]),o},numeric.dot=function(t,n){var r=numeric.dim;switch(r(t).length*1e3+r(n).length){case 2002:return n.length<10?numeric.dotMMsmall(t,n):numeric.dotMMbig(t,n);case 2001:return numeric.dotMV(t,n);case 1002:return numeric.dotVM(t,n);case 1001:return numeric.dotVV(t,n);case 1e3:return numeric.mulVS(t,n);case 1:return numeric.mulSV(t,n);case 0:return t*n;default:throw new Error("numeric.dot only works on vectors and matrices")}},numeric.diag=function(t){var n,r,i,s=t.length,o=Array(s),u;for(n=s-1;n>=0;n--){u=Array(s),r=n+2;for(i=s-1;i>=r;i-=2)u[i]=0,u[i-1]=0;i>n&&(u[i]=0),u[n]=t[n];for(i=n-1;i>=1;i-=2)u[i]=0,u[i-1]=0;i===0&&(u[0]=0),o[n]=u}return o},numeric.getDiag=function(e){var t=Math.min(e.length,e[0].length),n,r=Array(t);for(n=t-1;n>=1;--n)r[n]=e[n][n],--n,r[n]=e[n][n];return n===0&&(r[0]=e[0][0]),r},numeric.identity=function(t){return numeric.diag(numeric.rep([t],1))},numeric.pointwise=function(t,n,r){typeof r=="undefined"&&(r="");var i=[],s,o=/\[i\]$/,u,a="",f=!1;for(s=0;s=0;i--) ret[i] = arguments.callee("+t.join(",")+",_s,_k+1);\n"+" return ret;\n"+"}\n"+r+"\n"+"for(i=_n-1;i!==-1;--i) {\n"+" "+n+"\n"+"}\n"+"return ret;",Function.apply(null,i)},numeric.pointwise2=function(t,n,r){typeof r=="undefined"&&(r="");var i=[],s,o=/\[i\]$/,u,a="",f=!1;for(s=0;s=0;s--)_biforeach(typeof e=="object"?e[s]:e,typeof t=="object"?t[s]:t,n,r+1,i)},numeric._biforeach2=function _biforeach2(e,t,n,r,i){if(r===n.length-1)return i(e,t);var s,o=n[r],u=Array(o);for(s=o-1;s>=0;--s)u[s]=_biforeach2(typeof e=="object"?e[s]:e,typeof t=="object"?t[s]:t,n,r+1,i);return u},numeric._foreach=function _foreach(e,t,n,r){if(n===t.length-1){r(e);return}var i,s=t[n];for(i=s-1;i>=0;i--)_foreach(e[i],t,n+1,r)},numeric._foreach2=function _foreach2(e,t,n,r){if(n===t.length-1)return r(e);var i,s=t[n],o=Array(s);for(i=s-1;i>=0;i--)o[i]=_foreach2(e[i],t,n+1,r);return o},numeric.ops2={add:"+",sub:"-",mul:"*",div:"/",mod:"%",and:"&&",or:"||",eq:"===",neq:"!==",lt:"<",gt:">",leq:"<=",geq:">=",band:"&",bor:"|",bxor:"^",lshift:"<<",rshift:">>",rrshift:">>>"},numeric.opseq={addeq:"+=",subeq:"-=",muleq:"*=",diveq:"/=",modeq:"%=",lshifteq:"<<=",rshifteq:">>=",rrshifteq:">>>=",bandeq:"&=",boreq:"|=",bxoreq:"^="},numeric.mathfuns=["abs","acos","asin","atan","ceil","cos","exp","floor","log","round","sin","sqrt","tan","isNaN","isFinite"],numeric.mathfuns2=["atan2","pow","max","min"],numeric.ops1={neg:"-",not:"!",bnot:"~",clone:""},numeric.mapreducers={any:["if(xi) return true;","var accum = false;"],all:["if(!xi) return false;","var accum = true;"],sum:["accum += xi;","var accum = 0;"],prod:["accum *= xi;","var accum = 1;"],norm2Squared:["accum += xi*xi;","var accum = 0;"],norminf:["accum = max(accum,abs(xi));","var accum = 0, max = Math.max, abs = Math.abs;"],norm1:["accum += abs(xi)","var accum = 0, abs = Math.abs;"],sup:["accum = max(accum,xi);","var accum = -Infinity, max = Math.max;"],inf:["accum = min(accum,xi);","var accum = Infinity, min = Math.min;"]},function(){var e,t;for(e=0;em&&(v=h,m=d);a=o[v],o[v]=o[p],o[p]=a,c=f[v],f[v]=f[p],f[p]=c,t=a[p];for(d=p;d!==s;++d)a[d]/=t;for(d=s-1;d!==-1;--d)c[d]/=t;for(h=i-1;h!==-1;--h)if(h!==p){u=o[h],l=f[h],t=u[p];for(d=p+1;d!==s;++d)u[d]-=a[d]*t;for(d=s-1;d>0;--d)l[d]-=c[d]*t,--d,l[d]-=c[d]*t;d===0&&(l[0]-=c[0]*t)}}return f},numeric.det=function(t){var n=numeric.dim(t);if(n.length!==2||n[0]!==n[1])throw new Error("numeric: det() only works on square matrices");var r=n[0],i=1,s,o,u,a=numeric.clone(t),f,l,c,h,p,d,v;for(o=0;oMath.abs(a[u][o])&&(u=s);u!==o&&(h=a[u],a[u]=a[o],a[o]=h,i*=-1),f=a[o];for(s=o+1;s=1;n-=2){a=t[n],u=t[n-1];for(r=s-1;r>=1;--r)f=o[r],f[n]=a[r],f[n-1]=u[r],--r,f=o[r],f[n]=a[r],f[n-1]=u[r];r===0&&(f=o[0],f[n]=a[0],f[n-1]=u[0])}if(n===0){u=t[0];for(r=s-1;r>=1;--r)o[r][0]=u[r],--r,o[r][0]=u[r];r===0&&(o[0][0]=u[0])}return o},numeric.negtranspose=function(t){var n,r,i=t.length,s=t[0].length,o=Array(s),u,a,f;for(r=0;r=1;n-=2){a=t[n],u=t[n-1];for(r=s-1;r>=1;--r)f=o[r],f[n]=-a[r],f[n-1]=-u[r],--r,f=o[r],f[n]=-a[r],f[n-1]=-u[r];r===0&&(f=o[0],f[n]=-a[0],f[n-1]=-u[0])}if(n===0){u=t[0];for(r=s-1;r>=1;--r)o[r][0]=-u[r],--r,o[r][0]=-u[r];r===0&&(o[0][0]=-u[0])}return o},numeric._random=function _random(e,t){var n,r=e[t],i=Array(r),s;if(t===e.length-1){s=Math.random;for(n=r-1;n>=1;n-=2)i[n]=s(),i[n-1]=s();return n===0&&(i[0]=s()),i}for(n=r-1;n>=0;n--)i[n]=_random(e,t+1);return i},numeric.random=function(t){return numeric._random(t,0)},numeric.norm2=function(t){return Math.sqrt(numeric.norm2Squared(t))},numeric.linspace=function(t,n,r){typeof r=="undefined"&&(r=Math.max(Math.round(n-t)+1,1));if(r<2)return r===1?[t]:[];var i,s=Array(r);r--;for(i=r;i>=0;i--)s[i]=(i*n+(r-i)*t)/r;return s},numeric.getBlock=function(t,n,r){function s(e,t){var o,u=n[t],a=r[t]-u,f=Array(a);if(t===i.length-1){for(o=a;o>=0;o--)f[o]=e[o+u];return f}for(o=a;o>=0;o--)f[o]=s(e[o+u],t+1);return f}var i=numeric.dim(t);return s(t,0)},numeric.setBlock=function(t,n,r,i){function o(e,t,i){var u,a=n[i],f=r[i]-a;if(i===s.length-1)for(u=f;u>=0;u--)e[u+a]=t[u];for(u=f;u>=0;u--)o(e[u+a],t[u],i+1)}var s=numeric.dim(t);return o(t,i,0),t},numeric.getRange=function(t,n,r){var i=n.length,s=r.length,o,u,a=Array(i),f,l;for(o=i-1;o!==-1;--o){a[o]=Array(s),f=a[o],l=t[n[o]];for(u=s-1;u!==-1;--u)f[u]=l[r[u]]}return a},numeric.blockMatrix=function(t){var n=numeric.dim(t);if(n.length<4)return numeric.blockMatrix([t]);var r=n[0],i=n[1],s,o,u,a,f;s=0,o=0;for(u=0;u=0;f--){a=Array(o),c=t[f];for(l=o-1;l>=3;--l)a[l]=c*n[l],--l,a[l]=c*n[l],--l,a[l]=c*n[l],--l,a[l]=c*n[l];while(l>=0)a[l]=c*n[l],--l;u[f]=a}return u},numeric.T=function(t,n){this.x=t,this.y=n},numeric.t=function(t,n){return new numeric.T(t,n)},numeric.Tbinop=function(t,n,r,i,s){var o=numeric.indexOf;if(typeof s!="string"){var u;s="";for(u in numeric)numeric.hasOwnProperty(u)&&(t.indexOf(u)>=0||n.indexOf(u)>=0||r.indexOf(u)>=0||i.indexOf(u)>=0)&&u.length>1&&(s+="var "+u+" = numeric."+u+";\n")}return Function(["y"],"var x = this;\nif(!(y instanceof numeric.T)) { y = new numeric.T(y); }\n"+s+"\n"+"if(x.y) {"+" if(y.y) {"+" return new numeric.T("+i+");\n"+" }\n"+" return new numeric.T("+r+");\n"+"}\n"+"if(y.y) {\n"+" return new numeric.T("+n+");\n"+"}\n"+"return new numeric.T("+t+");\n")},numeric.T.prototype.add=numeric.Tbinop("add(x.x,y.x)","add(x.x,y.x),y.y","add(x.x,y.x),x.y","add(x.x,y.x),add(x.y,y.y)"),numeric.T.prototype.sub=numeric.Tbinop("sub(x.x,y.x)","sub(x.x,y.x),neg(y.y)","sub(x.x,y.x),x.y","sub(x.x,y.x),sub(x.y,y.y)"),numeric.T.prototype.mul=numeric.Tbinop("mul(x.x,y.x)","mul(x.x,y.x),mul(x.x,y.y)","mul(x.x,y.x),mul(x.y,y.x)","sub(mul(x.x,y.x),mul(x.y,y.y)),add(mul(x.x,y.y),mul(x.y,y.x))"),numeric.T.prototype.reciprocal=function(){var t=numeric.mul,n=numeric.div;if(this.y){var r=numeric.add(t(this.x,this.x),t(this.y,this.y));return new numeric.T(n(this.x,r),n(numeric.neg(this.y),r))}return new T(n(1,this.x))},numeric.T.prototype.div=function div(e){e instanceof numeric.T||(e=new numeric.T(e));if(e.y)return this.mul(e.reciprocal());var div=numeric.div;return this.y?new numeric.T(div(this.x,e.x),div(this.y,e.x)):new numeric.T(div(this.x,e.x))},numeric.T.prototype.dot=numeric.Tbinop("dot(x.x,y.x)","dot(x.x,y.x),dot(x.x,y.y)","dot(x.x,y.x),dot(x.y,y.x)","sub(dot(x.x,y.x),dot(x.y,y.y)),add(dot(x.x,y.y),dot(x.y,y.x))"),numeric.T.prototype.transpose=function(){var t=numeric.transpose,n=this.x,r=this.y;return r?new numeric.T(t(n),t(r)):new numeric.T(t(n))},numeric.T.prototype.transjugate=function(){var t=numeric.transpose,n=this.x,r=this.y;return r?new numeric.T(t(n),numeric.negtranspose(r)):new numeric.T(t(n))},numeric.Tunop=function(t,n,r){return typeof r!="string"&&(r=""),Function("var x = this;\n"+r+"\n"+"if(x.y) {"+" "+n+";\n"+"}\n"+t+";\n")},numeric.T.prototype.exp=numeric.Tunop("return new numeric.T(ex)","return new numeric.T(mul(cos(x.y),ex),mul(sin(x.y),ex))","var ex = numeric.exp(x.x), cos = numeric.cos, sin = numeric.sin, mul = numeric.mul;"),numeric.T.prototype.conj=numeric.Tunop("return new numeric.T(x.x);","return new numeric.T(x.x,numeric.neg(x.y));"),numeric.T.prototype.neg=numeric.Tunop("return new numeric.T(neg(x.x));","return new numeric.T(neg(x.x),neg(x.y));","var neg = numeric.neg;"),numeric.T.prototype.sin=numeric.Tunop("return new numeric.T(numeric.sin(x.x))","return x.exp().sub(x.neg().exp()).div(new numeric.T(0,2));"),numeric.T.prototype.cos=numeric.Tunop("return new numeric.T(numeric.cos(x.x))","return x.exp().add(x.neg().exp()).div(2);"),numeric.T.prototype.abs=numeric.Tunop("return new numeric.T(numeric.abs(x.x));","return new numeric.T(numeric.sqrt(numeric.add(mul(x.x,x.x),mul(x.y,x.y))));","var mul = numeric.mul;"),numeric.T.prototype.log=numeric.Tunop("return new numeric.T(numeric.log(x.x));","var theta = new numeric.T(numeric.atan2(x.y,x.x)), r = x.abs();\nreturn new numeric.T(numeric.log(r.x),theta.x);"),numeric.T.prototype.norm2=numeric.Tunop("return numeric.norm2(x.x);","var f = numeric.norm2Squared;\nreturn Math.sqrt(f(x.x)+f(x.y));"),numeric.T.prototype.inv=function(){var t=this;if(typeof t.y=="undefined")return new numeric.T(numeric.inv(t.x));var n=t.x.length,r,i,s,o=numeric.identity(n),u=numeric.rep([n,n],0),a=numeric.clone(t.x),f=numeric.clone(t.y),l,c,h,p,d,v,m,g,r,i,s,y,b,w,E,S,x,T;for(r=0;ry&&(s=i,y=b);s!==r&&(T=a[r],a[r]=a[s],a[s]=T,T=f[r],f[r]=f[s],f[s]=T,T=o[r],o[r]=o[s],o[s]=T,T=u[r],u[r]=u[s],u[s]=T),l=a[r],c=f[r],d=o[r],v=u[r],w=l[r],E=c[r];for(i=r+1;i0;r--){d=o[r],v=u[r];for(i=r-1;i>=0;i--){m=o[i],g=u[i],w=a[i][r],E=f[i][r];for(s=n-1;s>=0;s--)S=d[s],x=v[s],m[s]-=w*S-E*x,g[s]-=w*x+E*S}}return new numeric.T(o,u)},numeric.T.prototype.get=function(t){var n=this.x,r=this.y,i=0,s,o=t.length;if(r){while(i=0?1:-1,i=r*numeric.norm2(t);n[0]+=i;var s=numeric.norm2(n);if(s===0)throw new Error("eig: internal error");return numeric.div(n,s)},numeric.toUpperHessenberg=function(t){var n=numeric.dim(t);if(n.length!==2||n[0]!==n[1])throw new Error("numeric: toUpperHessenberg() only works on square matrices");var r=n[0],i,s,o,u,a,f=numeric.clone(t),l,c,h,p,d=numeric.identity(r),v;for(s=0;s0){a=numeric.house(u),l=numeric.getBlock(f,[s+1,s],[r-1,r-1]),c=numeric.tensor(a,numeric.dot(a,l));for(i=s+1;i=4*c){var k,L;k=.5*(h+Math.sqrt(h*h-4*c)),L=.5*(h-Math.sqrt(h*h-4*c)),p=numeric.add(numeric.sub(numeric.dot(p,p),numeric.mul(p,k+L)),numeric.diag(numeric.rep([3],k*L)))}else p=numeric.add(numeric.sub(numeric.dot(p,p),numeric.mul(p,h)),numeric.diag(numeric.rep([3],c)));s=[p[0][0],p[1][0],p[2][0]],o=numeric.house(s),g=[e[0],e[1],e[2]],y=numeric.tensor(o,numeric.dot(o,g));for(w=0;w<3;w++){m=e[w],b=y[w];for(S=0;S=0?(w<0?x=-0.5*(w-A(S)):x=-0.5*(w+A(S)),k=(m-x)*(m-x)+g*g,L=y*y+(b-x)*(b-x),k>L?(k=A(k),N=(m-x)/k,C=g/k):(L=A(L),N=y/L,C=(b-x)/L),p=new s([[C,-N],[N,C]]),h.setRows(u,v,p.dot(h.getRows(u,v)))):(x=-0.5*w,T=.5*A(-S),k=(m-x)*(m-x)+g*g,L=y*y+(b-x)*(b-x),k>L?(k=A(k+T*T),N=(m-x)/k,C=g/k,x=0,T/=k):(L=A(L+T*T),N=y/L,C=(b-x)/L,x=T/L,T=0),p=new s([[C,-N],[N,C]],[[x, +T],[T,-x]]),h.setRows(u,v,p.dot(h.getRows(u,v))))}}var O=h.dot(t).dot(h.transjugate()),o=t.length,M=numeric.T.identity(o);for(v=0;v0)for(a=v-1;a>=0;a--){var _=O.get([a,a]),D=O.get([v,v]);if(!numeric.neq(_.x,D.x)&&!numeric.neq(_.y,D.y)){M.setRow(v,M.getRow(a));continue}x=O.getRow(a).getBlock([a],[v-1]),T=M.getRow(v).getBlock([a],[v-1]),M.set([v,a],O.get([a,v]).neg().sub(x.dot(T)).div(_.sub(D)))}for(v=0;v=u.length)u[u.length]=0;i[o]!==0&&u[o]++}}var r=u.length,a=Array(r+1);a[0]=0;for(s=0;s=d){s[f]=h[u];if(u===0)return;++f,--u,p=l[u],d=c[u]}else a=o[r[p]],i[a]===0?(i[a]=1,l[u]=p,++u,h[u]=a,p=n[a],c[u]=d=n[a+1]):++p},numeric.ccsLPSolve=function(t,n,r,i,s,o,u){var a=t[0],f=t[1],l=t[2],c=a.length-1,h=0,p=n[0],d=n[1],v=n[2],m,g,y,b,w,E,S,x,T,N,C,k;g=p[s],y=p[s+1],i.length=0;for(m=g;mb&&(w=m,b=E)}C(h[d])=v){s[l]=o[p[a]];if(a===0)return;++l,--a,d=c[a],v=h[a]}else f=r[d],i[f]===0?(i[f]=1,c[a]=d,++a,p[a]=f,f=o[f],d=n[f],h[a]=v=n[f+1]):++d}},numeric.ccsLPSolve0=function(t,n,r,i,s,o,u,a){var f=t[0],l=t[1],c=t[2],h=f.length-1,p=0,d=n[0],v=n[1],m=n[2],g,y,b,w,E,S,x,T,N,C,k,L;y=d[s],b=d[s+1],i.length=0;for(g=y;gb&&(w=m,b=E)}C(h[k[d]])t[n]&&(t[n]=e.length);var r;for(r in e)e.hasOwnProperty(r)&&dim(e[r],t,n+1);return t},numeric.sclone=function clone(e,t,n){typeof t=="undefined"&&(t=0),typeof n=="undefined"&&(n=numeric.sdim(e).length);var r,i=Array(e.length);if(t===n-1){for(r in e)e.hasOwnProperty(r)&&(i[r]=e[r]);return i}for(r in e)e.hasOwnProperty(r)&&(i[r]=clone(e[r],t+1,n));return i},numeric.sdiag=function(t){var n=t.length,r,i=Array(n),s,o,u;for(r=n-1;r>=1;r-=2)s=r-1,i[r]=[],i[r][r]=t[r],i[s]=[],i[s][s]=t[s];return r===0&&(i[0]=[],i[0][0]=t[r]),i},numeric.sidentity=function(t){return numeric.sdiag(numeric.rep([t],1))},numeric.stranspose=function(t){var n=[],r=t.length,i,s,o;for(i in t){if(!t.hasOwnProperty(i))continue;o=t[i];for(s in o){if(!o.hasOwnProperty(s))continue;typeof n[s]!="object"&&(n[s]=[]),n[s][i]=o[s]}}return n},numeric.sLUP=function(t,n){throw new Error("The function numeric.sLUP had a bug in it and has been removed. Please use the new numeric.ccsLUP function instead.")},numeric.sdotMM=function(t,n){var r=t.length,i=n.length,s=numeric.stranspose(n),o=s.length,u,a,f,l,c,h,p=Array(r),d;for(f=r-1;f>=0;f--){d=[],u=t[f];for(c=o-1;c>=0;c--){h=0,a=s[c];for(l in u){if(!u.hasOwnProperty(l))continue;l in a&&(h+=u[l]*a[l])}h&&(d[c]=h)}p[f]=d}return p},numeric.sdotMV=function(t,n){var r=t.length,i,s,o,u=Array(r),a;for(s=r-1;s>=0;s--){i=t[s],a=0;for(o in i){if(!i.hasOwnProperty(o))continue;n[o]&&(a+=i[o]*n[o])}a&&(u[s]=a)}return u},numeric.sdotVM=function(t,n){var r,i,s,o,u=[],a;for(r in t){if(!t.hasOwnProperty(r))continue;s=n[r],o=t[r];for(i in s){if(!s.hasOwnProperty(i))continue;u[i]||(u[i]=0),u[i]+=o*s[i]}}return u},numeric.sdotVV=function(t,n){var r,i=0;for(r in t)t[r]&&n[r]&&(i+=t[r]*n[r]);return i},numeric.sdot=function(t,n){var r=numeric.sdim(t).length,i=numeric.sdim(n).length,s=r*1e3+i;switch(s){case 0:return t*n;case 1001:return numeric.sdotVV(t,n);case 2001:return numeric.sdotMV(t,n);case 1002:return numeric.sdotVM(t,n);case 2002:return numeric.sdotMM(t,n);default:throw new Error("numeric.sdot not implemented for tensors of order "+r+" and "+i)}},numeric.sscatter=function(t){var n=t[0].length,r,i,s,o=t.length,u=[],a;for(i=n-1;i>=0;--i){if(!t[o-1][i])continue;a=u;for(s=0;s=0;--i)t[i]=[];for(i=r;i>=0;--i)t[i].push(n[i]);t[r+1].push(s)}}else gather(s,t,n)}return n.length>r&&n.pop(),t},numeric.cLU=function(t){var n=t[0],r=t[1],i=t[2],s=n.length,o=0,u,a,f,l,c,h;for(u=0;uo&&(o=n[u]);o++;var p=Array(o),d=Array(o),v=numeric.rep([o],Infinity),m=numeric.rep([o],-Infinity),g,y,b;for(f=0;fm[u]&&(m[u]=a);for(u=0;um[u+1]&&(m[u+1]=m[u]);for(u=o-1;u>=1;u--)v[u]=0;v--){while(l[g]>v)s[v]-=c[g]*s[l[g]],g--;s[v]/=c[g],g--}return s},numeric.cgrid=function(t,n){typeof t=="number"&&(t=[t,t]);var r=numeric.rep(t,-1),i,s,o;if(typeof n!="function")switch(n){case"L":n=function(e,n){return e>=t[0]/2||nf&&(f=i[u]);f++,r=numeric.rep([f],0);for(u=0;u1)o=u((i+s)/2),n[o]<=t?i=o:s=o;return this._at(t,i)}var r=t.length,c,h=Array(r);for(c=r-1;c!==-1;--c)h[c]=this.at(t[c]);return h},numeric.Spline.prototype.diff=function(){var t=this.x,n=this.yl,r=this.yr,i=this.kl,s=this.kr,o=n.length,u,a,f,l=i,c=s,h=Array(o),p=Array(o),d=numeric.add,v=numeric.mul,m=numeric.div,g=numeric.sub;for(u=o-1;u!==-1;--u)a=t[u+1]-t[u],f=g(r[u+1],n[u]),h[u]=m(d(v(f,6),v(i[u],-4*a),v(s[u+1],-2*a)),a*a),p[u+1]=m(d(v(f,-6),v(i[u],2*a),v(s[u+1],4*a)),a*a);return new numeric.Spline(t,l,c,h,p)},numeric.Spline.prototype.roots=function(){function t(e){return e*e}function n(e,t,n,r,i){var s=n*2-(t-e),o=-r*2+(t-e),u=(i+1)*.5,a=u*(1-u);return(1-u)*e+u*t+s*a*(1-u)+o*a*u}var r=[],i=this.x,s=this.yl,o=this.yr,u=this.kl,a=this.kr;typeof s[0]=="number"&&(s=[s],o=[o],u=[u],a=[a]);var f=s.length,l=i.length-1,c,h,p,d,v,m,g,y,b,w,r=Array(f),E,S,x,T,N,C,k,L,A,O,M,_,D,P,H,B,j,F=Math.sqrt;for(c=0;c!==f;++c){g=s[c],y=o[c],b=u[c],w=a[c],E=[];for(h=0;h!==l;h++){h>0&&y[h]*g[h]<0&&E.push(i[h]),A=i[h+1]-i[h],O=i[h],T=g[h],N=y[h+1],S=b[h]/A,x=w[h+1]/A,L=t(S-x+3*(T-N))+12*x*T,C=x+3*T+2*S-3*N,k=3*(x+S+2*(T-N)),L<=0?(_=C/k,_>i[h]&&_i[h]&&_i[h]&&D0){H=B,_=D;continue}var I=0;for(;;){j=(_*B-D*H)/(_-D);if(j<=H||j>=B)break;P=this._at(j,h);if(P*D>0)B=j,D=P,I===-1&&(_*=.5),I=-1;else{if(!(P*_>0))break;H=j,_=P,I===1&&(D*=.5),I=1}}E.push(j),H=M[p+1],_=this._at(H,h)}D===0&&E.push(B)}r[c]=E}return typeof this.yl[0]=="number"?r[0]:r},numeric.spline=function(t,n,r,i){var s=t.length,o=[],u=[],a=[],f,l=numeric.sub,c=numeric.mul,h=numeric.add;for(f=s-2;f>=0;f--)u[f]=t[f+1]-t[f],a[f]=l(n[f+1],n[f]);if(typeof r=="string"||typeof i=="string")r=i="periodic";var p=[[],[],[]];switch(typeof r){case"undefined":o[0]=c(3/(u[0]*u[0]),a[0]),p[0].push(0,0),p[1].push(0,1),p[2].push(2/u[0],1/u[0]);break;case"string":o[0]=h(c(3/(u[s-2]*u[s-2]),a[s-2]),c(3/(u[0]*u[0]),a[0])),p[0].push(0,0,0),p[1].push(s-2,0,1),p[2].push(1/u[s-2],2/u[s-2]+2/u[0],1/u[0]);break;default:o[0]=r,p[0].push(0),p[1].push(0),p[2].push(1)}for(f=1;f20)throw new Error("Numerical gradient fails");u[o]=n[o]+N,a=t(u),u[o]=n[o]-N,f=t(u),u[o]=n[o];if(isNaN(a)||isNaN(f)){N/=16;continue}l[o]=(a-f)/(2*N),y=n[o]-N,b=n[o],w=n[o]+N,S=(a-i)/N,x=(i-f)/N,T=s(m(l[o]),m(i),m(a),m(f),m(y),m(b),m(w),1e-8),p=g(s(m(S-l[o]),m(x-l[o]),m(S-x))/T,N/T);if(!(p>v))break;N/=16}}return l},numeric.uncmin=function(t,n,r,i,s,o,u){var a=numeric.gradient;typeof u=="undefined"&&(u={}),typeof r=="undefined"&&(r=1e-8),typeof i=="undefined"&&(i=function(e){return a(t,e)}),typeof s=="undefined"&&(s=1e3),n=numeric.clone(n);var f=n.length,l=t(n),c,h;if(isNaN(l))throw new Error("uncmin: f(x0) is a NaN!");var p=Math.max,d=numeric.norm2;r=p(r,numeric.epsilon);var v,m,g,y=u.Hinv||numeric.identity(f),b=numeric.dot,w=numeric.inv,E=numeric.sub,S=numeric.add,x=numeric.tensor,T=numeric.div,N=numeric.mul,C=numeric.all,k=numeric.isFinite,L=numeric.neg,A=0,O,M,_,D,P,H,B,j,F,I,q,R,U="";m=i(n);while(A=.1*F*h||isNaN(c)){F*=.5,++A;continue}break}if(F*I1)i=s(.5*(n+r)),a[i]<=t?n=i:r=i;return this._at(t,n)},numeric.dopri=function(t,n,r,i,s,o,u){typeof s=="undefined"&&(s=1e-6),typeof o=="undefined"&&(o=1e3);var a=[t],f=[r],l=[i(t,r)],c,h,p,d,v,m,g=[],y=.2,b=[.075,.225],w=[44/45,-56/15,32/9],E=[19372/6561,-25360/2187,64448/6561,-212/729],S=[9017/3168,-355/33,46732/5247,49/176,-5103/18656],x=[35/384,0,500/1113,125/192,-2187/6784,11/84],T=[.10013431883002395,0,.3918321794184259,-0.02982460176594817,.05893268337240795,-0.04497888809104361,.023904308236133973],N=[.2,.3,.8,8/9,1,1],C=[-71/57600,0,71/16695,-71/1920,17253/339200,-22/525,.025],k=0,L,A,O=(n-t)/10,M=0,_=numeric.add,D=numeric.mul,P,H,B=Math.max,j=Math.min,F=Math.abs,I=numeric.norminf,q=Math.pow,R=numeric.any,U=numeric.lt,z=numeric.and,W=numeric.sub,X,V,$,J=new numeric.Dopri(a,f,l,g,-1,"");typeof u=="function"&&(X=u(t,r));while(tn&&(O=n-t),c=i(t+N[0]*O,_(r,D(y*O,l[k]))),h=i(t+N[1]*O,_(_(r,D(b[0]*O,l[k])),D(b[1]*O,c))),p=i(t+N[2]*O,_(_(_(r,D(w[0]*O,l[k])),D(w[1]*O,c)),D(w[2]*O,h))),d=i(t+N[3]*O,_(_(_(_(r,D(E[0]*O,l[k])),D(E[1]*O,c)),D(E[2]*O,h)),D(E[3]*O,p))),v=i(t+N[4]*O,_(_(_(_(_(r,D(S[0]*O,l[k])),D(S[1]*O,c)),D(S[2]*O,h)),D(S[3]*O,p)),D(S[4]*O,d))),P=_(_(_(_(_(r,D(l[k],O*x[0])),D(h,O*x[2])),D(p,O*x[3])),D(d,O*x[4])),D(v,O*x[5])),m=i(t+O,P),L=_(_(_(_(_(D(l[k],O*C[0]),D(h,O*C[2])),D(p,O*C[3])),D(d,O*C[4])),D(v,O*C[5])),D(m,O*C[6])),typeof L=="number"?H=F(L):H=I(L);if(H>s){O=.2*O*q(s/H,.25);if(t+O===t){J.msg="Step size became too small";break}continue}g[k]=_(_(_(_(_(_(r,D(l[k],O*T[0])),D(h,O*T[2])),D(p,O*T[3])),D(d,O*T[4])),D(v,O*T[5])),D(m,O*T[6])),++k,a[k]=t+O,f[k]=P,l[k]=m;if(typeof u=="function"){var K,Q=t,G=t+.5*O,Y;V=u(G,g[k-1]),$=z(U(X,0),U(0,V)),R($)||(Q=G,G=t+O,X=V,V=u(G,P),$=z(U(X,0),U(0,V)));if(R($)){var Z,et,tt,nt,rt=0,it=1,st=1;for(;;){if(typeof X=="number")Y=(st*V*Q-it*X*G)/(st*V-it*X);else{Y=G;for(A=X.length-1;A!==-1;--A)X[A]<0&&V[A]>0&&(Y=j(Y,(st*V[A]*Q-it*X[A]*G)/(st*V[A]-it*X[A])))}if(Y<=Q||Y>=G)break;K=J._at(Y,k-1),nt=u(Y,K),tt=z(U(X,0),U(0,nt)),R(tt)?(G=Y,V=nt,$=tt,st=1,rt===-1?it*=.5:it=1,rt=-1):(Q=Y,X=nt,it=1,rt===1?st*=.5:st=1,rt=1)}return P=J._at(.5*(t+Y),k-1),J.f[k]=i(Y,K),J.x[k]=Y,J.y[k]=K,J.ymid[k-1]=P,J.events=$,J.iterations=M,J}}t+=O,r=P,X=V,O=j(.8*O*q(s/H,.25),4*O)}return J.iterations=M,J},numeric.LU=function(e,t){t=t||!1;var n=Math.abs,r,i,s,o,u,a,f,l,c,h=e.length,p=h-1,d=new Array(h);t||(e=numeric.clone(e));for(s=0;s=0;--r){l=s[r];for(i=r+1;iK)E=K;W=d(t,l(E,B)),I=h(J,j);for(X=v-1;X!==-1;--X)I[X][X]+=1;$=q(I,p(W,E),!0);var Q=p(R,h(n,$)),G=1;for(X=m-1;X!==-1;--X)Q[X]<0&&(G=D(G,-0.999*Q[X]));g=c(o,l($,G)),R=c(r,h(n,g));if(!P(H(R,0)))return{solution:o,message:"",iterations:U};o=g;if(E=0?y=!1:y=!0;if(y)return{solution:g,message:"Unbounded",iterations:U}}return{solution:o,message:"maximum iteration count exceeded",iterations:U}},numeric._solveLP=function(t,n,r,i,s){var o=t.length,u=r.length,a,f=numeric.sum,l=numeric.log,c=numeric.mul,h=numeric.sub,p=numeric.dot,d=numeric.div,v=numeric.add,m=numeric.rep([o],0).concat([1]),g=numeric.rep([u,1],-1),y=numeric.blockMatrix([[n,g]]),b=r,a=numeric.rep([o],0).concat(Math.max(0,numeric.sup(numeric.neg(r)))+1),w=numeric.__solveLP(m,y,b,i,s,a,!1),E=numeric.clone(w.solution);E.length=o;var S=numeric.inf(h(r,p(n,E)));if(S<0)return{solution:NaN,message:"Infeasible",iterations:w.iterations};var x=numeric.__solveLP(t,n,r,i,s-w.iterations,E,!0);return x.iterations+=w.iterations,x},numeric.solveLP=function(t,n,r,i,s,o,u){typeof u=="undefined"&&(u=1e3),typeof o=="undefined"&&(o=numeric.epsilon);if(typeof i=="undefined")return numeric._solveLP(t,n,r,o,u);var a=i.length,f=i[0].length,l=n.length,c=numeric.echelonize(i),h=numeric.rep([f],0),p=c.P,d=[],v;for(v=p.length-1;v!==-1;--v)h[p[v]]=1;for(v=f-1;v!==-1;--v)h[v]===0&&d.push(v);var m=numeric.getRange,g=numeric.linspace(0,a-1),y=numeric.linspace(0,l-1),b=m(i,g,d),w=m(n,y,p),E=m(n,y,d),S=numeric.dot,x=numeric.sub,T=S(w,c.I),N=x(E,S(T,b)),C=x(r,S(T,s)),k=Array(p.length),L=Array(d.length);for(v=p.length-1;v!==-1;--v)k[v]=t[p[v]];for(v=d.length-1;v!==-1;--v)L[v]=t[d[v]];var A=x(L,S(k,S(c.I,b))),O=numeric._solveLP(A,N,C,o,u),M=O.solution;if(M!==M)return O;var _=S(c.I,x(s,S(b,M))),D=Array(t.length);for(v=p.length-1;v!==-1;--v)D[p[v]]=_[v];for(v=d.length-1;v!==-1;--v)D[d[v]]=M[v];return{solution:D,message:O.message,iterations:O.iterations}},numeric.MPStoLP=function(t){function y(e){throw new Error("MPStoLP: "+e+"\nLine "+s+": "+t[s]+"\nCurrent state: "+r[n]+"\n")}t instanceof String&&t.split("\n");var n=0,r=["Initial state","NAME","ROWS","COLUMNS","RHS","BOUNDS","ENDATA"],i=t.length,s,o,u,a=0,f={},l=[],c=0,h={},p=0,d,v=[],m=[],g=[];for(s=0;s=s)t/=2,u/=2,a>>>=1;return(t+a)/u},c},o=t.pow(n,r),i=t.pow(2,i),s=i*2,f(t.random(),e)}([],numeric.seedrandom,256,6,52),function(e){function t(e){if(typeof e!="object")return e;var n=[],r,i=e.length;for(r=0;rp)g[E]=P;else{g[E]=-Math.abs(P);if(P>0){for(w=1;w<=o;w+=1)f[w][b]=-f[w][b];l[b]=-l[b]}}}for(b=1;b<=v;b+=1)g[L+d[b]]=0;O=0,D=0;for(b=1;b<=h;b+=1)g[L+b]=1;b-=1){P=g[b],E=k+b*(b+3)/2,S=E-b;for(w=b+1;w<=v;w+=1)P-=g[E]*g[C+w],E+=w;P/=g[S],g[C+b]=P;if(d[b]p)g[L+O]=P;else{g[L+O]=-Math.abs(P);if(P>0){for(w=1;w<=o;w+=1)f[w][O]=-f[w][O];l[O]=-l[O]}}return 700}v+=1,d[v]=O,E=k+(v-1)*v/2+1;for(b=1;b<=v-1;b+=1)g[E]=g[b],E+=1;if(v===o)g[E]=g[o];else{for(b=o;b>=v+1;b-=1){if(g[b]===0)break;j=Math.max(Math.abs(g[b-1]),Math.abs(g[b])),F=Math.min(Math.abs(g[b-1]),Math.abs(g[b])),g[b-1]>=0?D=Math.abs(j*Math.sqrt(1+F*F/(j*j))):D=-Math.abs(j*Math.sqrt(1+F*F/(j*j))),j=g[b-1]/D,F=g[b]/D;if(j===1)break;if(j===0){g[b-1]=F*D;for(w=1;w<=o;w+=1)D=e[w][b-1],e[w][b-1]=e[w][b],e[w][b]=D}else{g[b-1]=D,I=F/(1+j);for(w=1;w<=o;w+=1)D=j*e[w][b-1]+F*e[w][b],e[w][b]=I*(e[w][b-1]+D)-e[w][b],e[w][b-1]=D}}g[E]=g[v]}return 0}function J(){E=k+T*(T+1)/2+1,S=E+T;if(g[S]===0)return 798;j=Math.max(Math.abs(g[S-1]),Math.abs(g[S])),F=Math.min(Math.abs(g[S-1]),Math.abs(g[S])),g[S-1]>=0?D=Math.abs(j*Math.sqrt(1+F*F/(j*j))):D=-Math.abs(j*Math.sqrt(1+F*F/(j*j))),j=g[S-1]/D,F=g[S]/D;if(j===1)return 798;if(j===0){for(b=T+1;b<=v;b+=1)D=g[S-1],g[S-1]=g[S],g[S]=D,S+=b;for(b=1;b<=o;b+=1)D=e[b][T],e[b][T]=e[b][T+1],e[b][T+1]=D}else{I=F/(1+j);for(b=T+1;b<=v;b+=1)D=j*g[S-1]+F*g[S],g[S]=I*(g[S-1]+D)-g[S],g[S-1]=D,S+=b;for(b=1;b<=o;b+=1)D=j*e[b][T]+F*e[b][T+1],e[b][T+1]=I*(e[b][T]+D)-e[b][T+1],e[b][T]=D}return 0}function K(){S=E-T;for(b=1;b<=T;b+=1)g[S]=g[E],E+=1,S+=1;return g[A+T]=g[A+T+1],d[T]=d[T+1],T+=1,Tt?e*Math.sqrt(1+t*t/e/e):t==0?e:t*Math.sqrt(1+e*e/t/t)}var n,r=numeric.epsilon,i=1e-64/r,s=50,o=0,u=0,a=0,f=0,l=0,c=numeric.clone(t),h=c.length,p=c[0].length;if(h=0&&(b=-b),w=y*b-T,c[u][u]=y-b;for(a=l;a=0&&(b=-b),w=y*b-T,c[u][u+1]=y-b;for(a=l;aE&&(E=S)}for(u=p-1;u!=-1;u+=-1){if(b!=0){w=b*c[u][u+1];for(a=l;a=s-1)throw"Error: no convergence.";E=v[l],S=v[f-1],b=d[f-1],w=d[f],y=((S-x)*(S+x)+(b-w)*(b+w))/(2*w*S),b=g(y,1),y<0?y=((E-x)*(E+x)+w*(S/(y-b)-w))/E:y=((E-x)*(E+x)+w*(S/(y+b)-w))/E,o=1,T=1;for(u=l+1;u=0;a--)if(v[a]