glen/dyna3
1
0
Fork 0
Code Issues 14 Pull Requests 1 Actions Projects Releases 2 Wiki Activity

Engine prototype #13

Merged
glen merged 133 commits from engine-proto into main 2024-10-21 03:18:48 +00:00
Conversation 9 Commits 133 Files Changed 20 +3977 -21
Vectornaut commented 2024-08-20 09:09:14 +00:00
Collaborator
Copy Link

The branch to be merged tracks our work on prototyping a geometric constraint-solving engine. This work is in the engine-proto folder. Most of the code is in Julia (version 1.10.0).

Approaches

We tried two broad approaches to the constraint geometry problem. Each one suggested various solution techniques. The Gram matrix approach, with the low-rank factorization technique, seems the most promising.

  • Algebraic (In the alg-test subfolder). Write the constraints as polynomials in the inversive coordinates of the elements, and use computational algebraic geometry techniques to solve the resulting system. We tried the following techniques.

    • Gröbner bases (Engine.Algebraic.jl). Symbolic. Find a Gröbner basis for the ideal generated by the constraint equations. Information about the solution variety, like its codimension, is then relatively easy to extract.
    • Homotopy continuation (Engine.Numerical.jl). Numerical. Cut the solution set along a random hyperplane to get a generic zero-dimensional slice, and then use a fancy homotopy technique to approximate the points in that slice.

    A few notes about our experiences can be found on the engine prototype wiki page.

  • Gram matrix (in the gram-test subfolder). A construction is described completely, up to conformal transformations, by the Gram matrix of the vectors representing its elements. Express the constraints as fixed entries of the Gram matrix, and use numerical linear algebra techniques to find a list of vectors whose Gram matrix fits the bill. We tried the following techniques.

    • LDL decomposition (gram-test.sage, gram-test.jl, overlap-test.jl). Find a cluster of up to five elements whose Gram matrix is completely filled in by the constraints. Use LDL decomposition to find a list of vectors with that Gram matrix. This technique can be made algebraic, as seen in overlap-test.jl.
    • Low-rank factorization (source files listed in findings section). Write down a quadratic loss function that says how far a set of vectors is from meeting the Gram matrix constraints. Use a smooth optimization technique like Newton's method or gradient descent to find a zero of the loss function. In addition to the polished prototype described in the results section, we have an early prototype using an off-the-shelf factorization package (low-rank-test.jl) and an visualization of the loss function landscape near global minima (basin-shapes.jl).

    The Gram matrix parameterization wiki page contains detailed notes on this approach.

Findings

With the algebraic approach, we hit a performance wall pretty quickly as our constructions grew. It was often hard to find real solutions of the polynomial system, since the techniques we use work most naturally in the complex world.

With the Gram matrix approach, on the other hand, we could solve interesting problems in acceptably short times using the low-rank factorization technique. We put the optimization routine in its own module (Engine.jl) and used it to solve five example problems:

  • overlapping-pyramids.jl
  • circles-in-triangle.jl
  • sphere-in-tetrahedron.jl
  • tetrahedron-radius-ratio.jl
  • irisawa-hexlet.jl

We plan to use low-rank factorization of the Gram matrix in our first app prototype.

Visualizations

We used the visualizer in the ganja-test folder to visually check our low-rank factorization results. The visualizer runs Ganja.js in an Electron app, made with Blink. Although Ganja.js makes beautiful pictures under most circumstances, we found two obstacles to using it in production.

  • It seems to have precision problems with low-curvature spheres.
  • We couldn't figure out how to customize its clipping and transparency settings, and the default settings often obscure construction details.
