Template:Team:Amsterdam/scripts/numeric.js

From 2012.igem.org

var numeric = (typeof exports === "undefined")?(function numeric() {}):(exports); if(typeof global !== "undefined") { global.numeric = numeric; }

numeric.version = "1.0.2";

// 1. Utility functions numeric.bench = function bench (f,interval) {

   var t1,t2,n,i;
   if(typeof interval === "undefined") { interval = 15; }
   n = 1;
   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);

}

/*

* IE doesn't have a global eval that works?
* This apparently returns nothing: window.execScript('function () {}').            
* Other things that don't work: window.execScript('foo = eval("function () {}");')
*/

if(1 || typeof window === "undefined" || typeof window.execScript !== "undefined") numeric.Function = Function; else numeric.Function = function() {

   var foo = '(function (';
   for(var k=0;k<arguments.length-1;++k) { if(k>0) foo+=','; foo += arguments[k]; }
   foo += ') { \n'+arguments[k]+'\n});';
   return window.eval(foo);

}

numeric.precision = 4; numeric.largeArray = 50;

numeric.prettyPrint = function prettyPrint(x) {

   function fmtnum(x) {
       if(x === 0) { return '0'; }
       if(isNaN(x)) { return 'NaN'; }
       if(x<0) { return '-'+fmtnum(-x); }
       if(isFinite(x)) {
           var scale = Math.floor(Math.log(x) / Math.log(10));
           var normalized = x / Math.pow(10,scale);
           var basic = normalized.toPrecision(numeric.precision);
           if(parseFloat(basic) === 10) { scale++; normalized = 1; basic = normalized.toPrecision(numeric.precision); }
           return parseFloat(basic).toString()+'e'+scale.toString();
       }
       return 'Infinity';
   }
   var ret = [];
   function foo(x) {
       var k;
       if(typeof x === "undefined") { ret.push(Array(numeric.precision+8).join(' ')); return false; }
       if(typeof x === "string") { ret.push('"'+x+'"'); return false; }
       if(typeof x === "boolean") { ret.push(x.toString()); return false; }
       if(typeof x === "number") {
           var a = fmtnum(x);
           var b = x.toPrecision(numeric.precision);
           var c = parseFloat(x.toString()).toString();
           var d = [a,b,c,parseFloat(b).toString(),parseFloat(c).toString()];
           for(k=1;k<d.length;k++) { if(d[k].length < a.length) a = d[k]; }
           ret.push(Array(numeric.precision+8-a.length).join(' ')+a);
           return false;
       }
       if(x === null) { ret.push("null"); return false; }
       if(typeof x === "function") { 
           ret.push(x.toString());
           var flag = false;
           for(k in x) { if(x.hasOwnProperty(k)) { 
               if(flag) ret.push(',\n');
               else ret.push('\n{');
               flag = true; 
               ret.push(k); 
               ret.push(': \n'); 
               foo(x[k]); 
           } }
           if(flag) ret.push('}\n');
           return true;
       }
       if(x instanceof Array) {
           if(x.length > numeric.largeArray) { ret.push('...Large Array...'); return true; }
           var flag = false;
           ret.push('[');
           for(k=0;k<x.length;k++) { if(k>0) { 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;k<d.length;k++) { ret[k] = foo(d[k]); }
       return ret;
   }
   return foo(d);

}

numeric.parseFloat = function parseFloat_(d) {

   function foo(d) {
       if(typeof d === 'string') { return parseFloat(d); }
       if(!(d instanceof Array)) { throw new Error("parseFloat: parameter must be arrays of strings"); }
       var ret = [],k;
       for(k=0;k<d.length;k++) { ret[k] = foo(d[k]); }
       return ret;
   }
   return foo(d);

}

numeric.parseCSV = function parseCSV(t) {

   var foo = t.split('\n');
   var j,k;
   var ret = [];
   var pat = /(([^'",]*)|('[^']*')|("[^"]*")),/g;
   var patnum = /^\s*(([+-]?[0-9]+(\.[0-9]*)?(e[+-]?[0-9]+)?)|([+-]?[0-9]*(\.[0-9]+)?(e[+-]?[0-9]+)?))\s*$/;
   var stripper = function(n) { return n.substr(0,n.length-1); }
   var count = 0;
   for(k=0;k<foo.length;k++) {
     var bar = (foo[k]+",").match(pat),baz;
     if(bar.length>0) {
         ret[count] = [];
         for(j=0;j<bar.length;j++) {
             baz = stripper(bar[j]);
             if(patnum.test(baz)) { ret[count][j] = parseFloat(baz); }
             else ret[count][j] = baz;
         }
         count++;
     }
   }
   return ret;

}

numeric.toCSV = function toCSV(A) {

   var s = numeric.dim(A);
   var i,j,m,n,row,ret;
   m = s[0];
   n = s[1];
   ret = [];
   for(i=0;i<m;i++) {
       row = [];
       for(j=0;j<m;j++) { row[j] = A[i][j].toString(); }
       ret[i] = row.join(', ');
   }
   return ret.join('\n')+'\n';

}

numeric.getURL = function getURL(url) {

   var client = new XMLHttpRequest();
   client.open("GET",url,false);
   client.send();
   return client;

}

numeric.imageURL = function imageURL(img) {

   function base64(A) {
       var n = A.length, i,x,y,z,p,q,r,s;
       var key = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=";
       var ret = "";
       for(i=0;i<n;i+=3) {
           x = A[i];
           y = A[i+1];
           z = A[i+2];
           p = x >> 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<h;i++) {
       if(i<h-1) { stream.push(0); }
       else { stream.push(1); }
       a = (3*w+1+(i===0))&255; b = ((3*w+1+(i===0))>>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;j<w;j++) {
           for(k=0;k<3;k++) {
               a = img[k][i][j];
               if(a>255) 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 numeric.Function('x','accum','_s','_k',
           'if(typeof accum === "undefined") accum = '+init+';\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.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<n;i++) {
       if(x[i] === y[i]) { continue; }
       if(typeof x[i] === "object") { if(!same(x[i],y[i])) return false; }
       else { return false; }
   }
   return true;

}

numeric.rep = function rep(s,v,k) {

   if(typeof k === "undefined") { k=0; }
   var n = s[k], ret = Array(n), i;
   if(k === s.length-1) {
       for(i=n-2;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.dotMMbig = function dotMMbig(x,y) {

   var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0,s1,s2,s3,baz,accum;
   var dotVV = numeric.dotVV,min = Math.min;
   p = x.length; q = y.length; r = y[0].length;
   ret = Array(p);
   woo = numeric.transpose(y);
   for(i0=0;i0<p;i0+=4) {
       p0 = min(i0+4,p);
       for(i=i0;i<p0;i++) { ret[i] = Array(r); }
       for(k0=0;k0<r;k0+=4) {
           r0 = min(k0+4,r);
           for(i=i0;i<p0;i++) {
               bar = x[i];
               foo = ret[i];
               for(k=k0;k<r0;k++) {
                   foo[k] = dotVV(bar,woo[k]);
               }
           }
       }
   }
   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.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 < 40) 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 = ;
   for(k=0;k<params.length;k++) {
       if(avec.test(params[k])) {
           p = params[k].substring(0,params[k].length-3);
           thevec = p;
       } else { p = params[k]; }
       fun.push(p);
   }
   fun[params.length] = '_s';
   fun[params.length+1] = '_k';
   fun[params.length+2] = (
           'if(typeof _s === "undefined") _s = numeric.dim('+thevec+');\n'+
           'if(typeof _k === "undefined") _k = 0;\n'+
           'var _n = _s[_k];\n'+
           'var i, ret = Array(_n);\n'+
           'if(_k < _s.length-1) {\n'+
           '    for(i=_n-1;i>=0;i--) ret[i] = arguments.callee('+params.join(',')+',_s,_k+1);\n'+
           '    return ret;\n'+
           '}\n'+
           setup+'\n'+
           'for(i=_n-1;i>=3;--i) { \n'+
           '    '+body+'\n'+
           '    --i;\n'+
           '    '+body+'\n'+
           '    --i;\n'+
           '    '+body+'\n'+
           '    --i;\n'+
           '    '+body+'\n'+
           '}\n'+
           'while(i>=0) {\n'+
           '    '+body+'\n'+
           '    --i;\n'+
           '}\n'+
           'return ret;'
           );
   return numeric.Function.apply(null,fun);

}

numeric._biforeach = (function _biforeach(x,y,s,k,f) {

   if(k === s.length-1) { f(x,y); return; }
   var i,n=s[k];
   for(i=n-1;i>=0;i--) { _biforeach(x[i],y[i],s,k+1,f); }

});

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: '>>>=',
       andeq: '&=',
       oreq: '|=',
       xoreq: '^='

}; numeric.mathfuns = ['abs','acos','asin','atan','ceil','cos',

                   'exp','floor','log','round','sin','sqrt','tan'];

numeric.ops1 = {

       neg: '-',
       not: '!',
       bnot: '~'

};

(function () {

   var i,o;
   for(i in numeric.ops2) {
       if(numeric.ops2.hasOwnProperty(i)) {
           o = numeric.ops2[i];
           numeric[i+'VV'] = numeric.pointwise(['x[i]','y[i]'],'ret[i] = x[i] '+o+' y[i];');
           numeric[i+'SV'] = numeric.pointwise(['x','y[i]'],'ret[i] = x '+o+' y[i];');
           numeric[i+'VS'] = numeric.pointwise(['x[i]','y'],'ret[i] = x[i] '+o+' y;');
           numeric[i] = numeric.Function(
                   'var n = arguments.length, i, x = arguments[0], y;\n'+
                   'var VV = numeric.'+i+'VV, VS = numeric.'+i+'VS, SV = numeric.'+i+'SV;\n'+
                   'for(i=1;i!==n;++i) { \n'+
                   '  y = arguments[i];'+
                   '  if(typeof x === "object") {\n'+
                   '      if(typeof y === "object") x = VV(x,y);\n'+
                   '      else x = VS(x,y);\n'+
                   '  } else if(typeof y === "object") x = SV(x,y);\n'+
                   '  else x = x '+o+' y;\n'+
                   '}\nreturn x;\n');
           numeric[o] = numeric[i];
       }
   }
   for(i in numeric.ops1) {
       if(numeric.ops1.hasOwnProperty(i)) {
           o = numeric.ops1[i];
           numeric[i+'V'] = numeric.pointwise(['x[i]'],'ret[i] = '+o+'x[i];');
           numeric[i] = numeric.Function('x','if(typeof x === "object") return numeric.'+i+'V(x);\nreturn '+o+'(x);');
       }
   }
   for(i=0;i<numeric.mathfuns.length;i++) {
       o = numeric.mathfuns[i];
       numeric[o+'V'] = numeric.pointwise(['x[i]'],'ret[i] = fun(x[i]);','var fun = Math.'+o+';');
       numeric[o] = numeric.Function('x','if(typeof x === "object") return numeric.'+o+'V(x);\nreturn Math.'+o+'(x);');
   }
   numeric.isNaNV = numeric.pointwise(['x[i]'],'ret[i] = isNaN(x[i]);');
   numeric.isNaN = function isNaN(x) { if(typeof x === "object") return numeric.isNaNV(x); return isNaN(x); }
   numeric.isFiniteV = numeric.pointwise(['x[i]'],'ret[i] = isFinite(x[i]);');
   numeric.isFinite = function isNaN(x) { if(typeof x === "object") return numeric.isFiniteV(x); return isFinite(x); }
   for(i in numeric.opseq) {
       if(numeric.opseq.hasOwnProperty(i)) {
           numeric[i+'S'] = numeric.Function('x','y',
                   'var n = x.length, i;\n'+
                   'for(i=n-1;i>=0;i--) x[i] '+numeric.opseq[i]+' y;');
           numeric[i+'V'] = numeric.Function('x','y',
                   'var n = x.length, i;\n'+
                   'for(i=n-1;i>=0;i--) x[i] '+numeric.opseq[i]+' y[i];');
           numeric[i] = numeric.Function('x','y',
                   'var s = numeric.dim(x);\n'+
                   'if(typeof y === "number") { numeric._biforeach(x,y,s,0,numeric.'+i+'S); return x; }\n'+
                   'numeric._biforeach(x,y,s,0,numeric.'+i+'V);\n'+
                   'return x;');
           numeric[numeric.opseq[i]] = numeric[i];
       }
   }

}());

numeric.atan2VV = numeric.pointwise(['x[i]','y[i]'],'ret[i] = atan2(x[i],y[i]);','var atan2 = Math.atan2;'); numeric.atan2VS = numeric.pointwise(['x[i]','y'],'ret[i] = atan2(x[i],y);','var atan2 = Math.atan2;'); numeric.atan2SV = numeric.pointwise(['x','y[i]'],'ret[i] = atan2(x,y[i]);','var atan2 = Math.atan2;'); numeric.atan2 = function atan2(x,y) {

   if(typeof x === "object") {
       if(typeof y === "object") return numeric.atan2VV(x,y);
       return numeric.atan2VS(x,y);
   }
   if (typeof y === "object") return numeric.atan2SV(x,y);
   return Math.atan2(x,y);

}

numeric.powVV = numeric.pointwise(['x[i]','y[i]'],'ret[i] = pow(x[i],y[i]);','var pow = Math.pow;'); numeric.powVS = numeric.pointwise(['x[i]','y'],'ret[i] = pow(x[i],y);','var pow = Math.pow;'); numeric.powSV = numeric.pointwise(['x','y[i]'],'ret[i] = pow(x,y[i]);','var pow = Math.pow;'); numeric.pow = function pow(x,y) {

   if(typeof x === "object") {
       if(typeof y === "object") return numeric.powVV(x,y);
       return numeric.powVS(x,y);
   }
   if (typeof y === "object") return numeric.powSV(x,y);
   return Math.pow(x,y);

}

numeric.clone = numeric.pointwise(['x[i]'],'ret[i] = x[i];');

numeric.inv = function inv(x) {

   var s = numeric.dim(x), abs = Math.abs;
   if(s.length !== 2 || s[0] !== s[1]) { throw new Error('numeric: inv() only works on square matrices'); }
   var n = s[0], ret = numeric.identity(n),i,j,k,A = numeric.clone(x),Aj,Ai,Ij,Ii,alpha,temp,k0,k1,k2,k3;
   var P = numeric.linspace(0,n-1), Q = numeric.rep([n],0);
   for(j=0;j<n-1;j++) {
       k=j;
       for(i=j+1;i<n;i++) { if(abs(A[i][j]) > abs(A[k][j])) { k = i; } }
       if(k!==j) {
           temp = A[k]; A[k] = A[j]; A[j] = temp;
           temp = ret[k]; ret[k] = ret[j]; ret[j] = temp;
           temp = P[k]; P[k] = P[j]; P[j] = temp;
       }
       Aj = A[j];
       Ij = ret[j];
       for(i=j+1;i<n;i++) {
           Ai = A[i];
           Ii = ret[i];
           alpha = Ai[j]/Aj[j];
           if(alpha === 0) continue;
           for(k=j+1;k<n-1;k+=2) {
               k1 = k+1;
               Ai[k] -= Aj[k]*alpha;
               Ai[k1] -= Aj[k1]*alpha;
           }
           if(k<n) { Ai[k] -= Aj[k]*alpha; }
           for(k=j;k>=1;k-=2) {
               k2 = P[k-1]; 
               k3 = P[k];
               Ii[k2] -= Ij[k2]*alpha;
               Ii[k3] -= Ij[k3]*alpha;
           }
           if(k>=0) { Ii[P[k]] -= Ij[P[k]]*alpha; }
       }
   }
   for(j=n-1;j>0;j--) {
       Aj = A[j];
       Ij = ret[j];
       for(i=0;i<j;i++) {
           Ii = ret[i];
           alpha = A[i][j]/Aj[j];
           if(alpha === 0) continue;
           for(k=0;k<n-1;k+=2) {
               k1 = k+1;
               Ii[k] -= Ij[k]*alpha;
               Ii[k1] -= Ij[k1]*alpha;
           }
           if(k!==n) { Ii[k] -= Ij[k]*alpha; }
       }
   }
   for(i=0;i<n;i++) {
       alpha = A[i][i];
       if(alpha === 1) continue;
       Ii = ret[i];
       for(j=0;j<n;j++) { Ii[j] /= alpha; }
   }
   return ret;

}

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<n-1;j++) {
       k=j;
       for(i=j+1;i<n;i++) { if(Math.abs(A[i][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<n;i++) {
           Ai = A[i];
           alpha = Ai[j]/Aj[j];
           for(k=j+1;k<n-1;k+=2) {
               k1 = k+1;
               Ai[k] -= Aj[k]*alpha;
               Ai[k1] -= Aj[k1]*alpha;
           }
           if(k!==n) { Ai[k] -= Aj[k]*alpha; }
       }
       if(Aj[j] === 0) { return 0; }
       ret *= Aj[j];
   }
   return ret*A[j][j];

}

numeric.transpose = function transpose(x) {

   var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj;
   for(j=0;j<n;j++) ret[j] = Array(m);
   for(i=m-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<n;j++) ret[j] = Array(m);
   for(i=m-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._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.norm2Squared = function norm2Squared(x) {} numeric.norm2Squared = numeric.mapreduce('accum += xi*xi;','0');

numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); } numeric.norminf = numeric.mapreduce('accum = max(abs(xi),accum);','0; var max = Math.max, abs = Math.abs;'); numeric.sum = numeric.mapreduce('accum += xi;','0');

numeric.linspace = function linspace(a,b,n) {

   if(typeof n === "undefined") n = Math.round(b-a)+1;
   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.tensor = function tensor(x,y) {

   if(typeof x === "number" || typeof y === "number") return numeric.mul(x,y);
   var s1 = numeric.dim(x), s2 = numeric.dim(y);
   if(s1.length !== 1 || s2.length !== 1) {
       throw new Error('numeric: tensor product is only defined for vectors');
   }
   var m = s1[0], n = s2[0], A = Array(m), Ai, i,j,xi;
   for(i=m-1;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 numeric.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 numeric.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<n;i++) {
       ax = Ax[i][i]; ay = Ay[i][i];
       d = ax*ax+ay*ay;
       k = i;
       for(j=i+1;j<n;j++) {
           ax = Ax[j][i]; ay = Ay[j][i];
           d1 = ax*ax+ay*ay;
           if(d1 > 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;j<n;j++) {
           bx = Aix[j]; by = Aiy[j];
           Aix[j] = (bx*ax+by*ay)/d;
           Aiy[j] = (by*ax-bx*ay)/d;
       }
       for(j=0;j<n;j++) {
           bx = Rix[j]; by = Riy[j];
           Rix[j] = (bx*ax+by*ay)/d;
           Riy[j] = (by*ax-bx*ay)/d;
       }
       for(j=i+1;j<n;j++) {
           Ajx = Ax[j]; Ajy = Ay[j];
           Rjx = Rx[j]; Rjy = Ry[j];
           ax = Ajx[i]; ay = Ajy[i];
           for(k=i+1;k<n;k++) {
               bx = Aix[k]; by = Aiy[k];
               Ajx[k] -= bx*ax-by*ay;
               Ajy[k] -= by*ax+bx*ay;
           }
           for(k=0;k<n;k++) {
               bx = Rix[k]; by = Riy[k];
               Rjx[k] -= bx*ax-by*ay;
               Rjy[k] -= by*ax+bx*ay;
           }
       }
   }
   for(i=n-1;i>0;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<n) {
           ik = i[k];
           x = x[ik];
           y = y[ik];
           k++;
       }
       return new numeric.T(x,y);
   }
   while(k<n) {
       ik = i[k];
       x = x[ik];
       k++;
   }
   return new numeric.T(x);

} numeric.T.prototype.set = function set(i,v) {

   var x = this.x, y = this.y, k = 0, ik, n = i.length, vx = v.x, vy = v.y;
   if(n===0) {
       if(vy) { this.y = vy; }
       else if(y) { this.y = undefined; }
       this.x = x;
       return this;
   }
   if(vy) {
       if(y) { /* ok */ }
       else {
           y = numeric.rep(numeric.dim(x),0);
           this.y = y;
       }
       while(k<n-1) {
           ik = i[k];
           x = x[ik];
           y = y[ik];
           k++;
       }
       ik = i[k];
       x[ik] = vx;
       y[ik] = vy;
       return this;
   }
   if(y) {
       while(k<n-1) {
           ik = i[k];
           x = x[ik];
           y = y[ik];
           k++;
       }
       ik = i[k];
       x[ik] = vx;
       if(vx instanceof Array) y[ik] = numeric.rep(numeric.dim(vx),0);
       else y[ik] = 0;
       return this;
   }
   while(k<n-1) {
       ik = i[k];
       x = x[ik];
       k++;
   }
   ik = i[k];
   x[ik] = vx;
   return this;

} numeric.T.prototype.getRows = function getRows(i0,i1) {

   var n = i1-i0+1, j;
   var rx = Array(n), ry, x = this.x, y = this.y;
   for(j=i0;j<=i1;j++) { rx[j-i0] = x[j]; }
   if(y) {
       ry = Array(n);
       for(j=i0;j<=i1;j++) { ry[j-i0] = y[j]; }
       return new numeric.T(rx,ry);
   }
   return new numeric.T(rx);

} numeric.T.prototype.setRows = function setRows(i0,i1,A) {

   var j;
   var rx = this.x, ry = this.y, x = A.x, y = A.y;
   for(j=i0;j<=i1;j++) { rx[j] = x[j-i0]; }
   if(y) {
       if(!ry) { ry = numeric.rep(numeric.dim(rx),0); this.y = ry; }
       for(j=i0;j<=i1;j++) { ry[j] = y[j-i0]; }
   } else if(ry) {
       for(j=i0;j<=i1;j++) { ry[j] = numeric.rep([x[j-i0].length],0); }
   }
   return this;

} numeric.T.prototype.getRow = function getRow(k) {

   var x = this.x, y = this.y;
   if(y) { return new numeric.T(x[k],y[k]); }
   return new numeric.T(x[k]);

} numeric.T.prototype.setRow = function setRow(i,v) {

   var rx = this.x, ry = this.y, x = v.x, y = v.y;
   rx[i] = x;
   if(y) {
       if(!ry) { ry = numeric.rep(numeric.dim(rx),0); this.y = ry; }
       ry[i] = y;
   } else if(ry) {
       ry = numeric.rep([x.length],0);
   }
   return this;

}

numeric.T.prototype.getBlock = function getBlock(from,to) {

   var x = this.x, y = this.y, b = numeric.getBlock;
   if(y) { return new numeric.T(b(x,from,to),b(y,from,to)); }
   return new numeric.T(b(x,from,to));

} numeric.T.prototype.setBlock = function setBlock(from,to,A) {

   if(!(A instanceof numeric.T)) A = new numeric.T(A);
   var x = this.x, y = this.y, b = numeric.setBlock, Ax = A.x, Ay = A.y;
   if(Ay) {
       if(!y) { this.y = numeric.rep(numeric.dim(this),0); y = this.y; }
       b(x,from,to,Ax);
       b(y,from,to,Ay);
       return this;
   }
   b(x,from,to,Ax);
   if(y) b(y,from,to,numeric.rep(numeric.dim(Ax),0));

} numeric.T.rep = function rep(s,v) {

   var T = numeric.T;
   if(!(v instanceof T)) v = new T(v);
   var x = v.x, y = v.y, r = numeric.rep;
   if(y) return new T(r(s,x),r(s,y));
   return new T(r(s,x));

} numeric.T.diag = function diag(d) {

   if(!(d instanceof numeric.T)) d = new numeric.T(d);
   var x = d.x, y = d.y, diag = numeric.diag;
   if(y) return new numeric.T(diag(x),diag(y));
   return new numeric.T(diag(x));

} numeric.T.eig = function eig() {

   if(this.y) { throw new Error('eig: not implemented for complex matrices.'); }
   return numeric.eig(this.x);

} numeric.T.identity = function identity(n) { return new numeric.T(numeric.identity(n)); } numeric.T.prototype.getDiag = function getDiag() {

   var n = numeric;
   var x = this.x, y = this.y;
   if(y) { return new n.T(n.getDiag(x),n.getDiag(y)); }
   return new n.T(n.getDiag(x));

}

// 4. Eigenvalues of real matrices

numeric.house = function house(x) {

   var v = numeric.clone(x);
   var s = x[0] >= 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;j<m-2;j++) {
       x = Array(m-j-1);
       for(i=j+1;i<m;i++) { x[i-j-1] = A[i][j]; }
       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<m;i++) { Ai = A[i]; Ci = C[i-j-1]; for(k=j;k<m;k++) Ai[k] -= 2*Ci[k-j]; }
       B = numeric.getBlock(A,[0,j+1],[m-1,m-1]);
       C = numeric.tensor(numeric.dot(B,v),v);
       for(i=0;i<m;i++) { Ai = A[i]; Ci = C[i]; for(k=j+1;k<m;k++) Ai[k] -= 2*Ci[k-j-1]; }
       B = Array(m-j-1);
       for(i=j+1;i<m;i++) B[i-j-1] = Q[i];
       C = numeric.tensor(v,numeric.dot(v,B));
       for(i=j+1;i<m;i++) { Qi = Q[i]; Ci = C[i-j-1]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; }
   }
   return {H:A, Q:Q};

}

numeric.epsilon = 2.220446049250313e-16;

numeric.QRFrancis = function(H,maxiter) {

   if(typeof maxiter === "undefined") { maxiter = 10000; }
   H = numeric.clone(H);
   var H0 = numeric.clone(H);
   var s = numeric.dim(H),m=s[0],x,v,a,b,c,d,det,tr, Hloc, Q = numeric.identity(m), Qi, Hi, B, C, Ci,i,j,k,iter;
   if(m<3) { return {Q:Q, B:[ [0,m-1] ]}; }
   var epsilon = numeric.epsilon;
   for(iter=0;iter<maxiter;iter++) {
       for(j=0;j<m-1;j++) {
           if(Math.abs(H[j+1][j]) < epsilon*(Math.abs(H[j][j])+Math.abs(H[j+1][j+1]))) {
               var QH1 = numeric.QRFrancis(numeric.getBlock(H,[0,0],[j,j]),maxiter);
               var QH2 = numeric.QRFrancis(numeric.getBlock(H,[j+1,j+1],[m-1,m-1]),maxiter);
               B = Array(j+1);
               for(i=0;i<=j;i++) { B[i] = Q[i]; }
               C = numeric.dot(QH1.Q,B);
               for(i=0;i<=j;i++) { Q[i] = C[i]; }
               B = Array(m-j-1);
               for(i=j+1;i<m;i++) { B[i-j-1] = Q[i]; }
               C = numeric.dot(QH2.Q,B);
               for(i=j+1;i<m;i++) { Q[i] = C[i-j-1]; }
               return {Q:Q,B:QH1.B.concat(numeric.add(QH2.B,j+1))};
           }
       }
       a = H[m-2][m-2]; b = H[m-2][m-1];
       c = H[m-1][m-2]; d = H[m-1][m-1];
       tr = a+d;
       det = (a*d-b*c);
       Hloc = numeric.getBlock(H, [0,0], [2,2]);
       if(tr*tr>=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<m;k++) Hi[k] -= 2*Ci[k]; }
       B = numeric.getBlock(H, [0,0],[m-1,2]);
       C = numeric.tensor(numeric.dot(B,v),v);
       for(i=0;i<m;i++) { Hi = H[i]; Ci = C[i]; for(k=0;k<3;k++) Hi[k] -= 2*Ci[k]; }
       B = [Q[0],Q[1],Q[2]];
       C = numeric.tensor(v,numeric.dot(v,B));
       for(i=0;i<3;i++) { Qi = Q[i]; Ci = C[i]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; }
       var J;
       for(j=0;j<m-2;j++) {
           for(k=j;k<=j+1;k++) {
               if(Math.abs(H[k+1][k]) < epsilon*(Math.abs(H[k][k])+Math.abs(H[k+1][k+1]))) {
                   var QH1 = numeric.QRFrancis(numeric.getBlock(H,[0,0],[k,k]),maxiter);
                   var QH2 = numeric.QRFrancis(numeric.getBlock(H,[k+1,k+1],[m-1,m-1]),maxiter);
                   B = Array(k+1);
                   for(i=0;i<=k;i++) { B[i] = Q[i]; }
                   C = numeric.dot(QH1.Q,B);
                   for(i=0;i<=k;i++) { Q[i] = C[i]; }
                   B = Array(m-k-1);
                   for(i=k+1;i<m;i++) { B[i-k-1] = Q[i]; }
                   C = numeric.dot(QH2.Q,B);
                   for(i=k+1;i<m;i++) { Q[i] = C[i-k-1]; }
                   return {Q:Q,B:QH1.B.concat(numeric.add(QH2.B,k+1))};
               }
           }
           J = Math.min(m-1,j+3);
           x = Array(J-j);
           for(i=j+1;i<=J;i++) { x[i-j-1] = H[i][j]; }
           v = numeric.house(x);
           B = numeric.getBlock(H, [j+1,j],[J,m-1]);
           C = numeric.tensor(v,numeric.dot(v,B));
           for(i=j+1;i<=J;i++) { Hi = H[i]; Ci = C[i-j-1]; for(k=j;k<m;k++) Hi[k] -= 2*Ci[k-j]; }
           B = numeric.getBlock(H, [0,j+1],[m-1,J]);
           C = numeric.tensor(numeric.dot(B,v),v);
           for(i=0;i<m;i++) { Hi = H[i]; Ci = C[i]; for(k=j+1;k<=J;k++) Hi[k] -= 2*Ci[k-j-1]; }
           B = Array(J-j);
           for(i=j+1;i<=J;i++) B[i-j-1] = Q[i];
           C = numeric.tensor(v,numeric.dot(v,B));
           for(i=j+1;i<=J;i++) { Qi = Q[i]; Ci = C[i-j-1]; for(k=0;k<m;k++) Qi[k] -= 2*Ci[k]; }
       }
   }
   throw new Error('numeric: eigenvalue iteration does not converge -- increase maxiter?');

}

numeric.eig = function eig(A,maxiter) {

   var QH = numeric.toUpperHessenberg(A);
   var QB = numeric.QRFrancis(QH.H,maxiter);
   var T = numeric.T;
   var n = A.length,i,k,flag = false,B = QB.B,H = numeric.dot(QB.Q,numeric.dot(QH.H,numeric.transpose(QB.Q)));
   var Q = new T(numeric.dot(QB.Q,QH.Q)),Q0;
   var m = B.length,j;
   var a,b,c,d,p1,p2,disc,x,y,p,q,n1,n2;
   var sqrt = Math.sqrt;
   for(k=0;k<m;k++) {
       i = B[k][0];
       if(i === B[k][1]) {
           // nothing
       } else {
           j = i+1;
           a = H[i][i];
           b = H[i][j];
           c = H[j][i];
           d = H[j][j];
           if(b === 0 && c === 0) continue;
           p1 = -a-d;
           p2 = a*d-b*c;
           disc = p1*p1-4*p2;
           if(disc>=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;j<n;j++) {
       if(j>0) {
           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<n;j++) {
       x = E.getRow(j);
       E.setRow(j,x.div(x.norm2()));
   }
   E = E.transpose();
   E = Q.transjugate().dot(E);
   return { lambda:R.getDiag(), E:E };

};

// 5. Real sparse linear algebra

numeric.sdim = function dim(A,ret,k) {

   if(typeof ret === "undefined") { ret = []; }
   if(typeof A !== "object") return ret;
   if(typeof k === "undefined") { k=0; }
   if(!(k in ret)) { ret[k] = 0; }
   if(A.length > 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) {

   if(typeof tol === "undefined") { tol = 1; }
   var n = A.length, i,j,k;
   var L = numeric.sidentity(n), U = numeric.sclone(A), UT = numeric.stranspose(U);
   var P = numeric.linspace(0,n-1),Q = numeric.linspace(0,n-1);
   var Ui, Uj, UTi, temp,alpha;
   var abs = Math.abs;
   for(i=0;i<n-1;i++) {
       UTi = UT[i];
       j = i;
       for(k in UTi) {
           if(!(Q.hasOwnProperty(k))) continue;
           k = Q[k];
           if(k<=i) continue;
           if(abs(U[k][i]) > abs(U[j][i])) { j = k; }
       }
       if(abs(U[i][i]) >= tol*abs(U[j][i])) { j = i; }
       if(j!==i) {
           temp = U[i]; U[i] = U[j]; U[j] = temp;
           temp = L[i]; L[i] = L[j]; L[j] = temp;
           temp = P[i]; P[i] = P[j]; P[j] = temp;
           Q = Array(n);
           for(j=0;j<n;j++) { Q[P[j]] = j; }
       }
       Ui = U[i];
       for(j in UTi) {
           if(!(Q.hasOwnProperty(j))) continue;
           j = Q[j];
           if(j<=i) continue;
           Uj = U[j];
           alpha = Uj[i]/Ui[i];
           L[j][i] = alpha;
           for(k in Ui) {
               if(!Ui.hasOwnProperty(k)) continue;
               if(k > i) {
                   if(!(k in Uj)) { Uj[k] = 0; UT[k][j] = 0; }
                   Uj[k] -= alpha*Ui[k];
               } else {
                   delete Uj[k];
               }
           }
       }
   }
   return {L:L, U:U, P:P, Pinv:Q};

};

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.sLUPsolve = function LUPsolve(lup,b) {

   var L = lup.L, U = lup.U, P = lup.P;
   var n = L.length, i,j, ret = Array(n), accum, Ai,foo;
   for(i = 0;i<n;i++) {
       if(b.hasOwnProperty(P[i])) accum = b[P[i]];
       else accum = 0;
       Ai = L[i];
       for(j in Ai) {
           if(!Ai.hasOwnProperty(j)) continue;
           if(j=0;i--) {
       accum = ret[i];
       Ai = U[i];
       for(j in Ai) {
           if(!Ai.hasOwnProperty(j)) continue;
           if(j>i) { accum -= Ai[j]*ret[j]; }
       }
       ret[i] = accum/Ai[i];
   }
   return ret;

   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<m-2;j++) {
           Vij = V[j][i];
           if(!Aj[Vij]) Aj[Vij] = [];
           Aj = Aj[Vij];
       }
       Aj[V[j][i]] = V[j+1][i];
   }
   return A;

   if(typeof ret === "undefined") ret = [];
   if(typeof k === "undefined") k = [];
   var n,i,Ai;
   n = k.length;
   for(i in A) {
       if(A.hasOwnProperty(i)) {
           k[n] = parseInt(i);
           Ai = A[i];
           if(typeof Ai === "number") {
               if(Ai) {
                   if(ret.length === 0) {
                       for(i=n+1;i>=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;

   var I = A[0], J = A[1], V = A[2];
   var p = I.length, m=0, i,j,k,a,b,c;