Using gomez for solving geometric constraints / intecepts

[twitchyliquid64]

Hi!

I'm wondering if gomez is appropriate to use for implementing a geometric solver (like this)?

Based on the description in the docs:

F(x) = 0,

where F(x) = { f1(x), ..., fn(x) }
and x = { x1, ..., xn }

I'm not sure if each Fn(x) can only be dependent on x, can it be dependent on other variables?

For instance if I have a system of equations that solve for the intersection of a line and a circle, I would need to reference the parameters of the line and circle across the same equation.

NB: Apologies if I'm totally missing the mark of understanding here, I'm learning :)

[pnevyk]

Yes, F(x) is represented by a struct which can have any fields that you can reference in your function evaluation. So, in your example, given a * x1 + b * x2 = 0 for line and (x1 - c)^2 + (x2 - d)^2 = r^2 for circle, you would create a struct

struct LineCircleIntersection {
    a: f64,
    b: f64,
    c: f64,
    d: f64,
    r: f64,
}

and implement all necessary traits:

impl Problem for LineCircleIntersection {
    // ...
}

impl System for LineCircleIntersection {
    fn eval<Sx, Sfx>(
        &self,
        x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
        fx: &mut na::Vector<Self::Scalar, Self::Dim, Sfx>,
    ) -> Result<(), ProblemError>
    where
        Sx: na::storage::Storage<Self::Scalar, Self::Dim> + IsContiguous,
        Sfx: na::storage::StorageMut<Self::Scalar, Self::Dim>,
    {
        // Fill `fx` which represents residuals of the system of equations.
        // Use fields in `self`.
        Ok(())
    }
}

All the trait bounds might feel overwhelming, that's because I use nalgebra. I am planning to rewrite everything and use faer instead, which should simplify the interface.

If you run into problems, feel free to ask here.

[twitchyliquid64]

Thanks so much for the pointers!

Is there a requirement that the variable at index n in the vector be specifically related to the function at index n?

I was originally thinking yes, but after just trying it without caring about which variables and functions go in which slots it just worked:

use gomez::nalgebra as na;
use gomez::prelude::*;
use gomez::solver::TrustRegion;
use na::{Dim, DimName, IsContiguous};

struct Point {
    x: f64,
    y: f64,
}

// Find a point equidistant from two other points
struct DistFromPoints {
    points: [Point; 2],
    dist: f64,
}

impl Problem for DistFromPoints {
    type Scalar = f64;
    type Dim = na::U2;

    fn dim(&self) -> Self::Dim {
        na::U2::name()
    }
}

impl System for DistFromPoints {
    fn eval<Sx, Sfx>(
        &self,
        x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
        fx: &mut na::Vector<Self::Scalar, Self::Dim, Sfx>,
    ) -> Result<(), ProblemError>
    where
        Sx: na::storage::Storage<Self::Scalar, Self::Dim> + IsContiguous,
        Sfx: na::storage::StorageMut<Self::Scalar, Self::Dim>,
    {
        // 0 = distance(point, current_guess) - desired_distance
        fx[0] = ((x[0] - self.points[0].x).powi(2) + (x[1] - self.points[0].y).powi(2))
            - self.dist.powi(2);
        fx[1] = ((x[0] - self.points[1].x).powi(2) + (x[1] - self.points[1].y).powi(2))
            - self.dist.powi(2);

        Ok(())
    }
}

// <snip func main() >

Output:

iter = 1	|| fx || = 286.35642126552705	x = [-19.0, 13.0]
iter = 2	|| fx || = 231.08517186240644	x = [-14.964180335249775, 21.01677717631165]
iter = 3	|| fx || = 231.08517186240644	x = [-14.964180335249775, 21.01677717631165]
iter = 4	|| fx || = 151.59562228173127	x = [-10.717027237388809, 22.077465495809374]
iter = 5	|| fx || = 151.59562228173127	x = [-10.717027237388809, 22.077465495809374]
iter = 6	|| fx || = 114.84592356742507	x = [-7.818162935273372, 22.735701708640608]
iter = 7	|| fx || = 100.66641542931053	x = [-5.39000586385896e-8, 25.907178638126684]
iter = 8	|| fx || = 2.6783437941030357	x = [-5.39000586385896e-8, 24.53352553260862]
iter = 9	|| fx || = 0.0022076401825397835	x = [-5.39000586385896e-8, 24.4949292923505]
iter = 10	|| fx || = 0.0000007622641508539823	x = [-5.39000586385896e-8, 24.494897427858568]
solved

