From a6da6f9925590cf9382fab0e6d20be470f638abe Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Thu, 15 Feb 2024 14:17:03 -0800
Subject: [PATCH] Investigate why witness sets aren't reproducible
All the random number generators seems to be seeded, so why aren't the results
reproducible?
---
engine-proto/Engine.Numerical.jl | 15 +++++++++------
engine-proto/Engine.jl | 24 ++++++++++++------------
2 files changed, 21 insertions(+), 18 deletions(-)
diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl
index 48fb682..c54e71c 100644
--- a/engine-proto/Engine.Numerical.jl
+++ b/engine-proto/Engine.Numerical.jl
@@ -1,5 +1,6 @@
module Numerical
+using Random: default_rng
using LinearAlgebra
using AbstractAlgebra
using HomotopyContinuation:
@@ -28,16 +29,18 @@ end
# --- sampling ---
-function real_samples(F::AbstractSystem, dim)
+function real_samples(F::AbstractSystem, dim; rng = default_rng())
# 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)))
+ ##normals = transpose(hcat(
+ ## (normalize(randn(rng, nvariables(F))) for _ in 1:dim)...
+ ##))
+ ##cut = LinearSubspace(normals, fill(0., dim))
+ ##filter(isreal, results(witness_set(F, cut, seed = 0x8af341df)))
+ ##filter(isreal, results(witness_set(F, seed = 0x8af341df)))
+ results(witness_set(F, seed = 0x8af341df))
end
AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index 0410f6d..be201cf 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -59,13 +59,13 @@ tangencies = [
ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
##small_eqns_tan_sph = eqns_tan_sph
-small_eqns_tan_sph = [
- eqns_tan_sph;
- spheres[2].coords - [1, 0, 0, 0, 1];
- spheres[3].coords - [1, 0, 0, 0, -1];
-]
-small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph)
-freedom = Engine.dimension(small_ideal_tan_sph)
+##small_eqns_tan_sph = [
+## eqns_tan_sph;
+## spheres[2].coords - [1, 0, 0, 0, 1];
+## spheres[3].coords - [1, 0, 0, 0, -1];
+##]
+##small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph)
+freedom = Engine.dimension(ideal_tan_sph)
println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom")
##points = [Engine.Point{CoeffType}() for _ in 1:3]
@@ -83,17 +83,17 @@ println("Three mutually tangent spheres, with two fixed: $freedom degrees of fre
# --- test rational cut ---
-coordring = base_ring(small_ideal_tan_sph)
+coordring = base_ring(ideal_tan_sph)
vbls = Variable.(symbols(coordring))
# test a random witness set
-system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls))
+system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
norm2 = vec -> real(dot(conj.(vec), vec))
-Random.seed!(6071)
-n_planes = 36
+rng = MersenneTwister(6701)
+n_planes = 6
samples = []
for _ in 1:n_planes
- real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
+ real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng))
for soln in real_solns
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
push!(samples, soln)