ickby wrote:ahhhh I did not yet work on the dogleg solver as promised. Sorry for that, but currently coding time is very limited for me. But I keep it in mind.

It is ok. Not all in life is coding for FC

When you have time

ickby wrote:Are you sure that Dogleg is used for dragging? As far as I remember the SQP solver is used for that. The reason behind this is that for normal nonlinear equation solvers like dogleg/bfsg/LM there is no notion of weak constraints. Weak constraints are such that are nice to be fullfilled, but are only considered if all others are fullfilled. That is what is used for dragging: the "follow the mouse" is a weak condition. So it is fullfilled if it is possible to fullfill all user constraints, but if it is not possible the user constraints are fullfilled and the point does not follow the mouse .

Well, I am not even sure what is a SQP solver and where is it implementation anymore...

I explain myself:

Current version of the code
DL seems to be used for dragging. If you activate the solver Debug (and my solvers, normal and redundant, are LM currently) and drag a point you get in python console:

Code: Select all

```
>>> App.ActiveDocument.Sketch.movePoint(5,1,App.Vector(-9.916887,-24.583879,0),0)
```

and in the report window:

Code: Select all

```
Sketcher::Solve()-DogLeg-T:0.011
Sketcher::Solve()-DogLeg-T:0.007
Sketcher::Solve()-DogLeg-T:0.006
```

It may be that this is not what is happening and that SQP is indeed executed and the code wrongly reporting that Dogleg is being executed...

Yes, there is/was even a comment in the code of the solver saying it so (SQP for dragging):

Old version
This is the version of the solver before I touched it (commit 280019ddf9):

Code: Select all

```
for (int soltype=0; soltype < (isInitMove ? 1 : 4); soltype++) {
std::string solvername;
switch (soltype) {
case 0: // solving with the default DogLeg solver
// (or with SQP if we are in moving mode)
solvername = isInitMove ? "SQP" : "DogLeg";
ret = GCSsys.solve(isFine, GCS::DogLeg);
break;
case 1: // solving with the LevenbergMarquardt solver
solvername = "LevenbergMarquardt";
ret = GCSsys.solve(isFine, GCS::LevenbergMarquardt);
break;
case 2: // solving with the BFGS solver
solvername = "BFGS";
ret = GCSsys.solve(isFine, GCS::BFGS);
break;
case 3: // last resort: augment the system with a second subsystem and use the SQP solver
solvername = "SQP(augmented system)";
InitParameters.resize(Parameters.size());
int i=0;
for (std::vector<double*>::iterator it = Parameters.begin(); it != Parameters.end(); ++it, i++) {
InitParameters[i] = **it;
GCSsys.addConstraintEqual(*it, &InitParameters[i], -1);
}
GCSsys.initSolution();
ret = GCSsys.solve(isFine);
break;
}
...
}
```

so if we are moving for (int soltype=0; soltype < 1; soltype++) , only soltype==0 is executed which finely defines the solver name to SQP, but it does not care because it executes DogLeg anyway: ret = GCSsys.solve(isFine, GCS::DogLeg);

For sure we are not doing what is under "case 3" that is labeled as SQP. Even though that "case 3", after decoupling the two subsystems calls GCSsys.solve(bool, solver=GCS::DogLeg). So it calls DogLeg... (or not... continue reading...)

However if we center ourselves in case 0 and initmove==true and go to the GCSsys call, we end up executing the same solve call as in case 3 (without decoupling, which might be fine because in case 3 we are augmenting the system, whereas while moving we already have two subsystems with different priorities,... but I do not fully understand it...) as forcing to GCS::DogLeg or not is the same as that is the "default value". This function is:

Code: Select all

