d7dbee4c05
We're using the Gram matrix engine for the next stage of development, so the algebraic engine shouldn't be at the top level anymore.
53 lines
No EOL
1.8 KiB
Julia
53 lines
No EOL
1.8 KiB
Julia
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 |