From 291d5c8ff685b642fd3535b6c46179c9f134c469 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 13:28:01 -0800 Subject: [PATCH 1/6] Study mutually tangent spheres with two fixed --- engine-proto/Engine.Algebraic.jl | 24 ++++----- engine-proto/Engine.jl | 93 +++++++++++++++++--------------- 2 files changed, 62 insertions(+), 55 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index 380cee1..b6fc7a7 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -186,19 +186,19 @@ function realize(ctx::Construction{T}) where T # add relations to center, orient, and scale the construction # [to do] the scaling constraint, as written, can be impossible to satisfy # when all of the spheres have to go through the origin - if !isempty(ctx.points) - append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) - end - if !isempty(ctx.spheres) - append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) - end - n_elts = length(ctx.points) + length(ctx.spheres) - if n_elts > 0 - push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) - end + ##if !isempty(ctx.points) + ## append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) + ##end + ##if !isempty(ctx.spheres) + ## append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) + ##end + ##n_elts = length(ctx.points) + length(ctx.spheres) + ##if n_elts > 0 + ## push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) + ##end - ## [test] (Generic.Ideal(coordring, eqns), eqns) - (nothing, eqns) + (Generic.Ideal(coordring, eqns), eqns) + ## [test] (nothing, eqns) end end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 49011c6..0410f6d 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -39,15 +39,15 @@ CoeffType = Rational{Int64} ##freedom = Engine.dimension(ideal_ab_s) ##println("Two points on a sphere: $freedom degrees of freedom") -##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] -##tangencies = [ -## Engine.AlignsWithBy{CoeffType}( -## spheres[n], -## spheres[mod1(n+1, length(spheres))], -## CoeffType(-1//1) -## ) -## for n in 1:3 -##] +spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] +tangencies = [ + Engine.AlignsWithBy{CoeffType}( + spheres[n], + spheres[mod1(n+1, length(spheres))], + CoeffType(-1)^n + ) + for n in 1:3 +] ##tangencies = [ ##Engine.LiesOn{CoeffType}(points[1], spheres[2]), ##Engine.LiesOn{CoeffType}(points[1], spheres[3]), @@ -56,34 +56,41 @@ CoeffType = Rational{Int64} ##Engine.LiesOn{CoeffType}(points[3], spheres[1]), ##Engine.LiesOn{CoeffType}(points[3], spheres[2]) ##] -##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) -##ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph) -##freedom = Engine.dimension(ideal_tan_sph) -##println("Three mutually tangent spheres: $freedom degrees of freedom") +ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) +ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph) +##small_eqns_tan_sph = eqns_tan_sph +small_eqns_tan_sph = [ + eqns_tan_sph; + spheres[2].coords - [1, 0, 0, 0, 1]; + spheres[3].coords - [1, 0, 0, 0, -1]; +] +small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph) +freedom = Engine.dimension(small_ideal_tan_sph) +println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom") -points = [Engine.Point{CoeffType}() for _ in 1:3] -spheres = [Engine.Sphere{CoeffType}() for _ in 1:2] -ctx_joined = Engine.Construction{CoeffType}( - elements = Set([points; spheres]), - relations= Set([ - Engine.LiesOn{CoeffType}(pt, sph) - for pt in points for sph in spheres - ]) -) -ideal_joined, eqns_joined = Engine.realize(ctx_joined) -freedom = Engine.dimension(ideal_joined) -println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom") +##points = [Engine.Point{CoeffType}() for _ in 1:3] +##spheres = [Engine.Sphere{CoeffType}() for _ in 1:2] +##ctx_joined = Engine.Construction{CoeffType}( +## elements = Set([points; spheres]), +## relations= Set([ +## Engine.LiesOn{CoeffType}(pt, sph) +## for pt in points for sph in spheres +## ]) +##) +##ideal_joined, eqns_joined = Engine.realize(ctx_joined) +##freedom = Engine.dimension(ideal_joined) +##println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom") # --- test rational cut --- -coordring = base_ring(ideal_joined) +coordring = base_ring(small_ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_joined, variables = vbls)) +system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) Random.seed!(6071) -n_planes = 3 +n_planes = 36 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) @@ -94,21 +101,21 @@ for _ in 1:n_planes end end println("$(length(samples)) sample solutions:") -for soln in samples - ## display([vbls round.(soln, digits = 6)]) ## [verbose] - k_sq = abs2(soln[1]) - if abs2(soln[end-2]) > 1e-12 - if k_sq < 1e-12 - println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))") - else - sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq - println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))") - end - else - sum_sq = sum(soln[[4, 7, 10]] .^ 2) - println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))") - end -end +##for soln in samples +## ## display([vbls round.(soln, digits = 6)]) ## [verbose] +## k_sq = abs2(soln[1]) +## if abs2(soln[end-2]) > 1e-12 +## if k_sq < 1e-12 +## println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))") +## else +## sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq +## println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))") +## end +## else +## sum_sq = sum(soln[[4, 7, 10]] .^ 2) +## println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))") +## end +##end # show a sample solution function show_solution(ctx, vals) From 8d8bc9162cd494798821b1d83c2f338a788e736d Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 14:27:41 -0800 Subject: [PATCH 2/6] Store elements in arrays to keep order stable This seems to restore reproducibility. --- engine-proto/Engine.Algebraic.jl | 8 ++++---- engine-proto/Engine.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index b6fc7a7..f14f58c 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -120,11 +120,11 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - # --- constructions --- mutable struct Construction{T} - points::Set{Point{T}} - spheres::Set{Sphere{T}} - relations::Set{Relation{T}} + points::Vector{Point{T}} + spheres::Vector{Sphere{T}} + relations::Vector{Relation{T}} - function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T + function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{Relation{T}}()) where T allelements = union(elements, (rel.elements for rel in relations)...) new{T}( filter(elt -> isa(elt, Point), allelements), diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 0410f6d..d283e72 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -56,7 +56,7 @@ tangencies = [ ##Engine.LiesOn{CoeffType}(points[3], spheres[1]), ##Engine.LiesOn{CoeffType}(points[3], spheres[2]) ##] -ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) +ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies) ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph) ##small_eqns_tan_sph = eqns_tan_sph small_eqns_tan_sph = [ From ae5db0f9eaa71a4f87b70e5d95efd43fa0fd9426 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 16:00:46 -0800 Subject: [PATCH 3/6] Make results reproducible --- engine-proto/Engine.Numerical.jl | 7 ++++--- engine-proto/Engine.jl | 6 +++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl index 48fb682..d1e14bd 100644 --- a/engine-proto/Engine.Numerical.jl +++ b/engine-proto/Engine.Numerical.jl @@ -1,5 +1,6 @@ module Numerical +using Random: default_rng using LinearAlgebra using AbstractAlgebra using HomotopyContinuation: @@ -28,16 +29,16 @@ end # --- sampling --- -function real_samples(F::AbstractSystem, dim) +function real_samples(F::AbstractSystem, dim; rng = default_rng()) # choose a random real hyperplane of codimension `dim` by intersecting # hyperplanes whose normal vectors are uniformly distributed over the unit # sphere # [to do] guard against the unlikely event that one of the normals is zero normals = transpose(hcat( - (normalize(randn(nvariables(F))) for _ in 1:dim)... + (normalize(randn(rng, nvariables(F))) for _ in 1:dim)... )) cut = LinearSubspace(normals, fill(0., dim)) - filter(isreal, results(witness_set(F, cut))) + filter(isreal, results(witness_set(F, cut, seed = 0x1974abba))) end AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) = diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index d283e72..7146e80 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -89,11 +89,11 @@ vbls = Variable.(symbols(coordring)) # test a random witness set system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) -Random.seed!(6071) -n_planes = 36 +rng = MersenneTwister(6071) +n_planes = 3 samples = [] for _ in 1:n_planes - real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) + real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng)) for soln in real_solns if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) push!(samples, soln) From ba365174d3a643c033a04178706fe04bfe2bc555 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 16:16:06 -0800 Subject: [PATCH 4/6] Find real solutions for three mutually tangent spheres I'm not sure why the solver wasn't working before. It might've been just an unlucky random number draw. --- engine-proto/Engine.Algebraic.jl | 20 ++++++++++---------- engine-proto/Engine.jl | 18 +++++++++--------- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index f14f58c..2d28017 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -186,16 +186,16 @@ function realize(ctx::Construction{T}) where T # add relations to center, orient, and scale the construction # [to do] the scaling constraint, as written, can be impossible to satisfy # when all of the spheres have to go through the origin - ##if !isempty(ctx.points) - ## append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) - ##end - ##if !isempty(ctx.spheres) - ## append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) - ##end - ##n_elts = length(ctx.points) + length(ctx.spheres) - ##if n_elts > 0 - ## push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) - ##end + if !isempty(ctx.points) + append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) + end + if !isempty(ctx.spheres) + append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) + end + n_elts = length(ctx.points) + length(ctx.spheres) + if n_elts > 0 + push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) + end (Generic.Ideal(coordring, eqns), eqns) ## [test] (nothing, eqns) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 7146e80..e92eed4 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -59,13 +59,13 @@ tangencies = [ ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies) ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph) ##small_eqns_tan_sph = eqns_tan_sph -small_eqns_tan_sph = [ - eqns_tan_sph; - spheres[2].coords - [1, 0, 0, 0, 1]; - spheres[3].coords - [1, 0, 0, 0, -1]; -] -small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph) -freedom = Engine.dimension(small_ideal_tan_sph) +##small_eqns_tan_sph = [ +## eqns_tan_sph; +## spheres[2].coords - [1, 0, 0, 0, 1]; +## spheres[3].coords - [1, 0, 0, 0, -1]; +##] +##small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph) +freedom = Engine.dimension(ideal_tan_sph) println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom") ##points = [Engine.Point{CoeffType}() for _ in 1:3] @@ -83,11 +83,11 @@ println("Three mutually tangent spheres, with two fixed: $freedom degrees of fre # --- test rational cut --- -coordring = base_ring(small_ideal_tan_sph) +coordring = base_ring(ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls)) +system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) rng = MersenneTwister(6071) n_planes = 3 From f2000e5731f52e2583d15619738acbd86c7c7974 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 16:25:09 -0800 Subject: [PATCH 5/6] Test different sign patterns for cosines It seems like there are real solutions if and only if the product of the cosines is positive. --- engine-proto/Engine.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index e92eed4..a517117 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -44,7 +44,7 @@ tangencies = [ Engine.AlignsWithBy{CoeffType}( spheres[n], spheres[mod1(n+1, length(spheres))], - CoeffType(-1)^n + CoeffType([1, 1, 1][n]) ) for n in 1:3 ] @@ -90,7 +90,7 @@ vbls = Variable.(symbols(coordring)) system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) rng = MersenneTwister(6071) -n_planes = 3 +n_planes = 36 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng)) From 3170a933e4594befb3395882dc453b3f96143ae7 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 17:16:37 -0800 Subject: [PATCH 6/6] Clean up example of three mutually tangent spheres --- engine-proto/Engine.Algebraic.jl | 1 - engine-proto/Engine.jl | 67 ++------------------------------ 2 files changed, 4 insertions(+), 64 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index 2d28017..a9b6667 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -198,7 +198,6 @@ function realize(ctx::Construction{T}) where T end (Generic.Ideal(coordring, eqns), eqns) - ## [test] (nothing, eqns) end end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index a517117..f6f92c5 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -23,63 +23,19 @@ using GLMakie CoeffType = Rational{Int64} -##a = Engine.Point{CoeffType}() -##s = Engine.Sphere{CoeffType}() -##a_on_s = Engine.LiesOn{CoeffType}(a, s) -##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s])) -##ideal_a_s = Engine.realize(ctx) -##println("A point on a sphere: $(Engine.dimension(ideal_a_s)) degrees of freedom") - -##b = Engine.Point{CoeffType}() -##b_on_s = Engine.LiesOn{CoeffType}(b, s) -##Engine.push!(ctx, b) -##Engine.push!(ctx, s) -##Engine.push!(ctx, b_on_s) -##ideal_ab_s, eqns_ab_s = Engine.realize(ctx) -##freedom = Engine.dimension(ideal_ab_s) -##println("Two points on a sphere: $freedom degrees of freedom") - spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] tangencies = [ Engine.AlignsWithBy{CoeffType}( spheres[n], spheres[mod1(n+1, length(spheres))], - CoeffType([1, 1, 1][n]) + CoeffType(1) ) for n in 1:3 ] -##tangencies = [ - ##Engine.LiesOn{CoeffType}(points[1], spheres[2]), - ##Engine.LiesOn{CoeffType}(points[1], spheres[3]), - ##Engine.LiesOn{CoeffType}(points[2], spheres[3]), - ##Engine.LiesOn{CoeffType}(points[2], spheres[1]), - ##Engine.LiesOn{CoeffType}(points[3], spheres[1]), - ##Engine.LiesOn{CoeffType}(points[3], spheres[2]) -##] ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies) ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph) -##small_eqns_tan_sph = eqns_tan_sph -##small_eqns_tan_sph = [ -## eqns_tan_sph; -## spheres[2].coords - [1, 0, 0, 0, 1]; -## spheres[3].coords - [1, 0, 0, 0, -1]; -##] -##small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph) freedom = Engine.dimension(ideal_tan_sph) -println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom") - -##points = [Engine.Point{CoeffType}() for _ in 1:3] -##spheres = [Engine.Sphere{CoeffType}() for _ in 1:2] -##ctx_joined = Engine.Construction{CoeffType}( -## elements = Set([points; spheres]), -## relations= Set([ -## Engine.LiesOn{CoeffType}(pt, sph) -## for pt in points for sph in spheres -## ]) -##) -##ideal_joined, eqns_joined = Engine.realize(ctx_joined) -##freedom = Engine.dimension(ideal_joined) -##println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom") +println("Three mutually tangent spheres: $freedom degrees of freedom") # --- test rational cut --- @@ -90,7 +46,7 @@ vbls = Variable.(symbols(coordring)) system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) rng = MersenneTwister(6071) -n_planes = 36 +n_planes = 6 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng)) @@ -100,22 +56,7 @@ for _ in 1:n_planes end end end -println("$(length(samples)) sample solutions:") -##for soln in samples -## ## display([vbls round.(soln, digits = 6)]) ## [verbose] -## k_sq = abs2(soln[1]) -## if abs2(soln[end-2]) > 1e-12 -## if k_sq < 1e-12 -## println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))") -## else -## sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq -## println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))") -## end -## else -## sum_sq = sum(soln[[4, 7, 10]] .^ 2) -## println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))") -## end -##end +println("Found $(length(samples)) sample solutions") # show a sample solution function show_solution(ctx, vals)