a7f9545a37
Since the self-product of the point at infinity is left unspecified, the first three components can vary without violating any constraints. To keep the point at infinity where it's supposed to be, we freeze all of its components.
109 lines
3.2 KiB
Julia
109 lines
3.2 KiB
Julia
include("Engine.jl")
|
|
|
|
using SparseArrays
|
|
|
|
# initialize the partial gram matrix for a sphere inscribed in a regular
|
|
# tetrahedron
|
|
J = Int64[]
|
|
K = Int64[]
|
|
values = BigFloat[]
|
|
for j in 1:9
|
|
for k in 1:9
|
|
filled = false
|
|
if j == 9
|
|
if (k <= 5 && k != 2)
|
|
push!(values, 0)
|
|
filled = true
|
|
end
|
|
elseif k == 9
|
|
if (j <= 5 && j != 2)
|
|
push!(values, 0)
|
|
filled = true
|
|
end
|
|
elseif 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))
|
|
#= make construction rigid
|
|
append!(J, [3, 4, 4, 5])
|
|
append!(K, [4, 3, 5, 4])
|
|
append!(values, fill(-0.5, 4))
|
|
=#
|
|
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)),
|
|
BigFloat[0, 0, 0, 0, 1]
|
|
)
|
|
frozen = [CartesianIndex(j, 9) for j in 1:5]
|
|
#=
|
|
guess = hcat(
|
|
Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
|
|
Engine.sphere(BigFloat[0, 0, 0], BigFloat(0.9)),
|
|
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(4//3*BigFloat[-1, 0, 0], BigFloat(1//3)),
|
|
Engine.sphere(4//3*BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//3)),
|
|
Engine.sphere(4//3*BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//3)),
|
|
BigFloat[0, 0, 0, 1, 1]
|
|
)
|
|
=#
|
|
|
|
# complete the gram matrix using gradient descent followed by Newton's method
|
|
#=
|
|
L, history = Engine.realize_gram_gradient(gram, guess, scaled_tol = 0.01)
|
|
L_pol, history_pol = Engine.realize_gram_newton(gram, L, rate = 0.3, scaled_tol = 1e-9)
|
|
L_pol2, history_pol2 = Engine.realize_gram_newton(gram, L_pol)
|
|
=#
|
|
L, success, history = Engine.realize_gram(gram, guess, frozen)
|
|
completed_gram = L'*Engine.Q*L
|
|
println("Completed Gram matrix:\n")
|
|
display(completed_gram)
|
|
#=
|
|
println(
|
|
"\nSteps: ",
|
|
size(history.scaled_loss, 1),
|
|
" + ", size(history_pol.scaled_loss, 1),
|
|
" + ", size(history_pol2.scaled_loss, 1)
|
|
)
|
|
println("Loss: ", history_pol2.scaled_loss[end], "\n")
|
|
=#
|
|
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") |