Move random cut routine into engine

Unit normals are uniformly distributed over the sphere.

engine-proto/Engine.Algebraic.jl

@@ -0,0 +1,201 @@
+module Algebraic
+
+export
+  codimension, dimension,
+  Construction, realize,
+  Element, Point, Sphere,
+  Relation, LiesOn, AlignsWithBy, mprod
+
+import Subscripts
+using LinearAlgebra
+using AbstractAlgebra
+using Groebner
+using ...HittingSet
+
+# --- commutative algebra ---
+
+# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate
+# polynomial rings when the coefficients are integers. we use Groebner to extend
+# support to rationals and to finite fields of prime order
+Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} =
+  Generic.Ideal{U}(base_ring(I), groebner(gens(I)))
+
+function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}}
+  leading = [exponent_vector(f, 1) for f in gens(I)]
+  targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading]
+  length(HittingSet.solve(HittingSetProblem(targets), maxdepth))
+end
+
+dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} =
+  length(gens(base_ring(I))) - codimension(I, maxdepth)
+
+# --- primitve elements ---
+
+abstract type Element{T} end
+
+mutable struct Point{T} <: Element{T}
+  coords::Vector{MPolyRingElem{T}}
+  vec::Union{Vector{MPolyRingElem{T}}, Nothing}
+  rel::Nothing
+  
+  ## [to do] constructor argument never needed?
+  Point{T}(
+    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
+    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
+  ) where T = new(coords, vec, nothing)
+end
+
+function buildvec!(pt::Point)
+  coordring = parent(pt.coords[1])
+  pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
+end
+
+mutable struct Sphere{T} <: Element{T}
+  coords::Vector{MPolyRingElem{T}}
+  vec::Union{Vector{MPolyRingElem{T}}, Nothing}
+  rel::Union{MPolyRingElem{T}, Nothing}
+  
+  ## [to do] constructor argument never needed?
+  Sphere{T}(
+    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
+    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
+    rel::Union{MPolyRingElem{T}, Nothing} = nothing
+  ) where T = new(coords, vec, rel)
+end
+
+function buildvec!(sph::Sphere)
+  coordring = parent(sph.coords[1])
+  sph.vec = sph.coords
+  sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
+end
+
+const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}(
+  nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ],
+  nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
+)
+
+coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index]
+
+function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex)
+  eltindex, elt = indexed_elt
+  name = coordname(elt, coordindex)
+  if !isnothing(name)
+    subscript = Subscripts.sub(string(eltindex))
+    push!(coordnamelist, Symbol(name, subscript))
+  end
+end
+
+function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex)
+  elt = indexed_elt[2]
+  if !isnothing(coordname(elt, coordindex))
+    push!(elt.coords, popfirst!(coordlist))
+  end
+end
+
+# --- primitive relations ---
+
+abstract type Relation{T} end
+
+mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end])
+
+# elements: point, sphere
+struct LiesOn{T} <: Relation{T}
+  elements::Vector{Element{T}}
+  
+  LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph])
+end
+
+equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec)
+
+# elements: sphere, sphere
+struct AlignsWithBy{T} <: Relation{T}
+  elements::Vector{Element{T}}
+  cos_angle::T
+  
+  AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle)
+end
+
+equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle
+
+# --- constructions ---
+
+mutable struct Construction{T}
+  points::Set{Point{T}}
+  spheres::Set{Sphere{T}}
+  relations::Set{Relation{T}}
+  
+  function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
+    allelements = union(elements, (rel.elements for rel in relations)...)
+    new{T}(
+      filter(elt -> isa(elt, Point), allelements),
+      filter(elt -> isa(elt, Sphere), allelements),
+      relations
+    )
+  end
+end
+
+function Base.push!(ctx::Construction{T}, elt::Point{T}) where T
+  push!(ctx.points, elt)
+end
+
+function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T
+  push!(ctx.spheres, elt)
+end
+
+function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
+  push!(ctx.relations, rel)
+  for elt in rel.elements
+    push!(ctx, elt)
+  end
+end
+
+function realize(ctx::Construction{T}) where T
+  # collect coordinate names
+  coordnamelist = Symbol[]
+  eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points)))
+  for coordindex in 1:5
+    for indexed_elt in eltenum
+      pushcoordname!(coordnamelist, indexed_elt, coordindex)
+    end
+  end
+  
+  # construct coordinate ring
+  coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
+  
+  # retrieve coordinates
+  for (_, elt) in eltenum
+    empty!(elt.coords)
+  end
+  for coordindex in 1:5
+    for indexed_elt in eltenum
+      takecoord!(coordqueue, indexed_elt, coordindex)
+    end
+  end
+  
+  # construct coordinate vectors
+  for (_, elt) in eltenum
+    buildvec!(elt)
+  end
+  
+  # turn relations into equations
+  eqns = vcat(
+    equation.(ctx.relations),
+    [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)]
+  )
+  
+  # add relations to center, orient, and scale the construction
+  if !isempty(ctx.points)
+    append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3])
+  end
+  if !isempty(ctx.spheres)
+    append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
+  end
+  n_elts = length(ctx.points) + length(ctx.spheres)
+  if n_elts > 0
+    push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
+  end
+  
+  (Generic.Ideal(coordring, eqns), eqns)
+end
+
+end