The branch to be merged tracks our work on prototyping a geometric constraint-solving engine. This work is in the `engine-proto` folder. Most of the code is in Julia (version 1.10.0). ### Approaches We tried two broad approaches to the constraint geometry problem. Each one suggested various solution techniques. The Gram matrix approach, with the low-rank factorization technique, seems the most promising. - **Algebraic** *(In the `alg-test` subfolder).* Write the constraints as polynomials in the inversive coordinates of the elements, and use computational algebraic geometry techniques to solve the resulting system. We tried the following techniques. - **Gröbner bases** *(`Engine.Algebraic.jl`).* Symbolic. Find a Gröbner basis for the ideal generated by the constraint equations. Information about the solution variety, like its codimension, is then relatively easy to extract. - **Homotopy continuation** *(`Engine.Numerical.jl`).* Numerical. Cut the solution set along a random hyperplane to get a generic zero-dimensional slice, and then use a fancy homotopy technique to approximate the points in that slice. A few notes about our experiences can be found on the [engine prototype](wiki/Engine-prototype) wiki page. - **Gram matrix** *(in the `gram-test` subfolder).* A construction is described completely, up to conformal transformations, by the Gram matrix of the vectors representing its elements. Express the constraints as fixed entries of the Gram matrix, and use numerical linear algebra techniques to find a list of vectors whose Gram matrix fits the bill. We tried the following techniques. - **LDL decomposition** *(`gram-test.sage`, `gram-test.jl`, `overlap-test.jl`).* Find a cluster of up to five elements whose Gram matrix is completely filled in by the constraints. Use LDL decomposition to find a list of vectors with that Gram matrix. This technique can be made algebraic, as seen in `overlap-test.jl`. - **Low-rank factorization** *(source files listed in findings section).* Write down a quadratic loss function that says how far a set of vectors is from meeting the Gram matrix constraints. Use a smooth optimization technique like Newton's method or gradient descent to find a zero of the loss function. In addition to the polished prototype described in the results section, we have an early prototype using an off-the-shelf factorization package (`low-rank-test.jl`) and an visualization of the loss function landscape near global minima (`basin-shapes.jl`). The [Gram matrix parameterization](wiki/Gram-matrix-parameterization) wiki page contains detailed notes on this approach. ### Findings With the algebraic approach, we hit a performance wall pretty quickly as our constructions grew. It was often hard to find real solutions of the polynomial system, since the techniques we use work most naturally in the complex world. With the Gram matrix approach, on the other hand, we could solve interesting problems in acceptably short times using the low-rank factorization technique. We put the optimization routine in its own module (`Engine.jl`) and used it to solve five example problems: - `overlapping-pyramids.jl` - `circles-in-triangle.jl` - `sphere-in-tetrahedron.jl` - `tetrahedron-radius-ratio.jl` - `irisawa-hexlet.jl` We plan to use low-rank factorization of the Gram matrix in our first app prototype. ### Visualizations We used the visualizer in the `ganja-test` folder to visually check our low-rank factorization results. The visualizer runs [Ganja.js](https://enkimute.github.io/ganja.js/) in an Electron app, made with [Blink](https://github.com/JuliaGizmos/Blink.jl). Although Ganja.js makes beautiful pictures under most circumstances, we found two obstacles to using it in production. - It seems to have precision problems with low-curvature spheres. - We couldn't figure out how to customize its clipping and transparency settings, and the default settings often obscure construction details.
Vectornaut added 113 commits 2024-08-20 09:09:16 +00:00 
In other words, order coordinates like
  (rₛ₁, rₛ₂, sₛ₁, sₛ₂, xₛ₁, xₛ₂, xₚ₃, yₛ₁, yₛ₂, yₚ₃, zₛ₁, zₛ₂, zₚ₃)
instead of like
  (rₛ₁, sₛ₁, xₛ₁, yₛ₁, zₛ₁, rₛ₂, sₛ₂, xₛ₂, yₛ₂, zₛ₂, xₚ₃, yₚ₃, zₚ₃).

