diff --git a/web/app/constr/solver.js b/web/app/constr/solver.js index ff38eb4b..00f74bd2 100644 --- a/web/app/constr/solver.js +++ b/web/app/constr/solver.js @@ -57,10 +57,9 @@ TCAD.parametric.System.prototype.makeJacobian = function() { return jacobi; }; -// ∇f(x) = E rj(x) ∇rj(x) = J(x)T r(x) -TCAD.parametric.System.prototype.calcGrad = function(out) { +TCAD.parametric.System.prototype.calcGrad_ = function(out) { var i; - for (i = 0; i < out.length; ++i) { + for (i = 0; i < out.length || i < this.params.length; ++i) { out[i][0] = 0; } @@ -74,10 +73,30 @@ TCAD.parametric.System.prototype.calcGrad = function(out) { for (var p = 0; p < cParams.length; p++) { var param = cParams[p]; var j = param.j; - out[j][0] += this.constraints[i].error() * grad[p]; // (10.4) + out[j][0] += this.constraints[i].error() * grad[p]; // (10.4) + } + } +}; + +TCAD.parametric.System.prototype.calcGrad = function(out) { + var i; + for (i = 0; i < out.length || i < this.params.length; ++i) { + out[i] = 0; + } + + for (i=0; i < this.constraints.length; i++) { + var c = this.constraints[i]; + + var cParams = c.params; + var grad = []; + c.gradient(grad); + + for (var p = 0; p < cParams.length; p++) { + var param = cParams[p]; + var j = param.j; + out[j] += this.constraints[i].error() * grad[p]; // (10.4) } } - }; TCAD.parametric.System.prototype.fillParams = function(out) { @@ -112,7 +131,7 @@ TCAD.parametric.System.prototype.errorSquare = function() { var t = this.constraints[i].error(); error += t * t; } - return error * .5; + return error * 0.5; }; TCAD.parametric.System.prototype.getValues = function() { @@ -152,7 +171,7 @@ TCAD.parametric.lock2 = function(constrs, locked) { } }; -TCAD.parametric.solve = function(constrs, locked, fineLevel) { +TCAD.parametric.solve = function(constrs, locked, fineLevel, alg) { if (constrs.length == 0) return; @@ -183,9 +202,24 @@ TCAD.parametric.solve = function(constrs, locked, fineLevel) { return sys.makeJacobian(); }; -// var res = TCAD.math.solve_BFGS(sys, 1e-4, 1e-4); -// console.log(res); -// return; + if (alg > 0) { + switch (alg) { + case 1: + var res = TCAD.math.solve_BFGS(sys, 1e-4, 1e-4); + console.log(res); + break; + case 2: + TCAD.math.solve_TR(sys); + break; + case 3: + TCAD.math.noptim(sys); + break; + case 4: + TCAD.math.solve_UNCMIN(sys); + break; + } + return sys; + } var opt = new LMOptimizer(sys.getParams(), arr(sys.constraints.length), model, jacobian); var eps; @@ -206,4 +240,5 @@ TCAD.parametric.solve = function(constrs, locked, fineLevel) { var res = opt.doOptimize(); sys.setParams(res[0]); // console.log("Solved with error: " + sys.error()); + return sys; }; diff --git a/web/app/main2d.js b/web/app/main2d.js index e1c32a08..9f92420f 100644 --- a/web/app/main2d.js +++ b/web/app/main2d.js @@ -120,7 +120,18 @@ TCAD.App2D = function() { solve : function() { app.viewer.parametricManager.solve([], 0); app.viewer.refresh(); + }, + + step : function() { + app.viewer.parametricManager.solve([], 0, 3); + app.viewer.refresh(); + }, + + stepUNCMIN : function() { + app.viewer.parametricManager.solve([], 0, 4); + app.viewer.refresh(); } + }; actionsF.add(actions, 'addSegment'); @@ -140,6 +151,8 @@ TCAD.App2D = function() { actionsF.add(actions, 'tangent'); actionsF.add(actions, 'lock'); actionsF.add(actions, 'solve'); + actionsF.add(actions, 'step'); + actionsF.add(actions, 'stepUNCMIN'); actionsF.add(actions, 'analyze'); actionsF.open(); diff --git a/web/app/math/math.js b/web/app/math/math.js index 3bbee53e..27ebea74 100644 --- a/web/app/math/math.js +++ b/web/app/math/math.js @@ -1,7 +1,16 @@ TCAD.math = {}; +TCAD.math._arr = function(size) { + var out = []; + out.length = size; + for (var i = 0; i < size; ++i) { + out[i] = 0; + } + return out; +}; + TCAD.math.distance = function(x1, y1, x2, y2) { var dx = x1 - x2; var dy = y1 - y2; return Math.sqrt(dx * dx + dy * dy); -} \ No newline at end of file +}; \ No newline at end of file diff --git a/web/app/math/matrix.js b/web/app/math/matrix.js new file mode 100644 index 00000000..f7356694 --- /dev/null +++ b/web/app/math/matrix.js @@ -0,0 +1,127 @@ + +TCAD.math.Matrix = function(r, c) { + this.data = []; + this.rSize = r; + this.cSize = c; + for (var i = 0; i < r; i++) { + this.data[i] = TCAD.math._arr(c) + } +}; + +TCAD.math.Matrix.prototype.identity = function() { + + for (var i = 0; i < this.rSize; i++) { + for (var j = 0; j < this.cSize; j++) { + this.data[i][j] = i === j ? 1 : 0; + } + } +}; + +TCAD.math.Matrix.prototype.subtract = function(m) { + var out = new TCAD.math.Matrix(this.rSize, this.cSize); + for (var i = 0; i < this.rSize; i++) { + for (var j = 0; j < this.cSize; j++) { + out.data[i][j] = this.data[i][j] - m.data[i][j]; + } + } + return out; +}; + +TCAD.math.Matrix.prototype.add = function(m) { + var out = new TCAD.math.Matrix(this.rSize, this.cSize); + for (var i = 0; i < this.rSize; i++) { + for (var j = 0; j < this.cSize; j++) { + out.data[i][j] = this.data[i][j] + m.data[i][j]; + } + } + return out; +}; + +TCAD.math.Matrix.prototype.multiply = function(m) { + + var nRows = this.rSize; + var nCols = m.cSize; + var nSum = this.cSize; + + var out = new TCAD.math.Matrix(nRows, nCols); + + var outData = out.data; + // Will hold a column of "m". + var mCol = TCAD.math._arr(nSum); + var mData = m.data; + + // Multiply. + for (var col = 0; col < nCols; col++) { + // Copy all elements of column "col" of "m" so that + // will be in contiguous memory. + for (var mRow = 0; mRow < nSum; mRow++) { + mCol[mRow] = mData[mRow][col]; + } + + for (var row = 0; row < nRows; row++) { + var dataRow = this.data[row]; + var sum = 0; + for (var i = 0; i < nSum; i++) { + sum += dataRow[i] * mCol[i]; + } + outData[row][col] = sum; + } + } + + return out; +}; + +TCAD.math.Matrix.prototype.scalarMultiply = function(s) { + var out = new TCAD.math.Matrix(this.rSize, this.cSize); + for (var i = 0; i < this.rSize; i++) { + for (var j = 0; j < this.cSize; j++) { + out.data[i][j] = this.data[i][j] * s; + } + } + return out; +}; + +TCAD.math.Matrix.prototype.transpose = function() { + var out = new TCAD.math.Matrix(this.cSize, this.rSize); + for (var i = 0; i < this.rSize; i++) { + for (var j = 0; j < this.cSize; j++) { + out.data[j][i] = this.data[i][j]; + } + } + return out; +}; + +TCAD.math.Matrix.prototype.copy = function() { + var out = new TCAD.math.Matrix(this.rSize, this.cSize); + for (var i = 0; i < this.rSize; i++) { + for (var j = 0; j < this.cSize; j++) { + out.data[i][j] = this.data[i][j]; + } + } + return out; +}; + +TCAD.math.Matrix.prototype.dot = function(v) { + var vData = v.data; + var dot = 0; + for (var i = 0; i < this.rSize; i++) { + dot += this.data[i][0] * vData[i][0]; + } + return dot; +}; + +TCAD.math.Matrix.prototype.norm = function(v) { + var sum = 0; + for (var i = 0; i < this.rSize; i++) { + var a = this.data[i][0]; + sum += a * a; + } + return Math.sqrt(sum); +}; + +TCAD.math.Vector = function(n) { + TCAD.math.Matrix.call(this, n, 1); +}; + +TCAD.TWO.utils.extend(TCAD.math.Vector, TCAD.math.Matrix); + diff --git a/web/app/math/noptim.js b/web/app/math/noptim.js new file mode 100644 index 00000000..30bafe87 --- /dev/null +++ b/web/app/math/noptim.js @@ -0,0 +1,59 @@ +numeric.grdec = 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}; +} \ No newline at end of file diff --git a/web/app/math/optim.js b/web/app/math/optim.js index 9efb529c..e916da58 100644 --- a/web/app/math/optim.js +++ b/web/app/math/optim.js @@ -10,7 +10,7 @@ TCAD.math.solve_BFGS = function(subsys, convergence, smallF) { if (xsize == 0) { return "Success"; } - + var xdir; //Vector var B = new TCAD.math.Matrix(xsize, xsize); B.identity(); @@ -21,7 +21,7 @@ TCAD.math.solve_BFGS = function(subsys, convergence, smallF) { // Initial unknowns vector and initial gradient vector TCAD.math.fillParams(subsys, x.data); - subsys.calcGrad(grad.data); + subsys.calcGrad_(grad.data); // Initial search direction oposed to gradient (steepest-descent) xdir = grad.scalarMultiply(-1); @@ -43,16 +43,16 @@ TCAD.math.solve_BFGS = function(subsys, convergence, smallF) { break; y = grad.copy(); - subsys.calcGrad(grad.data); + subsys.calcGrad_(grad.data); y = grad.subtract(y); // = grad - gradold //Now calculate the BFGS update on B -// TCAD.math.bfgsUpdate(B, h, y); - TCAD.math.bfgsUpdateInverse(B, y, h); - +// B = TCAD.math.bfgsUpdate(B, h, y); +// B = TCAD.math.bfgsUpdateInverse(B, y, h); xdir = B.multiply(grad).scalarMultiply(-1); // xdir = grad.scalarMultiply(-1); + TCAD.math.lineSearch(subsys, xdir); err = subsys.errorSquare(); @@ -68,6 +68,130 @@ TCAD.math.solve_BFGS = function(subsys, convergence, smallF) { return "Failed"; }; +TCAD.math.solve_UNCMIN = function(subsys) { + var x0 = []; + subsys.fillParams(x0); + + var f = function(x) { + subsys.setParams(x); + return subsys.errorSquare(); + }; + var gradient = function(x) { + subsys.setParams(x); + var grad = []; + subsys.calcGrad(grad); + return grad; + }; + numeric.uncmin(f,x0,0.01,gradient,1000); +}; + + +TCAD.math.solve_TR = function(subsys) { + + var xsize = subsys.params.length; + if (xsize == 0) { + return "Success"; + } + + var p; //Vector + var B = new TCAD.math.Matrix(xsize, xsize); + B.identity(); + var H = new TCAD.math.Matrix(xsize, xsize); + H.identity(); + + var x = new TCAD.math.Vector(xsize); + var grad = new TCAD.math.Vector(xsize); + var h = new TCAD.math.Vector(xsize); + var y = new TCAD.math.Vector(xsize); + var pZero = new TCAD.math.Vector(xsize); + + TCAD.math.fillParams(subsys, x.data); + subsys.calcGrad_(grad.data); + p = grad.scalarMultiply(-1); + + var err = subsys.errorSquare(); + + var delta = 0.01; + var deltaD = 1; + var nu = 0.1; + + + var maxIterNumber = 100 * xsize; + var divergingLim = 1e6*err + 1e12; + + + function m(fx, p, g, B) { + return fx + g.dot(p) + 0.5 * p.dot(B.multiply(p)); //4.3 + } + + function norm(x) { + var sum = 0; + for (var i = 0; i < x.rSize; i++) { + var a = x.data[i][0]; + sum += a * a; + } + return Math.sqrt(sum); + } + + var tolx = 0.01; + for (var iter=1; iter < maxIterNumber; iter++) { + + TCAD.math.setParams2(subsys, x.data); + var fx = subsys.errorSquare(); + + if (fx <= tolx) break; + //if (h.norm() <= tolx) break; +// if (delta <= tolx*(tolx + x.norm())) break; + if (fx > divergingLim || fx != fx) break; + + p = TCAD.math.cauchyPoint(delta, grad, B, norm); + + var xAddP = x.add(p); + + TCAD.math.setParams2(subsys, xAddP.data); + var fxAddP = subsys.errorSquare(); + + var r = ( fx - fxAddP ) / ( m(fx, pZero, grad, H) - m(fx, p, grad, H) ); + + if (r < 0.25) { + delta = 0.25 * norm(p); + } else { + if (r > 0.75) { + delta = Math.min(delta * 2, deltaD); + } + } + + if (r > nu) { + h = xAddP.subtract(x); // = x - xold + x = xAddP; + + TCAD.math.setParams2(subsys, x.data); + + y = grad.copy(); + subsys.calcGrad_(grad.data); + y = grad.subtract(y); + + //Now calculate the BFGS on B and H + B = TCAD.math.bfgsUpdateInverse(B, y, h); + H = TCAD.math.bfgsUpdate(H, y, h); + } + } + TCAD.math.setParams2(subsys, x.data); +}; + + +TCAD.math.cauchyPoint = function(delta, grad, B, normaF) { + var tau; + var tauCondition = grad.dot(B.multiply(grad)) + var norm = normaF(grad); + if (tauCondition <= 0) { + tau = 1; + } else { + tau = Math.min((norm*norm*norm)/(delta*tauCondition), 1); + } + return grad.scalarMultiply(- tau * delta / norm) ; +}; + TCAD.math.fillParams = function(sys, out) { for (var p = 0; p < sys.params.length; p++) { out[p][0] = sys.params[p].get(); @@ -113,13 +237,14 @@ TCAD.math.bfgsUpdate = function(B, y, h) { B = B.add( y_x_yT.scalarMultiply( 1 / yT_x_h ) ); B = B.subtract( ( B_x_h.multiply(hT_x_B) ).scalarMultiply( 1./hT_x_B_x_h ) ); + return B; }; TCAD.math.solve_SD = function(subsys) { var i = 0; var grad = new TCAD.math.Vector(subsys.params.length); while (subsys.errorSquare() > 0.1 ) { - subsys.calcGrad(grad.data); + subsys.calcGrad_(grad.data); var xdir = grad.scalarMultiply(-1); TCAD.math.lineSearch(subsys, xdir); if (i ++ > 100) { @@ -129,7 +254,7 @@ TCAD.math.solve_SD = function(subsys) { console.log(subsys.errorSquare()); }; -TCAD.math.lineSearch = function(subsys, xdir) { +TCAD.math.lineSearchOrig = function(subsys, xdir) { var f1,f2,f3,alpha1,alpha2,alpha3,alphaStar; @@ -199,11 +324,191 @@ TCAD.math.lineSearch = function(subsys, xdir) { //Take a final step to alphaStar x = x0 .add( xdir.scalarMultiply( alphaStar ) ); TCAD.math.setParams2(subsys, x.data); - return alphaStar; - + }; +TCAD.math.lineSearchWeight = function(subsys, xdir) { + + var f1,f2,f3,alpha1,alpha2,alpha3,alphaStar; + + + var alphaMax = 1e28; //maxStep(xdir); + + var x; + var x0 = new TCAD.math.Vector(subsys.params.length); + + var costs = []; + function updateCosts() { + var maxErr = -1; + var i; + var t; + for (i=0; i < subsys.constraints.length; i++) { + t = subsys.constraints[i].error(); + maxErr = Math.max(maxErr, t*t); + } + if (maxErr > 0) { + for (i=0; i < subsys.constraints.length; i++) { + t = subsys.constraints[i].error(); + costs[i] = t*t / maxErr; + } + } else { + TCAD.math.fill_array(costs, 0, subsys.constraints.length, 1) + } + } + updateCosts(); +// console.log(costs); + var xdir = new TCAD.math.Vector(subsys.params.length); + calcGrad = function(out) { + var i; + for (i = 0; i < out.length; ++i) { + out[i][0] = 0; + } + + for (i=0; i < subsys.constraints.length; i++) { + var c = subsys.constraints[i]; + + var cParams = c.params; + var grad = []; + c.gradient(grad); + + for (var p = 0; p < cParams.length; p++) { + var param = cParams[p]; + var j = param.j; + out[j][0] += costs[i] * grad[p]; // (10.4) + } + } + }; + calcGrad(xdir.data) + console.log(xdir.data); + + function errorSquare() { + var error = 0; + for (var i = 0; i < subsys.constraints.length; i++) { + var t = subsys.constraints[i].error(); + error += t * t * costs[i]; + } + return error * 0.5; + } + + + //Save initial values + TCAD.math.fillParams(subsys, x0.data); + + //Start at the initial position alpha1 = 0 + alpha1 = 0.; + f1 = errorSquare(); + + + //Take a step of alpha2 = 1 + alpha2 = 1.; + x = x0.add(xdir.scalarMultiply(alpha2)); + TCAD.math.setParams2(subsys, x.data); + f2 = errorSquare(); + + //Take a step of alpha3 = 2*alpha2 + alpha3 = alpha2*2; + x = x0.add(xdir.scalarMultiply(alpha3)); + TCAD.math.setParams2(subsys, x.data); + f3 = errorSquare(); + + //Now reduce or lengthen alpha2 and alpha3 until the minimum is + //Bracketed by the triplet f1>f2 f1 || f2 > f3) { + if (f2 > f1) { + //If f2 is greater than f1 then we shorten alpha2 and alpha3 closer to f1 + //Effectively both are shortened by a factor of two. + alpha3 = alpha2; + f3 = f2; + alpha2 = alpha2 / 2; + x = x0.add( xdir.scalarMultiply(alpha2 )); + TCAD.math.setParams2(subsys, x.data); + f2 = errorSquare(); + } + else if (f2 > f3) { + if (alpha3 >= alphaMax) + break; + //If f2 is greater than f3 then we increase alpha2 and alpha3 away from f1 + //Effectively both are lengthened by a factor of two. + alpha2 = alpha3; + f2 = f3; + alpha3 = alpha3 * 2; + x = x0.add( xdir.scalarMultiply(alpha3)); + TCAD.math.setParams2(subsys, x.data); + f3 = errorSquare(); + } + } + //Get the alpha for the minimum f of the quadratic approximation + alphaStar = alpha2 + ((alpha2-alpha1)*(f1-f3))/(3*(f1-2*f2+f3)); + + //Guarantee that the new alphaStar is within the bracket + if (alphaStar >= alpha3 || alphaStar <= alpha1) + alphaStar = alpha2; + + if (alphaStar > alphaMax) + alphaStar = alphaMax; + + if (alphaStar != alphaStar) + alphaStar = 0.; + + //Take a final step to alphaStar + x = x0 .add( xdir.scalarMultiply( alphaStar ) ); + TCAD.math.setParams2(subsys, x.data); +// console.log(alphaStar); + return alphaStar; + +}; + +TCAD.math.lineSearchWolfeCond = function(sys, d) { + + var c1 = 0.1; + var c2 = 0.9; + + var x0 = new TCAD.math.Vector(sys.params.length); + TCAD.math.fillParams(sys, x0.data); + + var alpha = 1; + var fx0 = sys.errorSquare(); + var grad = new TCAD.math.Vector(sys.params.length); + sys.calcGrad_(grad.data); + var gx0 = grad.dot(d); + + //bound the solution + var alphaL = 0; + var alphaR = 10000; + var maxit = 800; + + for (var iter = 1; iter <= maxit; iter++){ + var xp = x0.add(d.scalarMultiply(alpha)); + + TCAD.math.setParams2(sys, xp.data); + var erroralpha = sys.errorSquare(); //get the error at that point + + if (erroralpha >= fx0 + alpha * c1 * gx0) { // if error is not sufficiently reduced + alphaR = alpha;//move halfway between current alpha and lower alpha + alpha = (alphaL + alphaR)/2.0; + }else{//if error is sufficiently decreased + + TCAD.math.setParams2(sys, xp.data); + sys.calcGrad_(grad.data); + var slopealpha = grad.dot(d); // then get slope along search direction + if (slopealpha <= c2 * Math.abs(gx0)){ // if slope sufficiently closer to 0 + break;//then this is an acceptable point + }else if ( slopealpha >= c2 * gx0) { // if slope is too steep and positive then go to the left + alphaR = alpha;//move halfway between current alpha and lower alpha + alpha = (alphaL+ alphaR)/2; + }else{//if slope is too steep and negative then go to the right of this alpha + alphaL = alpha;//move halfway between current alpha and upper alpha + alpha = (alphaL+ alphaR)/2; + } + } + } + + //if ran out of iterations then return the best thing we got + var x = x0.add(d.scalarMultiply(alpha)); + TCAD.math.setParams2(sys, x.data); + return alpha; +}; TCAD.math.lineSearch3 = function(sys, xdir) { @@ -342,4 +647,12 @@ TCAD.math.lineSearch2 = function(subsys, xdir) { return alphaStar; +}; + + +TCAD.math.lineSearch = TCAD.math.lineSearchWeight; +//TCAD.math.lineSearch = TCAD.math.lineSearchOrig; + +TCAD.math.fill_array = function(a, fromIndex, toIndex,val) { + for (var i = fromIndex; i < toIndex; i++) a[i] = val; }; \ No newline at end of file diff --git a/web/app/math/test.js b/web/app/math/test.js new file mode 100644 index 00000000..86b52075 --- /dev/null +++ b/web/app/math/test.js @@ -0,0 +1,53 @@ + +function perpAndCoinc() { + var sys = []; + var ID = 0; + var l1 = { + p1 : { + x : new TCAD.parametric.Param(ID++, 10), + y : new TCAD.parametric.Param(ID++, 10) + }, + p2 : { + x : new TCAD.parametric.Param(ID++, 300), + y : new TCAD.parametric.Param(ID++, 300) + } + }; + + + var l2 = { + p1 : { + x : new TCAD.parametric.Param(ID++, 10 + 100), + y : new TCAD.parametric.Param(ID++, -10 + 100) + }, + p2 : { + x : new TCAD.parametric.Param(ID++, 300 + 100), + y : new TCAD.parametric.Param(ID++, -300 + 100) + } + }; + + + sys.push(new TCAD.constraints.Perpendicular([l1.p1.x, l1.p1.y, l1.p2.x, l1.p2.y, l2.p1.x, l2.p1.y, l2.p2.x, l2.p2.y])); + sys.push(new TCAD.constraints.Equal([l1.p1.x, l2.p1.x])); + sys.push(new TCAD.constraints.Equal([l1.p1.y, l2.p1.y])); + + return sys; + +} + + +function testCompare() { + var bfgs_ = perpAndCoinc(); + var lm_ = perpAndCoinc(); + var tr_ = perpAndCoinc(); + + var bfgs = TCAD.parametric.solve(bfgs_, [], 0, 1); + var lm = TCAD.parametric.solve(lm_, [], 0, 0); + var tr = TCAD.parametric.solve(tr_, [], 0, 2); + + console.log("bfgs: " + bfgs.errorSquare()); + console.log("lm: " + lm.errorSquare()); + console.log("trusted region: " + tr.errorSquare()); +} + + + diff --git a/web/app/parametric.js b/web/app/parametric.js index 12e7ccd8..a9a737fe 100644 --- a/web/app/parametric.js +++ b/web/app/parametric.js @@ -196,7 +196,7 @@ TCAD.TWO.ParametricManager.prototype.solve1 = function(locked, onSolved) { }; -TCAD.TWO.ParametricManager.prototype.solve = function(locked, fineLevel) { +TCAD.TWO.ParametricManager.prototype.solve = function(locked, fineLevel, alg) { var pdict = {}; var params; var _constrs = []; @@ -242,7 +242,7 @@ TCAD.TWO.ParametricManager.prototype.solve = function(locked, fineLevel) { } } - TCAD.parametric.solve(_constrs, _locked, fineLevel); + TCAD.parametric.solve(_constrs, _locked, fineLevel, alg); for (p in pdict) { _p = pdict[p]; diff --git a/web/canvas.html b/web/canvas.html index d574b32b..99ba9ccd 100644 --- a/web/canvas.html +++ b/web/canvas.html @@ -10,11 +10,15 @@ + + + + diff --git a/web/lib/numeric-1.2.6.js b/web/lib/numeric-1.2.6.js new file mode 100644 index 00000000..1714d6d8 --- /dev/null +++ b/web/lib/numeric-1.2.6.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;k + + TCAD + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file