diff --git a/web/app/math/noptim.js b/web/app/math/noptim.js index 09c1a5c1..4f3204cf 100644 --- a/web/app/math/noptim.js +++ b/web/app/math/noptim.js @@ -257,17 +257,17 @@ optim.inv = function inv(x) { return I; } -optim.dog_leg = function(subsys, rough) { +optim.dog_leg = function (subsys, rough) { rough = true - var tolg= rough ? 1e-4 : 1e-4; + var tolg = rough ? 1e-4 : 1e-4; - var tolx=1e-80, tolf=1e-10; + var tolx = 1e-80, tolf = 1e-10; var xsize = subsys.params.length; var csize = subsys.constraints.length; if (xsize == 0) - return 'Success'; + return 'Success'; var vec = TCAD.math._arr; var mx = TCAD.math._matrix; @@ -288,23 +288,22 @@ optim.dog_leg = function(subsys, rough) { var h_dl = vec(xsize); var r0 = vec(csize); - + var err; subsys.fillParams(x); - + // subsys.setParams(vec(xsize)); // subsys.calcResidual(r0); - + subsys.setParams(x); err = subsys.calcResidual(fx); - + subsys.fillJacobian(Jx); function lls(A, b) { var At = n.transpose(A); var J = n.dot(At, A); - var r = n.dot(At, b);; - + var r = n.dot(At, b); return nocadel_10_18(J, r); } @@ -334,7 +333,7 @@ optim.dog_leg = function(subsys, rough) { // return n.solve(A, b); // } var At = n.transpose(A); - var res = n.dot(n.dot(At, optim.inv(n.dot(A, At)) ), b); + var res = n.dot(n.dot(At, optim.inv(n.dot(A, At))), b); return res; } @@ -352,149 +351,148 @@ optim.dog_leg = function(subsys, rough) { var g_inf = n.norminf(gNg); var fx_inf = n.norminf(fx); - var maxIterNumber = xsize * 100; - var divergingLim = 1e6*err + 1e12; + var maxIterNumber = xsize * 100; + var divergingLim = 1e6 * err + 1e12; - var delta=0.1; - var alpha=0.; - var nu=2.; - var iter=0, stop=0, reduce=0; + var delta = 0.1; + var alpha = 0.; + var nu = 2.; + var iter = 0, stop = 0, reduce = 0; while (stop === 0) { - // check if finished - if (fx_inf <= tolf || (rough && err <= tolf)) // Success - stop = 1; - else if (g_inf <= tolg) - stop = 2; - else if (delta <= tolx*(tolx + n.norm2(x))) - stop = 2; - else if (iter >= maxIterNumber) - stop = 4; - else if (err > divergingLim || err != err) { // check for diverging and NaN - stop = 6; - } - else { - // get the steepest descent direction - var Jt = n.transpose(Jx); - var B = n.dot(Jt, Jx); - var gBg = n.dot(g, n.dot(B, g)); - alpha = n.norm2Squared(g) / gBg; - h_sd = n.mul(g, - alpha); + // check if finished + if (fx_inf <= tolf || (rough && err <= tolf)) // Success + stop = 1; + else if (g_inf <= tolg) + stop = 2; + else if (delta <= tolx * (tolx + n.norm2(x))) + stop = 2; + else if (iter >= maxIterNumber) + stop = 4; + else if (err > divergingLim || err != err) { // check for diverging and NaN + stop = 6; + } + else { + // get the steepest descent direction + var Jt = n.transpose(Jx); + var B = n.dot(Jt, Jx); + var gBg = n.dot(g, n.dot(B, g)); + alpha = n.norm2Squared(g) / gBg; + h_sd = n.mul(g, -alpha); - // get the gauss-newton step - //h_gn = n.solve(Jx, n.mul(fx, -1)); - h_gn = lsolve(Jx, n.mul(fx, -1)); + // get the gauss-newton step + //h_gn = n.solve(Jx, n.mul(fx, -1)); + h_gn = lsolve(Jx, n.mul(fx, -1)); - //LU-Decomposition + //LU-Decomposition // h_gn = lusolve(Jx, n.mul(fx, -1)); - //Conjugate gradient method - //h_gn = optim.cg(Jx, h_gn, n.mul(fx, -1), 1e-8, maxIterNumber); - - //solve linear problem using svd formula to get the gauss-newton step - //h_gn = lls(Jx, n.mul(fx, -1)); - - var rel_error = n.norm2(n.add(n.dot(Jx, h_gn), fx)) / n.norm2(fx); - if (rel_error > 1e15) - break; + //Conjugate gradient method + //h_gn = optim.cg(Jx, h_gn, n.mul(fx, -1), 1e-8, maxIterNumber); - var hitBoundary = false; + //solve linear problem using svd formula to get the gauss-newton step + //h_gn = lls(Jx, n.mul(fx, -1)); - // compute the dogleg step - var gnorm = n.norm2(gNg); - if (n.norm2(h_gn) < delta) { - h_dl = n.clone(h_gn); - if (n.norm2(h_dl) <= tolx*(tolx + n.norm2(x))) { - stop = 5; - break; - } - } - else if (alpha* gnorm >= delta) { - var normRadius = delta/gnorm; - if( alpha >= normRadius ) alpha = normRadius; - h_dl = n.mul(g, - alpha); - hitBoundary = true; - } - else { - //compute beta - var d = n.sub(h_gn, h_sd); + var rel_error = n.norm2(n.add(n.dot(Jx, h_gn), fx)) / n.norm2(fx); + if (rel_error > 1e15) + break; - var a = n.dot(d, d); - var b = 2 * n.dot(h_sd, d); - var c = n.dot(h_sd, h_sd) - delta*delta + var hitBoundary = false; - var sqrt_discriminant = Math.sqrt(b*b - 4*a*c) - - var beta = (-b + sqrt_discriminant) / (2*a) - - // and update h_dl and dL with beta - h_dl = n.add(h_sd, n.mul(beta,d)); - hitBoundary = true; - } - } - - - // see if we are already finished - if (stop) + // compute the dogleg step + var gnorm = n.norm2(gNg); + if (n.norm2(h_gn) < delta) { + h_dl = n.clone(h_gn); + if (n.norm2(h_dl) <= tolx * (tolx + n.norm2(x))) { + stop = 5; break; + } + } + else if (alpha * gnorm >= delta) { + var normRadius = delta / gnorm; + if (alpha >= normRadius) alpha = normRadius; + h_dl = n.mul(g, -alpha); + hitBoundary = true; + } + else { + //compute beta + var d = n.sub(h_gn, h_sd); - // get the new values - var err_new; - x_new = n.add(x, h_dl); - subsys.setParams(x_new); - err_new = subsys.calcResidual(fx_new); - subsys.fillJacobian(Jx_new); + var a = n.dot(d, d); + var b = 2 * n.dot(h_sd, d); + var c = n.dot(h_sd, h_sd) - delta * delta - // calculate the linear model and the update ratio + var sqrt_discriminant = Math.sqrt(b * b - 4 * a * c) - var dF = err - err_new; - var dL = - n.dot(g, h_dl) - 0.5 * n.dot(h_dl, n.dot(B, h_dl)); + var beta = (-b + sqrt_discriminant) / (2 * a) - var acceptCandidate; + // and update h_dl and dL with beta + h_dl = n.add(h_sd, n.mul(beta, d)); + hitBoundary = true; + } + } - if( dF == 0 || dL == 0 ) { - acceptCandidate = true; + + // see if we are already finished + if (stop) + break; + + // get the new values + var err_new; + x_new = n.add(x, h_dl); + subsys.setParams(x_new); + err_new = subsys.calcResidual(fx_new); + subsys.fillJacobian(Jx_new); + + // calculate the linear model and the update ratio + + var dF = err - err_new; + var dL = -n.dot(g, h_dl) - 0.5 * n.dot(h_dl, n.dot(B, h_dl)); + + var acceptCandidate; + + if (dF == 0 || dL == 0) { + acceptCandidate = true; + } else { + var rho = dF / dL; + // update delta + if (rho < 0.25) { + // if the model is a poor predictor reduce the size of the trust region + delta *= 0.5; } else { - var rho = dF/dL; - // update delta - if( rho < 0.25 ) { - // if the model is a poor predictor reduce the size of the trust region - delta *= 0.5; - } else { - // only increase the size of the trust region if it is taking a step of maximum size - // otherwise just assume it's doing good enough job - if( rho > 0.75 && hitBoundary) { - var r = n.norm2(h_dl); - delta *= 2; - } - } - acceptCandidate = rho > 0; // could be 0 .. 0.25 + // only increase the size of the trust region if it is taking a step of maximum size + // otherwise just assume it's doing good enough job + if (rho > 0.75 && hitBoundary) { + var r = n.norm2(h_dl); + delta *= 2; + } } + acceptCandidate = rho > 0; // could be 0 .. 0.25 + } - if (acceptCandidate) { - x = n.clone(x_new); - Jx = n.clone(Jx_new); - fx = n.clone(fx_new); - err = err_new; + if (acceptCandidate) { + x = n.clone(x_new); + Jx = n.clone(Jx_new); + fx = n.clone(fx_new); + err = err_new; - gNg = n.dot(n.transpose(Jx), n.mul(fx, -1)); - g = n.dot(n.transpose(Jx), fx); + gNg = n.dot(n.transpose(Jx), n.mul(fx, -1)); + g = n.dot(n.transpose(Jx), fx); - // get infinity norms - g_inf = n.norminf(gNg); - fx_inf = n.norminf(fx); - } + // get infinity norms + g_inf = n.norminf(gNg); + fx_inf = n.norminf(fx); + } - - // count this iteration and start again - iter++; + // count this iteration and start again + iter++; } return { - evalCount : iter, - error : err, - returnCode : stop + evalCount: iter, + error: err, + returnCode: stop }; };