Reimplementing constraint solver [all solvers implemented]

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

Reimplementing constraint solver [all solvers implemented]

Postby DeepSOIC » Thu Oct 31, 2019 9:10 pm

Dev branch: https://github.com/DeepSOIC/FreeCAD-ellipse (branch ConstraintSolver1, I made it default branch of the repo)

Gitter chat: https://gitter.im/FreeCAD-ellipse/community

So, as I wanted to fix sketcher solver, and scaling constraint error functions didn't give good results, I now want to scale the dofs themselves. It seems doable within the current situation. However, since I am thinking too hard about the solver, I decided to try initiate the effort of rewriting the solver.
The primary goal is to have a standalone solver decoupled from Sketcher, ready for 3d assembly solving, and with full python bindings.

So, here we go.

Design goals (i.e., improvements over current solver).
* full python support (one should be able to implement a 2d/3d sketcher fully in python using the solver)
* python-powered constraint support (with inevitable speed problems, of course)
* externally extendable (i.e. no hard-coding of constraint types)
* reentrant solver.
That is, allow create multiple solvers and run them in parallel. That implies no writing to ParameterStore while solving.
* built-in parameter scaling factor support, to eliminate the size-of-sketch magnitude problem
* use DualNumbers in constraint-related math to implicitly compute derivatives
(was DeriVectors, pretty much the same stuff)
* don't evaluate constraint derivatives for parameters the constraint does not depend on, for speed, yet avoid branching.
(by using list of parameters returned by constraints, and stored indexes in lookup tables, we can selectively fill matrix elements with very little overhead)
* support rigid groups of geometry

This is a rough architecture I came up with to fit the goals. I'm open for discussion.

Code: Select all

class Parameter
//parameter should always live in ParameterStore
//Py: nothing, accessed via ParameterRef
    HParameterStore host
    double value
    double magnitude = 1.0
    int par_index
    int masterindex //==par_index, unless the parameter was redirected
    std::string name

class ParameterRef
//Py: copy-construct
    HParameterStore host
    int par_index
    Parameter& operator*

class ParameterStore
    //py: hosted by PyObject, no constructor on stack
    vector<Parameter> parameters
    vector<int> redirectionmap (key = index of parameter, value = index of master parameter)
    PyObject* self

class HParameterStore
//handle to ParameterStore
    //Py: nothing, it's a C++ wrapper around python object
    Py::Object store_py
    ParameterStore* store
    ParameterStore& operator*

class ParameterSet
//list of parameters. e.g., list of unknowns for solver
    //py: copy-construct?
    smart_ptr vector<ParameterRef> parameters
    smart_ptr vector<int> lut //key = parameter index. value = index into the set
    operator[] (retrieves the parameter)
    PyObject* pyobj
    
class ParameterValueSet(ParameterSet)
    //py: construct around one on stack, with crash potential
    Eigen::vectorxD values
    Eigen::vector<double> duals //dual parts, for derivative computation
    operator[] (retrieves the alternate value, as values[lut[arg->masterindex]] if in set else arg->value())
    
    
class Solver:
    //py: hosted by PyObject, no constructor on stack
    ParameterSet unknowns
    Vector<HConstraint> constraints
    ...and heaps of other stuff
    PyObject* self

class HSolver:
    Py::Object solver_py
    Solver* solver
    Solver& operator*

class Constraint
    //py: hosted by PyObject, no constructor on stack
    PyObject* self
    int tag
  methods:
    vector<ParameterRef> params() //all params the constraint is using
    DualNumber error(ParameterValueSet &on)
    double maxStep(ParameterValueSet &on, ParameterValueSet &dir, double curstep)

class HConstraint:
    Py::Object constraint_py
    Constraint* solver
    Constraint& operator*




class ParaObject
//base class for all objects built on parameters
    //py: copy-construct
    ParameterSet pvec //to be filled in by constructors
    static makeWithParams(ParameterStore params)

2D:

class DualVector2

class ParaVector2(ParaObject)