In the test cases, this really cuts down the size of the Gröbner basis.
I think we need this to find the dimension of the solution variety.
The extension should also let us work over finite fields of prime order,
although we don't need to do that.
You're right: this naming convention seems to be standard for Julia
modules now.
In the cases I've tried so far, this leads to substantially smaller
Gröbner bases.
The example hyperplane yields a single solution, with multiplicity six. You can
find it analytically by hand, and homotopy continuation finds it numerically.
The previous cut was supposed to do this, but I was missing some parentheses.
Choose hyperplanes that go through the trivial solution.
In all but a few cases (for example, a single point on a plane), we
should be able to us the radius-coradius boost symmetry to make the
average co-radius—representing the "overall scale"—roughly one.
Unit normals are uniformly distributed over the sphere.
For three mutually tangent spheres, I couldn't find real solutions.
Guess conditions that make the scaling constraint impossible to satisfy.
Three points on two spheres is too much.
This seems to restore reproducibility.
I'm not sure why the solver wasn't working before. It might've been just
an unlucky random number draw.
It seems like there are real solutions if and only if the product of the
cosines is positive.
The signature of the Minkowski form on the subspace spanned by the Gram
matrix should tell us what the big Gram matrix has to look like
The best technique I've found so far is the homemade gradient descent
routine in `descent-test.jl`.
Declare JavaScript variables. Revise Julia comments to match new code.
Try configuration of five tangent spheres.
Visualize low-rank factorization results.
Get visibility controls.
The previos version accidentally returned steps in Float64.
The completed gram matrix from this commit matches the one from commit
e7dde58 to six decimal places.
This needs to be rewritten: it can fail at generating spacelike vectors.
Randomly perturb the pre-solved part of the guess, and randomly choose
the unsolved part.
Our formula for the improvement theshold works when the step size is
an absolute distance. However, in commit `4d5ea06`, the step size was
measured relative to the current gradient instead. This commit scales
the base step to unit length, so now the step size really is an absolute
distance.
This code is a mess, but I'm committing it to record a working state
before I start trying to clean up.
Drop experimental singularity handling strategies. Reduce the default
tolerance to within 64-bit floating point precision. Report success.
In previous commits, the `circles-in-triangle` example converged much
more slowly in BigFloat precision than in Float64 precision. This
turned out to be a sign of a bug in the Float64 computation: converting
the Gram matrix using `Float64.()` dropped the explicit zeros, removing
many constraints and making the problem much easier to solve. This
commit corrects the Gram matrix conversion. The Float64 search now
solves the same problem as the BigFloat search, with comparable
performance.
This performs much better than the trust region Newton's method for the
actual `circles-in-triangle` problem. (The trust region method performs
better for the simplified problem produced by the conversion bug.)
This changes the meaning of `indep_val` in the overlapping pyramids
example, so we adjust `indep_val` to get a nice-looking construction.
Add the vertices of the tetrahedron to the `sphere-in-tetrahedron`
example.
This reverts commit 4728959ae0, which
actually gave the spheres negative radii! I got confused by the sign
convention differences between the notes and the engine.
Jiggle the vertex guesses. Put the circumscribed sphere guess on-shell.
The balance between the light cone basis vectors was wrong, throwing the
point's coordinates off by a factor of two.
Clarify the relevant notes too.
Since the self-product of the point at infinity is left unspecified, the
first three components can vary without violating any constraints. To
keep the point at infinity where it's supposed to be, we freeze all of
its components.
Hat tip Romy, who sent me the article on sangaku that led me to this
problem.
The approach in the deleted file can't work, because the "sun" and
"moon" spheres can't be placed arbitrarily.
Abe uses the names "sun" and "moon" for what Wikipedia calls the nucleus
spheres.
We're using the Gram matrix engine for the next stage of development,
so the algebraic engine shouldn't be at the top level anymore.
glen added 1 commit 2024-09-19 02:50:49 +00:00
glen added 1 commit 2024-09-19 02:52:16 +00:00
glen added 1 commit 2024-09-19 03:12:29 +00:00
glen added 1 commit 2024-09-19 03:27:21 +00:00
glen commented 2024-09-19 03:39:35 +00:00
Owner
Copy Link

OK, by way of starting the review of this PR, I looked over the first six rows of notes/inversive.md, and had a few questions:

  1. I made a few changes, please double-check me, thanks.
  2. The right column of the third row had become garbled. Is there in fact a scaling choice in the point coordinates? I think so, I think they are like projective coordinates where we could represent (1,4,8) by (9,1,1,4,8) as described, but just as well by (18,2,2,4,8), etc. We simply pick the one with second coordinate 1 for canonicalness. Please fill in the ... as you feel appropriate.
  3. In the sixth row, you give a nice example about orientations matching. It seems to beg the question though about whether our sphere coordinates are similarly oriented, or whether in these coordinates they always have outward-pointing normals. Is there another representation of the unit sphere centered at (0,0,1) but with inward-pointing normals? Is that represented by negative radius, i.e. (0, -1, 0, 0, 1)? That vector has Q-norm -1, but sadly it still has Q-product -1 with the plane (0,0,0,0,1), which would seem to refute that now their normals are opposite at the point of tangency. So does (0, -1, 0, 0, 1) represent anything in our inversive coordinates, and if so what? Please put appropriate notes in the first few rows of these notes to clarify these points, thanks.

