From 6cf07dc6a10766026cf30365c3cd431d87f17967 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 12 Feb 2024 20:34:12 -0500 Subject: [PATCH 1/4] Evaluate and display elements --- engine-proto/Engine.Numerical.jl | 16 +++++- engine-proto/Engine.jl | 90 ++++++++++++++++---------------- 2 files changed, 58 insertions(+), 48 deletions(-) diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl index 1ae7b97..48fb682 100644 --- a/engine-proto/Engine.Numerical.jl +++ b/engine-proto/Engine.Numerical.jl @@ -2,7 +2,10 @@ module Numerical using LinearAlgebra using AbstractAlgebra -using HomotopyContinuation +using HomotopyContinuation: + Variable, Expression, AbstractSystem, System, LinearSubspace, + nvariables, isreal, witness_set, results +import GLMakie using ..Algebraic # --- polynomial conversion --- @@ -11,7 +14,7 @@ using ..Algebraic # https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/issues/520#issuecomment-1317681521 function Base.convert(::Type{Expression}, f::MPolyRingElem) variables = Variable.(symbols(parent(f))) - f_data = zip(AbstractAlgebra.coefficients(f), exponent_vectors(f)) + f_data = zip(coefficients(f), exponent_vectors(f)) sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data) end @@ -37,4 +40,13 @@ function real_samples(F::AbstractSystem, dim) filter(isreal, results(witness_set(F, cut))) end +AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) = + GLMakie.Point3f([evaluate(u, vals) for u in pt.coords]) + +function AbstractAlgebra.evaluate(sph::Sphere, vals::Vector{<:RingElement}) + radius = 1 / evaluate(sph.coords[1], vals) + center = radius * [evaluate(u, vals) for u in sph.coords[3:end]] + GLMakie.Sphere(GLMakie.Point3f(center), radius) +end + end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index f058085..099a5f0 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -19,6 +19,7 @@ using Distributions using LinearAlgebra using AbstractAlgebra using HomotopyContinuation +using GLMakie CoeffType = Rational{Int64} @@ -55,62 +56,59 @@ println("Two points on a sphere: ", freedom, " degrees of freedom") coordring = base_ring(ideal_ab_s) vbls = Variable.(symbols(coordring)) -##cut_system = CompiledSystem(System([eqns_ab_s; cut], variables = vbls)) -##cut_result = HomotopyContinuation.solve(cut_system) -##println("non-singular solutions:") -##for soln in solutions(cut_result) -## display(soln) -##end -##println("singular solutions:") -##for sing in singular(cut_result) -## display(sing.solution) -##end # test a random witness set system = CompiledSystem(System(eqns_ab_s, variables = vbls)) sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring)) println("sphere z variables: ", vbls[sph_z_ind]) -trivial_soln = fill(0, length(gens(coordring))) -trivial_soln[sph_z_ind] .= 1 -println("trivial solutions: $trivial_soln") +## [old] trivial_soln = fill(0, length(gens(coordring))) +## [old] trivial_soln[sph_z_ind] .= 1 +## [old] println("trivial solutions: $trivial_soln") norm2 = vec -> real(dot(conj.(vec), vec)) -is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring)) -##max_slope = 5 -##binom = Binomial(2max_slope, 1/2) +## [old] is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring)) Random.seed!(6071) -n_planes = 36 -for through_trivial in [false, true] - samples = [] - for _ in 1:n_planes - real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) - nontrivial_solns = filter(is_nontrivial, real_solns) - println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") - for soln in nontrivial_solns - ## [test] for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns)))) - if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) - push!(samples, soln) - end +n_planes = 3 +samples = [] +for _ in 1:n_planes + real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) + ## [old] nontrivial_solns = filter(is_nontrivial, real_solns) + ## [old] println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") + for soln in real_solns + if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) + push!(samples, soln) end end - if through_trivial - println("--- planes through trivial solution ---") - else - println("--- planes through origin ---") - end - println("$(length(samples)) sample solutions, not including the trivial one:") - 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 +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 = 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))") + 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(vals) + # evaluate elements + real_vals = real.(vals) + disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points] + disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres] + + # create scene + scene = Scene() + cam3d!(scene) + scatter!(scene, disp_points, color = :green) + for sph in disp_spheres + mesh!(scene, sph, color = :gray) + end + scene end \ No newline at end of file From a450f701fbfa1818c7795d654e6f2a65d2b75b15 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 12 Feb 2024 21:14:07 -0500 Subject: [PATCH 2/4] Try displaying a chain of spheres For three mutually tangent spheres, I couldn't find real solutions. --- engine-proto/Engine.jl | 56 ++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 29 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 099a5f0..113c6a3 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -23,23 +23,23 @@ 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])) +##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") +##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") +##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] +spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] ##tangencies = [ ## Engine.AlignsWithBy{CoeffType}( ## spheres[n], @@ -48,31 +48,29 @@ println("Two points on a sphere: ", freedom, " degrees of freedom") ## ) ## for n in 1:3 ##] -##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) -##ideal_tan_sph = Engine.realize(ctx_tan_sph) -##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom") +tangencies = [ + Engine.AlignsWithBy{CoeffType}(spheres[1], spheres[2], CoeffType(-1)), + Engine.AlignsWithBy{CoeffType}(spheres[2], spheres[3], CoeffType(-1//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") +println("Chain of three spheres: $freedom degrees of freedom") # --- test rational cut --- -coordring = base_ring(ideal_ab_s) +coordring = base_ring(ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_ab_s, variables = vbls)) -sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring)) -println("sphere z variables: ", vbls[sph_z_ind]) -## [old] trivial_soln = fill(0, length(gens(coordring))) -## [old] trivial_soln[sph_z_ind] .= 1 -## [old] println("trivial solutions: $trivial_soln") +system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) -## [old] is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring)) Random.seed!(6071) -n_planes = 3 +n_planes = 16 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) - ## [old] nontrivial_solns = filter(is_nontrivial, real_solns) - ## [old] println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") for soln in real_solns if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) push!(samples, soln) @@ -97,7 +95,7 @@ for soln in samples end # show a sample solution -function show_solution(vals) +function show_solution(ctx, vals) # evaluate elements real_vals = real.(vals) disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points] From 31d5e7e864704895dbf0273c46185476e53e40f1 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 12 Feb 2024 22:48:16 -0500 Subject: [PATCH 3/4] Play with two points on two spheres Guess conditions that make the scaling constraint impossible to satisfy. --- engine-proto/Engine.Algebraic.jl | 2 ++ engine-proto/Engine.jl | 43 +++++++++++++++++++++++--------- 2 files changed, 33 insertions(+), 12 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index b9b790a..ca39967 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -184,6 +184,8 @@ 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 diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 113c6a3..7ddf72d 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -39,7 +39,7 @@ 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] +##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] ##tangencies = [ ## Engine.AlignsWithBy{CoeffType}( ## spheres[n], @@ -48,26 +48,45 @@ spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] ## ) ## for n in 1:3 ##] -tangencies = [ - Engine.AlignsWithBy{CoeffType}(spheres[1], spheres[2], CoeffType(-1)), - Engine.AlignsWithBy{CoeffType}(spheres[2], spheres[3], CoeffType(-1//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) +##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 = 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") -println("Chain of three spheres: $freedom degrees of freedom") + +p = Engine.Point{CoeffType}() +q = Engine.Point{CoeffType}() +a = Engine.Sphere{CoeffType}() +b = Engine.Sphere{CoeffType}() +p_on_a = Engine.LiesOn{CoeffType}(p, a) +p_on_b = Engine.LiesOn{CoeffType}(p, b) +q_on_a = Engine.LiesOn{CoeffType}(q, a) +q_on_b = Engine.LiesOn{CoeffType}(q, b) +ctx_joined = Engine.Construction{CoeffType}( + elements = Set([p, q, a, b]), + relations= Set([p_on_a, p_on_b, q_on_a, q_on_b]) +) +ideal_joined, eqns_joined = Engine.realize(ctx_joined) +freedom = Engine.dimension(ideal_joined) +println("Two points on two spheres: $freedom degrees of freedom") # --- test rational cut --- -coordring = base_ring(ideal_tan_sph) +coordring = base_ring(ideal_joined) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) +system = CompiledSystem(System(eqns_joined, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) Random.seed!(6071) -n_planes = 16 +n_planes = 3 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) From e41bcc7e13d29217d51ed5e7810b2c8f4e904fed Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 13 Feb 2024 04:02:14 -0500 Subject: [PATCH 4/4] Explore the performance wall Three points on two spheres is too much. --- engine-proto/Engine.Algebraic.jl | 3 ++- engine-proto/Engine.jl | 19 ++++++++----------- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index ca39967..380cee1 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -197,7 +197,8 @@ function realize(ctx::Construction{T}) where T push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) end - (Generic.Ideal(coordring, eqns), eqns) + ## [test] (Generic.Ideal(coordring, eqns), eqns) + (nothing, eqns) end end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 7ddf72d..49011c6 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -61,21 +61,18 @@ CoeffType = Rational{Int64} ##freedom = Engine.dimension(ideal_tan_sph) ##println("Three mutually tangent spheres: $freedom degrees of freedom") -p = Engine.Point{CoeffType}() -q = Engine.Point{CoeffType}() -a = Engine.Sphere{CoeffType}() -b = Engine.Sphere{CoeffType}() -p_on_a = Engine.LiesOn{CoeffType}(p, a) -p_on_b = Engine.LiesOn{CoeffType}(p, b) -q_on_a = Engine.LiesOn{CoeffType}(q, a) -q_on_b = Engine.LiesOn{CoeffType}(q, b) +points = [Engine.Point{CoeffType}() for _ in 1:3] +spheres = [Engine.Sphere{CoeffType}() for _ in 1:2] ctx_joined = Engine.Construction{CoeffType}( - elements = Set([p, q, a, b]), - relations= Set([p_on_a, p_on_b, q_on_a, q_on_b]) + 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("Two points on two spheres: $freedom degrees of freedom") +println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom") # --- test rational cut ---