dyna3/engine-proto/gram-test/irisawa-hexlet.jl
Aaron Fenyes 9d69a900e2 Irisawa hexlet: use Abe's terminology in comments
Abe uses the names "sun" and "moon" for what Wikipedia calls the nucleus
spheres.
2024-07-18 03:39:41 -07:00

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