dyna3/engine-proto/gram-test/irisawa-hexlet.jl

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