class ParaGeom2(ParaObject)
    class ParaPoint2(ParaGeom2)
    class ParaEdge2(ParaGeom2)
        class ParaFullEdge2(ParaEdge)
        class ParaTrimmedEdge2(ParaEdge)
    class PlacedGeometry(ParaGeom2) //to support rigid blocks, think assembly
        class PlacedPoint(PlacedGeometry)
        class PlacedEdge(PlacedGeometry)
    class ParaCurve2(ParaGeom2)
    class ParaLine2(ParaCurve2)
    class ParaCircle2(ParaCurve2)
    class ParaEllipse2(ParaCurve2)
    ...

3D:

class DualVector3
class DualRotation
class DualPlacement

class ParaVector3(ParaObject)
class ParaRotation(ParaObject)
class ParaPlacement(ParaObject)

class ParaGeom3(ParaObject)
    class PlacedShape(ParaGeom3)
    class ParaEdge3(ParaGeom3)
        class ParaFullEdge3(ParaEdge3)
        class ParaTrimmedEdge3(ParaEdge3)
    class ParaCurve3(ParaGeom3)
    class ParaLine3(ParaCurve3)
    class ParaCircle3(ParaCurve3)
    class ParaEllipse3(ParaCurve3)
    ...
    class ParaSurface3(ParaGeom3)
        class ParaPlane3(ParaGeom3)
        class ParaSphere3(ParaGeom3)
        ...


Notes:
DualVector, DualRotation, ParaVector2,3, ParaRotation, ParaPlacement etc should support assignment from and to their Base::<> twins
Also nice to have:
* typesystem integration
* add some storage facilities in geometry, constraints and parameters. For example, add a Py::Object data member. And probably it's good to have something for C++ too.
* (super optional) interactive solving support (i.e. manual iteration advancement)
* template constraint implementations for coordinates (e.g., CurveValueX and CurveValueY constraints share the implementation, and just return differences in respective coordinates as error their function)


Tasks (same number = can be done in parallel)
1. Set up new FreeCAD module
2. implement Parameter, ParameterRef, ParameterStore, HParameterStore, ParameterSet, ParameterValueSet
2. implement DualNumber math functions, DualVector2
3. implement ParaVector2,
4. Implement ParaGeom2 and derived, and constraints
4. reimplement solver, taking care to account for parameter scaling
5. Test the solver with python
6. rewrite Sketch.cpp to use the new solver
7. some good testing!
------------
8. implement 3d placement, geometry and constraints
9. use in assembly workbench

Who is brave enough to join this? I feel like I can't lift it all on my own.
Last edited by DeepSOIC on Mon Dec 02, 2019 11:13 pm, edited 10 times in total.
Jee-Bee
Posts: 1966
Joined: Tue Jun 16, 2015 10:32 am
Location: Netherlands

Re: Call for team! Reimplementing constraint solver

Postby Jee-Bee » Fri Nov 01, 2019 9:48 am

Without wanting to stop this kind of efforts(i think it is great)!!

In my opinion the sketch solver and assembly solvers are two different solvers.
even within assembly solvers there are multiple methods.

I'm normaly a (PTC) Creo user, but last half year i need to work with Solid Works... The main issue i had with Solid works is the stabiliy of the solver. I think i figgured out why the solver of SW is kind of unstable. I think the solid works solver works this way:

Code: Select all

