Move Constraint solver into a separate crate and put it onto a solid mathematical foundation

[TobiasJacob]

Current approach

packages/cadmium/src/sketch/constraints.rs is currently one file containing the spring based constraint solver.

If I understand correctly you apply a virtual spring for every constraint that has a stiffness $K_p$ and a damping factor $K_d$. I assume that you are probably doing a full dynamic simulation (basically solving a second order differential equation, namely calculating accelerations $\ddot{q}_t$, then velocities $\dot{q}_t = \dot{q} _{t-1} + \ddot{q}_t \cdot dt$ and then positions $q_t=q _{t-1} + \dot{q} \cdot dt$) then, even though I didnt find this part in the code yet.

I think the idea is good in general and not to far away from how professional CAD systems do it. Modern PCs are powerful enough, and most sketches are simple enough that this will likely be the only approach we ever need. The proof of convergence for a small enough $dt$ is also tirival, because it is just a simulation of a mechanical system that due to the damping will always loose energy at every point in time, until it stops moving.

Limitations of current approach

Only the numerical stabililty is something that I worry about. Coming from robotics I know that these rigid body simulations tend to become unstable quickly and usually require very small $dt$ and manual tuning of the damping factors to get working nicely.

Mathematics behind current approach

Basically, what you do is setup a Loss function (or call it an error function) which is the mechanical energy of the system $L(q)$. $q$ are the parameters of all your geometry. The mechanical energy of the system is the energy stored in each spring $L(q) = \sum \frac{1}{2} D s^2$. On top of that is the kinetic energy of the moving particles and the damping factors, which lead to a more complex optimization problem touching Lagrange Mechanics.

But in the end, once the simulation converged, it will also have minimizied the energy stored in each spring, so $L(q) = \sum \frac{1}{2} D s^2$ is the only term we actually care about.

Removing the second order of the system

We can make our lives a lot easier here. First, instead of making it a second order dynamic simulation, we can remove the damping terms and make it a first order dynamic simulation by doing Gradient Descent:

$$ \min L(q) $$

$$ q_t = q _{t - 1} - \nabla L(q) \cdot dt $$

This will get rid of the problems with the damping factors and converge into the same solution. It will probably be enough for this program already.

Converging even faster and more numerically stable by using convex optimization algorithms

$L(q)$ is a convex quadratic function (since $\frac{1}{2} D s^2$ is convex). This is great, because it allows us to use a whole new suite of faster optimization algorithms specifically targeted at qudratic and convex functions. The most interesting one is the conjugate gradient methods, or methods that require the second order derivative as well: Newtons Method.

Making our lifes even easier by using a finished framework or building our own crate for it

Implementing a conjuage gradient algorithm is fairly straight forward (<300 lines of code), so defenitly a valid approach for us. Also, the algorithm is simple enough that we dont really have to intruduce an additional depenceny for it. However, we shouldn't reinvent the wheel. If we have $L(q)$, we can also just pass it directly into a general nonlinear or convex optimization framework, and have the mathematicans deals with the details such as preconditioning. E.g.

• Gomez
• Clarabel

Figuring out which constraints are in conflict with each other

This is relativly simple. If constraints are coherent, they will end up finding a configuration where each spring is relaxed, or the energy stored in each spring is $\frac{1}{2} D s^2 = 0$.

If two or more springs are conflicting with each other, they end up under stress $\frac{1}{2} D s^2 > 0$. The user interface can highlight these constraints.

Call to action

I am more than happy to help with this project. I work at Neura Robotics, and we the whole robotics industry is in desperate need of a Blender like free capable CAD program, e.g. Fusion 360 but free. I started working on my own CAD program in Rust geop, but there is no way for a single person to achieve this. This is why I am very exited to join forces on this project. I can probably contribute this part of the project, after we discussed some of the technical details. You can reach out to me on LinkedIn.

[MattFerraro]

Thank you Tobias, what a great writeup!

I have no attachment to the second-order nature of the current approach--it was just easy to think about during prototyping. I think the first order formulation is likely to be much better!

I strongly support the idea of breaking out the 2D constraint solver into its own standalone crate and giving it as few dependencies as possible. I also love the idea of trying out a few different approaches to the gradient descent step so we can weigh the pros and cons.

If you feel confident to lead this effort, I will support you however I can!

[TobiasJacob]

Hi @MattFerraro ,

Thanks, great to hear. I think I can start working on refactoring the code to use gradient descent. I can also document the maths behind it in a readme. This shouldnt be too much of a change.

Refactoring the code into a separate crate will probably be a big pullrequest leading to merge conflicts with other braches, so I'm not sure how to deal with this.

You wrote on Discord that you are working on settigung up a Github Organization anyways, which would make it easier to move the code around.

I guess suddenly coordinating 800+ people puts a bit of unexpected stress on you 😅 so yeah, just let me know what you think a good way forward

[MattFerraro]

I think it would be best to do a completely clean break here rather than keeping it inside the CADmium repo, so I've created https://github.com/CADmium-Co/ISOtope

I picked "ISOtope" because it's an Iterative SOlver and because it sorta fits the naming scheme of elemental things. Also "tope" like topography, and how it's all based on gradient descent? IDK, we'll workshop it.

Don't worry about matching the existing API too much--it was whipped up super fast as a proof of concept and it may well be bad. Build the right API for a standalone 2D Sketch Solver and CADmium will adapt to meet it.

[TobiasJacob]

I was able to make some progress in the repo https://github.com/CADmium-Co/ISOtope.

There is still room for improvement, but I think we can start working on integrating it into the software before the codebases diverge too far.

[MattFerraro]

This is fantastic work, Tobias! Unbelievable job

I pulled it locally and ran cargo test. On my m1 mac I'm seeing:

point_a: Point2 { data: [[0.5700834586643445, -0.06745530512562747]], gradient: [[-25.709630952671287, -3.5250742668599533]] }
point_b: Point2 { data: [[-1.3913236621093832, -0.33256804099923265]], gradient: [[-0.03432845891903502, -0.13300160678809736]] }
point_c: Point2 { data: [[-1.7858823110455389, 2.616758693265621]], gradient: [[0.04124606881310794, -0.04022156009781125]] }
point_d: Point2 { data: [[0.17054457073731163, 2.8818523972647134]], gradient: [[0.02043760840499649, -0.022759707622775023]] }
point_reference: Point2 { data: [[0.5592929825521086, -0.06193753782611579]], gradient: [[25.811662378731278, 3.5916709580860164]] }
thread 'examples::test_rectangle_rotated::tests::test_rectangle_rotated' panicked at src/examples/test_rectangle_rotated.rs:113:9:
assertion failed: (point_a.as_ref().borrow().data() - Vector2::new(0.0, 0.0)).norm() < 0.001


failures:
    examples::test_rectangle_rotated::tests::test_rectangle_rotated

test result: FAILED. 17 passed; 1 failed; 0 ignored; 0 measured; 0 filtered out; finished in 6.30s

Do you get the same?

[TobiasJacob]

Yes, one test case is not yet working, its written down in the Todos of the Readme.

[TobiasJacob]

PS: Some of the constraints could be improved still, but I think the overall structure shouldnt change too much anymore. Also having UI would help a lot for debugging, thats why integrating it into the code base would be helpful.

[hrgdavor]

have you considered a wasm build aswell ? for folks that would like to use this directly ?

[TobiasJacob]

Yeah, you can open up an issue for this in the new isotope repo

[MattFerraro]

Excellent work @TobiasJacob! I'm marking this done and we can have other tickets for integrating into the UI more fully