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

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