Name gradient descent test more specifically

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

@@ -0,0 +1,81 @@
+include("Engine.jl")
+
+using SparseArrays
+using AbstractAlgebra
+using PolynomialRoots
+
+# 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
+guess = sqrt(1/BigFloat(2)) * BigFloat[
+  1  1 -1 -1  0                 -0.1  0.3;
+  1 -1  1 -1  0                 -0.5  0.4;
+  1 -1 -1  1  0                  0.1 -0.2;
+  0  0  0  0 -sqrt(BigFloat(6))  0.3 -0.2;
+  1  1  1  1  2                  0.2  0.1;
+]
+
+# 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)