Move random cut routine into engine

engine-proto/Engine.Numerical.jl

@@ -1,7 +1,8 @@
 module Numerical
 
+using LinearAlgebra
 using AbstractAlgebra
-using HomotopyContinuation: Variable, Expression, System
+using HomotopyContinuation
 using ..Algebraic
 
 # --- polynomial conversion ---
@@ -10,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
 
@@ -22,4 +23,18 @@ function System(I::Generic.Ideal)
   System(eqns, variables = variables)
 end
 
+# --- sampling ---
+
+function real_samples(F::AbstractSystem, dim)
+  # choose a random real hyperplane of codimension `dim` by intersecting
+  # hyperplanes whose normal vectors are uniformly distributed over the unit
+  # sphere
+  # [to do] guard against the unlikely event that one of the normals is zero
+  normals = transpose(hcat(
+    (normalize(randn(nvariables(F))) for _ in 1:dim)...
+  ))
+  cut = LinearSubspace(normals, fill(0., dim))
+  filter(isreal, results(witness_set(F, cut)))
+end
+
 end

engine-proto/Engine.jl

@@ -82,25 +82,7 @@ n_planes = 36
 for through_trivial in [false, true]
   samples = []
   for _ in 1:n_planes
-    cut_matrix = transpose(hcat(
-      (normalize(randn(length(gens(coordring)))) for _ in 1:freedom)...
-    ))
-    ##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
-    ##]
-    ## display(cut_matrix) ## [verbose]
-    if through_trivial
-      cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
-      ## display(cut_offset) ## [verbose]
-      cut_subspace = LinearSubspace(cut_matrix, cut_offset)
-    else
-      cut_subspace = LinearSubspace(cut_matrix, fill(0., freedom))
-    end
-    wtns = witness_set(system, cut_subspace)
-    real_solns = solution.(filter(isreal, results(wtns)))
+    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