From ae5db0f9eaa71a4f87b70e5d95efd43fa0fd9426 Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Thu, 15 Feb 2024 16:00:46 -0800
Subject: [PATCH] Make results reproducible

---
 engine-proto/Engine.Numerical.jl | 7 ++++---
 engine-proto/Engine.jl           | 6 +++---
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl
index 48fb682..d1e14bd 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,16 @@ 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)...
+    (normalize(randn(rng, nvariables(F))) for _ in 1:dim)...
   ))
   cut = LinearSubspace(normals, fill(0., dim))
-  filter(isreal, results(witness_set(F, cut)))
+  filter(isreal, results(witness_set(F, cut, seed = 0x1974abba)))
 end
 
 AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index d283e72..7146e80 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -89,11 +89,11 @@ vbls = Variable.(symbols(coordring))
 # test a random witness set
 system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls))
 norm2 = vec -> real(dot(conj.(vec), vec))
-Random.seed!(6071)
-n_planes = 36
+rng = MersenneTwister(6071)
+n_planes = 3
 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)