engine-proto/Engine.Algebraic.jl

@@ -184,8 +184,6 @@ 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
@@ -197,8 +195,7 @@ 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
   
-  ## [test] (Generic.Ideal(coordring, eqns), eqns)
+  (Generic.Ideal(coordring, eqns), eqns)
-  (nothing, eqns)
 end
 
 end

engine-proto/Engine.Numerical.jl

@@ -2,10 +2,7 @@ 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 ---
@@ -14,7 +11,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(coefficients(f), exponent_vectors(f))
+  f_data = zip(AbstractAlgebra.coefficients(f), exponent_vectors(f))
   sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data)
 end
 
@@ -40,13 +37,4 @@ 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,25 +19,24 @@ using Distributions
 using LinearAlgebra
 using AbstractAlgebra
 using HomotopyContinuation
-using GLMakie
 
 CoeffType = Rational{Int64}
 
-##a = Engine.Point{CoeffType}()
+a = Engine.Point{CoeffType}()
-##s = Engine.Sphere{CoeffType}()
+s = Engine.Sphere{CoeffType}()
-##a_on_s = Engine.LiesOn{CoeffType}(a, s)
+a_on_s = Engine.LiesOn{CoeffType}(a, s)
-##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_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 = Engine.Point{CoeffType}()
-##b_on_s = Engine.LiesOn{CoeffType}(b, s)
+b_on_s = Engine.LiesOn{CoeffType}(b, s)
-##Engine.push!(ctx, b)
+Engine.push!(ctx, b)
-##Engine.push!(ctx, s)
+Engine.push!(ctx, s)
-##Engine.push!(ctx, b_on_s)
+Engine.push!(ctx, b_on_s)
-##ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
+ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
-##freedom = Engine.dimension(ideal_ab_s)
+freedom = Engine.dimension(ideal_ab_s)
-##println("Two points on a sphere: $freedom degrees of freedom")
+println("Two points on a sphere: ", freedom, " degrees of freedom")
 
 ##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
 ##tangencies = [
@@ -48,81 +47,70 @@ CoeffType = Rational{Int64}
 ##  )
 ##  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 = Set(spheres), relations = Set(tangencies))
-##ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
+##ideal_tan_sph = Engine.realize(ctx_tan_sph)
-##freedom = Engine.dimension(ideal_tan_sph)
+##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom")
-##println("Three mutually tangent 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(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_joined, variables = vbls))
+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")
 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)
 Random.seed!(6071)
-n_planes = 3
+n_planes = 36
-samples = []
+for through_trivial in [false, true]
-for _ in 1:n_planes
+  samples = []
-  real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
+  for _ in 1:n_planes
-  for soln in real_solns
+    real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
-    if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
+    nontrivial_solns = filter(is_nontrivial, real_solns)
-      push!(samples, soln)
+    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
     end
   end
-end
+  if through_trivial
-println("$(length(samples)) sample solutions:")
+    println("--- planes through trivial solution ---")
-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("--- planes through origin ---")
-    println("  center at origin:   r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
   end
-end
+  println("$(length(samples)) sample solutions, not including the trivial one:")
-
+  for soln in samples
-# show a sample solution
+    ## display([vbls round.(soln, digits = 6)]) ## [verbose]
-function show_solution(ctx, vals)
+    k_sq = abs2(soln[1])
-  # evaluate elements
+    if abs2(soln[end-2]) > 1e-12
-  real_vals = real.(vals)
+      if k_sq < 1e-12
-  disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
+        println("  center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
-  disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
+      else
-  
+        sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
-  # create scene
+        println("  center on z axis:   r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
-  scene = Scene()
+      end
-  cam3d!(scene)
+    else
-  scatter!(scene, disp_points, color = :green)
+      sum_sq = sum(soln[[4, 7, 10]] .^ 2)
-  for sph in disp_spheres
+      println("  center at origin:   r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
-    mesh!(scene, sph, color = :gray)
+    end
   end
-  scene
 end