Store elements in arrays to keep order stable

This seems to restore reproducibility.

engine-proto/Engine.Algebraic.jl

@@ -0,0 +1,204 @@
+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::Vector{Point{T}}
+  spheres::Vector{Sphere{T}}
+  relations::Vector{Relation{T}}
+  
+  function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{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
+  # [to do] the scaling constraint, as written, can be impossible to satisfy
+  #         when all of the spheres have to go through the origin
+  ##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)
+  ## [test] (nothing, eqns)
+end
+
+end

engine-proto/Engine.Numerical.jl

@@ -0,0 +1,55 @@
+module Numerical
+
+using Random: default_rng
+using LinearAlgebra
+using AbstractAlgebra
+using HomotopyContinuation:
+  Variable, Expression, AbstractSystem, System, LinearSubspace,
+  nvariables, isreal, witness_set, results
+import GLMakie
+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(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; 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(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}) =
+  GLMakie.Point3f([evaluate(u, vals) for u in pt.coords])
+
+function AbstractAlgebra.evaluate(sph::Sphere, vals::Vector{<:RingElement})
+  radius = 1 / evaluate(sph.coords[1], vals)
+  center = radius * [evaluate(u, vals) for u in sph.coords[3:end]]
+  GLMakie.Sphere(GLMakie.Point3f(center), radius)
+end
+
+end

engine-proto/Engine.jl

@@ -0,0 +1,135 @@
+include("HittingSet.jl")
+
+module Engine
+
+include("Engine.Algebraic.jl")
+include("Engine.Numerical.jl")
+
+using .Algebraic
+using .Numerical
+
+export Construction, mprod, codimension, dimension
+
+end
+
+# ~~~ sandbox setup ~~~
+
+using Random
+using Distributions
+using LinearAlgebra
+using AbstractAlgebra
+using HomotopyContinuation
+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]))
+##ideal_a_s = Engine.realize(ctx)
+##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")
+
+spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
+tangencies = [
+  Engine.AlignsWithBy{CoeffType}(
+    spheres[n],
+    spheres[mod1(n+1, length(spheres))],
+    CoeffType(-1)^n
+  )
+  for n in 1:3
+]
+##tangencies = [
+  ##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
+  ##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
+  ##Engine.LiesOn{CoeffType}(points[2], spheres[3]),
+  ##Engine.LiesOn{CoeffType}(points[2], spheres[1]),
+  ##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
+  ##Engine.LiesOn{CoeffType}(points[3], spheres[2])
+##]
+ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = 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(ideal_tan_sph)
+println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom")
+
+##points = [Engine.Point{CoeffType}() for _ in 1:3]
+##spheres = [Engine.Sphere{CoeffType}() for _ in 1:2]
+##ctx_joined = Engine.Construction{CoeffType}(
+##  elements = Set([points; spheres]),
+##  relations= Set([
+##    Engine.LiesOn{CoeffType}(pt, sph)
+##    for pt in points for sph in spheres
+##  ])
+##)
+##ideal_joined, eqns_joined = Engine.realize(ctx_joined)
+##freedom = Engine.dimension(ideal_joined)
+##println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom")
+
+# --- test rational cut ---
+
+coordring = base_ring(ideal_tan_sph)
+vbls = Variable.(symbols(coordring))
+
+# test a random witness set
+system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
+norm2 = vec -> real(dot(conj.(vec), vec))
+rng = MersenneTwister(6701)
+n_planes = 6
+samples = []
+for _ in 1:n_planes
+  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)
+    end
+  end
+end
+println("$(length(samples)) sample solutions:")
+##for soln in samples
+##  ## display([vbls round.(soln, digits = 6)]) ## [verbose]
+##  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
+
+# show a sample solution
+function show_solution(ctx, vals)
+  # evaluate elements
+  real_vals = real.(vals)
+  disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
+  disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
+  
+  # create scene
+  scene = Scene()
+  cam3d!(scene)
+  scatter!(scene, disp_points, color = :green)
+  for sph in disp_spheres
+    mesh!(scene, sph, color = :gray)
+  end
+  scene
+end

engine-proto/HittingSet.jl

@@ -0,0 +1,111 @@
+module HittingSet
+
+export HittingSetProblem, solve
+
+HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}}
+
+# `targets` should be a collection of Set objects
+function HittingSetProblem(targets, chosen = Set())
+  wholeset = union(targets...)
+  T = eltype(wholeset)
+  unsorted_moves = [
+    elt => Set(filter(s -> elt ∉ s, targets))
+    for elt in wholeset
+  ]
+  moves = sort(unsorted_moves, by = pair -> length(pair.second))
+  Set{T}(chosen) => moves
+end
+
+function Base.display(problem::HittingSetProblem{T}) where T
+  println("HittingSetProblem{$T}")
+  
+  chosen = problem.first
+  println(" {", join(string.(chosen), ", "), "}")
+  
+  moves = problem.second
+  for (choice, missed) in moves
+    println(" | ", choice)
+    for s in missed
+      println(" | | {", join(string.(s), ", "), "}")
+    end
+  end
+  println()
+end
+
+function solve(pblm::HittingSetProblem{T}, maxdepth = Inf) where T
+  problems = Dict(pblm)
+  while length(first(problems).first) < maxdepth
+    subproblems = typeof(problems)()
+    for (chosen, moves) in problems
+      if isempty(moves)
+        return chosen
+      else
+        for (choice, missed) in moves
+          to_be_chosen = union(chosen, Set([choice]))
+          if isempty(missed)
+            return to_be_chosen
+          elseif !haskey(subproblems, to_be_chosen)
+            push!(subproblems, HittingSetProblem(missed, to_be_chosen))
+          end
+        end
+      end
+    end
+    problems = subproblems
+  end
+  problems
+end
+
+function test(n = 1)
+  T = [Int64, Int64, Symbol, Symbol][n]
+  targets = Set{T}.([
+    [
+      [1, 3, 5],
+      [2, 3, 4],
+      [1, 4],
+      [2, 3, 4, 5],
+      [4, 5]
+    ],
+    # example from Amit Chakrabarti's graduate-level algorithms class (CS 105)
+    # notes by Valika K. Wan and Khanh Do Ba, Winter 2005
+    # https://www.cs.dartmouth.edu/~ac/Teach/CS105-Winter05/
+    [
+      [1, 3], [1,    4], [1,    5],
+      [1, 3], [1, 2, 4], [1, 2, 5],
+      [4, 3], [   2, 4], [   2, 5],
+      [6, 3], [6,    4], [      5]
+    ],
+    [
+      [:w, :x, :y],
+      [:x, :y, :z],
+      [:w, :z],
+      [:x, :y]
+    ],
+    # Wikipedia showcases this as an example of a problem where the greedy
+    # algorithm performs especially poorly
+    [
+      [:a, :x, :t1],
+      [:a, :y, :t2],
+      [:a, :y, :t3],
+      [:a, :z, :t4],
+      [:a, :z, :t5],
+      [:a, :z, :t6],
+      [:a, :z, :t7],
+      [:b, :x, :t8],
+      [:b, :y, :t9],
+      [:b, :y, :t10],
+      [:b, :z, :t11],
+      [:b, :z, :t12],
+      [:b, :z, :t13],
+      [:b, :z, :t14]
+    ]
+  ][n])
+  problem = HittingSetProblem(targets)
+  if isa(problem, HittingSetProblem{T})
+    println("Correct type")
+  else
+    println("Wrong type: ", typeof(problem))
+  end
+  problem
+end
+
+end