[Discussion, user input] Sketcher stability improvements

Here's the place for discussion related to coding in FreeCAD, C++ or Python. Design, interfaces and structures.
abdullah
Posts: 3174
Joined: Sun May 04, 2014 3:16 pm

[Discussion, user input] Sketcher stability improvements

Postby abdullah » Thu Aug 20, 2015 8:47 pm

Hi folks,

My next FC task will be to dedicate some time to improve the stability of the Sketcher.

I have identified some problems (including sometimes crashes) when the solver fails. Additionally I would like to make LM solver be used (instead of/in addition to/as a fallback of) DL for dragging. Currently only DL is used for dragging, so if you arrived to an Sketcher where DogLeg fails, you won't be able to drag a single point (but you should be able to drag it using LM if LM succeeds when DL failed).

If any of you has identified/identifies a wrong behaviour of the Sketcher leading to a crash of FC or wrong behaviour/representation as a consequence of a solver failure, please let me know how to reproduce it in this thread (if you filed a Mantis ticket you can just paste the link).
User avatar
DeepSOIC
Posts: 6639
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: [Discussion, user input] Sketcher stability improvements

Postby DeepSOIC » Thu Aug 20, 2015 8:57 pm

Hi, Abdullah! I've been using your new sketcher for a little bit already, and I find it fantastic. I haven't got even a single bug so far, I think :D
abdullah
Posts: 3174
Joined: Sun May 04, 2014 3:16 pm

Re: [Discussion, user input] Sketcher stability improvements

Postby abdullah » Thu Aug 20, 2015 9:05 pm

DeepSOIC wrote:Hi, Abdullah! I've been using your new sketcher for a little bit already, and I find it fantastic. I haven't got even a single bug so far, I think :D
Thanks!! That is probably because you use it in a smart way and you do not get into troubles. I am quite "unexperienced" to put it mildly and manage to push its current boundaries from time to time... ;)
ickby
Posts: 2903
Joined: Wed Oct 05, 2011 7:36 am

Re: [Discussion, user input] Sketcher stability improvements

Postby ickby » Fri Aug 21, 2015 6:04 am

Hello,

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.

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 .
abdullah
Posts: 3174
Joined: Sun May 04, 2014 3:16 pm

Re: [Discussion, user input] Sketcher stability improvements

Postby abdullah » Fri Aug 21, 2015 12:59 pm

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 :lol:

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... :lol: 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. :)
triplus
Posts: 8415
Joined: Mon Dec 12, 2011 4:45 pm

Re: [Discussion, user input] Sketcher stability improvements

Postby triplus » Sat Aug 22, 2015 12:19 am

abdullah wrote:If any of you has identified/identifies a wrong behaviour of the Sketcher leading to a crash of FC or wrong behaviour/representation as a consequence of a solver failure, please let me know how to reproduce it in this thread (if you filed a Mantis ticket you can just paste the link).
I am not sure if it is related to solver failure but sure why not. Here is one crash:

issue #0002225
Jee-Bee
Posts: 1911
Joined: Tue Jun 16, 2015 10:32 am
Location: Netherlands

Re: [Discussion, user input] Sketcher stability improvements

Postby Jee-Bee » Sat Aug 22, 2015 8:18 am

triplus wrote:issue #0002225
On my mac it doesn't crash i just can't remove the last external reference...
triplus
Posts: 8415
Joined: Mon Dec 12, 2011 4:45 pm

Re: [Discussion, user input] Sketcher stability improvements

Postby triplus » Sat Aug 22, 2015 12:28 pm

Jee-Bee wrote:
triplus wrote:issue #0002225
On my mac it doesn't crash i just can't remove the last external reference...
Yes you are correct as that is one more of the possible symptoms i get some times (not able to delete some external links).
abdullah
Posts: 3174
Joined: Sun May 04, 2014 3:16 pm

Re: [Discussion, user input] Sketcher stability improvements

Postby abdullah » Sat Aug 22, 2015 5:35 pm

I have assigned the bug to myself. I have now some bugs assigned to myself, so somewhere in the near future I will probably try to handle those tickets. ;)
triplus
Posts: 8415
Joined: Mon Dec 12, 2011 4:45 pm

Re: [Discussion, user input] Sketcher stability improvements

Postby triplus » Sat Aug 22, 2015 11:29 pm

abdullah wrote:I have assigned the bug to myself. I have now some bugs assigned to myself, so somewhere in the near future I will probably try to handle those tickets. ;)
Great to hear that and its nice to see that level of commitment.