53 lines
1.8 KiB
Julia
53 lines
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
|