Compare commits
29 Commits
b185fd4b83
...
9d69a900e2
| Author | SHA1 | Date | |
|---|---|---|---|
|
9d69a900e2 | ||
|
|
8a77cd7484 | ||
|
|
a26f1e3927 | ||
|
|
19a4d49497 | ||
|
|
71c10adbdd | ||
|
|
33c09917d0 | ||
|
|
b24dcc9af8 | ||
|
|
b040bbb7fe | ||
|
|
9007c8bc7c | ||
|
|
a7f9545a37 | ||
|
|
3764fde2f6 | ||
|
|
24dae6807b | ||
|
|
74c7f64b0c | ||
|
|
d0340c0b65 | ||
|
|
69a704d414 | ||
|
|
01f44324c1 | ||
|
|
96ffc59642 | ||
|
|
a02b76544a | ||
|
|
6e719f9943 | ||
|
|
d51d43f481 | ||
|
|
6d233b5ee9 | ||
|
|
5abd4ca6e1 | ||
|
|
ea640f4861 | ||
|
|
4728959ae0 | ||
|
|
2038103d80 | ||
|
|
bde42ebac0 | ||
|
|
e6cf08a9b3 | ||
|
|
7c77481f5e | ||
|
|
1ce609836b |
@ -4,6 +4,8 @@ using Blink
|
|||||||
using Colors
|
using Colors
|
||||||
using Printf
|
using Printf
|
||||||
|
|
||||||
|
using Main.Engine
|
||||||
|
|
||||||
export ConstructionViewer, display!, opentools!, closetools!
|
export ConstructionViewer, display!, opentools!, closetools!
|
||||||
|
|
||||||
# === Blink utilities ===
|
# === Blink utilities ===
|
||||||
@ -133,7 +135,7 @@ mprod(v, w) =
|
|||||||
function display!(viewer::ConstructionViewer, elements::Matrix)
|
function display!(viewer::ConstructionViewer, elements::Matrix)
|
||||||
# load elements
|
# load elements
|
||||||
elements_full = []
|
elements_full = []
|
||||||
for elt in eachcol(elements)
|
for elt in eachcol(Engine.unmix * elements)
|
||||||
if mprod(elt, elt) < 0.5
|
if mprod(elt, elt) < 0.5
|
||||||
elt_full = [0; elt; fill(0, 26)]
|
elt_full = [0; elt; fill(0, 26)]
|
||||||
else
|
else
|
||||||
@ -205,14 +207,16 @@ end
|
|||||||
|
|
||||||
# ~~~ sandbox setup ~~~
|
# ~~~ sandbox setup ~~~
|
||||||
|
|
||||||
# in the default view, e4 + e5 is the point at infinity
|
elements = let
|
||||||
elements = sqrt(0.5) * BigFloat[
|
a = sqrt(BigFloat(3)/2)
|
||||||
1 1 -1 -1 0;
|
sqrt(0.5) * BigFloat[
|
||||||
1 -1 1 -1 0;
|
1 1 -1 -1 0
|
||||||
1 -1 -1 1 0;
|
1 -1 1 -1 0
|
||||||
0 0 0 0 -sqrt(6);
|
1 -1 -1 1 0
|
||||||
1 1 1 1 2
|
0.5 0.5 0.5 0.5 1+a
|
||||||
|
0.5 0.5 0.5 0.5 1-a
|
||||||
]
|
]
|
||||||
|
end
|
||||||
|
|
||||||
# show construction
|
# show construction
|
||||||
viewer = Viewer.ConstructionViewer()
|
viewer = Viewer.ConstructionViewer()
|
||||||
|
|||||||
@ -29,7 +29,7 @@ function rand_on_shell(rng::AbstractRNG, shell::T) where T <: Number
|
|||||||
space_part = rand_on_sphere(rng, T, 4)
|
space_part = rand_on_sphere(rng, T, 4)
|
||||||
rapidity = randn(rng, T)
|
rapidity = randn(rng, T)
|
||||||
sig = sign(shell)
|
sig = sign(shell)
|
||||||
[sconh(rapidity, sig)*space_part; sconh(rapidity, -sig)]
|
nullmix * [sconh(rapidity, sig)*space_part; sconh(rapidity, -sig)]
|
||||||
end
|
end
|
||||||
|
|
||||||
rand_on_shell(rng::AbstractRNG, shells::Array{T}) where T <: Number =
|
rand_on_shell(rng::AbstractRNG, shells::Array{T}) where T <: Number =
|
||||||
@ -39,21 +39,28 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she
|
|||||||
|
|
||||||
# === elements ===
|
# === elements ===
|
||||||
|
|
||||||
plane(normal, offset) = [normal; offset; offset]
|
point(pos) = [pos; 0.5; 0.5 * dot(pos, pos)]
|
||||||
|
|
||||||
|
plane(normal, offset) = [-normal; 0; -offset]
|
||||||
|
|
||||||
function sphere(center, radius)
|
function sphere(center, radius)
|
||||||
dist_sq = dot(center, center)
|
dist_sq = dot(center, center)
|
||||||
return [
|
[
|
||||||
center / radius;
|
center / radius;
|
||||||
0.5 * ((dist_sq - 1) / radius - radius);
|
0.5 / radius;
|
||||||
0.5 * ((dist_sq + 1) / radius - radius)
|
0.5 * (dist_sq / radius - radius)
|
||||||
]
|
]
|
||||||
end
|
end
|
||||||
|
|
||||||
# === Gram matrix realization ===
|
# === Gram matrix realization ===
|
||||||
|
|
||||||
|
# basis changes
|
||||||
|
nullmix = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [-1 1; 1 1]//2]
|
||||||
|
unmix = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [-1 1; 1 1]]
|
||||||
|
|
||||||
# the Lorentz form
|
# the Lorentz form
|
||||||
Q = diagm([1, 1, 1, 1, -1])
|
## [old] Q = diagm([1, 1, 1, 1, -1])
|
||||||
|
Q = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [0 -2; -2 0]]
|
||||||
|
|
||||||
# project a matrix onto the subspace of matrices whose entries vanish at the
|
# project a matrix onto the subspace of matrices whose entries vanish at the
|
||||||
# given indices
|
# given indices
|
||||||
@ -311,7 +318,8 @@ end
|
|||||||
# explicit entry of `gram`. use gradient descent starting from `guess`
|
# explicit entry of `gram`. use gradient descent starting from `guess`
|
||||||
function realize_gram(
|
function realize_gram(
|
||||||
gram::SparseMatrixCSC{T, <:Any},
|
gram::SparseMatrixCSC{T, <:Any},
|
||||||
guess::Matrix{T};
|
guess::Matrix{T},
|
||||||
|
frozen = nothing;
|
||||||
scaled_tol = 1e-30,
|
scaled_tol = 1e-30,
|
||||||
min_efficiency = 0.5,
|
min_efficiency = 0.5,
|
||||||
init_rate = 1.0,
|
init_rate = 1.0,
|
||||||
@ -336,6 +344,15 @@ function realize_gram(
|
|||||||
scale_adjustment = sqrt(T(length(constrained)))
|
scale_adjustment = sqrt(T(length(constrained)))
|
||||||
tol = scale_adjustment * scaled_tol
|
tol = scale_adjustment * scaled_tol
|
||||||
|
|
||||||
|
# list the un-frozen indices
|
||||||
|
has_frozen = !isnothing(frozen)
|
||||||
|
if has_frozen
|
||||||
|
is_unfrozen = fill(true, size(guess))
|
||||||
|
is_unfrozen[frozen] .= false
|
||||||
|
unfrozen = findall(is_unfrozen)
|
||||||
|
unfrozen_stacked = reshape(is_unfrozen, total_dim)
|
||||||
|
end
|
||||||
|
|
||||||
# initialize variables
|
# initialize variables
|
||||||
grad_rate = init_rate
|
grad_rate = init_rate
|
||||||
L = copy(guess)
|
L = copy(guess)
|
||||||
@ -371,7 +388,23 @@ function realize_gram(
|
|||||||
if min_eigval <= 0
|
if min_eigval <= 0
|
||||||
hess -= reg_scale * min_eigval * I
|
hess -= reg_scale * min_eigval * I
|
||||||
end
|
end
|
||||||
base_step = reshape(hess \ reshape(neg_grad, total_dim), dims)
|
|
||||||
|
# compute the Newton step
|
||||||
|
neg_grad_stacked = reshape(neg_grad, total_dim)
|
||||||
|
if has_frozen
|
||||||
|
hess = hess[unfrozen_stacked, unfrozen_stacked]
|
||||||
|
neg_grad_compressed = neg_grad_stacked[unfrozen_stacked]
|
||||||
|
else
|
||||||
|
neg_grad_compressed = neg_grad_stacked
|
||||||
|
end
|
||||||
|
base_step_compressed = hess \ neg_grad_compressed
|
||||||
|
if has_frozen
|
||||||
|
base_step_stacked = zeros(total_dim)
|
||||||
|
base_step_stacked[unfrozen_stacked] .= base_step_compressed
|
||||||
|
else
|
||||||
|
base_step_stacked = base_step_compressed
|
||||||
|
end
|
||||||
|
base_step = reshape(base_step_stacked, dims)
|
||||||
push!(history.base_step, base_step)
|
push!(history.base_step, base_step)
|
||||||
|
|
||||||
# store the current position, loss, and slope
|
# store the current position, loss, and slope
|
||||||
@ -412,7 +445,7 @@ function realize_gram(
|
|||||||
|
|
||||||
# return the factorization and its history
|
# return the factorization and its history
|
||||||
push!(history.scaled_loss, loss / scale_adjustment)
|
push!(history.scaled_loss, loss / scale_adjustment)
|
||||||
L, true, history
|
L, loss < tol, history
|
||||||
end
|
end
|
||||||
|
|
||||||
end
|
end
|
||||||
@ -1,6 +1,7 @@
|
|||||||
include("Engine.jl")
|
include("Engine.jl")
|
||||||
|
|
||||||
using SparseArrays
|
using SparseArrays
|
||||||
|
using Random
|
||||||
|
|
||||||
# initialize the partial gram matrix for a sphere inscribed in a regular
|
# initialize the partial gram matrix for a sphere inscribed in a regular
|
||||||
# tetrahedron
|
# tetrahedron
|
||||||
@ -10,23 +11,23 @@ values = BigFloat[]
|
|||||||
for j in 1:9
|
for j in 1:9
|
||||||
for k in 1:9
|
for k in 1:9
|
||||||
filled = false
|
filled = false
|
||||||
if j == k
|
if j == 9
|
||||||
push!(values, j < 9 ? 1 : 0)
|
if k <= 5 && k != 2
|
||||||
filled = true
|
|
||||||
elseif (j == 9)
|
|
||||||
if (k <= 5 && k != 2)
|
|
||||||
push!(values, 0)
|
push!(values, 0)
|
||||||
filled = true
|
filled = true
|
||||||
end
|
end
|
||||||
elseif (k == 9)
|
elseif k == 9
|
||||||
if (j <= 5 && j != 2)
|
if j <= 5 && j != 2
|
||||||
push!(values, 0)
|
push!(values, 0)
|
||||||
filled = true
|
filled = true
|
||||||
end
|
end
|
||||||
elseif (j == 1 || k == 1)
|
elseif j == k
|
||||||
|
push!(values, 1)
|
||||||
|
filled = true
|
||||||
|
elseif j == 1 || k == 1
|
||||||
push!(values, 0)
|
push!(values, 0)
|
||||||
filled = true
|
filled = true
|
||||||
elseif (j == 2 || k == 2)
|
elseif j == 2 || k == 2
|
||||||
push!(values, -1)
|
push!(values, -1)
|
||||||
filled = true
|
filled = true
|
||||||
end
|
end
|
||||||
@ -46,59 +47,26 @@ append!(values, fill(-0.5, 4))
|
|||||||
=#
|
=#
|
||||||
gram = sparse(J, K, values)
|
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
|
# set initial guess
|
||||||
|
Random.seed!(58271)
|
||||||
guess = hcat(
|
guess = hcat(
|
||||||
Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
|
Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
|
||||||
Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)),
|
Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
|
||||||
Engine.plane(BigFloat[1, 0, 0], BigFloat(1)),
|
Engine.plane(-BigFloat[1, 0, 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([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)) + 0.1*Engine.rand_on_shell([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)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
|
||||||
Engine.sphere(BigFloat[-1, 0, 0], BigFloat(1//5)),
|
Engine.sphere(BigFloat[-1, 0, 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
|
||||||
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)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
|
||||||
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)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
|
||||||
BigFloat[0, 0, 0, 1, 1]
|
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
|
# complete the gram matrix using Newton's method with backtracking
|
||||||
#=
|
L, success, history = Engine.realize_gram(gram, guess, frozen)
|
||||||
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)
|
|
||||||
completed_gram = L'*Engine.Q*L
|
completed_gram = L'*Engine.Q*L
|
||||||
println("Completed Gram matrix:\n")
|
println("Completed Gram matrix:\n")
|
||||||
display(completed_gram)
|
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
|
if success
|
||||||
println("\nTarget accuracy achieved!")
|
println("\nTarget accuracy achieved!")
|
||||||
else
|
else
|
||||||
|
|||||||
77
engine-proto/gram-test/irisawa-hexlet.jl
Normal file
77
engine-proto/gram-test/irisawa-hexlet.jl
Normal file
@ -0,0 +1,77 @@
|
|||||||
|
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
|
||||||
@ -23,12 +23,6 @@ end
|
|||||||
gram = sparse(J, K, values)
|
gram = sparse(J, K, values)
|
||||||
|
|
||||||
# set the independent variable
|
# 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
|
indep_val = -9//5
|
||||||
gram[6, 1] = BigFloat(indep_val)
|
gram[6, 1] = BigFloat(indep_val)
|
||||||
gram[1, 6] = gram[6, 1]
|
gram[1, 6] = gram[6, 1]
|
||||||
@ -36,30 +30,25 @@ gram[1, 6] = gram[6, 1]
|
|||||||
# in this initial guess, the mutual tangency condition is satisfied for spheres
|
# in this initial guess, the mutual tangency condition is satisfied for spheres
|
||||||
# 1 through 5
|
# 1 through 5
|
||||||
Random.seed!(50793)
|
Random.seed!(50793)
|
||||||
guess = hcat(
|
guess = let
|
||||||
|
a = sqrt(BigFloat(3)/2)
|
||||||
|
hcat(
|
||||||
sqrt(1/BigFloat(2)) * BigFloat[
|
sqrt(1/BigFloat(2)) * BigFloat[
|
||||||
1 1 -1 -1 0;
|
1 1 -1 -1 0
|
||||||
1 -1 1 -1 0;
|
1 -1 1 -1 0
|
||||||
1 -1 -1 1 0;
|
1 -1 -1 1 0
|
||||||
0 0 0 0 -sqrt(BigFloat(6));
|
0.5 0.5 0.5 0.5 1+a
|
||||||
1 1 1 1 2;
|
0.5 0.5 0.5 0.5 1-a
|
||||||
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
|
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
|
||||||
Engine.rand_on_shell(fill(BigFloat(-1), 2))
|
Engine.rand_on_shell(fill(BigFloat(-1), 2))
|
||||||
)
|
)
|
||||||
|
end
|
||||||
|
|
||||||
# complete the gram matrix
|
# complete the gram matrix using Newton's method with backtracking
|
||||||
#=
|
|
||||||
L, history = Engine.realize_gram_gradient(gram, guess, scaled_tol = 0.01)
|
|
||||||
L_pol, history_pol = Engine.realize_gram_newton(gram, L)
|
|
||||||
=#
|
|
||||||
L, success, history = Engine.realize_gram(gram, guess)
|
L, success, history = Engine.realize_gram(gram, guess)
|
||||||
completed_gram = L'*Engine.Q*L
|
completed_gram = L'*Engine.Q*L
|
||||||
println("Completed Gram matrix:\n")
|
println("Completed Gram matrix:\n")
|
||||||
display(completed_gram)
|
display(completed_gram)
|
||||||
#=
|
|
||||||
println("\nSteps: ", size(history.scaled_loss, 1), " + ", size(history_pol.scaled_loss, 1))
|
|
||||||
println("Loss: ", history_pol.scaled_loss[end], "\n")
|
|
||||||
=#
|
|
||||||
if success
|
if success
|
||||||
println("\nTarget accuracy achieved!")
|
println("\nTarget accuracy achieved!")
|
||||||
else
|
else
|
||||||
|
|||||||
@ -8,16 +8,32 @@ using Random
|
|||||||
J = Int64[]
|
J = Int64[]
|
||||||
K = Int64[]
|
K = Int64[]
|
||||||
values = BigFloat[]
|
values = BigFloat[]
|
||||||
for j in 1:5
|
for j in 1:6
|
||||||
for k in 1:5
|
for k in 1:6
|
||||||
push!(J, j)
|
filled = false
|
||||||
push!(K, k)
|
if j == 6
|
||||||
if j == k
|
if k <= 4
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
elseif k == 6
|
||||||
|
if j <= 4
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
elseif j == k
|
||||||
push!(values, 1)
|
push!(values, 1)
|
||||||
elseif (j <= 4 && k <= 4)
|
filled = true
|
||||||
|
elseif j <= 4 && k <= 4
|
||||||
push!(values, -1/BigFloat(3))
|
push!(values, -1/BigFloat(3))
|
||||||
|
filled = true
|
||||||
else
|
else
|
||||||
push!(values, -1)
|
push!(values, -1)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
if filled
|
||||||
|
push!(J, j)
|
||||||
|
push!(K, k)
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
@ -25,19 +41,20 @@ gram = sparse(J, K, values)
|
|||||||
|
|
||||||
# set initial guess
|
# set initial guess
|
||||||
Random.seed!(99230)
|
Random.seed!(99230)
|
||||||
guess = sqrt(1/BigFloat(3)) * BigFloat[
|
guess = hcat(
|
||||||
|
sqrt(1/BigFloat(3)) * BigFloat[
|
||||||
1 1 -1 -1 0
|
1 1 -1 -1 0
|
||||||
1 -1 1 -1 0
|
1 -1 1 -1 0
|
||||||
1 -1 -1 1 0
|
1 -1 -1 1 0
|
||||||
1 1 1 1 -2
|
0 0 0 0 1.5
|
||||||
1 1 1 1 1
|
1 1 1 1 -0.5
|
||||||
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5))
|
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
|
||||||
|
BigFloat[0, 0, 0, 0, 1]
|
||||||
|
)
|
||||||
|
frozen = [CartesianIndex(j, 6) for j in 1:5]
|
||||||
|
|
||||||
# complete the gram matrix
|
# complete the gram matrix using Newton's method with backtracking
|
||||||
#=
|
L, success, history = Engine.realize_gram(gram, guess, frozen)
|
||||||
L, history = Engine.realize_gram_newton(gram, guess)
|
|
||||||
=#
|
|
||||||
L, success, history = Engine.realize_gram(gram, guess)
|
|
||||||
completed_gram = L'*Engine.Q*L
|
completed_gram = L'*Engine.Q*L
|
||||||
println("Completed Gram matrix:\n")
|
println("Completed Gram matrix:\n")
|
||||||
display(completed_gram)
|
display(completed_gram)
|
||||||
|
|||||||
96
engine-proto/gram-test/tetrahedron-radius-ratio.jl
Normal file
96
engine-proto/gram-test/tetrahedron-radius-ratio.jl
Normal file
@ -0,0 +1,96 @@
|
|||||||
|
include("Engine.jl")
|
||||||
|
|
||||||
|
using LinearAlgebra
|
||||||
|
using SparseArrays
|
||||||
|
using Random
|
||||||
|
|
||||||
|
# initialize the partial gram matrix for a sphere inscribed in a regular
|
||||||
|
# tetrahedron
|
||||||
|
J = Int64[]
|
||||||
|
K = Int64[]
|
||||||
|
values = BigFloat[]
|
||||||
|
for j in 1:11
|
||||||
|
for k in 1:11
|
||||||
|
filled = false
|
||||||
|
if j == 11
|
||||||
|
if k <= 4
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
elseif k == 11
|
||||||
|
if j <= 4
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
elseif j == k
|
||||||
|
push!(values, j <= 6 ? 1 : 0)
|
||||||
|
filled = true
|
||||||
|
elseif j <= 4
|
||||||
|
if k <= 4
|
||||||
|
push!(values, -1/BigFloat(3))
|
||||||
|
filled = true
|
||||||
|
elseif k == 5
|
||||||
|
push!(values, -1)
|
||||||
|
filled = true
|
||||||
|
elseif 7 <= k <= 10 && k - j != 6
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
elseif k <= 4
|
||||||
|
if j == 5
|
||||||
|
push!(values, -1)
|
||||||
|
filled = true
|
||||||
|
elseif 7 <= j <= 10 && j - k != 6
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
elseif j == 6 && 7 <= k <= 10 || k == 6 && 7 <= j <= 10
|
||||||
|
push!(values, 0)
|
||||||
|
filled = true
|
||||||
|
end
|
||||||
|
if filled
|
||||||
|
push!(J, j)
|
||||||
|
push!(K, k)
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
gram = sparse(J, K, values)
|
||||||
|
|
||||||
|
# set initial guess
|
||||||
|
Random.seed!(99230)
|
||||||
|
guess = hcat(
|
||||||
|
sqrt(1/BigFloat(3)) * BigFloat[
|
||||||
|
1 1 -1 -1 0 0
|
||||||
|
1 -1 1 -1 0 0
|
||||||
|
1 -1 -1 1 0 0
|
||||||
|
0 0 0 0 1.5 0.5
|
||||||
|
1 1 1 1 -0.5 -1.5
|
||||||
|
] + 0.0*Engine.rand_on_shell(fill(BigFloat(-1), 6)),
|
||||||
|
Engine.point([-0.5, -0.5, -0.5] + 0.3*randn(3)),
|
||||||
|
Engine.point([-0.5, 0.5, 0.5] + 0.3*randn(3)),
|
||||||
|
Engine.point([ 0.5, -0.5, 0.5] + 0.3*randn(3)),
|
||||||
|
Engine.point([ 0.5, 0.5, -0.5] + 0.3*randn(3)),
|
||||||
|
BigFloat[0, 0, 0, 0, 1]
|
||||||
|
)
|
||||||
|
frozen = vcat(
|
||||||
|
[CartesianIndex(4, k) for k in 7:10],
|
||||||
|
[CartesianIndex(j, 11) for j in 1:5]
|
||||||
|
)
|
||||||
|
|
||||||
|
# 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])
|
||||||
|
if success
|
||||||
|
infty = BigFloat[0, 0, 0, 0, 1]
|
||||||
|
radius_ratio = dot(infty, Engine.Q * L[:,5]) / dot(infty, Engine.Q * L[:,6])
|
||||||
|
println("\nCircumradius / inradius: ", radius_ratio)
|
||||||
|
end
|
||||||
@ -6,22 +6,22 @@ These coordinates are of form $I=(c, r, x, y, z)$ where we think of $c$ as the c
|
|||||||
|
|
||||||
| Entity or Relationship | Representation | Comments/questions |
|
| Entity or Relationship | Representation | Comments/questions |
|
||||||
| ------------------------------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
| ------------------------------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
|
||||||
| Sphere s with radius r>0 centered on P = (x,y,z) | $I_s = (1/c, 1/r, x/r, y/r, z/r)$ satisfying $Q(I_s,I_s) = -1$, i.e., $c = r/(\|P\|^2 - r^2)$. | Can also write $I_s = (\|P\|^2/r - r, 1/r, x/r. y/r, z/r)$ -- so there is no trouble if $\|E_{I_s}\| = r$, just get first coordinate to be 0. |
|
| Sphere s with radius r>0 centered on P = (x,y,z) | $I_s = (1/c, 1/r, x/r, y/r, z/r)$ satisfying $Q(I_s,I_s) = -1$, i.e., $c = r/(\|P\|^2 - r^2)$. | Can also write $I_s = (\|P\|^2/r - r, 1/r, x/r. y/r, z/r)$—so there is no trouble if $\|E_{I_s}\| = r$, just get first coordinate to be 0. |
|
||||||
| Plane p with unit normal (x,y,z), a distance s from origin | $I_p = (2s, 0, x, y, z)$ | Note $Q(I_p, I_p)$ is still -1. Also, there are two representations for each plane through the origin, namely $(0,0,x,y,z)$ and $(0,0,-x,-y,-z)$ |
|
| Plane p with unit normal (x,y,z) through the point s(x,y,z) | $I_p = (-2s, 0, -x, -y, -z)$ | Note $Q(I_p, I_p)$ is still −1. |
|
||||||
| Point P with Euclidean coordinates (x,y,z) | $I_P = (\|P\|^2, 1, x, y, z)$ | Note $Q(I_P,I_P) = 0$. Because of this we might choose some other scaling of the inversive coordinates, say $(\||P\||,1/\||P\||,x/\||P\||,y/\||P\||,z/\||P\||)$ instead, but that fails at the origin, and likely won't have some of the other nice properties listed below. Note that scaling just the co-radius by $s$ and the radius by $1/s$ (which still preserves $Q=0$) dilates by a factor of $s$ about the origin, so that $(\|P\|, \|P\|, x, y, z)$, which might look symmetric, would actually have to represent the Euclidean point $(x/\||P\||, y/\||P\||, z/\||P\||)$ . |
|
| Point P with Euclidean coordinates (x,y,z) | $I_P = (\|P\|^2, 1, x, y, z)$ | Note $Q(I_P,I_P) = 0$. Because of this we might choose some other scaling of the inversive coordinates, say $(\||P\||,1/\||P\||,x/\||P\||,y/\||P\||,z/\||P\||)$ instead, but that fails at the origin, and likely won't have some of the other nice properties listed below. Note that scaling just the co-radius by $s$ and the radius by $1/s$ (which still preserves $Q=0$) dilates by a factor of $s$ about the origin, so that $(\|P\|, \|P\|, x, y, z)$, which might look symmetric, would actually have to represent the Euclidean point $(x/\||P\||, y/\||P\||, z/\||P\||)$ . |
|
||||||
| ∞, the "point at infinity" | $I_\infty = (1,0,0,0,0)$ | The only solution to $Q(I,I) = 0$ not covered by the above case. |
|
| ∞, the "point at infinity" | $I_\infty = (1,0,0,0,0)$ | The only solution to $Q(I,I) = 0$ not covered by the above case. |
|
||||||
| P lies on sphere or plane given by I | $Q(I_P, I) = 0$ | |
|
| P lies on sphere or plane given by I | $Q(I_P, I) = 0$ | |
|
||||||
| Sphere/planes represented by I and J are tangent | $Q(I,J) = 1$ (??, see note at right) | Seems as though this must be $Q(I,J) = \pm1$ ? For example, the $xy$ plane represented by (0,0,0,0,1) is tangent to the unit circle centered at (0,0,1) rep'd by (0,1,0,0,1), but their Q-product is -1. And in general you can reflect any sphere tangent to any plane through the plane and it should flip the sign of $Q(I,J)$, if I am not mistaken. |
|
| Sphere/planes represented by I and J are tangent | If $I$ and $J$ have the same orientation where they touch, $Q(I,J) = -1$. If they have opposing orientations, $Q(I,J) = 1$. | For example, the $xy$ plane with normal $-e_z$, represented by $(0,0,0,0,1)$, is tangent with matching orientation to the unit sphere centered at $(0,0,1)$ with outward normals, represented by $(0,1,0,0,1)$. Accordingly, their $Q$-product is −1. |
|
||||||
| Sphere/planes represented by I and J intersect (respectively, don't intersect) | $\|Q(I,J)\| < (\text{resp. }>)\; 1$ | Follows from the angle formula, at least conceptually. |
|
| Sphere/planes represented by I and J intersect (respectively, don't intersect) | $\|Q(I,J)\| < (\text{resp. }>)\; 1$ | Follows from the angle formula, at least conceptually. |
|
||||||
| P is center of sphere represented by I | Well, $Q(I_P, I)$ comes out to be $(\|P\|^2/r - r + \|P\|^2/r)/2 - \|P\|^2/r$ or just $-r/2$ . | Is it if and only if ? No this probably doesn't work because center is not conformal quantity. |
|
| P is center of sphere represented by I | Well, $Q(I_P, I)$ comes out to be $(\|P\|^2/r - r + \|P\|^2/r)/2 - \|P\|^2/r$ or just $-r/2$ . | Is it if and only if ? No this probably doesn't work because center is not conformal quantity. |
|
||||||
| Distance between P and R is d | $Q(I_P, I_R) = d^2/2$ | |
|
| Distance between P and R is d | $Q(I_P, I_R) = d^2/2$ | |
|
||||||
| Distance between P and sphere/plane rep by I | | In the very simple case of a plane $I$ rep'd by $(2s, 0, x, y, z)$ and a point $P$ that lies on its perpendicular through the origin, rep'd by $(r^2, 1, rx, ry, rz)$ we get $Q(I, I_p) = s-r$, which is indeed the signed distance between $I$ and $P$. Not sure if this generalizes to other combinations? |
|
| Distance between P and sphere/plane rep by I | | In the very simple case of a plane $I$ rep'd by $(2s, 0, x, y, z)$ and a point $P$ that lies on its perpendicular through the origin, rep'd by $(r^2, 1, rx, ry, rz)$ we get $Q(I, I_p) = s-r$, which is indeed the signed distance between $I$ and $P$. Not sure if this generalizes to other combinations? |
|
||||||
| Distance between sphere/planes rep by I and J | Note that for any two Euclidean-concentric spheres rep by $I$ and $J$ with radius $r$ and $s,$ $Q(I,J) = -\frac12\left(\frac rs + \frac sr\right)$ depends only on the ratio of $r$ and $s$. So this can't give something that determines the Euclidean distance between the two spheres, which presumably grows as the two spheres are blown up proportionally. For another example, for any two parallel planes, $Q(I,J) = \pm1$. | Alex had said: Q(I,J)=cosh^2 (d/2) maybe where d is distance in usual hyperbolic metric. Or maybe cosh d. That may be right depending on what's meant by the hyperbolic metric there, but it seems like it won't determine a reasonable Euclidean distance between planes, which should differ between different pairs of parallel planes. |
|
| Distance between sphere/planes rep by I and J | Note that for any two Euclidean-concentric spheres rep by $I$ and $J$ with radius $r$ and $s,$ $Q(I,J) = -\frac12\left(\frac rs + \frac sr\right)$ depends only on the ratio of $r$ and $s$. So this can't give something that determines the Euclidean distance between the two spheres, which presumably grows as the two spheres are blown up proportionally. For another example, for any two parallel planes, $Q(I,J) = \pm1$. | Alex had said: $Q(I,J)=\cosh(d/2)^2$ maybe where d is distance in usual hyperbolic metric. Or maybe $\cosh(d)$. That may be right depending on what's meant by the hyperbolic metric there, but it seems like it won't determine a reasonable Euclidean distance between planes, which should differ between different pairs of parallel planes. |
|
||||||
| Sphere centered on P through R | | Probably just calculate distance etc. |
|
| Sphere centered on P through R | | Probably just calculate distance etc. |
|
||||||
| Plane rep'd by I goes through center of sphere rep'd by J | I think this is equivalent to the plane being perpendicular to the sphere, i.e. $Q(I,J) = 0$. | |
|
| Plane rep'd by I goes through center of sphere rep'd by J | I think this is equivalent to the plane being perpendicular to the sphere, i.e. $Q(I,J) = 0$. | |
|
||||||
| Dihedral angle between planes (or spheres?) rep by I and J | $\theta = \arccos(Q(I,J))$ | Aaron Fenyes points out: The angle between spheres in $S^3$ matches the angle between the planes they bound in $R^{(1,4)}$, which matches the angle between the spacelike vectors perpendicular to those planes. So we should have $Q(I,J) = \cos\theta$. Note that when the spheres do not intersect, we can interpret this as the "imaginary angle" between them, via $\cosh t = \cos it$. |
|
| Dihedral angle between planes (or spheres?) rep by I and J | $\theta = \arccos(Q(I,J))$ | Aaron Fenyes points out: The angle between spheres in $S^3$ matches the angle between the planes they bound in $R^{(1,4)}$, which matches the angle between the spacelike vectors perpendicular to those planes. So we should have $Q(I,J) = \cos(\theta)$. Note that when the spheres do not intersect, we can interpret this as the "imaginary angle" between them, via $\cosh(t) = \cos(it)$. |
|
||||||
| R, P, S are collinear | Maybe just cross product of two differences is 0. Or, $R,P,S,\infty$ lie on a circle, or equivalently, $I_R,I_P,I_S,I_\infty$ span a plane (rather than a three-space). | Not a conformal property, but $R,P,S,\infty$ lying on a circle _is_. |
|
| R, P, S are collinear | Maybe just cross product of two differences is 0. Or, $R,P,S,\infty$ lie on a circle, or equivalently, $I_R,I_P,I_S,I_\infty$ span a plane (rather than a three-space). | $R,P,S$ lying on a line isn't a conformal property, but $R,P,S,\infty$ lying on a circle is. |
|
||||||
| Plane through noncollinear R, P, S | Should be, just solve Q(I, I_R) = 0 etc. | |
|
| Plane through noncollinear R, P, S | Should be, just solve $Q(I, I_R) = 0$ etc. | |
|
||||||
| circle | Maybe concentric sphere and the containing plane? Note it is easy to constrain the relationship between those two: they must be perpendicular. | Defn: circle is intersection of two spheres. That does cover lines. But you lose the canonicalness |
|
| circle | Maybe concentric sphere and the containing plane? Note it is easy to constrain the relationship between those two: they must be perpendicular. | Defn: circle is intersection of two spheres. That does cover lines. But you lose the canonicalness |
|
||||||
| line | Maybe two perpendicular containing planes? Maybe the plane perpendicular to the line and through origin, together with the point of the line on that plane? Or maybe just as a bag of collinear points? | The first is the limiting case of the possible circle rep, but it is not canonical. The second appears to be canonical, but I don't see a circle rep that corresponds to it. |
|
| line | Maybe two perpendicular containing planes? Maybe the plane perpendicular to the line and through origin, together with the point of the line on that plane? Or maybe just as a bag of collinear points? | The first is the limiting case of the possible circle rep, but it is not canonical. The second appears to be canonical, but I don't see a circle rep that corresponds to it. |
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user