From a450f701fbfa1818c7795d654e6f2a65d2b75b15 Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Mon, 12 Feb 2024 21:14:07 -0500
Subject: [PATCH] Try displaying a chain of spheres

For three mutually tangent spheres, I couldn't find real solutions.
---
 engine-proto/Engine.jl | 56 ++++++++++++++++++++----------------------
 1 file changed, 27 insertions(+), 29 deletions(-)

diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index 099a5f0..113c6a3 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -23,23 +23,23 @@ using GLMakie
 
 CoeffType = Rational{Int64}
 
-a = Engine.Point{CoeffType}()
-s = Engine.Sphere{CoeffType}()
-a_on_s = Engine.LiesOn{CoeffType}(a, s)
-ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
+##a = Engine.Point{CoeffType}()
+##s = Engine.Sphere{CoeffType}()
+##a_on_s = Engine.LiesOn{CoeffType}(a, 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_on_s = Engine.LiesOn{CoeffType}(b, s)
-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")
+##b = Engine.Point{CoeffType}()
+##b_on_s = Engine.LiesOn{CoeffType}(b, s)
+##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")
 
-##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
+spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
 ##tangencies = [
 ##  Engine.AlignsWithBy{CoeffType}(
 ##    spheres[n],
@@ -48,31 +48,29 @@ println("Two points on a sphere: ", freedom, " degrees of freedom")
 ##  )
 ##  for n in 1:3
 ##]
-##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
-##ideal_tan_sph = Engine.realize(ctx_tan_sph)
-##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom")
+tangencies = [
+  Engine.AlignsWithBy{CoeffType}(spheres[1], spheres[2], CoeffType(-1)),
+  Engine.AlignsWithBy{CoeffType}(spheres[2], spheres[3], CoeffType(-1//2))
+]
+ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
+ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
+freedom = Engine.dimension(ideal_tan_sph)
+##println("Three mutually tangent spheres: $freedom degrees of freedom")
+println("Chain of three spheres: $freedom degrees of freedom")
 
 # --- test rational cut ---
 
-coordring = base_ring(ideal_ab_s)
+coordring = base_ring(ideal_tan_sph)
 vbls = Variable.(symbols(coordring))
 
 # 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])
-## [old] trivial_soln = fill(0, length(gens(coordring)))
-## [old] trivial_soln[sph_z_ind] .= 1
-## [old] println("trivial solutions: $trivial_soln")
+system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
 norm2 = vec -> real(dot(conj.(vec), vec))
-## [old] is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring))
 Random.seed!(6071)
-n_planes = 3
+n_planes = 16
 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)
@@ -97,7 +95,7 @@ for soln in samples
 end
 
 # show a sample solution
-function show_solution(vals)
+function show_solution(ctx, vals)
   # evaluate elements
   real_vals = real.(vals)
   disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]