Compare commits
No commits in common. "macaulay2" and "main" have entirely different histories.
|
@ -1,227 +0,0 @@
|
||||||
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)
|
|
||||||
|
|
||||||
m2_ordering(R::MPolyRing) = Dict(
|
|
||||||
:lex => :Lex,
|
|
||||||
:deglex => :GLex,
|
|
||||||
:degrevlex => :GRevLex
|
|
||||||
)[ordering(R)]
|
|
||||||
|
|
||||||
string_m2(ring::MPolyRing) =
|
|
||||||
"QQ[$(join(symbols(ring), ", ")), MonomialOrder => $(m2_ordering(ring))]"
|
|
||||||
|
|
||||||
string_m2(f::MPolyRingElem) =
|
|
||||||
replace(string(f), "//" => "/")
|
|
||||||
|
|
||||||
# --- 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
|
|
||||||
|
|
||||||
# output options:
|
|
||||||
# nothing - find a Gröbner basis
|
|
||||||
# :m2 - write a system of polynomials to a Macaulay2 file
|
|
||||||
function realize(ctx::Construction{T}; output = nothing) 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
|
|
||||||
|
|
||||||
if output == :m2
|
|
||||||
file = open("macaulay2/construction.m2", "w")
|
|
||||||
write(file, string(
|
|
||||||
"coordring = $(string_m2(coordring))\n",
|
|
||||||
"eqns = {\n $(join(string_m2.(eqns), ",\n "))\n}"
|
|
||||||
))
|
|
||||||
close(file)
|
|
||||||
else
|
|
||||||
return (Generic.Ideal(coordring, eqns), eqns)
|
|
||||||
end
|
|
||||||
end
|
|
||||||
|
|
||||||
end
|
|
|
@ -1,53 +0,0 @@
|
||||||
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 = 0x1974abba)))
|
|
||||||
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
|
|
|
@ -1,77 +0,0 @@
|
||||||
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}
|
|
||||||
|
|
||||||
spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
|
|
||||||
tangencies = [
|
|
||||||
Engine.AlignsWithBy{CoeffType}(
|
|
||||||
spheres[n],
|
|
||||||
spheres[mod1(n+1, length(spheres))],
|
|
||||||
CoeffType(1)
|
|
||||||
)
|
|
||||||
for n in 1:3
|
|
||||||
]
|
|
||||||
ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies)
|
|
||||||
##ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
|
|
||||||
Engine.realize(ctx_tan_sph, output = :m2)
|
|
||||||
##freedom = Engine.dimension(ideal_tan_sph)
|
|
||||||
##println("Three mutually tangent 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(6071)
|
|
||||||
##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("Found $(length(samples)) sample solutions")
|
|
||||||
|
|
||||||
# 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
|
|
|
@ -1,111 +0,0 @@
|
||||||
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
|
|
|
@ -1,3 +0,0 @@
|
||||||
needsPackage "TriangularSets"
|
|
||||||
|
|
||||||
mprod = (v, w) -> (v#0*w#1 + w#0*v#1) / 2 - v#2*w#2 - v#3*w#3 - v#4*w#4
|
|
Loading…
Reference in New Issue