Vectornaut
b92be312e8
This PR adds code for a Julia-language prototype of a configuration solver, in the `engine-proto` folder. It uses Julia version 1.10.0. ### Approaches Development of this PR 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. Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo> Co-authored-by: Glen Whitney <glen@studioinfinity.org> Reviewed-on: #13 Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net> Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
77 lines
2.1 KiB
Julia
77 lines
2.1 KiB
Julia
include("Engine.jl")
|
|
|
|
using SparseArrays
|
|
|
|
# this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article below
|
|
# includes a nice translation of the problem statement, which was recorded in
|
|
# Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and Present_)
|
|
#
|
|
# "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki
|
|
# https://www.nippon.com/en/japan-topics/c12801/
|
|
#
|
|
|
|
# initialize the partial gram matrix
|
|
J = Int64[]
|
|
K = Int64[]
|
|
values = BigFloat[]
|
|
for s in 1:9
|
|
# each sphere is represented by a spacelike vector
|
|
push!(J, s)
|
|
push!(K, s)
|
|
push!(values, 1)
|
|
|
|
# the circumscribing sphere is internally tangent to all of the other spheres
|
|
if s > 1
|
|
append!(J, [1, s])
|
|
append!(K, [s, 1])
|
|
append!(values, [1, 1])
|
|
end
|
|
|
|
if s > 3
|
|
# each chain sphere is externally tangent to the "sun" and "moon" spheres
|
|
for n in 2:3
|
|
append!(J, [s, n])
|
|
append!(K, [n, s])
|
|
append!(values, [-1, -1])
|
|
end
|
|
|
|
# each chain sphere is externally tangent to the next chain sphere
|
|
s_next = 4 + mod(s-3, 6)
|
|
append!(J, [s, s_next])
|
|
append!(K, [s_next, s])
|
|
append!(values, [-1, -1])
|
|
end
|
|
end
|
|
gram = sparse(J, K, values)
|
|
|
|
# make an initial guess
|
|
guess = hcat(
|
|
Engine.sphere(BigFloat[0, 0, 0], BigFloat(15)),
|
|
Engine.sphere(BigFloat[0, 0, -9], BigFloat(5)),
|
|
Engine.sphere(BigFloat[0, 0, 11], BigFloat(3)),
|
|
(
|
|
Engine.sphere(9*BigFloat[cos(k*π/3), sin(k*π/3), 0], BigFloat(2.5))
|
|
for k in 1:6
|
|
)...
|
|
)
|
|
frozen = [CartesianIndex(4, k) for k in 1:4]
|
|
|
|
# complete the gram matrix using Newton's method with backtracking
|
|
L, success, history = Engine.realize_gram(gram, guess, frozen)
|
|
completed_gram = L'*Engine.Q*L
|
|
println("Completed Gram matrix:\n")
|
|
display(completed_gram)
|
|
if success
|
|
println("\nTarget accuracy achieved!")
|
|
else
|
|
println("\nFailed to reach target accuracy")
|
|
end
|
|
println("Steps: ", size(history.scaled_loss, 1))
|
|
println("Loss: ", history.scaled_loss[end], "\n")
|
|
if success
|
|
println("Chain diameters:")
|
|
println(" ", 1 / L[4,4], " sun (given)")
|
|
for k in 5:9
|
|
println(" ", 1 / L[4,k], " sun")
|
|
end
|
|
end |