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