Reimplementing constraint solver

Discussion about the development of the Assembly workbench.
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by DeepSOIC » Wed Nov 27, 2019 4:11 pm

Testing if the scaling works. It does!
Running this

Code: Select all

sc = 1e4
import ConstraintSolver as CS
ps = CS.ParameterStore()
p1 = CS.G2D.ParaPoint(ps)
p2 = CS.G2D.ParaPoint(ps)
p2.x.Value = 3 * sc
p2.y.Value = 4 * sc
p1.update()
p2.update()
for p in p1.Parameters + p2.Parameters:
    p.OwnScale = sc

c = CS.G2D.ConstraintDistance(
    p1 = p1,
    p2 = p2,
    store = ps,
    Label = "Constraint1"
)
c.dist.Value = sc
c.Weight = 1.0/sc

c.update()
c.NetError

sys = CS.SubSystem()
sys.addUnknown(p1.Parameters)
sys.addUnknown(p2.Parameters)
sys.addConstraint(c)

vs = CS.ValueSet(CS.ParameterSubset(ps.allFree()))
slv = CS.SketchSolver()
slv.solveDogLeg(sys, vs)
, setting different values to sc doesn't change solver output at all.
Which means, the scaling is working as expected. So hopes are high, those faults related to meter and larger sized sketches can be finally solved. Now, fingers crossed, I ever get that far.
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by DeepSOIC » Wed Nov 27, 2019 5:32 pm

Now onto LM solver, I immediately stumbled on the reason it has accuracy problems. Its termination condition is:
e1^2 + e2^2 + .. < eps, where default eps = 1e-10
In dogleg, it is:
max(abs(e1), abs(e2), ...) < tolf , where tolf default is again 1e-10
Since DL is comparing linearly and LM is comparing squared, there is a 5 order of magnitude difference in accuracy.
toralf
Posts: 44
Joined: Fri May 03, 2019 3:54 am

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by toralf » Wed Nov 27, 2019 6:40 pm

Great progress. Sound interesting. Could you make use of the dual numbers?
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by DeepSOIC » Wed Nov 27, 2019 6:44 pm

toralf wrote:
Wed Nov 27, 2019 6:40 pm
Could you make use of the dual numbers?
:shock: :? I am, see first post of the thread :roll:
toralf
Posts: 44
Joined: Fri May 03, 2019 3:54 am

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by toralf » Wed Nov 27, 2019 9:23 pm

I know that you wrote it in the first post, but these were your goals. Did they prove to be useful for this reimplementation?
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by DeepSOIC » Wed Nov 27, 2019 10:07 pm

They already proved useful in sketcher constraints, there is just no question about it. In the implementation so far, they aren't very useful, because I have not written geometry and constraints yet. The two I wrote are too trivial to say "dual numbers were useful", but for more complex stuff, dual numbers are just a lifesaver IMO.
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by DeepSOIC » Thu Nov 28, 2019 12:27 am

LM solver is up!
Here's a log where it solves the same super-simple test problem:

Code: Select all

Begin Levenberg-Marquardt solving
  iteration 0
    err = 2.000000
    max(abs(gradient[i])) = 1.600000
    damping factor = 0.000640, computing step...
    sq(step) = 1.998721
    step changed the error by 4.000000 (expected: 4.000000)
    step accepted
    changing damping by factor 0.333300
  iteration 1
    err = 0.000640
    max(abs(gradient[i])) = 0.000512
    damping factor = 0.000213, computing step...
    sq(step) = 0.000000
    step changed the error by 0.000000 (expected: 0.000000)
    step accepted
    changing damping by factor 0.333300
  iteration 2
    err = 0.000000
    max(abs(gradient[i])) = 0.000000
    damping factor = 0.000071, computing step...
    sq(step) = 0.000000
    step changed the error by 0.000000 (expected: 0.000000)
    step accepted
    changing damping by factor 0.333300
  iteration 3
    err = 0.000000
error is small, solved
ickby
Posts: 2936
Joined: Wed Oct 05, 2011 7:36 am

Re: Reimplementing constraint solver [dogleg solver is up!]

Post by ickby » Thu Nov 28, 2019 6:19 am

Nice progress!
DeepSOIC wrote:
Wed Nov 27, 2019 11:56 am
It seems there is an error in the implementation.
double dL = err - 0.5*(fx + Jx*h_dl).squaredNorm();
double dF = err - err_new;
double rho = dL/dF;
Hm, back in the days when I implemented dogleg and LM solvers I closely followed a paper (or excerpt of a bock, not sure anymore). I will have a look this evening if I still can find it on my other machine. There is a good chance I switched some variables and made a error, as I experimented a lot with those solvers.
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [LM solver is up!]

Post by DeepSOIC » Thu Nov 28, 2019 11:17 am

ickby wrote:
Thu Nov 28, 2019 6:19 am
I closely followed a paper (or excerpt of a bock, not sure anymore). I will have a look this evening if I still can find it on my other machine.
That'd be very nice to have, as I'm struggling to find a good source for dogleg for example. LM is decently explained on Wikipedia (although I'm still struggling to wrap my head around that J^T * J thing).
If I just drop all that "understand the solver algo" part, I feel like I can copy-paste and translate the implementation to the new system a lot faster. But it's interesting :geek: , and potentially useful for the future.
User avatar
DeepSOIC
Posts: 7291
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Reimplementing constraint solver [LM solver is up!]

Post by DeepSOIC » Thu Nov 28, 2019 11:29 am

LM's damping factor changer function

Code: Select all

                double dF = e.squaredNorm() - e_new.squaredNorm();
                double dL = h.dot(mu*h+g);

                if (dF>0. && dL>0.) { // reduction in error, increment is accepted
                    double tmp=2*dF/dL-1.;
                    mu *= std::max(1./3., 1.-tmp*tmp*tmp); //   <- this multiplier is plotted
Attachments
LM damping changer func.png
LM damping changer func.png (14.27 KiB) Viewed 330 times
Post Reply

Who is online

Users browsing this forum: No registered users and 2 guests