glen/dyna3
1
0
Fork 0
Code Issues 11 Pull Requests 1 Actions Projects Releases 2 Wiki Activity

Try displaying a chain of spheres

Browse Source 
For three mutually tangent spheres, I couldn't find real solutions.
This commit is contained in:
Aaron Fenyes 2024-02-12 21:14:07 -05:00
parent 6cf07dc6a1
commit a450f701fb
1 changed files with 27 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
engine-proto/Engine.jl
View File

@ -23,23 +23,23 @@ using GLMakie
CoeffType = Rational{Int64} CoeffType = Rational{Int64}
a = Engine.Point{CoeffType}() ##a = Engine.Point{CoeffType}()
s = Engine.Sphere{CoeffType}() ##s = Engine.Sphere{CoeffType}()
a_on_s = Engine.LiesOn{CoeffType}(a, s) ##a_on_s = Engine.LiesOn{CoeffType}(a, s)
ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s])) ##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
##ideal_a_s = Engine.realize(ctx) ##ideal_a_s = Engine.realize(ctx)
##println("A point on a sphere: ", Engine.dimension(ideal_a_s), " degrees of freedom") ##println("A point on a sphere: $(Engine.dimension(ideal_a_s)) degrees of freedom")
b = Engine.Point{CoeffType}() ##b = Engine.Point{CoeffType}()
b_on_s = Engine.LiesOn{CoeffType}(b, s) ##b_on_s = Engine.LiesOn{CoeffType}(b, s)
Engine.push!(ctx, b) ##Engine.push!(ctx, b)
Engine.push!(ctx, s) ##Engine.push!(ctx, s)
Engine.push!(ctx, b_on_s) ##Engine.push!(ctx, b_on_s)
ideal_ab_s, eqns_ab_s = Engine.realize(ctx) ##ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
freedom = Engine.dimension(ideal_ab_s) ##freedom = Engine.dimension(ideal_ab_s)
println("Two points on a sphere: ", freedom, " degrees of freedom") ##println("Two points on a sphere: $freedom degrees of freedom")
##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
##tangencies = [ ##tangencies = [
## Engine.AlignsWithBy{CoeffType}( ## Engine.AlignsWithBy{CoeffType}(
## spheres[n], ## spheres[n],
@ -48,31 +48,29 @@ println("Two points on a sphere: ", freedom, " degrees of freedom")
## ) ## )
## for n in 1:3 ## for n in 1:3
##] ##]
##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) tangencies = [
##ideal_tan_sph = Engine.realize(ctx_tan_sph) Engine.AlignsWithBy{CoeffType}(spheres[1], spheres[2], CoeffType(-1)),
##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom") Engine.AlignsWithBy{CoeffType}(spheres[2], spheres[3], CoeffType(-1//2))
]
ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
freedom = Engine.dimension(ideal_tan_sph)
##println("Three mutually tangent spheres: $freedom degrees of freedom")
println("Chain of three spheres: $freedom degrees of freedom")
# --- test rational cut --- # --- test rational cut ---
coordring = base_ring(ideal_ab_s) coordring = base_ring(ideal_tan_sph)
vbls = Variable.(symbols(coordring)) vbls = Variable.(symbols(coordring))
# test a random witness set # test a random witness set
system = CompiledSystem(System(eqns_ab_s, variables = vbls)) system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
println("sphere z variables: ", vbls[sph_z_ind])
## [old] trivial_soln = fill(0, length(gens(coordring)))
## [old] trivial_soln[sph_z_ind] .= 1
## [old] println("trivial solutions: $trivial_soln")
norm2 = vec -> real(dot(conj.(vec), vec)) norm2 = vec -> real(dot(conj.(vec), vec))
## [old] is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring))
Random.seed!(6071) Random.seed!(6071)
n_planes = 3 n_planes = 16
samples = [] samples = []
for _ in 1:n_planes for _ in 1:n_planes
real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
## [old] nontrivial_solns = filter(is_nontrivial, real_solns)
## [old] println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found")
for soln in real_solns for soln in real_solns
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
push!(samples, soln) push!(samples, soln)
@ -97,7 +95,7 @@ for soln in samples
end end
# show a sample solution # show a sample solution
function show_solution(vals) function show_solution(ctx, vals)
# evaluate elements # evaluate elements
real_vals = real.(vals) real_vals = real.(vals)
disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points] disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]