If they aren't related at all, can I have a different number of functions in the system to variables to guess?

[pnevyk]

Is there a requirement that the variable at index n in the vector be specifically related to the function at index n?

There isn't a requirement on the order, it's basically a naming (like if you would use y, x instead of x, y).

If they aren't related at all, can I have a different number of functions in the system to variables to guess?

Unfortunately, this is not supported. Having less equations than variables is not possible in general, the system would have infinite number of solutions (or zero). Having more equations than variables could be possible, but it requires special support by the algorithm, and it is not implemented in this library.

I recommend reading this Wikipedia article. This library supports only "exactly determined and consistent" systems.

[twitchyliquid64]

Thanks! That all makes sense.

Why is it that this library won't work for underdetermined, consistent systems? Indeed in my DistFromPoints system above, there were infinite solutions, but the TrustRegion solver had no issue converging on one of the solutions, heavily dependent on the initial guess.

Solvespace relies heavily on being able to solve underdetermined systems, to power its feature where it lets you drag a point and it will remain constrained to valid positions for that point. It does this with a numerical solver, initial guesses based on the current drag, and by scaling the variables such that the solver impacts the ones not being dragged the most. Would something like that not work here?

[pnevyk]

Interesting. I have to admit that I am not that strong in the theory behind, I just needed to solve systems of equations in Rust and so I implemented something. If you believe that your use case should work and it also seems like it works, then you might be fine. I can't say whether it will or will not work, because I don't know :) But I hope it will! And if not, but we can do anything about it, let me know.

[twitchyliquid64]

No worries at all, thank you for making this library! Fwiw, this is definitely one of the more thought-out rust libraries in terms of algorithms used, tests, API interface etc.

I'm going to go ahead and try and use it for my underdetermined-systems use-cases, will let you know how it goes!!

[pnevyk]

this is definitely one of the more thought-out rust libraries in terms of algorithms used, tests, API interface etc

Oh thanks for the kind words. There is a big room for improvement regarding the API, but your interest motivated me to give this library some love again. I made some changes to the API, you can read the list in #9. Any feedback would be appreciated. I will now give it a few days and if I can't think of any other changes I would like to make, I will release a new version. There will be breaking changes, I may write a short upgrade guide once I release the version.

[twitchyliquid64]

Sounds great! Thanks for your work on this library :)

Circling back on how its going with my use case: I had a lot of trouble getting gomez to work for specifically geometric CAD / underconstrained systems. Whereas a slow solver like gradient descent will tend to land on solutions that are close to the initial values, this is rarely the case for TrustedRegion. As a result, I'll be dragging a point thats very underconstrained (usually by distance or something to another point), and it will be zooming all over the screen, along with any other underconstrained points.

The other issue I had was when my constraint functions had multiple solutions, such as the distance function. Unless the residual bounced about zero, having multiple solutions which could be bounced back and forth made the solver unstable and usually unable to progress.

This is all on me, because as you said gomez is meant for exactly constrained systems. But thought you would be curious how it went.

[pnevyk]

This is all on me, because as you said gomez is meant for exactly constrained systems. But thought you would be curious how it went.

Yes, I am definitely curious. Thanks for the update.

In order to debug the problems, you can install a log-compatible logger (e.g., env_logger) and set it to output debug level logs. It will print a lot of things, but it should give a better idea about what is going on under the hood. If it is possible, you could also paste the logs here so I can try to make sense out of them.

Whereas a slow solver like gradient descent will tend to land on solutions that are close to the initial values, this is rarely the case for TrustedRegion.

This is really interesting. Seeing the logs could help us to understand what is going on, maybe it's a bug in gomez or it can be fixed by tweaking the settings. By the way, do you use your own implementation of gradient descent or a library?

Unless the residual bounced about zero, having multiple solutions which could be bounced back and forth made the solver unstable and usually unable to progress.

You could try to set the maximum size and initial size of the trust region to a small value, so the solver makes much smaller, but hopefully less chaotic steps.

