Compare commits
7 Commits
e41bcc7e13
...
74529048de
Author | SHA1 | Date | |
---|---|---|---|
|
74529048de | ||
|
3170a933e4 | ||
|
f2000e5731 | ||
|
ba365174d3 | ||
|
ae5db0f9ea | ||
|
8d8bc9162c | ||
|
291d5c8ff6 |
engine-proto
@ -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
|
||||
@ -120,11 +132,11 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) -
|
||||
# --- constructions ---
|
||||
|
||||
mutable struct Construction{T}
|
||||
points::Set{Point{T}}
|
||||
spheres::Set{Sphere{T}}
|
||||
relations::Set{Relation{T}}
|
||||
points::Vector{Point{T}}
|
||||
spheres::Vector{Sphere{T}}
|
||||
relations::Vector{Relation{T}}
|
||||
|
||||
function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where 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),
|
||||
@ -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,8 +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
|
||||
|
||||
## [test] (Generic.Ideal(coordring, eqns), eqns)
|
||||
(nothing, 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
|
@ -1,5 +1,6 @@
|
||||
module Numerical
|
||||
|
||||
using Random: default_rng
|
||||
using LinearAlgebra
|
||||
using AbstractAlgebra
|
||||
using HomotopyContinuation:
|
||||
@ -28,16 +29,16 @@ end
|
||||
|
||||
# --- sampling ---
|
||||
|
||||
function real_samples(F::AbstractSystem, dim)
|
||||
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(nvariables(F))) for _ in 1:dim)...
|
||||
(normalize(randn(rng, nvariables(F))) for _ in 1:dim)...
|
||||
))
|
||||
cut = LinearSubspace(normals, fill(0., dim))
|
||||
filter(isreal, results(witness_set(F, cut)))
|
||||
filter(isreal, results(witness_set(F, cut, seed = 0x1974abba)))
|
||||
end
|
||||
|
||||
AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
|
||||
|
@ -23,106 +23,55 @@ using GLMakie
|
||||
|
||||
CoeffType = Rational{Int64}
|
||||
|
||||
##a = Engine.Point{CoeffType}()
|
||||
##s = Engine.Sphere{CoeffType}()
|
||||
##a_on_s = Engine.LiesOn{CoeffType}(a, s)
|
||||
##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
|
||||
##ideal_a_s = Engine.realize(ctx)
|
||||
##println("A point on a sphere: $(Engine.dimension(ideal_a_s)) degrees of freedom")
|
||||
|
||||
##b = Engine.Point{CoeffType}()
|
||||
##b_on_s = Engine.LiesOn{CoeffType}(b, s)
|
||||
##Engine.push!(ctx, b)
|
||||
##Engine.push!(ctx, s)
|
||||
##Engine.push!(ctx, b_on_s)
|
||||
##ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
|
||||
##freedom = Engine.dimension(ideal_ab_s)
|
||||
##println("Two points on a sphere: $freedom degrees of freedom")
|
||||
|
||||
##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
|
||||
##tangencies = [
|
||||
## Engine.AlignsWithBy{CoeffType}(
|
||||
## spheres[n],
|
||||
## spheres[mod1(n+1, length(spheres))],
|
||||
## CoeffType(-1//1)
|
||||
## )
|
||||
## for n in 1:3
|
||||
##]
|
||||
##tangencies = [
|
||||
##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
|
||||
##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
|
||||
##Engine.LiesOn{CoeffType}(points[2], spheres[3]),
|
||||
##Engine.LiesOn{CoeffType}(points[2], spheres[1]),
|
||||
##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
|
||||
##Engine.LiesOn{CoeffType}(points[3], spheres[2])
|
||||
##]
|
||||
##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
|
||||
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")
|
||||
|
||||
points = [Engine.Point{CoeffType}() for _ in 1:3]
|
||||
spheres = [Engine.Sphere{CoeffType}() for _ in 1:2]
|
||||
ctx_joined = Engine.Construction{CoeffType}(
|
||||
elements = Set([points; spheres]),
|
||||
relations= Set([
|
||||
Engine.LiesOn{CoeffType}(pt, sph)
|
||||
for pt in points for sph in spheres
|
||||
])
|
||||
)
|
||||
ideal_joined, eqns_joined = Engine.realize(ctx_joined)
|
||||
freedom = Engine.dimension(ideal_joined)
|
||||
println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom")
|
||||
|
||||
# --- test rational cut ---
|
||||
|
||||
coordring = base_ring(ideal_joined)
|
||||
vbls = Variable.(symbols(coordring))
|
||||
##coordring = base_ring(ideal_tan_sph)
|
||||
##vbls = Variable.(symbols(coordring))
|
||||
|
||||
# test a random witness set
|
||||
system = CompiledSystem(System(eqns_joined, variables = vbls))
|
||||
norm2 = vec -> real(dot(conj.(vec), vec))
|
||||
Random.seed!(6071)
|
||||
n_planes = 3
|
||||
samples = []
|
||||
for _ in 1:n_planes
|
||||
real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
|
||||
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("$(length(samples)) sample solutions:")
|
||||
for soln in samples
|
||||
## display([vbls round.(soln, digits = 6)]) ## [verbose]
|
||||
k_sq = abs2(soln[1])
|
||||
if abs2(soln[end-2]) > 1e-12
|
||||
if k_sq < 1e-12
|
||||
println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
|
||||
else
|
||||
sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
|
||||
println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
|
||||
end
|
||||
else
|
||||
sum_sq = sum(soln[[4, 7, 10]] .^ 2)
|
||||
println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
|
||||
end
|
||||
end
|
||||
##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
|
3
engine-proto/macaulay2/engine.m2
Normal file
3
engine-proto/macaulay2/engine.m2
Normal 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
|
Loading…
Reference in New Issue
Block a user