Compare commits
No commits in common. "1f173708eb57fc2305b9a17ab88bdf1375cc26fd" and "af1d31f6e61552b0320c7127ef3c720f75e40dea" have entirely different histories.
1f173708eb
...
af1d31f6e6
@ -1,201 +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)
|
|
||||||
|
|
||||||
# --- 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::Set{Point{T}}
|
|
||||||
spheres::Set{Sphere{T}}
|
|
||||||
relations::Set{Relation{T}}
|
|
||||||
|
|
||||||
function Construction{T}(; elements = Set{Element{T}}(), relations = Set{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
|
|
||||||
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
|
|
||||||
@ -1,40 +0,0 @@
|
|||||||
module Numerical
|
|
||||||
|
|
||||||
using LinearAlgebra
|
|
||||||
using AbstractAlgebra
|
|
||||||
using HomotopyContinuation
|
|
||||||
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(AbstractAlgebra.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)
|
|
||||||
# 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)...
|
|
||||||
))
|
|
||||||
cut = LinearSubspace(normals, fill(0., dim))
|
|
||||||
filter(isreal, results(witness_set(F, cut)))
|
|
||||||
end
|
|
||||||
|
|
||||||
end
|
|
||||||
@ -2,14 +2,229 @@ include("HittingSet.jl")
|
|||||||
|
|
||||||
module Engine
|
module Engine
|
||||||
|
|
||||||
include("Engine.Algebraic.jl")
|
|
||||||
include("Engine.Numerical.jl")
|
|
||||||
|
|
||||||
using .Algebraic
|
|
||||||
using .Numerical
|
|
||||||
|
|
||||||
export Construction, mprod, codimension, dimension
|
export Construction, mprod, codimension, dimension
|
||||||
|
|
||||||
|
import Subscripts
|
||||||
|
using LinearAlgebra
|
||||||
|
using AbstractAlgebra
|
||||||
|
using Groebner
|
||||||
|
using HomotopyContinuation: Variable, Expression, System
|
||||||
|
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)
|
||||||
|
|
||||||
|
# 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
|
||||||
|
|
||||||
|
## [to do] not needed right now
|
||||||
|
# create a ModelKit.System from a list of elements of a multivariate polynomial
|
||||||
|
# ring. the variable ordering is taken from the polynomial ring
|
||||||
|
##function System(eqns::AbstractVector{MPolyRingElem})
|
||||||
|
## if isempty(eqns)
|
||||||
|
## return System([])
|
||||||
|
## else
|
||||||
|
## variables = Variable.(symbols(parent(f)))
|
||||||
|
## return System(Expression.(eqns), variables = variables)
|
||||||
|
## end
|
||||||
|
##end
|
||||||
|
|
||||||
|
# --- 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::Set{Point{T}}
|
||||||
|
spheres::Set{Sphere{T}}
|
||||||
|
relations::Set{Relation{T}}
|
||||||
|
|
||||||
|
function Construction{T}(; elements = Set{Element{T}}(), relations = Set{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
|
||||||
|
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
|
end
|
||||||
|
|
||||||
# ~~~ sandbox setup ~~~
|
# ~~~ sandbox setup ~~~
|
||||||
@ -75,18 +290,33 @@ trivial_soln[sph_z_ind] .= 1
|
|||||||
println("trivial solutions: $trivial_soln")
|
println("trivial solutions: $trivial_soln")
|
||||||
norm2 = vec -> real(dot(conj.(vec), vec))
|
norm2 = vec -> real(dot(conj.(vec), vec))
|
||||||
is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring))
|
is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring))
|
||||||
##max_slope = 5
|
max_slope = 5
|
||||||
##binom = Binomial(2max_slope, 1/2)
|
binom = Binomial(2max_slope, 1/2)
|
||||||
Random.seed!(6071)
|
Random.seed!(6071)
|
||||||
n_planes = 36
|
n_planes = 36
|
||||||
for through_trivial in [false, true]
|
for through_trivial in [false, true]
|
||||||
samples = []
|
samples = []
|
||||||
for _ in 1:n_planes
|
for _ in 1:n_planes
|
||||||
real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
|
cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope
|
||||||
|
##cut_matrix = [
|
||||||
|
## 1 1 1 1 0 1 1 0 1 1 0;
|
||||||
|
## 1 2 1 2 0 1 1 0 1 1 0;
|
||||||
|
## 1 1 0 1 0 1 2 0 2 0 0
|
||||||
|
##]
|
||||||
|
## [verbose] display(cut_matrix)
|
||||||
|
if through_trivial
|
||||||
|
cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
|
||||||
|
## [verbose] display(cut_offset)
|
||||||
|
cut_subspace = LinearSubspace(cut_matrix, cut_offset)
|
||||||
|
else
|
||||||
|
cut_subspace = LinearSubspace(cut_matrix, fill(0, freedom))
|
||||||
|
end
|
||||||
|
wtns = witness_set(system, cut_subspace)
|
||||||
|
real_solns = solution.(filter(isreal, results(wtns)))
|
||||||
nontrivial_solns = filter(is_nontrivial, real_solns)
|
nontrivial_solns = filter(is_nontrivial, real_solns)
|
||||||
println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found")
|
println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found")
|
||||||
for soln in nontrivial_solns
|
for soln in nontrivial_solns
|
||||||
## [test] for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns))))
|
##[test] for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns))))
|
||||||
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
|
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
|
||||||
push!(samples, soln)
|
push!(samples, soln)
|
||||||
end
|
end
|
||||||
@ -99,7 +329,7 @@ for through_trivial in [false, true]
|
|||||||
end
|
end
|
||||||
println("$(length(samples)) sample solutions, not including the trivial one:")
|
println("$(length(samples)) sample solutions, not including the trivial one:")
|
||||||
for soln in samples
|
for soln in samples
|
||||||
## display([vbls round.(soln, digits = 6)]) ## [verbose]
|
## [verbose] display([vbls round.(soln, digits = 6)])
|
||||||
k_sq = abs2(soln[1])
|
k_sq = abs2(soln[1])
|
||||||
if abs2(soln[end-2]) > 1e-12
|
if abs2(soln[end-2]) > 1e-12
|
||||||
if k_sq < 1e-12
|
if k_sq < 1e-12
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user