You could also try to use Nelder-Mead algorithm instead of trust region. It's a simple function optimization method, so it's not really powerful for solving systems of equations, but I would be curious how that performs for your case (probably badly).

[twitchyliquid64]

I'll give Nelder-Mead a try!

On the gradient-descent side, I made my own, but its pretty terrible & I don't know what I'm doing: https://gist.github.com/twitchyliquid64/65bb3cddf4287d58f56c24035febd775

But it mostly works lol

[pnevyk]

When I started with this library, I had no idea what I was doing either :) But in time, with trial and error, it gets better.

[papyDoctor]

Hi @pnevyk @twitchyliquid64
I'm working since a few months on an open-source 2D CAD/CAM design for my own using Rust on WebAssembly. Since a few weeks I'm working on implementing geometric constrains. After web searching some informations (state of the art CGS,...) I fell in this crate and I'm very happy to find something in Rust that address the numerical solving of non linear equation. But as twitchyliquid64 wrote, our problem of geometric constrains are mainly underdetermined and consistent.
Nevertheless I'll give a try to this crate to see what happens and let you know :)

[twitchyliquid64]

@papyDoctor I ended up with my own, very strange solver for https://github.com/twitchyliquid64/liquid-cad (online demo: https://liquid-cad.twitchyliquid.net/)

But lmk what you come up with, always excited to swap notes!

@MattFerraro is also working in this space, his approach is a spring-based model with an RK4 solver: https://github.com/MattFerraro/CADmium

[papyDoctor]

@twitchyliquid64 Ah, you'r using egui :)
I don't use it since it's too heavy for me (all the wasm code must stay <2MB to be embedded in this device or this one.
I've an (obsolete) online demo too.
I'll check also your solver as it is related to my constraint problem :)
It will be interesting to have a kind of place where we discuss the various problems we faced for this kind of app.

[pnevyk]

It's good to know that there is a demand for underdetermined systems. I am definitely open to adjusting this crate to support this type of problems, at least API-wise. Then it would be possible to contribute an implementation of such solver to this crate.

@papyDoctor @twitchyliquid64 I would appreciate if you keep sharing your progress, code and literature here. You can also use this issue for discussion although there are probably better tools for that. I could enable GitHub Discussions for this repo if you want.

[papyDoctor]

@pnevyk Yes please enable a discussion because I would like to share some use case I have. And indeed I'm no more sure about the absolute necessity of an underdetermined system. Also let me some times to implement your crate in my project.
Then wait before working on it :)

[pnevyk]

Discussions are now enabled.

I now realized that the underdetermined systems should already be possible (theoretically). One just needs to calculate fx[i] for i up to the number of equations, and (to be safe) fill the rest with zeros. It's overdetermined systems that would cause the trouble, because the fx vector is of size n and if the number of equations is higher than n then there is (at the moment) no way how to pass the values of extra equations.

If you will run into issues with the API and documentation during your implementation, let me know.

[twitchyliquid64]

@papyDoctor My use-case is defs PC & web-oriented so defs enjoying egui :) Im curious as to what your use-case is with trying to do this embedded? Ive always thought of embedded stuff to be the CNC controller, which would be working on G-code or something simpler like that.

@pnevyk lovely, thanks!

[papyDoctor]

@twitchyliquid64 The little boards I've linked will be the CNC controller :) but specifically for cutting machine (plasma and laser). G-code is an old good format but it has some limitation (for ex: no bezier path) that I want to overcome.
The purpose of the application/system is to have CAD with integrated CAM with paths in mind. This application is served by the embedded device to a browser, at the end it receives the file from browser and translate it to work with the motors, switches,…

[twitchyliquid64]

Gotcha! I definitely want to add basic CAM support to liquid-cad, my approach was going to be to generate an offset curve with the tool radius subtracted, and covert that to G-code: https://docs.rs/kurbo/latest/x86_64-apple-darwin/kurbo/offset/struct.CubicOffset.html

The G5 code supports bezier paths in the X-Y plane.

[papyDoctor]

Ohhh, I wasn't aware of this crate. I'll check and use it if not too heavy.
Thks for G5 command, I've never see it (it's a shame).
Let's get in touch, we're working nearly with the same business logic!

[papyDoctor]

@twitchyliquid64 @pnevyk I've made a use case on the discussion