[X   O...  ]
[ X   ... O]
[  X  ...  ]
[O  X ...  ]
....O..X....
[  O  ... X]
Basically it create a big matrix with ALL constraints(the O's). if one constraint is broken or needs to be flipped the solver don't know which one is or broken or which one to flip.
Creo have i think a itarative solver(master-slave) so it runs verry fast through a loop... more or less this way:

Code: Select all

[...[[X11]22
     [1X1]22
     [11X]22]
     [   2X2]
     [    2X]
[              ...]
What i meant with this this one the number one is first iteration (part/ assy 1 on main coordinate system/ planes /whatever). the number 2 is next iteration of second part on main planes/ coordinate system + previous part, next part is on all previos part and so on...
When now now a constraint have to be swapped. only the new part can flip.
if the solver loop through the hole tree it can exactly point to ALL failed constraints!
Yes i beleve the creo method is slower but much less error prone and stability is much better!

what the best method is for handling a sketch i realy can't say anything about it so i don't exactly know if both solvers are based on same pricinciples...
kbwbe
Posts: 881
Joined: Tue Apr 10, 2018 3:12 pm
Location: Germany, near Köln (Cologne)

Re: Call for team! Reimplementing constraint solver

Postby kbwbe » Sat Nov 02, 2019 1:15 pm

DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
The primary goal is to have a standalone solver decoupled from Sketcher, ready for 3d assembly solving, and with full python bindings.
Hi @DeepSOIC,
thank you for taking the initative and doing this call. My special respect to you...

I am also not quite sure, whether a new stand alone solver will be able to handle 2/3D sketches and 3D assemblies. Sketches usually are working with edges and points, whereas assemblies additionaly have to deal with faces for example. So the required data structures could get very complicated.
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* full python support (one should be able to implement a 2d/3d sketcher fully in python using the solver)
* python-powered constraint support (with inevitable speed problems, of course)
ok
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* externally extendable (i.e. no hard-coding of constraint types)
How do you want to do this ? Whats your idea behind ?
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* reentrant solver.
ok
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
That is, allow create multiple solvers and run them in parallel. That implies no writing to ParameterStore while solving.
Just speaking for assembly side 3D solving. Is that possible ? I could imagine this works with one solver per subassembly, but not for one assembly level.
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* built-in parameter scaling factor support, to eliminate the size-of-sketch magnitude problem
ok, i saw your other thread. For assemblies the problem will remain if we are assembling very small parts to very big ones, i fear.
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* use DualNumbers in constraint-related math to implicitly compute derivatives
(was DeriVectors, pretty much the same stuff)
I do not really know the math behind dual numbers. If there are for example 3 different constraints necessary to fix to2 parts together: Is it possible to calculate a placement operation for each constraint and then calculate an resulting average of the calculated placements? (That's the way A2plus solver is working)
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* don't evaluate constraint derivatives for parameters the constraint does not depend on, for speed, yet avoid branching.
(by using list of parameters returned by constraints, and stored indexes in lookup tables, we can selectively fill matrix elements with very little overhead)
I do not understand this ATM, (my personal lack of math perhaps). Assembly side again: A2plus is using some structure and dof analysis for solving a reduced number of constraints.
DeepSOIC wrote:
Thu Oct 31, 2019 9:10 pm
* support rigid groups of geometry
For the assembly side, that is clear.


Sorry if my comments mostly touch the assembly side. (But your post is here in the assembly forum, so i hope this is ok). Up to now i have been very happy with the sketcher and its GCS solver. Is it really necessary to redo this ?

Some remarks to the A2plus 3D solver:
- It's solving algoithm seems to be mostly very stable (except perhaps some small bugs)
- ATM it is not a standalone solver and it is missing a clear external interface. It is deeply integrated into A2plus, what should be changed.
- It is completely written in Python and therefore slow. I guess, reimplementing it in C could give a speed increase of factor 50-100. This could be enough for many applications.
- It does not make use of LCS constraints. They should be implemented. Assembly4 shows the potential of them.

So my question is: Isn't it better to have a solver for 2/3D sketches and a different one for assemblies ?
KBWBE

https://github.com/kbwbe/A2plus
latest release: v0.4.30, installable via FreeCAD's addon manager
Tutorials:
Paul Randall: https://youtu.be/mnkecA9S7kc
anisim (deutsch): https://www.youtube.com/watch?v=vDcaFq6IEJM
User avatar
DeepSOIC
Posts: 7145
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Call for team! Reimplementing constraint solver

Postby DeepSOIC » Sat Nov 02, 2019 2:57 pm

kbwbe wrote:
Sat Nov 02, 2019 1:15 pm
I am also not quite sure, whether a new stand alone solver will be able to handle 2/3D sketches and 3D assemblies. Sketches usually are working with edges and points, whereas assemblies additionaly have to deal with faces for example. So the required data structures could get very complicated.
I'm not sure either. That's what I want to be achieved, ultimately.
kbwbe wrote:
Sat Nov 02, 2019 1:15 pm
* externally extendable (i.e. no hard-coding of constraint types)
How do you want to do this ? Whats your idea behind ?
Well, actually, the current sketcher solver is extendable as well. The main thing I'd like to do is to make it easy, and that's what this ParaGeometry and DualNumbers are all about. Constraint is essentially just an error function. Currently, partial derivatives are programmed into an independent function. And with all that vector and placement math, and lots of involved parameters, programming all these derivatives quickly becomes impractical and difficult to debug and fix. Reimplementing all the math with dual numbers allows to calculate derivatives implicitly. Some constraints are already using dual numbers a little bit. I want to step up this game.

kbwbe wrote:
Sat Nov 02, 2019 1:15 pm
Is it possible to calculate a placement operation for each constraint and then calculate an resulting average of the calculated placements? (That's the way A2plus solver is working)
No. That is not how sketcher solver is working. Sketcher solver is a completely algebraic construct. Its job is to figure out a values for parameters, given a set of error functions, by zeroing out the error functions. It has no idea what is a placement, what is a vector, and so on.

It is possible to launch a solver for a simplified system of just one constraint, for each constraint. With this proposed design, it should be possible to run as many solvers as you want in threads, around the same set of parameters and constraints.
kbwbe wrote:
Sat Nov 02, 2019 1:15 pm
ok, i saw your other thread. For assemblies the problem will remain if we are assembling very small parts to very big ones, i fear.
Probably, yes. maybe not. This is one more thing I'm interested to find out. Requires implementing the thing.

kbwbe wrote:
Sat Nov 02, 2019 1:15 pm
* don't evaluate constraint derivatives for parameters the constraint does not depend on, for speed, yet avoid branching.
(by using list of parameters returned by constraints, and stored indexes in lookup tables, we can selectively fill matrix elements with very little overhead)

I do not understand this ATM,
No problem, it's a thing of current implementation: grad() is called for every parameter in the system. Then a constraint checks if the parameter is in its dependency list. That is a little inefficient.

kbwbe wrote:
Sat Nov 02, 2019 1:15 pm
Sorry if my comments mostly touch the assembly side. (But your post is here in the assembly forum, so i hope this is ok). Up to now i have been very happy with the sketcher and its GCS solver. Is it really necessary to redo this ?
To make a good python api for it, yes, I think it is necessary. I think that the solver can be useful for a wide variety of problems. Making specialized gearbox assembly workbench (2d or 3d), for example.
kbwbe
Posts: 881
Joined: Tue Apr 10, 2018 3:12 pm
Location: Germany, near Köln (Cologne)

Re: Call for team! Reimplementing constraint solver

Postby kbwbe » Sat Nov 02, 2019 4:20 pm

.
Thank you for all the information.
DeepSOIC wrote:
Sat Nov 02, 2019 2:57 pm
That is not how sketcher solver is working. Sketcher solver is a completely algebraic construct. Its job is to figure out a values for parameters, given a set of error functions, by zeroing out the error functions. It has no idea what is a placement, what is a vector, and so on.
I understand. (Also it is quite the opposite of what my A2plus solver is doing)
DeepSOIC wrote:
Sat Nov 02, 2019 2:57 pm
It is possible to launch a solver for a simplified system of just one constraint, for each constraint. With this proposed design, it should be possible to run as many solvers as you want in threads, around the same set of parameters and constraints.
Sounds very promising. Perhaps then it is possible to keep a solver class less complex. But how are all the results of the multiple solvers joined together, later? Please excuse if i am asking something stupid.

P.S.
If you are thinking about some iteration. If the job per solver is very small, do you think it is worth of running so many threads with the introduced overhead ?
Last edited by kbwbe on Sat Nov 02, 2019 4:42 pm, edited 1 time in total.
KBWBE

https://github.com/kbwbe/A2plus
latest release: v0.4.30, installable via FreeCAD's addon manager
Tutorials:
Paul Randall: https://youtu.be/mnkecA9S7kc
anisim (deutsch): https://www.youtube.com/watch?v=vDcaFq6IEJM
User avatar
DeepSOIC
Posts: 7145
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Call for team! Reimplementing constraint solver

Postby DeepSOIC » Sat Nov 02, 2019 4:40 pm

kbwbe wrote:
Sat Nov 02, 2019 4:20 pm
DeepSOIC wrote:
Sat Nov 02, 2019 2:57 pm
It is possible to launch a solver for a simplified system of just one constraint, for each constraint. With this proposed design, it should be possible to run as many solvers as you want in threads, around the same set of parameters and constraints.
Sounds very promising. Perhaps then it is possible to keep a solver class less complex. But how are all the results of the multiple solvers joined together, later? Please excuse if i am asking something stupid.
Well, I don't think they can be joined, if the systems aren't decoupled. So, while it is possible, the usefullness is questionable. I assumed, that is what you wanted when you asked "Is it possible to calculate a placement operation for each constraint and then calculate an resulting average of the calculated placements?".


Sketcher solver is capable of splitting the whole problem into subproblems in some situations. Then, the subproblems are solved independently. That can be paralleled (I think, that is even quite doable with the current architecture). With my proposed architecture, it should also be possible to run different solver algorithms in parallel (i.e., launching Dogleg, BFGS and LM in parallel, and picking the fastest or the best solution). It may also be possible to split up gradient matrix computation, although the overhead of creating threads may instead slow it down for reasonable-sized problems, and speed up only very large problems.

For sketches, the current solver has good-enough speed, and there is no need for multithreading. For 3d sketches, that may be a totally different story, as dofs pile up much faster.
kbwbe
Posts: 881
Joined: Tue Apr 10, 2018 3:12 pm
Location: Germany, near Köln (Cologne)

Re: Call for team! Reimplementing constraint solver

Postby kbwbe » Sat Nov 02, 2019 5:13 pm

.
Crosspost...
DeepSOIC wrote:
Sat Nov 02, 2019 4:40 pm
Well, I don't think they can be joined, if the systems aren't decoupled. So, while it is possible, the usefullness is questionable. I assumed, that is what you wanted when you asked "Is it possible to calculate a placement operation for each constraint and then calculate an resulting average of the calculated placements?".
Yes, that's more correct. If having several placement operations, averaging is possible for sure.
DeepSOIC wrote:
Sat Nov 02, 2019 4:40 pm
Sketcher solver is capable of splitting the whole problem into subproblems in some situations. Then, the subproblems are solved independently.
ok. Recognition of these subproblems in 3D is a complex job for it's own. We had to fight with that long time at A2plus. (DOF analysis...). And it was not done for multiple solver instances.
DeepSOIC wrote:
Sat Nov 02, 2019 4:40 pm
.although the overhead of creating threads may instead slow it down for reasonable-sized problems, and speed up only very large problems.
This answers the question of my P.S.
KBWBE

https://github.com/kbwbe/A2plus
latest release: v0.4.30, installable via FreeCAD's addon manager
Tutorials:
Paul Randall: https://youtu.be/mnkecA9S7kc
anisim (deutsch): https://www.youtube.com/watch?v=vDcaFq6IEJM
User avatar
DeepSOIC
Posts: 7145
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Call for team! Reimplementing constraint solver

Postby DeepSOIC » Mon Nov 04, 2019 10:03 pm

Started the development, although I am still a one-man-band.
https://github.com/DeepSOIC/FreeCAD-ellipse (default branch, ConstraintSolver1, is the development branch for now)

I have not settled the basics, so it may not be ready for contributions just yet. Particularly, I'm not too happy about how PyHandle turned out (see also https://forum.freecadweb.org/viewtopic.php?f=10&t=40629)
Mark Szlazak
Posts: 407
Joined: Tue Apr 04, 2017 6:06 pm
Location: Edmonton, Canada

Re: Call for team! Reimplementing constraint solver

Postby Mark Szlazak » Tue Nov 05, 2019 9:37 pm

In case you missed it, Ben Kenwright, the author of the Dual Quaternion papers people have posted, has authored a paper on programming rigid body constraint and physics solvers. Maybe the approaches are similar for constraints solvers alone or they could be made compatible with a future solvers that do physics simulation. Here is the link:

https://pdfs.semanticscholar.org/e2db/ ... 276005.pdf
User avatar
DeepSOIC
Posts: 7145
Joined: Fri Aug 29, 2014 12:45 am
Location: Saint-Petersburg, Russia

Re: Call for team! Reimplementing constraint solver

Postby DeepSOIC » Tue Nov 05, 2019 11:44 pm

Mark Szlazak wrote:
Tue Nov 05, 2019 9:37 pm
Thanks, I'm not sure if it's useful for me right now... since I am just reimplementing our solver, not really modifying the algorithm yet. But that may happen too some day.