From 8e33987f596ae9572b53885f5b4642105d45165c Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Sat, 10 Feb 2024 13:46:01 -0500
Subject: [PATCH] Systematically try out different cut planes

---
 engine-proto/Engine.jl | 80 +++++++++++++++++++++++++++---------------
 1 file changed, 51 insertions(+), 29 deletions(-)

diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index 9841223..ba4300c 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -227,6 +227,7 @@ end
 
 using Random
 using Distributions
+using LinearAlgebra
 using AbstractAlgebra
 using HomotopyContinuation
 
@@ -278,39 +279,60 @@ vbls = Variable.(symbols(coordring))
 
 # test a random witness set
 system = CompiledSystem(System(eqns_ab_s, variables = vbls))
-max_slope = 2
+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)
-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))}")
+n_planes = 36
+for through_trivial in [false, true]
+  samples = []
+  for _ in 1:n_planes
+    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
+    ##]
+    ## [verbose] display(cut_matrix)
+    if through_trivial
+      cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
+      ## [verbose] display(cut_offset)
+      cut_subspace = LinearSubspace(cut_matrix, cut_offset)
     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))")
+      cut_subspace = LinearSubspace(cut_matrix, fill(0, 3))
     end
+    wtns = witness_set(system, cut_subspace)
+    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
+  if through_trivial
+    println("--- planes through trivial solution ---")
   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))")
+    println("--- planes through origin ---")
+  end
+  println("$(length(samples)) sample solutions, not including the trivial one:")
+  for soln in samples
+    ## [verbose] 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))")
+    end
   end
 end
\ No newline at end of file