Find dimension of solution variety

engine-proto/engine.jl

@@ -1,3 +1,5 @@
+include("hitting-set.jl")
+
 module Engine
 
 export Construction, mprod
@@ -6,27 +8,25 @@ 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 coefficients are integers. in `reduce_gens`, the
-# `lmnode` constructor requires < to be defined on the coefficients, and the
-# `reducer_size` heuristic requires `ndigits` to be defined on the coefficients.
-# this patch for `reducer_size` removes the `ndigits` dependency
-##function Generic.reducer_size(f::T) where {U <: MPolyRingElem{<:FieldElement}, V, N, T <: Generic.lmnode{U, V, N}}
-##   if f.size != 0.0
-##      return f.size
-##   end
-##   return 0.0 + sum(j^2 for j in 1:length(f.poly))
-##end
-
 # 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
@@ -183,6 +183,7 @@ 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 freeom")
 
 b = Engine.Point{CoeffType}()
 b_on_s = Engine.LiesOn{CoeffType}(b, s)
@@ -190,6 +191,7 @@ Engine.push!(ctx, b)
 Engine.push!(ctx, s)
 Engine.push!(ctx, b_on_s)
 ideal_ab_s = Engine.realize(ctx)
+println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freeom")
 
 spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
 tangencies = [
@@ -200,5 +202,6 @@ tangencies = [
   )
   for n in 1:3
 ]
-ctx_chain = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
-ideal_chain = Engine.realize(ctx_chain)
+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 freeom")

engine-proto/hitting-set.jl

@@ -1,13 +1,15 @@
 module HittingSet
 
+export HittingSetProblem, solve
+
 HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}}
 
-# `subsets` should be a collection of Set objects
-function HittingSetProblem(subsets, chosen = Set())
-  wholeset = union(subsets...)
+# `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, subsets))
+    elt => Set(filter(s -> elt ∉ s, targets))
     for elt in wholeset
   ]
   moves = sort(unsorted_moves, by = pair -> length(pair.second))
@@ -32,7 +34,6 @@ end
 
 function solve(pblm::HittingSetProblem{T}, maxdepth = Inf) where T
   problems = Dict(pblm)
-  println(typeof(problems))
   while length(first(problems).first) < maxdepth
     subproblems = typeof(problems)()
     for (chosen, moves) in problems
@@ -56,7 +57,7 @@ end
 
 function test(n = 1)
   T = [Int64, Int64, Symbol, Symbol][n]
-  subsets = Set{T}.([
+  targets = Set{T}.([
     [
       [1, 3, 5],
       [2, 3, 4],
@@ -98,7 +99,7 @@ function test(n = 1)
       [:b, :z, :t14]
     ]
   ][n])
-  problem = HittingSetProblem(subsets)
+  problem = HittingSetProblem(targets)
   if isa(problem, HittingSetProblem{T})
     println("Correct type")
   else