From 86dbd9ea45ce9908c6895355c936a29ffe00629c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 27 Jan 2024 14:21:03 -0500 Subject: [PATCH 1/2] Order variables by coordinate and then element MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In other words, order coordinates like (rₛ₁, rₛ₂, sₛ₁, sₛ₂, xₛ₁, xₛ₂, xₚ₃, yₛ₁, yₛ₂, yₚ₃, zₛ₁, zₛ₂, zₚ₃) instead of like (rₛ₁, sₛ₁, xₛ₁, yₛ₁, zₛ₁, rₛ₂, sₛ₂, xₛ₂, yₛ₂, zₛ₂, xₚ₃, yₚ₃, zₚ₃). In the test cases, this really cuts down the size of the Gröbner basis. --- engine-proto/engine.jl | 76 +++++++++++++++++++++++++++++++----------- 1 file changed, 57 insertions(+), 19 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index f672745..7e68fe4 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -12,46 +12,68 @@ using Groebner abstract type Element{T} end mutable struct Point{T} <: Element{T} - coords::Union{Vector{MPolyRingElem{T}}, Nothing} + coords::Vector{MPolyRingElem{T}} vec::Union{Vector{MPolyRingElem{T}}, Nothing} rel::Nothing ## [to do] constructor argument never needed? Point{T}( - coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing ) where T = new(coords, vec, nothing) end -coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] +##coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] -function buildvec(pt::Point, coordqueue) - coordring = parent(coordqueue[1]) - pt.coords = splice!(coordqueue, 1:3) +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::Union{Vector{MPolyRingElem{T}}, Nothing} + coords::Vector{MPolyRingElem{T}} vec::Union{Vector{MPolyRingElem{T}}, Nothing} rel::Union{MPolyRingElem{T}, Nothing} + ## [to do] constructor argument never needed? Sphere{T}( - coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + 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 -coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] +##coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] -function buildvec(sph::Sphere, coordqueue) - coordring = parent(coordqueue[1]) - sph.coords = splice!(coordqueue, 1:5) +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(elem::Element, index) = coordnames[nameof(typeof(elem))][index] + +function pushcoordname!(coordnamelist, indexed_elem::Tuple{Any, Element}, coordindex) + elemindex, elem = indexed_elem + name = coordname(elem, coordindex) + if !isnothing(name) + subscript = Subscripts.sub(string(elemindex)) + push!(coordnamelist, Symbol(name, subscript)) + end +end + +function takecoord!(coordlist, indexed_elem::Tuple{Any, Element}, coordindex) + elem = indexed_elem[2] + if !isnothing(coordname(elem, coordindex)) + push!(elem.coords, popfirst!(coordlist)) + end +end + # --- primitive relations --- abstract type Relation{T} end @@ -99,22 +121,38 @@ function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T end function realize(ctx::Construction{T}) where T - # collect variable names + # collect coordinate names coordnamelist = Symbol[] elemenum = enumerate(ctx.elements) - for (index, elem) in elemenum - subscript = Subscripts.sub(string(index)) - append!(coordnamelist, - [Symbol(name, subscript) for name in coordnames(elem)] - ) + for coordindex in 1:5 + for indexed_elem in elemenum + pushcoordname!(coordnamelist, indexed_elem, coordindex) + end end + display(collect(elemenum)) + display(coordnamelist) + println() + # construct coordinate ring coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) + # retrieve coordinates + for (_, elem) in elemenum + empty!(elem.coords) + end + for coordindex in 1:5 + for indexed_elem in elemenum + takecoord!(coordqueue, indexed_elem, coordindex) + end + end + # construct coordinate vectors for (_, elem) in elemenum - buildvec(elem, coordqueue) + buildvec!(elem) + display(elem.coords) + display(elem.vec) + println() end # turn relations into equations From c29000d912a15896f98c33af5aefd691899721ad Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 28 Jan 2024 01:34:13 -0500 Subject: [PATCH 2/2] Write a simple solver for the hitting set problem I think we need this to find the dimension of the solution variety. --- engine-proto/hitting-set.jl | 110 ++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 engine-proto/hitting-set.jl diff --git a/engine-proto/hitting-set.jl b/engine-proto/hitting-set.jl new file mode 100644 index 0000000..e9aacf4 --- /dev/null +++ b/engine-proto/hitting-set.jl @@ -0,0 +1,110 @@ +module HittingSet + +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...) + T = eltype(wholeset) + unsorted_moves = [ + elt => Set(filter(s -> elt ∉ s, subsets)) + 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) + println(typeof(problems)) + 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] + subsets = 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(subsets) + if isa(problem, HittingSetProblem{T}) + println("Correct type") + else + println("Wrong type: ", typeof(problem)) + end + problem +end + +end \ No newline at end of file