I will finish the review when you've had a chance to do these things, thanks.

OK, by way of starting the review of this PR, I looked over the first six rows of notes/inversive.md, and had a few questions: 1) I made a few changes, please double-check me, thanks. 2) The right column of the third row had become garbled. Is there in fact a scaling choice in the point coordinates? I think so, I think they are like projective coordinates where we could represent (1,4,8) by (9,1,1,4,8) as described, but just as well by (18,2,2,4,8), etc. We simply pick the one with second coordinate 1 for canonicalness. Please fill in the ... as you feel appropriate. 3) In the sixth row, you give a nice example about orientations matching. It seems to beg the question though about whether our sphere coordinates are similarly oriented, or whether in these coordinates they always have outward-pointing normals. Is there another representation of the unit sphere centered at (0,0,1) but with inward-pointing normals? Is that represented by negative radius, i.e. (0, -1, 0, 0, 1)? That vector has Q-norm -1, but sadly it still has Q-product -1 with the plane (0,0,0,0,1), which would seem to refute that now their normals are opposite at the point of tangency. So does (0, -1, 0, 0, 1) represent anything in our inversive coordinates, and if so what? Please put appropriate notes in the first few rows of these notes to clarify these points, thanks. I will finish the review when you've had a chance to do these things, thanks.
glen commented 2024-09-20 02:48:06 +00:00
Owner
Copy Link

P.S. I have now understood from today's discussion that the coordinates in the current prototype are not quite exactly these. I think we should probably leave these notes alone, but just clearly document wherever the coordinate entities are declared exactly what they are are and how they relate to these coordinates. Thanks!

P.S. I have now understood from today's discussion that the coordinates in the current prototype are not quite exactly these. I think we should probably leave these notes alone, but just clearly document wherever the coordinate entities are declared exactly what they are are and how they relate to these coordinates. Thanks!
glen commented 2024-09-21 06:04:17 +00:00
Owner
Copy Link

P.P.S By "leave alone," I didn't mean that we shouldn't do the three things I requested in the review above; I just meant that we shouldn't change the inversive coordinate system described in the notes, but rather document in the code how the one there differs from what's described in the notes.

P.P.S By "leave alone," I didn't mean that we shouldn't do the three things I requested in the review above; I just meant that we shouldn't change the inversive coordinate system described in the notes, but rather document in the code how the one there differs from what's described in the notes.
Vectornaut added 2 commits 2024-09-26 08:07:44 +00:00 
In the process, clarify the signed distance from a point to a sphere and
add inversion across a sphere.
Vectornaut commented 2024-09-26 08:34:12 +00:00
Author
Collaborator
Copy Link
  1. I agree with the changes you made in 8084fde, bd3e350, a182b66, 23ecca3.
  2. I added an example of an alternate normalization for points. The Q(I_\infty, I_P) = 1/2 normalization that we use is useful for working with distances, because the direction of I_\infty says where the point at infinity is and the normalization of I_\infty encodes a choice of scale.
  3. I wrote more about orientation in the sphere and plane rows. The spheres represented by (0, 1, 0, 0, 1) and (0, -1, 0, 0, 1) are related by a reflection across the origin followed by a reversal of orientation. The pair you're looking for is (0, 1, 0, 0, 1) and (0, -1, 0, 0, -1), which coincide and have opposite orientations.
1. I agree with the changes you made in 8084fde, bd3e350, a182b66, 23ecca3. 2. I added an example of an alternate normalization for points. The $Q(I_\infty, I_P) = 1/2$ normalization that we use is useful for working with distances, because the direction of $I_\infty$ says where the point at infinity is and the normalization of $I_\infty$ encodes a choice of scale. 3. I wrote more about orientation in the sphere and plane rows. The spheres represented by $(0, 1, 0, 0, 1)$ and $(0, -1, 0, 0, 1)$ are related by a reflection across the origin followed by a reversal of orientation. The pair you're looking for is $(0, 1, 0, 0, 1)$ and $(0, -1, 0, 0, -1)$, which coincide and have opposite orientations.
glen added 1 commit 2024-10-04 03:39:54 +00:00
glen added 1 commit 2024-10-04 04:20:46 +00:00
glen commented 2024-10-04 04:53:10 +00:00
Owner
Copy Link