engine-proto/Engine.Numerical.jl

@@ -0,0 +1,40 @@
+module Numerical
+
+using LinearAlgebra
+using AbstractAlgebra
+using HomotopyContinuation
+using ..Algebraic
+
+# --- polynomial conversion ---
+
+# hat tip Sascha Timme
+#   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(AbstractAlgebra.coefficients(f), exponent_vectors(f))
+  sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data)
+end
+
+# create a ModelKit.System from an ideal in a multivariate polynomial ring. the
+# variable ordering is taken from the polynomial ring
+function System(I::Generic.Ideal)
+  eqns = Expression.(gens(I))
+  variables = Variable.(symbols(base_ring(I)))
+  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

@@ -2,229 +2,14 @@ include("HittingSet.jl")
 
 module Engine
 
+include("Engine.Algebraic.jl")
+include("Engine.Numerical.jl")
+
+using .Algebraic
+using .Numerical
+
 export Construction, mprod, codimension, dimension
 
-import Subscripts
-using LinearAlgebra
-using AbstractAlgebra
-using Groebner
-using HomotopyContinuation: Variable, Expression, System
-using ..HittingSet
-
-# --- commutative algebra ---
-
-# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate
-# polynomial rings when the coefficients are integers. we use Groebner to extend
-# support to rationals and to finite fields of prime order
-Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} =
-  Generic.Ideal{U}(base_ring(I), groebner(gens(I)))
-
-function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}}
-  leading = [exponent_vector(f, 1) for f in gens(I)]
-  targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading]
-  length(HittingSet.solve(HittingSetProblem(targets), maxdepth))
-end
-
-dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} =
-  length(gens(base_ring(I))) - codimension(I, maxdepth)
-
-# hat tip Sascha Timme
-#   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))
-  sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data)
-end
-
-# create a ModelKit.System from an ideal in a multivariate polynomial ring. the
-# variable ordering is taken from the polynomial ring
-function System(I::Generic.Ideal)
-  eqns = Expression.(gens(I))
-  variables = Variable.(symbols(base_ring(I)))
-  System(eqns, variables = variables)
-end
-
-## [to do] not needed right now
-# create a ModelKit.System from a list of elements of a multivariate polynomial
-# ring. the variable ordering is taken from the polynomial ring
-##function System(eqns::AbstractVector{MPolyRingElem})
-##  if isempty(eqns)
-##    return System([])
-##  else
-##    variables = Variable.(symbols(parent(f)))
-##    return System(Expression.(eqns), variables = variables)
-##  end
-##end
-
-# --- primitve elements ---
-
-abstract type Element{T} end
-
-mutable struct Point{T} <: Element{T}
-  coords::Vector{MPolyRingElem{T}}
-  vec::Union{Vector{MPolyRingElem{T}}, Nothing}
-  rel::Nothing
-  
-  ## [to do] constructor argument never needed?
-  Point{T}(
-    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
-    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
-  ) where T = new(coords, vec, nothing)
-end
-
-function buildvec!(pt::Point)
-  coordring = parent(pt.coords[1])
-  pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
-end
-
-mutable struct Sphere{T} <: Element{T}
-  coords::Vector{MPolyRingElem{T}}
-  vec::Union{Vector{MPolyRingElem{T}}, Nothing}
-  rel::Union{MPolyRingElem{T}, Nothing}
-  
-  ## [to do] constructor argument never needed?
-  Sphere{T}(
-    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
-    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
-    rel::Union{MPolyRingElem{T}, Nothing} = nothing
-  ) where T = new(coords, vec, rel)
-end
-
-function buildvec!(sph::Sphere)
-  coordring = parent(sph.coords[1])
-  sph.vec = sph.coords
-  sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
-end
-
-const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}(
-  nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ],
-  nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
-)
-
-coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index]
-
-function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex)
-  eltindex, elt = indexed_elt
-  name = coordname(elt, coordindex)
-  if !isnothing(name)
-    subscript = Subscripts.sub(string(eltindex))
-    push!(coordnamelist, Symbol(name, subscript))
-  end
-end
-
-function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex)
-  elt = indexed_elt[2]
-  if !isnothing(coordname(elt, coordindex))
-    push!(elt.coords, popfirst!(coordlist))
-  end
-end
-
-# --- primitive relations ---
-
-abstract type Relation{T} end
-
-mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end])
-
-# elements: point, sphere
-struct LiesOn{T} <: Relation{T}
-  elements::Vector{Element{T}}
-  
-  LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph])
-end
-
-equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec)
-
-# elements: sphere, sphere
-struct AlignsWithBy{T} <: Relation{T}
-  elements::Vector{Element{T}}
-  cos_angle::T
-  
-  AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle)
-end
-
-equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle
-
-# --- constructions ---
-
-mutable struct Construction{T}
-  points::Set{Point{T}}
-  spheres::Set{Sphere{T}}
-  relations::Set{Relation{T}}
-  
-  function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
-    allelements = union(elements, (rel.elements for rel in relations)...)
-    new{T}(
-      filter(elt -> isa(elt, Point), allelements),
-      filter(elt -> isa(elt, Sphere), allelements),
-      relations
-    )
-  end
-end
-
-function Base.push!(ctx::Construction{T}, elt::Point{T}) where T
-  push!(ctx.points, elt)
-end
-
-function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T
-  push!(ctx.spheres, elt)
-end
-
-function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
-  push!(ctx.relations, rel)
-  for elt in rel.elements
-    push!(ctx, elt)
-  end
-end
-
-function realize(ctx::Construction{T}) where T
-  # collect coordinate names
-  coordnamelist = Symbol[]
-  eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points)))
-  for coordindex in 1:5
-    for indexed_elt in eltenum
-      pushcoordname!(coordnamelist, indexed_elt, coordindex)
-    end
-  end
-  
-  # construct coordinate ring
-  coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
-  
-  # retrieve coordinates
-  for (_, elt) in eltenum
-    empty!(elt.coords)
-  end
-  for coordindex in 1:5
-    for indexed_elt in eltenum
-      takecoord!(coordqueue, indexed_elt, coordindex)
-    end
-  end
-  
-  # construct coordinate vectors
-  for (_, elt) in eltenum
-    buildvec!(elt)
-  end
-  
-  # turn relations into equations
-  eqns = vcat(
-    equation.(ctx.relations),
-    [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)]
-  )
-  
-  # add relations to center, orient, and scale the construction
-  if !isempty(ctx.points)
-    append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3])
-  end
-  if !isempty(ctx.spheres)
-    append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
-  end
-  n_elts = length(ctx.points) + length(ctx.spheres)
-  if n_elts > 0
-    push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
-  end
-  
-  (Generic.Ideal(coordring, eqns), eqns)
-end
-
 end
 
 # ~~~ sandbox setup ~~~
@@ -290,29 +75,14 @@ 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
+##max_slope = 5
-binom = Binomial(2max_slope, 1/2)
+##binom = Binomial(2max_slope, 1/2)
 Random.seed!(6071)
 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
+    real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
-    ##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
-      cut_subspace = LinearSubspace(cut_matrix, fill(0, freedom))
-    end
-    wtns = witness_set(system, cut_subspace)
-    real_solns = solution.(filter(isreal, results(wtns)))
     nontrivial_solns = filter(is_nontrivial, real_solns)
     println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found")
     for soln in nontrivial_solns
@@ -329,7 +99,7 @@ for through_trivial in [false, true]
   end
   println("$(length(samples)) sample solutions, not including the trivial one:")
   for soln in samples
-    ## [verbose] display([vbls round.(soln, digits = 6)])
+    ## display([vbls round.(soln, digits = 6)]) ## [verbose]
     k_sq = abs2(soln[1])
     if abs2(soln[end-2]) > 1e-12
       if k_sq < 1e-12