forked from glen/dyna3
a6da6f9925
All the random number generators seems to be seeded, so why aren't the results reproducible?
55 lines
No EOL
1.9 KiB
Julia
55 lines
No EOL
1.9 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 = 0x8af341df)))
|
|
##filter(isreal, results(witness_set(F, seed = 0x8af341df)))
|
|
results(witness_set(F, seed = 0x8af341df))
|
|
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 |