dyna3/engine-proto/gram-test/overlapping-pyramids.jl

87 lines
2.3 KiB
Julia
Raw Normal View History

include("Engine.jl")
using SparseArrays
using AbstractAlgebra
using PolynomialRoots
using Random
# initialize the partial gram matrix for an arrangement of seven spheres in
# which spheres 1 through 5 are mutually tangent, and spheres 3 through 7 are
# also mutually tangent
J = Int64[]
K = Int64[]
values = BigFloat[]
for j in 1:7
for k in 1:7
if (j <= 5 && k <= 5) || (j >= 3 && k >= 3)
push!(J, j)
push!(K, k)
push!(values, j == k ? 1 : -1)
end
end
end
gram = sparse(J, K, values)
# set the independent variable
#
# using gram[6, 2] or gram[7, 1] as the independent variable seems to stall
# convergence, even if its value comes from a known solution, like
#
# gram[6, 2] = 0.9936131705272925
#
indep_val = -9//5
gram[6, 1] = BigFloat(indep_val)
gram[1, 6] = gram[6, 1]
# in this initial guess, the mutual tangency condition is satisfied for spheres
# 1 through 5
Random.seed!(50793)
guess = hcat(
sqrt(1/BigFloat(2)) * BigFloat[
1 1 -1 -1 0;
1 -1 1 -1 0;
1 -1 -1 1 0;
0 0 0 0 -sqrt(BigFloat(6));
1 1 1 1 2;
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
Engine.rand_on_shell(fill(BigFloat(-1), 2))
)
# complete the gram matrix using gradient descent
L, history = Engine.realize_gram(gram, guess)
completed_gram = L'*Engine.Q*L
println("Completed Gram matrix:\n")
display(completed_gram)
println("\nSteps: ", size(history.stepsize, 1))
println("Loss: ", history.scaled_loss[end], "\n")
# === algebraic check ===
R, gens = polynomial_ring(AbstractAlgebra.Rationals{BigInt}(), ["x", "t₁", "t₂", "t₃"])
x = gens[1]
t = gens[2:4]
S, u = polynomial_ring(AbstractAlgebra.Rationals{BigInt}(), "u")
M = matrix_space(R, 7, 7)
gram_symb = M(R[
1 -1 -1 -1 -1 t[1] t[2];
-1 1 -1 -1 -1 x t[3]
-1 -1 1 -1 -1 -1 -1;
-1 -1 -1 1 -1 -1 -1;
-1 -1 -1 -1 1 -1 -1;
t[1] x -1 -1 -1 1 -1;
t[2] t[3] -1 -1 -1 -1 1
])
rank_constraints = det.([
gram_symb[1:6, 1:6],
gram_symb[2:7, 2:7],
gram_symb[[1, 3, 4, 5, 6, 7], [1, 3, 4, 5, 6, 7]]
])
# solve for x and t
x_constraint = 25//16 * to_univariate(S, evaluate(rank_constraints[1], [2], [indep_val]))
t₂_constraint = 25//16 * to_univariate(S, evaluate(rank_constraints[3], [2], [indep_val]))
x_vals = PolynomialRoots.roots(x_constraint.coeffs)
t₂_vals = PolynomialRoots.roots(t₂_constraint.coeffs)