```
int System::solve(bool isFine, Algorithm alg)
{
if (!isInit)
return Failed;
bool isReset = false;
// return success by default in order to permit coincidence constraints to be applied
// even if no other system has to be solved
int res = Success;
for (int cid=0; cid < int(subSystems.size()); cid++) {
if ((subSystems[cid] || subSystemsAux[cid]) && !isReset) {
resetToReference();
isReset = true;
}
if (subSystems[cid] && subSystemsAux[cid])
res = std::max(res, solve(subSystems[cid], subSystemsAux[cid], isFine));
else if (subSystems[cid])
res = std::max(res, solve(subSystems[cid], isFine, alg));
else if (subSystemsAux[cid])
res = std::max(res, solve(subSystemsAux[cid], isFine, alg));
}
if (res == Success) {
for (std::set<Constraint *>::const_iterator constr=redundant.begin();
constr != redundant.end(); constr++){
//DeepSOIC: there used to be a comparison of signed error value to
//convergence, which makes no sense. Potentially I fixed bug, and
//chances are low I've broken anything.
double err = (*constr)->error();
if (err*err > XconvergenceFine) {
res = Converged;
return res;
}
}
}
return res;
}
```

All the options end up calling DogLeg, except for: "res = std::max(res, solve(subSystems[cid], subSystemsAux[cid], isFine));", which calls:

Code: Select all

```
// The following solver variant solves a system compound of two subsystems
// treating the first of them as of higher priority than the second
int System::solve(SubSystem *subsysA, SubSystem *subsysB, bool isFine)
{
int xsizeA = subsysA->pSize();
int xsizeB = subsysB->pSize();
int csizeA = subsysA->cSize();
VEC_pD plistAB(xsizeA+xsizeB);
{
VEC_pD plistA, plistB;
subsysA->getParamList(plistA);
subsysB->getParamList(plistB);
std::sort(plistA.begin(),plistA.end());
std::sort(plistB.begin(),plistB.end());
VEC_pD::const_iterator it;
it = std::set_union(plistA.begin(),plistA.end(),
plistB.begin(),plistB.end(),plistAB.begin());
plistAB.resize(it-plistAB.begin());
}
int xsize = plistAB.size();
Eigen::MatrixXd B = Eigen::MatrixXd::Identity(xsize, xsize);
Eigen::MatrixXd JA(csizeA, xsize);
Eigen::MatrixXd Y,Z;
Eigen::VectorXd resA(csizeA);
Eigen::VectorXd lambda(csizeA), lambda0(csizeA), lambdadir(csizeA);
Eigen::VectorXd x(xsize), x0(xsize), xdir(xsize), xdir1(xsize);
Eigen::VectorXd grad(xsize);
Eigen::VectorXd h(xsize);
Eigen::VectorXd y(xsize);
Eigen::VectorXd Bh(xsize);
// We assume that there are no common constraints in subsysA and subsysB
subsysA->redirectParams();
subsysB->redirectParams();
subsysB->getParams(plistAB,x);
subsysA->getParams(plistAB,x);
subsysB->setParams(plistAB,x); // just to ensure that A and B are synchronized
subsysB->calcGrad(plistAB,grad);
subsysA->calcJacobi(plistAB,JA);
subsysA->calcResidual(resA);
double convergence = isFine ? XconvergenceFine : XconvergenceRough;
int maxIterNumber = MaxIterations;
double divergingLim = 1e6*subsysA->error() + 1e12;
double mu = 0;
lambda.setZero();
for (int iter=1; iter < maxIterNumber; iter++) {
int status = qp_eq(B, grad, JA, resA, xdir, Y, Z);
if (status)
break;
x0 = x;
lambda0 = lambda;
lambda = Y.transpose() * (B * xdir + grad);
lambdadir = lambda - lambda0;
// line search
{
double eta=0.25;
double tau=0.5;
double rho=0.5;
double alpha=1;
alpha = std::min(alpha, subsysA->maxStep(plistAB,xdir));
// Eq. 18.32
// double mu = lambda.lpNorm<Eigen::Infinity>() + 0.01;
// Eq. 18.33
// double mu = grad.dot(xdir) / ( (1.-rho) * resA.lpNorm<1>());
// Eq. 18.36
mu = std::max(mu,
(grad.dot(xdir) + std::max(0., 0.5*xdir.dot(B*xdir))) /
( (1. - rho) * resA.lpNorm<1>() ) );
// Eq. 18.27
double f0 = subsysB->error() + mu * resA.lpNorm<1>();
// Eq. 18.29
double deriv = grad.dot(xdir) - mu * resA.lpNorm<1>();
x = x0 + alpha * xdir;
subsysA->setParams(plistAB,x);
subsysB->setParams(plistAB,x);
subsysA->calcResidual(resA);
double f = subsysB->error() + mu * resA.lpNorm<1>();
// line search, Eq. 18.28
bool first = true;
while (f > f0 + eta * alpha * deriv) {
if (first) { // try a second order step
// xdir1 = JA.jacobiSvd(Eigen::ComputeThinU |
// Eigen::ComputeThinV).solve(-resA);
xdir1 = -Y*resA;
x += xdir1; // = x0 + alpha * xdir + xdir1
subsysA->setParams(plistAB,x);
subsysB->setParams(plistAB,x);
subsysA->calcResidual(resA);
f = subsysB->error() + mu * resA.lpNorm<1>();
if (f < f0 + eta * alpha * deriv)
break;
}
alpha = tau * alpha;
if (alpha < 1e-8) // let the linesearch fail
alpha = 0.;
x = x0 + alpha * xdir;
subsysA->setParams(plistAB,x);
subsysB->setParams(plistAB,x);
subsysA->calcResidual(resA);
f = subsysB->error() + mu * resA.lpNorm<1>();
if (alpha < 1e-8) // let the linesearch fail
break;
}
lambda = lambda0 + alpha * lambdadir;
}
h = x - x0;
y = grad - JA.transpose() * lambda;
{
subsysB->calcGrad(plistAB,grad);
subsysA->calcJacobi(plistAB,JA);
subsysA->calcResidual(resA);
}
y = grad - JA.transpose() * lambda - y; // Eq. 18.13
if (iter > 1) {
double yTh = y.dot(h);
if (yTh != 0) {
Bh = B * h;
//Now calculate the BFGS update on B
B += 1./yTh * y * y.transpose();
B -= 1./h.dot(Bh) * (Bh * Bh.transpose());
}
}
double err = subsysA->error();
if (h.norm() <= convergence && err <= smallF)
break;
if (err > divergingLim || err != err) // check for diverging and NaN
break;
}
int ret;
if (subsysA->error() <= smallF)
ret = Success;
else if (h.norm() <= convergence)
ret = Converged;
else
ret = Failed;
subsysA->revertParams();
subsysB->revertParams();
return ret;
}
```

