From cc3e0ffe8cca969669a015fefa757187dbe6b534 Mon Sep 17 00:00:00 2001 From: Martin von Gagern Date: Tue, 18 Nov 2014 16:36:35 +0100 Subject: [PATCH 1/2] Change CRLF to newlines only. --- src/svd.js | 592 ++++++++++++++++++++++++++--------------------------- 1 file changed, 296 insertions(+), 296 deletions(-) diff --git a/src/svd.js b/src/svd.js index 24a07d3..1db88e6 100644 --- a/src/svd.js +++ b/src/svd.js @@ -1,296 +1,296 @@ -/* -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;k 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;k Date: Tue, 18 Nov 2014 16:42:46 +0100 Subject: [PATCH 2/2] Explicitely state all semicolons. This avoids relying on automatic semicolon insertion, to make the code cleaner and avoid pitfalls. The svd module in particular had many unterminated lines, likely due to its Python heritage. The insertion places were identified by JSHint. --- src/numeric.js | 250 ++++++++++++++++++++++++------------------------- src/sparse2.js | 18 ++-- src/svd.js | 214 +++++++++++++++++++++--------------------- 3 files changed, 241 insertions(+), 241 deletions(-) diff --git a/src/numeric.js b/src/numeric.js index 537b68f..1a115c6 100644 --- a/src/numeric.js +++ b/src/numeric.js @@ -22,7 +22,7 @@ numeric.bench = function bench (f,interval) { 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; @@ -95,7 +95,7 @@ numeric.prettyPrint = function prettyPrint(x) { } foo(x); return ret.join(''); -} +}; numeric.parseDate = function parseDate(d) { function foo(d) { @@ -106,7 +106,7 @@ numeric.parseDate = function parseDate(d) { return ret; } return foo(d); -} +}; numeric.parseFloat = function parseFloat_(d) { function foo(d) { @@ -117,7 +117,7 @@ numeric.parseFloat = function parseFloat_(d) { return ret; } return foo(d); -} +}; numeric.parseCSV = function parseCSV(t) { var foo = t.split('\n'); @@ -125,7 +125,7 @@ numeric.parseCSV = function parseCSV(t) { 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 stripper = function(n) { return n.substr(0,n.length-1); }; var count = 0; for(k=0;k=0;i--) { ret[i] = numeric.rep(s,v,k+1); } return ret; -} +}; numeric.dotMMsmall = function dotMMsmall(x,y) { @@ -413,7 +413,7 @@ numeric.dotMMsmall = function dotMMsmall(x,y) { ret[i] = foo; } return ret; -} +}; numeric._getCol = function _getCol(A,j,x) { var n = A.length, i; for(i=n-1;i>0;--i) { @@ -422,7 +422,7 @@ numeric._getCol = function _getCol(A,j,x) { 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; @@ -441,14 +441,14 @@ numeric.dotMMbig = function dotMMbig(x,y){ } } 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; @@ -464,7 +464,7 @@ numeric.dotVM = function dotVM(x,y) { ret[k] = woo; } return ret; -} +}; numeric.dotVV = function dotVV(x,y) { var i,n=x.length,i1,ret = x[n-1]*y[n-1]; @@ -474,7 +474,7 @@ numeric.dotVV = function dotVV(x,y) { } if(i===0) { ret += x[0]*y[0]; } return ret; -} +}; numeric.dot = function dot(x,y) { var d = numeric.dim; @@ -490,7 +490,7 @@ numeric.dot = function dot(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; @@ -511,7 +511,7 @@ numeric.diag = function diag(d) { 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) { @@ -523,9 +523,9 @@ numeric.getDiag = function(A) { ret[0] = A[0][0]; } return ret; -} +}; -numeric.identity = function identity(n) { return numeric.diag(numeric.rep([n],1)); } +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 = []; @@ -558,7 +558,7 @@ numeric.pointwise = function pointwise(params,body,setup) { 'return ret;' ); return Function.apply(null,fun); -} +}; numeric.pointwise2 = function pointwise2(params,body,setup) { if(typeof setup === "undefined") { setup = ""; } var fun = []; @@ -583,7 +583,7 @@ numeric.pointwise2 = function pointwise2(params,body,setup) { 'return ret;' ); return 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]; @@ -790,7 +790,7 @@ numeric.trunc = function trunc(x,y) { } if (typeof y === "object") return numeric.truncSV(x,y); return Math.round(x/y)*y; -} +}; numeric.inv = function inv(x) { var s = numeric.dim(x), abs = Math.abs, m = s[0], n = s[1]; @@ -818,7 +818,7 @@ numeric.inv = function inv(x) { } } return I; -} +}; numeric.det = function det(x) { var s = numeric.dim(x); @@ -846,7 +846,7 @@ numeric.det = function det(x) { 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; @@ -873,7 +873,7 @@ numeric.transpose = function transpose(x) { 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=0;i--) ret[i] = _random(s,k+1); return ret; -} -numeric.random = function random(s) { return numeric._random(s,0); } +}; +numeric.random = function random(s) { return numeric._random(s,0); }; -numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); } +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); @@ -926,7 +926,7 @@ numeric.linspace = function linspace(a,b,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); @@ -940,7 +940,7 @@ numeric.getBlock = function getBlock(x,from,to) { return ret; } return foo(x,0); -} +}; numeric.setBlock = function setBlock(x,from,to,B) { var s = numeric.dim(x); @@ -951,7 +951,7 @@ numeric.setBlock = function setBlock(x,from,to,B) { } foo(x,B,0); return x; -} +}; numeric.getRange = function getRange(A,I,J) { var m = I.length, n = J.length; @@ -964,7 +964,7 @@ numeric.getRange = function getRange(A,I,J) { for(j=n-1;j!==-1;--j) Bi[j] = AI[J[j]]; } return B; -} +}; numeric.blockMatrix = function blockMatrix(X) { var s = numeric.dim(X); @@ -990,7 +990,7 @@ numeric.blockMatrix = function blockMatrix(X) { I += X[i][0].length; } return Z; -} +}; numeric.tensor = function tensor(x,y) { if(typeof x === "number" || typeof y === "number") return numeric.mul(x,y); @@ -1015,11 +1015,11 @@ numeric.tensor = function tensor(x,y) { 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.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; @@ -1047,7 +1047,7 @@ numeric.Tbinop = function Tbinop(rr,rc,cr,cc,setup) { '}\n'+ 'return new numeric.T('+rr+');\n' ); -} +}; numeric.T.prototype.add = numeric.Tbinop( 'add(x.x,y.x)', @@ -1072,14 +1072,14 @@ numeric.T.prototype.reciprocal = function reciprocal() { 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)', @@ -1090,12 +1090,12 @@ 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( @@ -1106,7 +1106,7 @@ numeric.Tunop = function Tunop(r,c,s) { '}\n'+ r+';\n' ); -} +}; numeric.T.prototype.exp = numeric.Tunop( 'return new numeric.T(ex)', @@ -1202,7 +1202,7 @@ numeric.T.prototype.inv = function inv() { } } 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) { @@ -1220,7 +1220,7 @@ numeric.T.prototype.get = function get(i) { 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) { @@ -1267,7 +1267,7 @@ numeric.T.prototype.set = function set(i,v) { 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; @@ -1278,7 +1278,7 @@ numeric.T.prototype.getRows = function getRows(i0,i1) { 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; @@ -1290,12 +1290,12 @@ numeric.T.prototype.setRows = function setRows(i0,i1,A) { 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; @@ -1306,13 +1306,13 @@ numeric.T.prototype.setRow = function setRow(i,v) { 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; @@ -1324,31 +1324,31 @@ numeric.T.prototype.setBlock = function setBlock(from,to,A) { } 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.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 @@ -1360,7 +1360,7 @@ numeric.house = function house(x) { 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); @@ -1384,7 +1384,7 @@ numeric.toUpperHessenberg = function toUpperHessenberg(me) { } } return {H:A, Q:Q}; -} +}; numeric.epsilon = 2.220446049250313e-16; @@ -1473,7 +1473,7 @@ numeric.QRFrancis = function(H,maxiter) { } } throw new Error('numeric: eigenvalue iteration does not converge -- increase maxiter?'); -} +}; numeric.eig = function eig(A,maxiter) { var QH = numeric.toUpperHessenberg(A); @@ -1589,7 +1589,7 @@ numeric.ccsSparse = function ccsSparse(A) { } } return [Ai,Aj,Av]; -} +}; numeric.ccsFull = function ccsFull(A) { var Ai = A[0], Aj = A[1], Av = A[2], s = numeric.ccsDim(A), m = s[0], n = s[1], i,j,j0,j1,k; var B = numeric.rep([m,n],0); @@ -1599,7 +1599,7 @@ numeric.ccsFull = function ccsFull(A) { for(j=j0;jn) k.pop(); return ret; -} +}; // 6. Coordinate matrices numeric.cLU = function LU(A) { @@ -2338,7 +2338,7 @@ numeric.cgrid = function grid(n,shape) { if(typeof shape !== "function") { switch(shape) { case 'L': - shape = function(i,j) { return (i>=n[0]/2 || j=n[0]/2 || j b) - return a*Math.sqrt(1.0+(b*b/a/a)) + 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)) + return a; + return b*Math.sqrt(1.0+(a*a/b/b)); } //Householder's reduction to bidiagonal form @@ -64,44 +64,44 @@ numeric.svd= function svd(A) { f= u[i][i]; g= Math.sqrt(s); if (f >= 0.0) g= -g; - h= f*g-s + h= f*g-s; u[i][i]=f-g; for (j=l; j < n; j++) { - s= 0.0 + s= 0.0; for (k=i; k < m; k++) - s += u[k][i]*u[k][j] - f= s/h + s += u[k][i]*u[k][j]; + f= s/h; for (k=i; k < m; k++) - u[k][j]+=f*u[k][i] + u[k][j]+=f*u[k][i]; } } - q[i]= g - s= 0.0 + q[i]= g; + s= 0.0; for (j=l; j < n; j++) - s= s + u[i][j]*u[i][j] + s= s + u[i][j]*u[i][j]; if (s <= tolerance) - g= 0.0 + g= 0.0; else { - f= u[i][i+1] - g= Math.sqrt(s) - if (f >= 0.0) g= -g - h= f*g - s + 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 < n; j++) e[j]= u[i][j]/h; for (j=l; j < m; j++) { - s=0.0 + s=0.0; for (k=l; k < n; k++) - s += (u[j][k]*u[i][k]) + s += (u[j][k]*u[i][k]); for (k=l; k < n; k++) - u[j][k]+=s*e[k] + u[j][k]+=s*e[k]; } } - y= Math.abs(q[i])+Math.abs(e[i]) + y= Math.abs(q[i])+Math.abs(e[i]); if (y>x) - x=y + x=y; } // accumulation of right hand gtransformations @@ -109,16 +109,16 @@ numeric.svd= function svd(A) { { if (g != 0.0) { - h= g*u[i][i+1] + h= g*u[i][i+1]; for (j=l; j < n; j++) - v[j][i]=u[i][j]/h + v[j][i]=u[i][j]/h; for (j=l; j < n; j++) { - s=0.0 + s=0.0; for (k=l; k < n; k++) - s += u[i][k]*v[k][j] + s += u[i][k]*v[k][j]; for (k=l; k < n; k++) - v[k][j]+=(s*v[k][i]) + v[k][j]+=(s*v[k][i]); } } for (j=l; j < n; j++) @@ -127,25 +127,25 @@ numeric.svd= function svd(A) { v[j][i] = 0; } v[i][i] = 1; - g= e[i] - l= i + g= e[i]; + l= i; } // accumulation of left hand transformations for (i=n-1; i != -1; i+= -1) { - l= i+1 - g= q[i] + 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 + h= u[i][i]*g; for (j=l; j < n; j++) { - s=0.0 + s=0.0; for (k=l; k < m; k++) s += u[k][i]*u[k][j]; - f= s/h + 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; @@ -156,119 +156,119 @@ numeric.svd= function svd(A) { } // diagonalization of the bidiagonal form - prec= prec*x + 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 + var test_convergence = false; for (l=k; l != -1; l+= -1) { if (Math.abs(e[l]) <= prec) - { test_convergence= true - break + { test_convergence= true; + break; } if (Math.abs(q[l-1]) <= prec) - break + break; } if (!test_convergence) { // cancellation of e[l] if l>0 - c= 0.0 - s= 1.0 - var l1= l-1 + c= 0.0; + s= 1.0; + var l1= l-1; for (i =l; i= itmax-1) - throw 'Error: no convergence.' + 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) + 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 + f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x; else - f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x + f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x; // next QR transformation - c= 1.0 - s= 1.0 + 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 + 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 + 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 + 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 + 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 + e[l]= 0.0; + e[k]= f; + q[k]= x; } } //vt= transpose(v) //return (u,q,vt) for (i=0;i