mirror of
https://github.com/xibyte/jsketcher
synced 2025-12-11 02:45:52 +01:00
format
This commit is contained in:
Val Erastov
parent
b977dc96dc
commit
573fc98684
1 changed files with 125 additions and 127 deletions
|
|
@ -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
|
||||
};
|
||||
|
||||
};
|
||||
|
|
|
|||
Loading…
Reference in a new issue