engine-proto/Engine.jl

@@ -225,8 +225,6 @@ end
 
 # ~~~ sandbox setup ~~~
 
-using Random
-using Distributions
 using AbstractAlgebra
 using HomotopyContinuation
 
@@ -245,8 +243,7 @@ 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")
+println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freedom")
 
 ##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
 ##tangencies = [
@@ -263,54 +260,33 @@ println("Two points on a sphere: ", freedom, " degrees of freedom")
 
 # --- test rational cut ---
 
-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))
-max_slope = 2
-binom = Binomial(2max_slope, 1/2)
-Random.seed!(6071)
-samples = []
-for _ in 1:3
-  cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope
-  ##cut_matrix = [
-  ##  1 1 1 1 0 1 1 0 1 1 0;
-  ##  1 2 1 2 0 1 1 0 1 1 0;
-  ##  1 1 0 1 0 1 2 0 2 0 0
-  ##]
-  sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
-  cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
-  println("sphere z variables: ", vbls[sph_z_ind])
-  display(cut_matrix)
-  display(cut_offset)
-  cut_subspace = LinearSubspace(cut_matrix, cut_offset)
-  wtns = witness_set(system, cut_subspace)
-  append!(samples, solution.(filter(isreal, results(wtns))))
-end
-println("$(length(samples)) sample solutions:")
-for soln in samples
-  display([vbls round.(soln, digits = 6)])
-  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))")
+cut = [
+  sum(vcat([1, 1, 1] .* a.coords, [1, 1, 1, 1, 1] .* (s.coords - [0, 0, 0, 0, 1])))
+  sum(vcat([2, 1, 1] .* a.coords, [1, 2, 1, 1, 1] .* (s.coords - [0, 0, 0, 0, 1])))
+  sum(vcat([1, 2, 0] .* a.coords, [1, 1, 0, 1, 2] .* (s.coords - [0, 0, 0, 0, 1])))
+]
+cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut])
+cut_dim = Engine.dimension(cut_ideal_ab_s)
+println("Two points on a sphere, after cut: ", cut_dim, " degrees of freedom")
+if cut_dim == 0
+  vbls = Variable.(symbols(base_ring(ideal_ab_s)))
+  cut_system = 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 corresponding witness set
+  cut_matrix = [1 1 1 1 0 1 1 0 1 1 0; 1 2 1 2 0 1 1 0 1 1 0; 1 1 0 1 0 1 2 0 2 0 0]
+  cut_subspace = LinearSubspace(cut_matrix, [1, 1, 2])
+  witness = witness_set(System(eqns_ab_s, variables = vbls), cut_subspace)
+  println("witness solutions:")
+  for wtns in solutions(witness)
+    display(wtns)
   end
 end