forked from glen/dyna3
ba365174d3
I'm not sure why the solver wasn't working before. It might've been just an unlucky random number draw.
135 lines
No EOL
4.2 KiB
Julia
135 lines
No EOL
4.2 KiB
Julia
include("HittingSet.jl")
|
|
|
|
module Engine
|
|
|
|
include("Engine.Algebraic.jl")
|
|
include("Engine.Numerical.jl")
|
|
|
|
using .Algebraic
|
|
using .Numerical
|
|
|
|
export Construction, mprod, codimension, dimension
|
|
|
|
end
|
|
|
|
# ~~~ sandbox setup ~~~
|
|
|
|
using Random
|
|
using Distributions
|
|
using LinearAlgebra
|
|
using AbstractAlgebra
|
|
using HomotopyContinuation
|
|
using GLMakie
|
|
|
|
CoeffType = Rational{Int64}
|
|
|
|
##a = Engine.Point{CoeffType}()
|
|
##s = Engine.Sphere{CoeffType}()
|
|
##a_on_s = Engine.LiesOn{CoeffType}(a, s)
|
|
##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
|
|
##ideal_a_s = Engine.realize(ctx)
|
|
##println("A point on a sphere: $(Engine.dimension(ideal_a_s)) degrees of freedom")
|
|
|
|
##b = Engine.Point{CoeffType}()
|
|
##b_on_s = Engine.LiesOn{CoeffType}(b, s)
|
|
##Engine.push!(ctx, b)
|
|
##Engine.push!(ctx, s)
|
|
##Engine.push!(ctx, b_on_s)
|
|
##ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
|
|
##freedom = Engine.dimension(ideal_ab_s)
|
|
##println("Two points on a sphere: $freedom degrees of freedom")
|
|
|
|
spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
|
|
tangencies = [
|
|
Engine.AlignsWithBy{CoeffType}(
|
|
spheres[n],
|
|
spheres[mod1(n+1, length(spheres))],
|
|
CoeffType(-1)^n
|
|
)
|
|
for n in 1:3
|
|
]
|
|
##tangencies = [
|
|
##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
|
|
##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
|
|
##Engine.LiesOn{CoeffType}(points[2], spheres[3]),
|
|
##Engine.LiesOn{CoeffType}(points[2], spheres[1]),
|
|
##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
|
|
##Engine.LiesOn{CoeffType}(points[3], spheres[2])
|
|
##]
|
|
ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies)
|
|
ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
|
|
##small_eqns_tan_sph = eqns_tan_sph
|
|
##small_eqns_tan_sph = [
|
|
## eqns_tan_sph;
|
|
## spheres[2].coords - [1, 0, 0, 0, 1];
|
|
## spheres[3].coords - [1, 0, 0, 0, -1];
|
|
##]
|
|
##small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph)
|
|
freedom = Engine.dimension(ideal_tan_sph)
|
|
println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom")
|
|
|
|
##points = [Engine.Point{CoeffType}() for _ in 1:3]
|
|
##spheres = [Engine.Sphere{CoeffType}() for _ in 1:2]
|
|
##ctx_joined = Engine.Construction{CoeffType}(
|
|
## elements = Set([points; spheres]),
|
|
## relations= Set([
|
|
## Engine.LiesOn{CoeffType}(pt, sph)
|
|
## for pt in points for sph in spheres
|
|
## ])
|
|
##)
|
|
##ideal_joined, eqns_joined = Engine.realize(ctx_joined)
|
|
##freedom = Engine.dimension(ideal_joined)
|
|
##println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom")
|
|
|
|
# --- test rational cut ---
|
|
|
|
coordring = base_ring(ideal_tan_sph)
|
|
vbls = Variable.(symbols(coordring))
|
|
|
|
# test a random witness set
|
|
system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
|
|
norm2 = vec -> real(dot(conj.(vec), vec))
|
|
rng = MersenneTwister(6071)
|
|
n_planes = 3
|
|
samples = []
|
|
for _ in 1:n_planes
|
|
real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng))
|
|
for soln in real_solns
|
|
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
|
|
push!(samples, soln)
|
|
end
|
|
end
|
|
end
|
|
println("$(length(samples)) sample solutions:")
|
|
##for soln in samples
|
|
## ## display([vbls round.(soln, digits = 6)]) ## [verbose]
|
|
## k_sq = abs2(soln[1])
|
|
## if abs2(soln[end-2]) > 1e-12
|
|
## if k_sq < 1e-12
|
|
## println(" center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
|
|
## else
|
|
## sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
|
|
## println(" center on z axis: r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
|
|
## end
|
|
## else
|
|
## sum_sq = sum(soln[[4, 7, 10]] .^ 2)
|
|
## println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
|
|
## end
|
|
##end
|
|
|
|
# show a sample solution
|
|
function show_solution(ctx, vals)
|
|
# evaluate elements
|
|
real_vals = real.(vals)
|
|
disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
|
|
disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
|
|
|
|
# create scene
|
|
scene = Scene()
|
|
cam3d!(scene)
|
|
scatter!(scene, disp_points, color = :green)
|
|
for sph in disp_spheres
|
|
mesh!(scene, sph, color = :gray)
|
|
end
|
|
scene
|
|
end |