Compare commits

...

36 Commits

Author SHA1 Message Date
Aaron Fenyes
74529048de 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.
2024-02-16 12:47:06 -08:00
Aaron Fenyes
3170a933e4 Clean up example of three mutually tangent spheres 2024-02-15 17:16:37 -08:00
Aaron Fenyes
f2000e5731 Test different sign patterns for cosines
It seems like there are real solutions if and only if the product of the
cosines is positive.
2024-02-15 16:25:09 -08:00
Aaron Fenyes
ba365174d3 Find real solutions for three mutually tangent spheres
I'm not sure why the solver wasn't working before. It might've been just
an unlucky random number draw.
2024-02-15 16:16:06 -08:00
Aaron Fenyes
ae5db0f9ea Make results reproducible 2024-02-15 16:00:46 -08:00
Aaron Fenyes
8d8bc9162c Store elements in arrays to keep order stable
This seems to restore reproducibility.
2024-02-15 15:42:26 -08:00
Aaron Fenyes
291d5c8ff6 Study mutually tangent spheres with two fixed 2024-02-15 13:28:01 -08:00
Aaron Fenyes
e41bcc7e13 Explore the performance wall
Three points on two spheres is too much.
2024-02-13 04:02:14 -05:00
Aaron Fenyes
31d5e7e864 Play with two points on two spheres
Guess conditions that make the scaling constraint impossible to satisfy.
2024-02-12 22:48:16 -05:00
Aaron Fenyes
a450f701fb Try displaying a chain of spheres
For three mutually tangent spheres, I couldn't find real solutions.
2024-02-12 21:14:07 -05:00
Aaron Fenyes
6cf07dc6a1 Evaluate and display elements 2024-02-12 20:34:12 -05:00
Aaron Fenyes
1f173708eb Move random cut routine into engine 2024-02-10 17:39:26 -05:00
Aaron Fenyes
6f18d4efcc Test lots of uniformly distributed hyperplanes 2024-02-10 15:10:48 -05:00
Aaron Fenyes
621c4c5776 Try uniformly distributed hyperplane orientations
Unit normals are uniformly distributed over the sphere.
2024-02-10 15:02:26 -05:00
Aaron Fenyes
b3b7c2026d Separate the algebraic and numerical parts of the engine 2024-02-10 14:50:50 -05:00
Aaron Fenyes
af1d31f6e6 Test a scale constraint
In all but a few cases (for example, a single point on a plane), we
should be able to us the radius-coradius boost symmetry to make the
average co-radius—representing the "overall scale"—roughly one.
2024-02-10 14:21:52 -05:00
Aaron Fenyes
8e33987f59 Systematically try out different cut planes 2024-02-10 13:46:01 -05:00
Aaron Fenyes
06872a04af Say how many sample solutions we found 2024-02-10 01:06:06 -05:00
Aaron Fenyes
becefe0c47 Try switching to compiled system 2024-02-10 00:59:50 -05:00
Aaron Fenyes
34358a8728 Find witnesses on random rational hyperplanes
Choose hyperplanes that go through the trivial solution.
2024-02-09 23:44:10 -05:00
Aaron Fenyes
95c0ff14b2 Show explicitly that all coefficients are 1 in first cut equation 2024-02-09 17:09:43 -05:00
Aaron Fenyes
f97090c997 Try a cut that goes through the trivial solution
The previous cut was supposed to do this, but I was missing some parentheses.
2024-02-08 01:58:12 -05:00
Aaron Fenyes
45aaaafc8f Seek sample solutions by cutting with a hyperplane
The example hyperplane yields a single solution, with multiplicity six. You can
find it analytically by hand, and homotopy continuation finds it numerically.
2024-02-08 01:53:55 -05:00
Aaron Fenyes
43cbf8a3a0 Add relations to center and orient the construction 2024-02-05 00:10:13 -05:00
Aaron Fenyes
21f09c4a4d Switch element abbreviation from "elem" to "elt" 2024-02-04 16:08:13 -05:00
Aaron Fenyes
a3f3f6a31b Order spheres before points within each coordinate block
In the cases I've tried so far, this leads to substantially smaller
Gröbner bases.
2024-02-01 16:13:22 -05:00
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
Aaron Fenyes
0731c7aac1 Correct relation equations 2024-01-29 12:41:07 -05:00
Aaron Fenyes
59a527af43 Correct Minkowski product; build chain of three spheres 2024-01-29 12:28:57 -05:00
Aaron Fenyes
c29000d912 Write a simple solver for the hitting set problem
I think we need this to find the dimension of the solution variety.
2024-01-28 01:34:13 -05:00
Aaron Fenyes
86dbd9ea45 Order variables by coordinate and then element
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.
2024-01-27 14:21:03 -05:00
Aaron Fenyes
463a3b21e1 Realize relations as equations 2024-01-27 12:28:29 -05:00
Aaron Fenyes
4d5aa3b327 Realize geometric elements as symbolic vectors 2024-01-26 11:14:32 -05:00
Aaron Fenyes
b864cf7866 Start drafting engine prototype 2024-01-24 11:16:24 -05:00
5 changed files with 471 additions and 0 deletions

View File

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

View File

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

77
engine-proto/Engine.jl Normal file
View File

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

111
engine-proto/HittingSet.jl Normal file
View File

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

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