More about the inversive coordinates, sorry.

  1. I slightly amended the condition that two sphere/planes (sphlanes?) intersect: seems it should be |Q(I,J)| \leq 1, with ordinary (not double) absolute value bars and \leq rather than <. Please confirm you agree.

  2. I am now confused by the row "P is center of sphere rep'd by I":
    a. Would it be clearer in the formula cell instead of "..., and I_P is normalized to have $Q(I_\infty, I_P) = 1/2$" to say "..., as long as I_P is normalized in the standard way in which Q(I_\infty, I_P) = 1/2?
    b. In the notes cell, shouldn't it be "(1) the point P has signed distance -r from the sphere", rather than -r/2?

  3. I am also confused by the row "Signed distance between point rep'd by V and sphere/plane rep'd by I":
    a. In the notes cell, you give the example that the point rep'd by V=(d^2, 1, 0, 0, d) and the sphere by I=(0, 1/r, 0, 0, -1) are obviously at signed distance d from one another, which I agree with. But given this, it seems to me from the formula that in this case Q(V,I) = d^2/2r + d from the formula for Q (not d^2/2r - d).
    b. Presuming that's the case, and that the formula then works in general, then the formula cell in this row should have Q(I_\infty, I)d^2 + d on the right-hand side. Would it be helpful to also include the following in the formula cell?

So if V is normalized in the standard way and the radius of the sphere is r, this is just Q(V,I) = d^2/2r + d.

And possibly also:

So that in the case I represents a plane and V is normalized in the standard way, this reduces to Q(V,I) = d.

