Compare commits

...

29 Commits

Author SHA1 Message Date
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
Aaron Fenyes
8a77cd7484 Irisawa hexlet: drop unviable approach
The approach in the deleted file can't work, because the "sun" and
"moon" spheres can't be placed arbitrarily.
2024-07-18 03:21:46 -07:00
Aaron Fenyes
a26f1e3927 Add Irisawa hexlet example
Hat tip Romy, who sent me the article on sangaku that led me to this
problem.
2024-07-18 03:16:57 -07:00
Aaron Fenyes
19a4d49497 Clean up example formatting 2024-07-18 01:48:05 -07:00
Aaron Fenyes
71c10adbdd Overlapping pyramids: drop outdated comment 2024-07-18 01:12:49 -07:00
Aaron Fenyes
33c09917d0 Correct scope of guess constants 2024-07-18 01:05:13 -07:00
Aaron Fenyes
b24dcc9af8 Report success correctly when step limit is reached 2024-07-18 01:04:40 -07:00
Aaron Fenyes
b040bbb7fe Drop old code from examples 2024-07-18 00:50:48 -07:00
Aaron Fenyes
9007c8bc7c Circles in triangle: jiggle the guess 2024-07-18 00:49:09 -07:00
Aaron Fenyes
a7f9545a37 Circles in triangle: correct frozen variables
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.
2024-07-18 00:43:00 -07:00
Aaron Fenyes
3764fde2f6 Clean up formatting of notes 2024-07-18 00:27:10 -07:00
Aaron Fenyes
24dae6807b Clarify notes on tangency 2024-07-18 00:16:23 -07:00
Aaron Fenyes
74c7f64b0c Correct sign of normal in plane utility
Clarify the relevant notes too.
2024-07-18 00:03:12 -07:00
Aaron Fenyes
d0340c0b65 Correct point utility again
The balance between the light cone basis vectors was wrong, throwing the
point's coordinates off by a factor of two.
2024-07-17 23:37:28 -07:00
Aaron Fenyes
69a704d414 Use notes' sign convention for light cone basis 2024-07-17 23:07:34 -07:00
Aaron Fenyes
01f44324c1 Tetrahedron radius ratio: find radius ratio 2024-07-17 22:45:17 -07:00
Aaron Fenyes
96ffc59642 Tetrahedron radius ratio: tweak guess
Jiggle the vertex guesses. Put the circumscribed sphere guess on-shell.
2024-07-17 19:01:34 -07:00
Aaron Fenyes
a02b76544a Tetrahedron radius ratio: add circumscribed sphere 2024-07-17 18:55:36 -07:00
Aaron Fenyes
6e719f9943 Tetrahedron radius ratio: correct vertex guesses 2024-07-17 18:27:58 -07:00
Aaron Fenyes
d51d43f481 Correct point utility 2024-07-17 18:27:22 -07:00
Aaron Fenyes
6d233b5ee9 Tetrahedron radius ratio: correct signs 2024-07-17 18:08:36 -07:00
Aaron Fenyes
5abd4ca6e1 Revert "Give spheres positive radii in examples"
This reverts commit 4728959ae0, which
actually gave the spheres negative radii! I got confused by the sign
convention differences between the notes and the engine.
2024-07-17 17:49:43 -07:00
Aaron Fenyes
ea640f4861 Start tetrahedron radius ratio example
Add the vertices of the tetrahedron to the `sphere-in-tetrahedron`
example.
2024-07-17 17:33:32 -07:00
Aaron Fenyes
4728959ae0 Give spheres positive radii in examples
This changes the meaning of `indep_val` in the overlapping pyramids
example, so we adjust `indep_val` to get a nice-looking construction.
2024-07-17 17:22:33 -07:00
Aaron Fenyes
2038103d80 Write examples directly in light cone basis 2024-07-17 15:37:14 -07:00
Aaron Fenyes
bde42ebac0 Switch engine to light cone basis 2024-07-17 14:30:43 -07:00
Aaron Fenyes
e6cf08a9b3 Make tetrahedron faces planar 2024-07-15 23:54:59 -07:00
Aaron Fenyes
7c77481f5e Don't constrain self-product of frozen vector 2024-07-15 23:39:05 -07:00
Aaron Fenyes
1ce609836b Implement frozen variables 2024-07-15 22:11:54 -07:00
8 changed files with 312 additions and 128 deletions

View File

@ -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()

View File

@ -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

View File

@ -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
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 filled = true
elseif (j == 9) elseif j == 1 || k == 1
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 == 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

View 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

View File

@ -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
sqrt(1/BigFloat(2)) * BigFloat[ a = sqrt(BigFloat(3)/2)
1 1 -1 -1 0; hcat(
1 -1 1 -1 0; sqrt(1/BigFloat(2)) * BigFloat[
1 -1 -1 1 0; 1 1 -1 -1 0
0 0 0 0 -sqrt(BigFloat(6)); 1 -1 1 -1 0
1 1 1 1 2; 1 -1 -1 1 0
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), 0.5 0.5 0.5 0.5 1+a
Engine.rand_on_shell(fill(BigFloat(-1), 2)) 0.5 0.5 0.5 0.5 1-a
) ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
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

View File

@ -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(
1 1 -1 -1 0 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 -2 1 -1 -1 1 0
1 1 1 1 1 0 0 0 0 1.5
] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)) 1 1 1 1 -0.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)

View 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

View File

@ -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. |