Compare commits

...

3 Commits

Author SHA1 Message Date
Aaron Fenyes
65d23fb667 Use module names as filenames
You're right: this naming convention seems to be standard for Julia
modules now.
2024-01-30 02:49:33 -05:00
Aaron Fenyes
4e02ee16fc Find dimension of solution variety 2024-01-30 02:45:14 -05:00
Aaron Fenyes
6349f298ae Extend AbstractAlgebra ideals to rational coefficients
The extension should also let us work over finite fields of prime order,
although we don't need to do that.
2024-01-29 19:11:21 -05:00
2 changed files with 38 additions and 23 deletions

View File

@ -1,3 +1,5 @@
include("HittingSet.jl")
module Engine module Engine
export Construction, mprod export Construction, mprod
@ -6,6 +8,24 @@ import Subscripts
using LinearAlgebra using LinearAlgebra
using AbstractAlgebra using AbstractAlgebra
using Groebner 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 --- # --- primitve elements ---
@ -23,8 +43,6 @@ mutable struct Point{T} <: Element{T}
) where T = new(coords, vec, nothing) ) where T = new(coords, vec, nothing)
end end
##coordnames(_::Point) = [:xₚ, :yₚ, :zₚ]
function buildvec!(pt::Point) function buildvec!(pt::Point)
coordring = parent(pt.coords[1]) coordring = parent(pt.coords[1])
pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
@ -43,8 +61,6 @@ mutable struct Sphere{T} <: Element{T}
) where T = new(coords, vec, rel) ) where T = new(coords, vec, rel)
end end
##coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
function buildvec!(sph::Sphere) function buildvec!(sph::Sphere)
coordring = parent(sph.coords[1]) coordring = parent(sph.coords[1])
sph.vec = sph.coords sph.vec = sph.coords
@ -130,10 +146,6 @@ function realize(ctx::Construction{T}) where T
end end
end end
display(collect(elemenum))
display(coordnamelist)
println()
# construct coordinate ring # construct coordinate ring
coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
@ -150,16 +162,14 @@ function realize(ctx::Construction{T}) where T
# construct coordinate vectors # construct coordinate vectors
for (_, elem) in elemenum for (_, elem) in elemenum
buildvec!(elem) buildvec!(elem)
display(elem.coords)
display(elem.vec)
println()
end end
# turn relations into equations # turn relations into equations
vcat( eqns = vcat(
equation.(ctx.relations), equation.(ctx.relations),
[elem.rel for elem in ctx.elements if !isnothing(elem.rel)] [elem.rel for elem in ctx.elements if !isnothing(elem.rel)]
) )
Generic.Ideal(coordring, eqns)
end end
end end
@ -172,22 +182,26 @@ a = Engine.Point{CoeffType}()
s = Engine.Sphere{CoeffType}() s = Engine.Sphere{CoeffType}()
a_on_s = Engine.LiesOn{CoeffType}(a, s) a_on_s = Engine.LiesOn{CoeffType}(a, s)
ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s])) ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
eqns_a_s = Engine.realize(ctx) 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 = Engine.Point{CoeffType}()
b_on_s = Engine.LiesOn{CoeffType}(b, s) b_on_s = Engine.LiesOn{CoeffType}(b, s)
Engine.push!(ctx, b) Engine.push!(ctx, b)
Engine.push!(ctx, s) Engine.push!(ctx, s)
Engine.push!(ctx, b_on_s) Engine.push!(ctx, b_on_s)
eqns_ab_s = Engine.realize(ctx) 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] spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
tangencies = [ tangencies = [
Engine.AlignsWithBy{CoeffType}( Engine.AlignsWithBy{CoeffType}(
spheres[n], spheres[n],
spheres[mod1(n+1, length(spheres))], spheres[mod1(n+1, length(spheres))],
-1//1 CoeffType(-1//1)
) )
for n in 1:3 for n in 1:3
] ]
ctx_chain = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) 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")

View File

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