dyna3/engine-proto/Engine.Numerical.jl

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