dyna3/engine-proto/Engine.Algebraic.jl

203 lines
5.9 KiB
Julia
Raw Normal View History

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)
end
end