From 1ce609836bf45c16131bdec879c6c3fb613385af Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 15 Jul 2024 22:11:54 -0700 Subject: [PATCH 01/29] Implement frozen variables --- engine-proto/gram-test/Engine.jl | 30 +++++++++++++++++-- engine-proto/gram-test/circles-in-triangle.jl | 3 +- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 10f4b02..dccb514 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -311,7 +311,8 @@ end # explicit entry of `gram`. use gradient descent starting from `guess` function realize_gram( gram::SparseMatrixCSC{T, <:Any}, - guess::Matrix{T}; + guess::Matrix{T}, + frozen = nothing; scaled_tol = 1e-30, min_efficiency = 0.5, init_rate = 1.0, @@ -336,6 +337,15 @@ function realize_gram( scale_adjustment = sqrt(T(length(constrained))) 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 grad_rate = init_rate L = copy(guess) @@ -371,7 +381,23 @@ function realize_gram( if min_eigval <= 0 hess -= reg_scale * min_eigval * I 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) # store the current position, loss, and slope diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index 12975c2..b173ed9 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -66,6 +66,7 @@ guess = hcat( Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5)), BigFloat[0, 0, 0, 1, 1] ) +frozen = [CartesianIndex(j, 9) for j in 4:5] #= guess = hcat( Engine.plane(BigFloat[0, 0, 1], BigFloat(0)), @@ -86,7 +87,7 @@ 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) +L, success, history = Engine.realize_gram(gram, guess, frozen) completed_gram = L'*Engine.Q*L println("Completed Gram matrix:\n") display(completed_gram) From 7c77481f5e9990bc5a9723b64a7bd7532943ce4d Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 15 Jul 2024 23:39:05 -0700 Subject: [PATCH 02/29] Don't constrain self-product of frozen vector --- engine-proto/gram-test/circles-in-triangle.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index b173ed9..dd01ebf 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -10,19 +10,19 @@ values = BigFloat[] for j in 1:9 for k in 1:9 filled = false - if j == k - push!(values, j < 9 ? 1 : 0) - filled = true - elseif (j == 9) + if j == 9 if (k <= 5 && k != 2) push!(values, 0) filled = true end - elseif (k == 9) + 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 From e6cf08a9b3535bbb2c2b364f669e76534a123dce Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 15 Jul 2024 23:54:59 -0700 Subject: [PATCH 03/29] Make tetrahedron faces planar --- .../gram-test/sphere-in-tetrahedron.jl | 46 +++++++++++++------ 1 file changed, 33 insertions(+), 13 deletions(-) diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 273d3ad..1d02d45 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -8,16 +8,32 @@ using Random J = Int64[] K = Int64[] values = BigFloat[] -for j in 1:5 - for k in 1:5 - push!(J, j) - push!(K, k) - if j == k +for j in 1:6 + for k in 1:6 + filled = false + if j == 6 + 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) + filled = true elseif (j <= 4 && k <= 4) push!(values, -1/BigFloat(3)) + filled = true else push!(values, -1) + filled = true + end + if filled + push!(J, j) + push!(K, k) end end end @@ -25,19 +41,23 @@ gram = sparse(J, K, values) # set initial guess Random.seed!(99230) -guess = sqrt(1/BigFloat(3)) * BigFloat[ - 1 1 -1 -1 0 - 1 -1 1 -1 0 - 1 -1 -1 1 0 - 1 1 1 1 -2 - 1 1 1 1 1 -] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)) +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 -2 + 1 1 1 1 1 + ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), + BigFloat[0, 0, 0, 1, 1] +) +frozen = [CartesianIndex(j, 6) for j in 1:5] # complete the gram matrix #= L, history = Engine.realize_gram_newton(gram, guess) =# -L, success, history = Engine.realize_gram(gram, guess) +L, success, history = Engine.realize_gram(gram, guess, frozen) completed_gram = L'*Engine.Q*L println("Completed Gram matrix:\n") display(completed_gram) From bde42ebac0548d1665de75af7722301229f95d24 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 14:30:43 -0700 Subject: [PATCH 04/29] Switch engine to light cone basis --- engine-proto/ConstructionViewer.jl | 6 ++++-- engine-proto/gram-test/Engine.jl | 14 ++++++++++++-- engine-proto/gram-test/circles-in-triangle.jl | 2 +- engine-proto/gram-test/overlapping-pyramids.jl | 2 +- engine-proto/gram-test/sphere-in-tetrahedron.jl | 2 +- 5 files changed, 19 insertions(+), 7 deletions(-) diff --git a/engine-proto/ConstructionViewer.jl b/engine-proto/ConstructionViewer.jl index 29af212..c4845fa 100644 --- a/engine-proto/ConstructionViewer.jl +++ b/engine-proto/ConstructionViewer.jl @@ -4,6 +4,8 @@ using Blink using Colors using Printf +using Main.Engine + export ConstructionViewer, display!, opentools!, closetools! # === Blink utilities === @@ -133,7 +135,7 @@ mprod(v, w) = function display!(viewer::ConstructionViewer, elements::Matrix) # load elements elements_full = [] - for elt in eachcol(elements) + for elt in eachcol(Engine.unmix * elements) if mprod(elt, elt) < 0.5 elt_full = [0; elt; fill(0, 26)] else @@ -206,7 +208,7 @@ end # ~~~ sandbox setup ~~~ # in the default view, e4 + e5 is the point at infinity -elements = sqrt(0.5) * BigFloat[ +elements = Engine.nullmix * sqrt(0.5) * BigFloat[ 1 1 -1 -1 0; 1 -1 1 -1 0; 1 -1 -1 1 0; diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index dccb514..93102f8 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -39,11 +39,16 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she # === elements === +## [temp] in light cone coordinates +point(pos) = [pos; 1; dot(pos, pos)] + +## [temp] in standard coordinates plane(normal, offset) = [normal; offset; offset] +## [temp] in standard coordinates function sphere(center, radius) dist_sq = dot(center, center) - return [ + [ center / radius; 0.5 * ((dist_sq - 1) / radius - radius); 0.5 * ((dist_sq + 1) / radius - radius) @@ -52,8 +57,13 @@ end # === 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 -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 # given indices diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index dd01ebf..ec6e1e0 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -55,7 +55,7 @@ gram = sparse(J, K, values) ## guess = Engine.rand_on_shell(fill(BigFloat(-1), 8)) # set initial guess -guess = hcat( +guess = Engine.nullmix * 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)), diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index ee4c1fc..706a334 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -36,7 +36,7 @@ gram[1, 6] = gram[6, 1] # in this initial guess, the mutual tangency condition is satisfied for spheres # 1 through 5 Random.seed!(50793) -guess = hcat( +guess = Engine.nullmix * hcat( sqrt(1/BigFloat(2)) * BigFloat[ 1 1 -1 -1 0; 1 -1 1 -1 0; diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 1d02d45..6bfb2c0 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -41,7 +41,7 @@ gram = sparse(J, K, values) # set initial guess Random.seed!(99230) -guess = hcat( +guess = Engine.nullmix * hcat( sqrt(1/BigFloat(3)) * BigFloat[ 1 1 -1 -1 0 1 -1 1 -1 0 From 2038103d805dce9e39f61bac3b2e3cfbf1b11c2c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 15:37:14 -0700 Subject: [PATCH 05/29] Write examples directly in light cone basis --- engine-proto/ConstructionViewer.jl | 18 ++++++++------- engine-proto/gram-test/Engine.jl | 11 ++++----- engine-proto/gram-test/circles-in-triangle.jl | 4 ++-- .../gram-test/overlapping-pyramids.jl | 23 +++++++++++-------- .../gram-test/sphere-in-tetrahedron.jl | 8 +++---- 5 files changed, 33 insertions(+), 31 deletions(-) diff --git a/engine-proto/ConstructionViewer.jl b/engine-proto/ConstructionViewer.jl index c4845fa..8cfa632 100644 --- a/engine-proto/ConstructionViewer.jl +++ b/engine-proto/ConstructionViewer.jl @@ -207,14 +207,16 @@ end # ~~~ sandbox setup ~~~ -# in the default view, e4 + e5 is the point at infinity -elements = Engine.nullmix * sqrt(0.5) * BigFloat[ - 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 2 -] +elements = begin + const a = sqrt(BigFloat(3)/2) + sqrt(0.5) * BigFloat[ + 1 1 -1 -1 0 + 1 -1 1 -1 0 + 1 -1 -1 1 0 + -0.5 -0.5 -0.5 -0.5 -a-1 + 0.5 0.5 0.5 0.5 -a+1 + ] +end # show construction viewer = Viewer.ConstructionViewer() diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 93102f8..920b043 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -29,7 +29,7 @@ function rand_on_shell(rng::AbstractRNG, shell::T) where T <: Number space_part = rand_on_sphere(rng, T, 4) rapidity = randn(rng, T) sig = sign(shell) - [sconh(rapidity, sig)*space_part; sconh(rapidity, -sig)] + nullmix * [sconh(rapidity, sig)*space_part; sconh(rapidity, -sig)] end rand_on_shell(rng::AbstractRNG, shells::Array{T}) where T <: Number = @@ -39,19 +39,16 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she # === elements === -## [temp] in light cone coordinates point(pos) = [pos; 1; dot(pos, pos)] -## [temp] in standard coordinates -plane(normal, offset) = [normal; offset; offset] +plane(normal, offset) = [normal; 0; offset] -## [temp] in standard coordinates function sphere(center, radius) dist_sq = dot(center, center) [ center / radius; - 0.5 * ((dist_sq - 1) / radius - radius); - 0.5 * ((dist_sq + 1) / radius - radius) + -0.5 / radius; + 0.5 * (dist_sq / radius - radius) ] end diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index ec6e1e0..fc5e13d 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -55,7 +55,7 @@ gram = sparse(J, K, values) ## guess = Engine.rand_on_shell(fill(BigFloat(-1), 8)) # set initial guess -guess = Engine.nullmix * hcat( +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)), @@ -64,7 +64,7 @@ guess = Engine.nullmix * hcat( 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, 1, 1] + BigFloat[0, 0, 0, 0, 1] ) frozen = [CartesianIndex(j, 9) for j in 4:5] #= diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index 706a334..8edb981 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -36,16 +36,19 @@ gram[1, 6] = gram[6, 1] # in this initial guess, the mutual tangency condition is satisfied for spheres # 1 through 5 Random.seed!(50793) -guess = Engine.nullmix * hcat( - sqrt(1/BigFloat(2)) * BigFloat[ - 1 1 -1 -1 0; - 1 -1 1 -1 0; - 1 -1 -1 1 0; - 0 0 0 0 -sqrt(BigFloat(6)); - 1 1 1 1 2; - ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), - Engine.rand_on_shell(fill(BigFloat(-1), 2)) -) +guess = begin + const a = sqrt(BigFloat(3)/2) + hcat( + sqrt(1/BigFloat(2)) * BigFloat[ + 1 1 -1 -1 0 + 1 -1 1 -1 0 + 1 -1 -1 1 0 + -0.5 -0.5 -0.5 -0.5 -a-1 + 0.5 0.5 0.5 0.5 -a+1 + ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), + Engine.rand_on_shell(fill(BigFloat(-1), 2)) + ) +end # complete the gram matrix #= diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 6bfb2c0..1c0dda8 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -41,15 +41,15 @@ gram = sparse(J, K, values) # set initial guess Random.seed!(99230) -guess = Engine.nullmix * hcat( +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 -2 - 1 1 1 1 1 + 0 0 0 0 -1.5 + 1 1 1 1 -0.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), - BigFloat[0, 0, 0, 1, 1] + BigFloat[0, 0, 0, 0, 1] ) frozen = [CartesianIndex(j, 6) for j in 1:5] From 4728959ae049f728c204e206236865427bad762c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 17:22:33 -0700 Subject: [PATCH 06/29] 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. --- engine-proto/ConstructionViewer.jl | 4 ++-- engine-proto/gram-test/overlapping-pyramids.jl | 12 +++--------- engine-proto/gram-test/sphere-in-tetrahedron.jl | 6 +++--- 3 files changed, 8 insertions(+), 14 deletions(-) diff --git a/engine-proto/ConstructionViewer.jl b/engine-proto/ConstructionViewer.jl index 8cfa632..0ce6a6d 100644 --- a/engine-proto/ConstructionViewer.jl +++ b/engine-proto/ConstructionViewer.jl @@ -213,8 +213,8 @@ elements = begin 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - -0.5 -0.5 -0.5 -0.5 -a-1 - 0.5 0.5 0.5 0.5 -a+1 + 0.5 0.5 0.5 0.5 a+1 + -0.5 -0.5 -0.5 -0.5 a-1 ] end diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index 8edb981..757a18c 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -23,13 +23,7 @@ end gram = sparse(J, K, values) # 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 = 2//5 gram[6, 1] = BigFloat(indep_val) gram[1, 6] = gram[6, 1] @@ -43,8 +37,8 @@ guess = begin 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - -0.5 -0.5 -0.5 -0.5 -a-1 - 0.5 0.5 0.5 0.5 -a+1 + 0.5 0.5 0.5 0.5 a+1 + -0.5 -0.5 -0.5 -0.5 a-1 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), Engine.rand_on_shell(fill(BigFloat(-1), 2)) ) diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 1c0dda8..d703321 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -28,7 +28,7 @@ for j in 1:6 push!(values, -1/BigFloat(3)) filled = true else - push!(values, -1) + push!(values, 1) filled = true end if filled @@ -46,8 +46,8 @@ guess = hcat( 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - 0 0 0 0 -1.5 - 1 1 1 1 -0.5 + 0 0 0 0 1.5 + 1 1 1 1 0.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), BigFloat[0, 0, 0, 0, 1] ) From ea640f48617b3423edb70344cf504fbdb2b16c1e Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 17:33:32 -0700 Subject: [PATCH 07/29] Start tetrahedron radius ratio example Add the vertices of the tetrahedron to the `sphere-in-tetrahedron` example. --- .../gram-test/tetrahedron-radius-ratio.jl | 87 +++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 engine-proto/gram-test/tetrahedron-radius-ratio.jl diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl new file mode 100644 index 0000000..bde7272 --- /dev/null +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -0,0 +1,87 @@ +include("Engine.jl") + +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:10 + for k in 1:10 + filled = false + if j == 10 + if k <= 4 + push!(values, 0) + filled = true + end + elseif k == 10 + if j <= 4 + push!(values, 0) + filled = true + end + elseif j == k + push!(values, j <= 5 ? 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 k <= 9 && k - j != 5 + push!(values, 0) + filled = true + end + elseif k <= 4 + if j == 5 + push!(values, 1) + filled = true + elseif j <= 9 && j - k != 5 + push!(values, 0) + filled = true + end + 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 + 1 -1 1 -1 0 + 1 -1 -1 1 0 + 0 0 0 0 1.5 + 1 1 1 1 0.5 + ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), + Engine.point([-1, -1, -1]), + Engine.point([ 1, -1, -1]), + Engine.point([-1, 1, -1]), + Engine.point([ 1, -1, 1]), + BigFloat[0, 0, 0, 0, 1] +) +frozen = vcat( + [CartesianIndex(4, k) for k in 6:9], + [CartesianIndex(j, 10) for j in 1:5] +) + +# complete the gram matrix +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") \ No newline at end of file From 5abd4ca6e16b1efc2c1ca8ad3d1610c920c41c5a Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 17:49:43 -0700 Subject: [PATCH 08/29] Revert "Give spheres positive radii in examples" This reverts commit 4728959ae049f728c204e206236865427bad762c, which actually gave the spheres negative radii! I got confused by the sign convention differences between the notes and the engine. --- engine-proto/ConstructionViewer.jl | 4 ++-- engine-proto/gram-test/overlapping-pyramids.jl | 12 +++++++++--- engine-proto/gram-test/sphere-in-tetrahedron.jl | 6 +++--- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/engine-proto/ConstructionViewer.jl b/engine-proto/ConstructionViewer.jl index 0ce6a6d..8cfa632 100644 --- a/engine-proto/ConstructionViewer.jl +++ b/engine-proto/ConstructionViewer.jl @@ -213,8 +213,8 @@ elements = begin 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - 0.5 0.5 0.5 0.5 a+1 - -0.5 -0.5 -0.5 -0.5 a-1 + -0.5 -0.5 -0.5 -0.5 -a-1 + 0.5 0.5 0.5 0.5 -a+1 ] end diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index 757a18c..8edb981 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -23,7 +23,13 @@ end gram = sparse(J, K, values) # set the independent variable -indep_val = 2//5 +# +# 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 gram[6, 1] = BigFloat(indep_val) gram[1, 6] = gram[6, 1] @@ -37,8 +43,8 @@ guess = begin 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - 0.5 0.5 0.5 0.5 a+1 - -0.5 -0.5 -0.5 -0.5 a-1 + -0.5 -0.5 -0.5 -0.5 -a-1 + 0.5 0.5 0.5 0.5 -a+1 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), Engine.rand_on_shell(fill(BigFloat(-1), 2)) ) diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index d703321..1c0dda8 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -28,7 +28,7 @@ for j in 1:6 push!(values, -1/BigFloat(3)) filled = true else - push!(values, 1) + push!(values, -1) filled = true end if filled @@ -46,8 +46,8 @@ guess = hcat( 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - 0 0 0 0 1.5 - 1 1 1 1 0.5 + 0 0 0 0 -1.5 + 1 1 1 1 -0.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), BigFloat[0, 0, 0, 0, 1] ) From 6d233b5ee90c32cbc83f8e242ddbdc0288abec9f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 18:08:36 -0700 Subject: [PATCH 09/29] Tetrahedron radius ratio: correct signs --- engine-proto/gram-test/tetrahedron-radius-ratio.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index bde7272..6172aa3 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -29,7 +29,7 @@ for j in 1:10 push!(values, -1/BigFloat(3)) filled = true elseif k == 5 - push!(values, 1) + push!(values, -1) filled = true elseif k <= 9 && k - j != 5 push!(values, 0) @@ -37,7 +37,7 @@ for j in 1:10 end elseif k <= 4 if j == 5 - push!(values, 1) + push!(values, -1) filled = true elseif j <= 9 && j - k != 5 push!(values, 0) @@ -59,8 +59,8 @@ guess = hcat( 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - 0 0 0 0 1.5 - 1 1 1 1 0.5 + 0 0 0 0 -1.5 + 1 1 1 1 -0.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), Engine.point([-1, -1, -1]), Engine.point([ 1, -1, -1]), From d51d43f481e9819b2a312e82da13235a6dd638ad Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 18:27:22 -0700 Subject: [PATCH 10/29] Correct point utility --- engine-proto/gram-test/Engine.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 920b043..bedda00 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -39,7 +39,7 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she # === elements === -point(pos) = [pos; 1; dot(pos, pos)] +point(pos) = [pos; -1; 0.25 * dot(pos, pos)] plane(normal, offset) = [normal; 0; offset] From 6e719f9943754b1b120e728738390b26752c54c7 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 18:27:58 -0700 Subject: [PATCH 11/29] Tetrahedron radius ratio: correct vertex guesses --- engine-proto/gram-test/tetrahedron-radius-ratio.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index 6172aa3..bb89268 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -63,9 +63,9 @@ guess = hcat( 1 1 1 1 -0.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), Engine.point([-1, -1, -1]), - Engine.point([ 1, -1, -1]), - Engine.point([-1, 1, -1]), + Engine.point([-1, 1, 1]), Engine.point([ 1, -1, 1]), + Engine.point([ 1, 1, -1]), BigFloat[0, 0, 0, 0, 1] ) frozen = vcat( From a02b76544a9930d09a86d36e0ea4560c4f339c45 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 18:55:36 -0700 Subject: [PATCH 12/29] Tetrahedron radius ratio: add circumscribed sphere --- .../gram-test/tetrahedron-radius-ratio.jl | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index bb89268..957f031 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -8,21 +8,21 @@ using Random J = Int64[] K = Int64[] values = BigFloat[] -for j in 1:10 - for k in 1:10 +for j in 1:11 + for k in 1:11 filled = false - if j == 10 + if j == 11 if k <= 4 push!(values, 0) filled = true end - elseif k == 10 + elseif k == 11 if j <= 4 push!(values, 0) filled = true end elseif j == k - push!(values, j <= 5 ? 1 : 0) + push!(values, j <= 6 ? 1 : 0) filled = true elseif j <= 4 if k <= 4 @@ -31,7 +31,7 @@ for j in 1:10 elseif k == 5 push!(values, -1) filled = true - elseif k <= 9 && k - j != 5 + elseif 7 <= k <= 10 && k - j != 6 push!(values, 0) filled = true end @@ -39,10 +39,13 @@ for j in 1:10 if j == 5 push!(values, -1) filled = true - elseif j <= 9 && j - k != 5 + 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) @@ -56,12 +59,12 @@ gram = sparse(J, K, values) Random.seed!(99230) guess = hcat( sqrt(1/BigFloat(3)) * BigFloat[ - 1 1 -1 -1 0 - 1 -1 1 -1 0 - 1 -1 -1 1 0 - 0 0 0 0 -1.5 - 1 1 1 1 -0.5 - ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), + 1 1 -1 -1 0 0 + 1 -1 1 -1 0 0 + 1 -1 -1 1 0 0 + 0 0 0 0 -1.5 -3 + 1 1 1 1 -0.5 -1 + ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 6)), Engine.point([-1, -1, -1]), Engine.point([-1, 1, 1]), Engine.point([ 1, -1, 1]), @@ -69,8 +72,8 @@ guess = hcat( BigFloat[0, 0, 0, 0, 1] ) frozen = vcat( - [CartesianIndex(4, k) for k in 6:9], - [CartesianIndex(j, 10) for j in 1:5] + [CartesianIndex(4, k) for k in 7:10], + [CartesianIndex(j, 11) for j in 1:5] ) # complete the gram matrix From 96ffc59642ca656193193ba6ec9d87e12cb0036f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 19:01:34 -0700 Subject: [PATCH 13/29] Tetrahedron radius ratio: tweak guess Jiggle the vertex guesses. Put the circumscribed sphere guess on-shell. --- engine-proto/gram-test/tetrahedron-radius-ratio.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index 957f031..9e79c05 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -62,13 +62,13 @@ guess = hcat( 1 1 -1 -1 0 0 1 -1 1 -1 0 0 1 -1 -1 1 0 0 - 0 0 0 0 -1.5 -3 - 1 1 1 1 -0.5 -1 + 0 0 0 0 -1.5 -0.5 + 1 1 1 1 -0.5 -1.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 6)), - Engine.point([-1, -1, -1]), - Engine.point([-1, 1, 1]), - Engine.point([ 1, -1, 1]), - Engine.point([ 1, 1, -1]), + Engine.point([-1, -1, -1] + 0.3*randn(3)), + Engine.point([-1, 1, 1] + 0.3*randn(3)), + Engine.point([ 1, -1, 1] + 0.3*randn(3)), + Engine.point([ 1, 1, -1] + 0.3*randn(3)), BigFloat[0, 0, 0, 0, 1] ) frozen = vcat( From 01f44324c12108fd04ff33577a32eb3d770dc4a5 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 22:45:17 -0700 Subject: [PATCH 14/29] Tetrahedron radius ratio: find radius ratio --- engine-proto/gram-test/tetrahedron-radius-ratio.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index 9e79c05..4218cb7 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -1,5 +1,6 @@ include("Engine.jl") +using LinearAlgebra using SparseArrays using Random @@ -87,4 +88,9 @@ else println("\nFailed to reach target accuracy") end println("Steps: ", size(history.scaled_loss, 1)) -println("Loss: ", history.scaled_loss[end], "\n") \ No newline at end of file +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 \ No newline at end of file From 69a704d4145d581ec0c3599a05a1ae88b42698af Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 23:07:34 -0700 Subject: [PATCH 15/29] Use notes' sign convention for light cone basis --- engine-proto/ConstructionViewer.jl | 4 ++-- engine-proto/gram-test/Engine.jl | 10 +++++----- engine-proto/gram-test/overlapping-pyramids.jl | 4 ++-- engine-proto/gram-test/sphere-in-tetrahedron.jl | 2 +- engine-proto/gram-test/tetrahedron-radius-ratio.jl | 2 +- 5 files changed, 11 insertions(+), 11 deletions(-) diff --git a/engine-proto/ConstructionViewer.jl b/engine-proto/ConstructionViewer.jl index 8cfa632..c9b0b7a 100644 --- a/engine-proto/ConstructionViewer.jl +++ b/engine-proto/ConstructionViewer.jl @@ -213,8 +213,8 @@ elements = begin 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - -0.5 -0.5 -0.5 -0.5 -a-1 - 0.5 0.5 0.5 0.5 -a+1 + 0.5 0.5 0.5 0.5 1+a + 0.5 0.5 0.5 0.5 1-a ] end diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index bedda00..73c5c31 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -39,7 +39,7 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she # === elements === -point(pos) = [pos; -1; 0.25 * dot(pos, pos)] +point(pos) = [pos; 1; 0.25 * dot(pos, pos)] plane(normal, offset) = [normal; 0; offset] @@ -47,7 +47,7 @@ function sphere(center, radius) dist_sq = dot(center, center) [ center / radius; - -0.5 / radius; + 0.5 / radius; 0.5 * (dist_sq / radius - radius) ] end @@ -55,12 +55,12 @@ end # === 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]] +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 ## [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]] +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 # given indices diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index 8edb981..0d1f018 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -43,8 +43,8 @@ guess = begin 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - -0.5 -0.5 -0.5 -0.5 -a-1 - 0.5 0.5 0.5 0.5 -a+1 + 0.5 0.5 0.5 0.5 1+a + 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)) ) diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 1c0dda8..631f0e5 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -46,7 +46,7 @@ guess = hcat( 1 1 -1 -1 0 1 -1 1 -1 0 1 -1 -1 1 0 - 0 0 0 0 -1.5 + 0 0 0 0 1.5 1 1 1 1 -0.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)), BigFloat[0, 0, 0, 0, 1] diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index 4218cb7..ed3ceb0 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -63,7 +63,7 @@ guess = hcat( 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 + 0 0 0 0 1.5 0.5 1 1 1 1 -0.5 -1.5 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 6)), Engine.point([-1, -1, -1] + 0.3*randn(3)), From d0340c0b658b428eeb6293edf1006c7ce9d7f093 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 17 Jul 2024 23:37:28 -0700 Subject: [PATCH 16/29] 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. --- engine-proto/gram-test/Engine.jl | 2 +- engine-proto/gram-test/tetrahedron-radius-ratio.jl | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 73c5c31..2662a17 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -39,7 +39,7 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she # === elements === -point(pos) = [pos; 1; 0.25 * dot(pos, pos)] +point(pos) = [pos; 0.5; 0.5 * dot(pos, pos)] plane(normal, offset) = [normal; 0; offset] diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index ed3ceb0..c284078 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -65,11 +65,11 @@ guess = hcat( 1 -1 -1 1 0 0 0 0 0 0 1.5 0.5 1 1 1 1 -0.5 -1.5 - ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 6)), - Engine.point([-1, -1, -1] + 0.3*randn(3)), - Engine.point([-1, 1, 1] + 0.3*randn(3)), - Engine.point([ 1, -1, 1] + 0.3*randn(3)), - Engine.point([ 1, 1, -1] + 0.3*randn(3)), + ] + 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( From 74c7f64b0c708b013db0175f0e729c8f3fc2bbc8 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 00:03:12 -0700 Subject: [PATCH 17/29] Correct sign of normal in plane utility Clarify the relevant notes too. --- engine-proto/gram-test/Engine.jl | 2 +- engine-proto/gram-test/circles-in-triangle.jl | 6 +++--- notes/inversive.md | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 2662a17..40b77b0 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -41,7 +41,7 @@ rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), she point(pos) = [pos; 0.5; 0.5 * dot(pos, pos)] -plane(normal, offset) = [normal; 0; offset] +plane(normal, offset) = [-normal; 0; -offset] function sphere(center, radius) dist_sq = dot(center, center) diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index fc5e13d..ca49574 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -58,9 +58,9 @@ gram = sparse(J, K, values) 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.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)), diff --git a/notes/inversive.md b/notes/inversive.md index 6de7ef2..9cb2e19 100644 --- a/notes/inversive.md +++ b/notes/inversive.md @@ -7,7 +7,7 @@ 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 | | ------------------------------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | | 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\||)$ . | | ∞, 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$ | | From 24dae6807bf9799d39a94088d02f7980d826d057 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 00:16:23 -0700 Subject: [PATCH 18/29] Clarify notes on tangency --- notes/inversive.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notes/inversive.md b/notes/inversive.md index 9cb2e19..acac4cd 100644 --- a/notes/inversive.md +++ b/notes/inversive.md @@ -11,7 +11,7 @@ These coordinates are of form $I=(c, r, x, y, z)$ where we think of $c$ as the c | 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. | | 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. | | 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$ | | From 3764fde2f6ca75cd56686030d019b11c44182e12 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 00:27:10 -0700 Subject: [PATCH 19/29] Clean up formatting of notes --- notes/inversive.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/notes/inversive.md b/notes/inversive.md index acac4cd..c845e3d 100644 --- a/notes/inversive.md +++ b/notes/inversive.md @@ -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 | | ------------------------------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| 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) 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\||)$ . | | ∞, 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$ | | -| 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 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. | | 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 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. | -| 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$. | -| 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_. | -| Plane through noncollinear R, P, S | Should be, just solve Q(I, I_R) = 0 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$. | | +| 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). | $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. | | | 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. | From a7f9545a3704002c328e71a5ac1b601a14e74ca5 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 00:43:00 -0700 Subject: [PATCH 20/29] 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. --- engine-proto/gram-test/circles-in-triangle.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index ca49574..457ac0d 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -66,7 +66,7 @@ guess = hcat( 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 4:5] +frozen = [CartesianIndex(j, 9) for j in 1:5] #= guess = hcat( Engine.plane(BigFloat[0, 0, 1], BigFloat(0)), From 9007c8bc7c906067e03f040180a2bc48b179053a Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 00:43:44 -0700 Subject: [PATCH 21/29] Circles in triangle: jiggle the guess --- engine-proto/gram-test/circles-in-triangle.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index 457ac0d..c1f8bf2 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -1,6 +1,7 @@ include("Engine.jl") using SparseArrays +using Random # initialize the partial gram matrix for a sphere inscribed in a regular # tetrahedron @@ -55,15 +56,16 @@ gram = sparse(J, K, values) ## guess = Engine.rand_on_shell(fill(BigFloat(-1), 8)) # set initial guess +Random.seed!(58271) 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)), + 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)) + 0.1*Engine.rand_on_shell([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)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]), + 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)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]), + 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, 0, 1] ) frozen = [CartesianIndex(j, 9) for j in 1:5] From b040bbb7feff00f3399eb1a85f938d32c814608c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 00:50:48 -0700 Subject: [PATCH 22/29] Drop old code from examples --- engine-proto/gram-test/circles-in-triangle.jl | 37 +------------------ .../gram-test/overlapping-pyramids.jl | 10 +---- .../gram-test/sphere-in-tetrahedron.jl | 5 +-- .../gram-test/tetrahedron-radius-ratio.jl | 2 +- 4 files changed, 4 insertions(+), 50 deletions(-) diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index c1f8bf2..f7248df 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -47,14 +47,6 @@ 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 Random.seed!(58271) guess = hcat( @@ -69,39 +61,12 @@ guess = hcat( 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) -=# +# 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) -#= -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 diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index 0d1f018..c530296 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -50,19 +50,11 @@ guess = begin ) end -# complete the gram matrix -#= -L, history = Engine.realize_gram_gradient(gram, guess, scaled_tol = 0.01) -L_pol, history_pol = Engine.realize_gram_newton(gram, L) -=# +# complete the gram matrix using Newton's method with backtracking L, success, history = Engine.realize_gram(gram, guess) 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)) -println("Loss: ", history_pol.scaled_loss[end], "\n") -=# if success println("\nTarget accuracy achieved!") else diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 631f0e5..1c36e99 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -53,10 +53,7 @@ guess = hcat( ) frozen = [CartesianIndex(j, 6) for j in 1:5] -# complete the gram matrix -#= -L, history = Engine.realize_gram_newton(gram, guess) -=# +# 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") diff --git a/engine-proto/gram-test/tetrahedron-radius-ratio.jl b/engine-proto/gram-test/tetrahedron-radius-ratio.jl index c284078..7ceb794 100644 --- a/engine-proto/gram-test/tetrahedron-radius-ratio.jl +++ b/engine-proto/gram-test/tetrahedron-radius-ratio.jl @@ -77,7 +77,7 @@ frozen = vcat( [CartesianIndex(j, 11) 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) completed_gram = L'*Engine.Q*L println("Completed Gram matrix:\n") From b24dcc9af8c26241825093762a8e8b984ece5f23 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 01:04:40 -0700 Subject: [PATCH 23/29] Report success correctly when step limit is reached --- engine-proto/gram-test/Engine.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index 40b77b0..ac5fe54 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -445,7 +445,7 @@ function realize_gram( # return the factorization and its history push!(history.scaled_loss, loss / scale_adjustment) - L, true, history + L, loss < tol, history end end \ No newline at end of file From 33c09917d0da99ad8d973254b37272ae0b4ef47e Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 01:05:13 -0700 Subject: [PATCH 24/29] Correct scope of guess constants --- engine-proto/ConstructionViewer.jl | 4 ++-- engine-proto/gram-test/overlapping-pyramids.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/engine-proto/ConstructionViewer.jl b/engine-proto/ConstructionViewer.jl index c9b0b7a..b9c8ffb 100644 --- a/engine-proto/ConstructionViewer.jl +++ b/engine-proto/ConstructionViewer.jl @@ -207,8 +207,8 @@ end # ~~~ sandbox setup ~~~ -elements = begin - const a = sqrt(BigFloat(3)/2) +elements = let + a = sqrt(BigFloat(3)/2) sqrt(0.5) * BigFloat[ 1 1 -1 -1 0 1 -1 1 -1 0 diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index c530296..cf4b88d 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -36,8 +36,8 @@ gram[1, 6] = gram[6, 1] # in this initial guess, the mutual tangency condition is satisfied for spheres # 1 through 5 Random.seed!(50793) -guess = begin - const a = sqrt(BigFloat(3)/2) +guess = let + a = sqrt(BigFloat(3)/2) hcat( sqrt(1/BigFloat(2)) * BigFloat[ 1 1 -1 -1 0 From 71c10adbdd2f82d9979b3fae32be3da2720e1efe Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 01:12:49 -0700 Subject: [PATCH 25/29] Overlapping pyramids: drop outdated comment --- engine-proto/gram-test/overlapping-pyramids.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index cf4b88d..a4ae01a 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -23,12 +23,6 @@ end gram = sparse(J, K, values) # 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 gram[6, 1] = BigFloat(indep_val) gram[1, 6] = gram[6, 1] From 19a4d49497061b4c0523f32e0c5d7ed5af045af3 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 01:48:05 -0700 Subject: [PATCH 26/29] Clean up example formatting --- engine-proto/gram-test/circles-in-triangle.jl | 8 ++++---- engine-proto/gram-test/sphere-in-tetrahedron.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index f7248df..1bd22a7 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -12,22 +12,22 @@ for j in 1:9 for k in 1:9 filled = false if j == 9 - if (k <= 5 && k != 2) + if k <= 5 && k != 2 push!(values, 0) filled = true end elseif k == 9 - if (j <= 5 && j != 2) + if j <= 5 && j != 2 push!(values, 0) filled = true end elseif j == k push!(values, 1) filled = true - elseif (j == 1 || k == 1) + elseif j == 1 || k == 1 push!(values, 0) filled = true - elseif (j == 2 || k == 2) + elseif j == 2 || k == 2 push!(values, -1) filled = true end diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 1c36e99..97f0720 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -24,7 +24,7 @@ for j in 1:6 elseif j == k push!(values, 1) filled = true - elseif (j <= 4 && k <= 4) + elseif j <= 4 && k <= 4 push!(values, -1/BigFloat(3)) filled = true else From a26f1e3927472cdd1b2453509ba0b29a1d4ae702 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 03:16:57 -0700 Subject: [PATCH 27/29] Add Irisawa hexlet example Hat tip Romy, who sent me the article on sangaku that led me to this problem. --- engine-proto/gram-test/irisawa-hexlet.jl | 77 ++++++++++++++ engine-proto/gram-test/irisawa-hexlet_bad.jl | 105 +++++++++++++++++++ 2 files changed, 182 insertions(+) create mode 100644 engine-proto/gram-test/irisawa-hexlet.jl create mode 100644 engine-proto/gram-test/irisawa-hexlet_bad.jl diff --git a/engine-proto/gram-test/irisawa-hexlet.jl b/engine-proto/gram-test/irisawa-hexlet.jl new file mode 100644 index 0000000..8b43ad4 --- /dev/null +++ b/engine-proto/gram-test/irisawa-hexlet.jl @@ -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 two nucleus 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 sphere in the chain + 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 \ No newline at end of file diff --git a/engine-proto/gram-test/irisawa-hexlet_bad.jl b/engine-proto/gram-test/irisawa-hexlet_bad.jl new file mode 100644 index 0000000..8786778 --- /dev/null +++ b/engine-proto/gram-test/irisawa-hexlet_bad.jl @@ -0,0 +1,105 @@ +include("Engine.jl") + +using SparseArrays + +# --- construct the nucleus spheres --- + +println("--- Nucleus spheres ---\n") + +# initialize the partial gram matrix for the circumscribing and nucleus spheres +J = Int64[] +K = Int64[] +values = BigFloat[] +for n in 1:3 + push!(J, n) + push!(K, n) + push!(values, 1) + if n > 1 + append!(J, [1, n]) + append!(K, [n, 1]) + append!(values, [1, 1]) + end +end +gram_nuc = sparse(J, K, values) + +# make an initial guess +guess_nuc = hcat( + Engine.sphere(BigFloat[0, 0, 0], BigFloat(15)), + Engine.sphere(BigFloat[0, 0, -10], BigFloat(5)), + Engine.sphere(BigFloat[0, 0, 11], BigFloat(3)), +) +frozen_nuc = [CartesianIndex(4, k) for k in 1:3] + +# complete the gram matrix using Newton's method with backtracking +L_nuc, success_nuc, history_nuc = Engine.realize_gram(gram_nuc, guess_nuc, frozen_nuc) +completed_gram_nuc = L_nuc'*Engine.Q*L_nuc +println("Completed Gram matrix:\n") +display(completed_gram_nuc) +if success_nuc + println("\nTarget accuracy achieved!") +else + println("\nFailed to reach target accuracy") +end +println("Steps: ", size(history_nuc.scaled_loss, 1)) +println("Loss: ", history_nuc.scaled_loss[end], "\n") + +# --- construct the chain of spheres --- + +# initialize the partial gram matrix for the chain of spheres +J = Int64[] +K = Int64[] +values = BigFloat[] +for a in 4:9 + push!(J, a) + push!(K, a) + push!(values, 1) + + # each chain sphere is internally tangent to the circumscribing sphere + append!(J, [a, 1]) + append!(K, [1, a]) + append!(values, [1, 1]) + + # each chain sphere is externally tangent to the nucleus spheres + for n in 2:3 + append!(J, [a, n]) + append!(K, [n, a]) + append!(values, [-1, -1]) + end + + # each chain sphere is externally tangent to the next sphere in the chain + #= + a_next = 4 + mod(a-3, 6) + append!(J, [a, a_next]) + append!(K, [a_next, a]) + append!(values, [-1, -1]) + =# +end +gram_chain = sparse(J, K, values) + +if success_nuc + println("--- Chain spheres ---\n") + + # make an initial guess, with the circumscribing and nucleus spheres included + # as frozen elements + guess_chain = hcat( + L_nuc, + ( + Engine.sphere(10*BigFloat[cos(k*π/3), sin(k*π/3), 0], BigFloat(2.5)) + for k in 1:6 + )... + ) + frozen_chain = [CartesianIndex(j, k) for k in 1:3 for j in 1:5] + + # complete the gram matrix using Newton's method with backtracking + L_chain, success_chain, history_chain = Engine.realize_gram(gram_chain, guess_chain, frozen_chain) + completed_gram_chain = L_chain'*Engine.Q*L_chain + println("Completed Gram matrix:\n") + display(completed_gram_chain) + if success_chain + println("\nTarget accuracy achieved!") + else + println("\nFailed to reach target accuracy") + end + println("Steps: ", size(history_chain.scaled_loss, 1)) + println("Loss: ", history_chain.scaled_loss[end], "\n") +end \ No newline at end of file From 8a77cd74846f972cd1721fbefa130618ac03e048 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 03:21:46 -0700 Subject: [PATCH 28/29] 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. --- engine-proto/gram-test/irisawa-hexlet_bad.jl | 105 ------------------- 1 file changed, 105 deletions(-) delete mode 100644 engine-proto/gram-test/irisawa-hexlet_bad.jl diff --git a/engine-proto/gram-test/irisawa-hexlet_bad.jl b/engine-proto/gram-test/irisawa-hexlet_bad.jl deleted file mode 100644 index 8786778..0000000 --- a/engine-proto/gram-test/irisawa-hexlet_bad.jl +++ /dev/null @@ -1,105 +0,0 @@ -include("Engine.jl") - -using SparseArrays - -# --- construct the nucleus spheres --- - -println("--- Nucleus spheres ---\n") - -# initialize the partial gram matrix for the circumscribing and nucleus spheres -J = Int64[] -K = Int64[] -values = BigFloat[] -for n in 1:3 - push!(J, n) - push!(K, n) - push!(values, 1) - if n > 1 - append!(J, [1, n]) - append!(K, [n, 1]) - append!(values, [1, 1]) - end -end -gram_nuc = sparse(J, K, values) - -# make an initial guess -guess_nuc = hcat( - Engine.sphere(BigFloat[0, 0, 0], BigFloat(15)), - Engine.sphere(BigFloat[0, 0, -10], BigFloat(5)), - Engine.sphere(BigFloat[0, 0, 11], BigFloat(3)), -) -frozen_nuc = [CartesianIndex(4, k) for k in 1:3] - -# complete the gram matrix using Newton's method with backtracking -L_nuc, success_nuc, history_nuc = Engine.realize_gram(gram_nuc, guess_nuc, frozen_nuc) -completed_gram_nuc = L_nuc'*Engine.Q*L_nuc -println("Completed Gram matrix:\n") -display(completed_gram_nuc) -if success_nuc - println("\nTarget accuracy achieved!") -else - println("\nFailed to reach target accuracy") -end -println("Steps: ", size(history_nuc.scaled_loss, 1)) -println("Loss: ", history_nuc.scaled_loss[end], "\n") - -# --- construct the chain of spheres --- - -# initialize the partial gram matrix for the chain of spheres -J = Int64[] -K = Int64[] -values = BigFloat[] -for a in 4:9 - push!(J, a) - push!(K, a) - push!(values, 1) - - # each chain sphere is internally tangent to the circumscribing sphere - append!(J, [a, 1]) - append!(K, [1, a]) - append!(values, [1, 1]) - - # each chain sphere is externally tangent to the nucleus spheres - for n in 2:3 - append!(J, [a, n]) - append!(K, [n, a]) - append!(values, [-1, -1]) - end - - # each chain sphere is externally tangent to the next sphere in the chain - #= - a_next = 4 + mod(a-3, 6) - append!(J, [a, a_next]) - append!(K, [a_next, a]) - append!(values, [-1, -1]) - =# -end -gram_chain = sparse(J, K, values) - -if success_nuc - println("--- Chain spheres ---\n") - - # make an initial guess, with the circumscribing and nucleus spheres included - # as frozen elements - guess_chain = hcat( - L_nuc, - ( - Engine.sphere(10*BigFloat[cos(k*π/3), sin(k*π/3), 0], BigFloat(2.5)) - for k in 1:6 - )... - ) - frozen_chain = [CartesianIndex(j, k) for k in 1:3 for j in 1:5] - - # complete the gram matrix using Newton's method with backtracking - L_chain, success_chain, history_chain = Engine.realize_gram(gram_chain, guess_chain, frozen_chain) - completed_gram_chain = L_chain'*Engine.Q*L_chain - println("Completed Gram matrix:\n") - display(completed_gram_chain) - if success_chain - println("\nTarget accuracy achieved!") - else - println("\nFailed to reach target accuracy") - end - println("Steps: ", size(history_chain.scaled_loss, 1)) - println("Loss: ", history_chain.scaled_loss[end], "\n") -end \ No newline at end of file From 9d69a900e28e0694b984ed326d3c9b6d6edca668 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 18 Jul 2024 03:39:41 -0700 Subject: [PATCH 29/29] Irisawa hexlet: use Abe's terminology in comments Abe uses the names "sun" and "moon" for what Wikipedia calls the nucleus spheres. --- engine-proto/gram-test/irisawa-hexlet.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/engine-proto/gram-test/irisawa-hexlet.jl b/engine-proto/gram-test/irisawa-hexlet.jl index 8b43ad4..67def8c 100644 --- a/engine-proto/gram-test/irisawa-hexlet.jl +++ b/engine-proto/gram-test/irisawa-hexlet.jl @@ -28,14 +28,14 @@ for s in 1:9 end if s > 3 - # each chain sphere is externally tangent to the two nucleus spheres + # 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 sphere in the chain + # 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])