Start interface to Macaulay2

I did this to try out Macaulay2's "triangularize" function, but that
turns out to use Maple for rings with more than three variables.
This commit is contained in:
Aaron Fenyes 2024-02-16 12:47:06 -08:00
parent 3170a933e4
commit 74529048de
3 changed files with 64 additions and 36 deletions

View File

@ -29,6 +29,18 @@ 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
@ -149,7 +161,10 @@ function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
end
end
function realize(ctx::Construction{T}) where T
# 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)))
@ -197,7 +212,16 @@ function realize(ctx::Construction{T}) where T
push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
end
(Generic.Ideal(coordring, eqns), eqns)
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

View File

@ -33,44 +33,45 @@ tangencies = [
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)
freedom = Engine.dimension(ideal_tan_sph)
println("Three mutually tangent spheres: $freedom degrees of freedom")
##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))
##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")
##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
##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

View File

@ -0,0 +1,3 @@
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