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