so maybe this last function is SQP, and the previous code is counter-intuitive showing that we are calling DogLeg when we are clearly not if there are two subsystems (if the is one DogLeg is called)...

Current version
The "end result" is the same in the current code, as the same function is called when moving (I just do not "lie" to the user saying it is SQP, or maybe it is SQP and that one is a bug and solvername should be "SQP", I mean the "notification" to the user, the actual functional code is the same):

Code: Select all

```
if(isInitMove){
solvername = "DogLeg"; // DogLeg is used for dragging (same as before)
ret = GCSsys.solve(isFine, GCS::DogLeg);
}
else{
switch (defaultSolver) {
case 0:
solvername = "BFGS";
ret = GCSsys.solve(isFine, GCS::BFGS);
defaultsoltype=2;
break;
case 1: // solving with the LevenbergMarquardt solver
solvername = "LevenbergMarquardt";
ret = GCSsys.solve(isFine, GCS::LevenbergMarquardt);
defaultsoltype=1;
break;
case 2: // solving with the BFGS solver
solvername = "DogLeg";
ret = GCSsys.solve(isFine, GCS::DogLeg);
defaultsoltype=0;
break;
}
}
```

Conclusions:

1. The code

is very obscure on what algorithm is being executed when moving and should be improved even if it is only for readibility. At the very least a comment should be introduced to clarify what is being executed.

2. There is one clear observation. I have a sketch where I can work perfectly with LM (add new lines and constraints) and I can not with DL. I can not move any point in this sketch by dragging:

http://www.freecadweb.org/tracker/view.php?id=2172
I do not know myself if what is executed currently when dragging is what it should be executed (is SQP activated?) or not. Or if SQP is executed (although it seems that DL was executed) but SQP also fails (but LM would not fail when not having priority constraints or weak constraints, that I assume is the implementation of weak constraints and what gives rise to the multiple subsystems with different priority).

I will probably take a closer look when doing the things said above. When you have time, I guess you could give a look to this "dragging thing solver" as you are the algorithm expert.