(Assuming I am correct in all of this, that's a very pleasant reduction (and means we can easily add such a constraint for plane-point distance in the current engine prototype.)

I note that 3a and 2b check with each other, in that when P has signed distance -r from a sphere of radius r rep'd by I, then by the formula, Q(I_P,I) = r^2/2r - r = r/2.

If you agree with these proposed changes please either just push a commit with them or let me know that you would like me to. Then I can proceed with the review.

More about the inversive coordinates, sorry. 1. I slightly amended the condition that two sphere/planes (sphlanes?) intersect: seems it should be $|Q(I,J)| \leq 1$, with ordinary (not double) absolute value bars and $\leq$ rather than $<$. Please confirm you agree. 2. I am now confused by the row "P is center of sphere rep'd by I": a. Would it be clearer in the formula cell instead of "..., and $I_P$ is normalized to have $Q(I_\infty, I_P) = 1/2$" to say "..., as long as $I_P$ is normalized in the standard way in which $Q(I_\infty, I_P) = 1/2$? b. In the notes cell, shouldn't it be "(1) the point $P$ has signed distance $-r$ from the sphere", rather than $-r/2$? 3. I am also confused by the row "Signed distance between point rep'd by V and sphere/plane rep'd by I": a. In the notes cell, you give the example that the point rep'd by $V=(d^2, 1, 0, 0, d)$ and the sphere by $I=(0, 1/r, 0, 0, -1)$ are obviously at signed distance $d$ from one another, which I agree with. But given this, it seems to me from the formula that in this case $Q(V,I) = d^2/2r + d$ from the formula for $Q$ (not $d^2/2r - d$). b. Presuming that's the case, and that the formula then works in general, then the formula cell in this row should have $Q(I_\infty, I)d^2 + d$ on the right-hand side. Would it be helpful to also include the following in the formula cell? So if $V$ is normalized in the standard way and the radius of the sphere is $r$, this is just $Q(V,I) = d^2/2r + d$. And possibly also: So that in the case $I$ represents a plane and $V$ is normalized in the standard way, this reduces to $Q(V,I) = d$. (Assuming I am correct in all of this, that's a very pleasant reduction (and means we can easily add such a constraint for plane-point distance in the current engine prototype.) I note that 3a and 2b check with each other, in that when $P$ has signed distance $-r$ from a sphere of radius $r$ rep'd by I, then by the formula, $Q(I_P,I) = r^2/2r - r = r/2$. If you agree with these proposed changes please either just push a commit with them or let me know that you would like me to. Then I can proceed with the review.
Vectornaut added 2 commits 2024-10-15 21:42:28 +00:00
Vectornaut commented 2024-10-15 21:45:48 +00:00
Author
Collaborator
Copy Link

I've checked and incorporated the changes you suggested, so you can start reviewing again!

I've checked and incorporated the changes you suggested, so you can start reviewing again!
Vectornaut added 1 commit 2024-10-16 23:01:25 +00:00
Vectornaut commented 2024-10-16 23:02:55 +00:00
Author
Collaborator
Copy Link

I just pushed a minor change to Engine.jl, which fixes a mistake in a comment.

I just pushed a minor change to `Engine.jl`, which fixes a mistake in a comment.
glen added 1 commit 2024-10-20 23:06:38 +00:00
glen added 2 commits 2024-10-20 23:15:43 +00:00
glen added 1 commit 2024-10-20 23:17:31 +00:00
glen added 1 commit 2024-10-21 02:27:22 +00:00
glen added 1 commit 2024-10-21 02:33:47 +00:00
glen added 1 commit 2024-10-21 02:41:30 +00:00
glen added 1 commit 2024-10-21 02:53:50 +00:00
glen added 1 commit 2024-10-21 03:03:05 +00:00
glen commented 2024-10-21 03:11:11 +00:00
Owner
Copy Link

  1. OK, I just made a bunch more edits to notes/inversive.md; mostly they are typographical, but there is some expanded content. I am not going to wait for your review to merge this PR, because I am pretty confident, but I would like you to read through the latest version of this file wherever it has ended up when you get a chance. Thanks.

  2. I see you've checked in entire 3rd-party files (ganja.js). Going forward, we don't want to do that, but instead have a build process that downloads and installs 3rd party files without their needing to be in our repo. However, since the engine-proto code is destined eventually to disappear from the repo, I am not going to hold up this PR for that sake.

1. OK, I just made a bunch more edits to notes/inversive.md; mostly they are typographical, but there is some expanded content. I am not going to wait for your review to merge this PR, because I am pretty confident, but I would like you to read through the latest version of this file wherever it has ended up when you get a chance. Thanks. 2. I see you've checked in entire 3rd-party files (ganja.js). Going forward, we don't want to do that, but instead have a build process that downloads and installs 3rd party files without their needing to be in our repo. However, since the engine-proto code is destined eventually to disappear from the repo, I am not going to hold up this PR for that sake.
glen commented 2024-10-21 03:13:38 +00:00
Owner
Copy Link
  1. Similarly, the code seems to be organized with implementation code and tests interspersed, as opposed to clearly delineated. Not appropriate for long-term code, but I am not going to hold this prototype up on that basis.
3. Similarly, the code seems to be organized with implementation code and tests interspersed, as opposed to clearly delineated. Not appropriate for long-term code, but I am not going to hold this prototype up on that basis.
glen merged commit b92be312e8 into main 2024-10-21 03:18:48 +00:00
glen referenced this issue from a commit 2024-10-21 03:18:49 +00:00
Engine prototype (#13)
Sign in to join this conversation.
Reviewers
No reviewers
Labels
Clear labels   
bug

Something is not working

  
duplicate

This issue or pull request already exists

  
enhancement

New feature

  
question

More information is needed

  
todo

Planned step in implementation

  
wontfix

This won't be fixed

No Label
bug
duplicate
enhancement
question
todo
wontfix
Milestone
Clear milestone
No items
No Milestone
Projects
Clear projects
No project
Assignees
Clear assignees
No Assignees
2 Participants
Notifications
Due Date
The due date is invalid or out of range. Please use the format 'yyyy-mm-dd'.

No due date set.

Dependencies

No dependencies set.

Reference: glen/dyna3#13
No description provided.
Delete Branch "engine-proto"

Deleting a branch is permanent. Although the deleted branch may continue to exist for a short time before it actually gets removed, it CANNOT be undone in most cases. Continue?