engine-proto/gram-test/circles-in-triangle.jl

@@ -0,0 +1,62 @@
+include("Engine.jl")
+
+using SparseArrays
+using AbstractAlgebra
+using PolynomialRoots
+
+# initialize the partial gram matrix for a sphere inscribed in a regular
+# tetrahedron
+J = Int64[]
+K = Int64[]
+values = BigFloat[]
+for j in 1:8
+  for k in 1:8
+    filled = false
+    if j == k
+      push!(values, 1)
+      filled = true
+    elseif (j == 1 || k == 1)
+      push!(values, 0)
+      filled = true
+    elseif (j == 2 || k == 2)
+      push!(values, -1)
+      filled = true
+    end
+    if filled
+      push!(J, j)
+      push!(K, k)
+    end
+  end
+end
+append!(J, [6, 4, 6, 5, 7, 5, 7, 3, 8, 3, 8, 4])
+append!(K, [4, 6, 5, 6, 5, 7, 3, 7, 3, 8, 4, 8])
+append!(values, fill(-1, 12))
+gram = sparse(J, K, values)
+
+# set initial guess (random)
+## Random.seed!(58271) # stuck; step size collapses on step 48
+## Random.seed!(58272) # good convergence
+## Random.seed!(58273) # stuck; step size collapses on step 18
+## Random.seed!(58274) # stuck
+## Random.seed!(58275) #
+## guess = Engine.rand_on_shell(fill(BigFloat(-1), 8))
+
+# set initial guess
+guess = hcat(
+  Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
+  Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)),
+  Engine.plane(BigFloat[1, 0, 0], BigFloat(1)),
+  Engine.plane(BigFloat[cos(2pi/3), sin(2pi/3), 0], BigFloat(1)),
+  Engine.plane(BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(1)),
+  Engine.sphere(BigFloat[-1, 0, 0], BigFloat(1//5)),
+  Engine.sphere(BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//5)),
+  Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5))
+)
+
+# complete the gram matrix using gradient descent
+L, history = Engine.realize_gram(gram, guess, max_descent_steps = 200)
+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")