Compare commits
3 Commits
95c0ff14b2
...
06872a04af
Author | SHA1 | Date | |
---|---|---|---|
|
06872a04af | ||
|
becefe0c47 | ||
|
34358a8728 |
@ -225,6 +225,8 @@ end
|
||||
|
||||
# ~~~ sandbox setup ~~~
|
||||
|
||||
using Random
|
||||
using Distributions
|
||||
using AbstractAlgebra
|
||||
using HomotopyContinuation
|
||||
|
||||
@ -243,7 +245,8 @@ Engine.push!(ctx, b)
|
||||
Engine.push!(ctx, s)
|
||||
Engine.push!(ctx, b_on_s)
|
||||
ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
|
||||
println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freedom")
|
||||
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 = [
|
||||
@ -260,33 +263,54 @@ println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of f
|
||||
|
||||
# --- test rational cut ---
|
||||
|
||||
cut = [
|
||||
sum(vcat([1, 1, 1] .* a.coords, [1, 1, 1, 1, 1] .* (s.coords - [0, 0, 0, 0, 1])))
|
||||
sum(vcat([2, 1, 1] .* a.coords, [1, 2, 1, 1, 1] .* (s.coords - [0, 0, 0, 0, 1])))
|
||||
sum(vcat([1, 2, 0] .* a.coords, [1, 1, 0, 1, 2] .* (s.coords - [0, 0, 0, 0, 1])))
|
||||
]
|
||||
cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut])
|
||||
cut_dim = Engine.dimension(cut_ideal_ab_s)
|
||||
println("Two points on a sphere, after cut: ", cut_dim, " degrees of freedom")
|
||||
if cut_dim == 0
|
||||
vbls = Variable.(symbols(base_ring(ideal_ab_s)))
|
||||
cut_system = System([eqns_ab_s; cut], variables = vbls)
|
||||
cut_result = HomotopyContinuation.solve(cut_system)
|
||||
println("non-singular solutions:")
|
||||
for soln in solutions(cut_result)
|
||||
display(soln)
|
||||
end
|
||||
println("singular solutions:")
|
||||
for sing in singular(cut_result)
|
||||
display(sing.solution)
|
||||
end
|
||||
coordring = base_ring(ideal_ab_s)
|
||||
vbls = Variable.(symbols(coordring))
|
||||
##cut_system = CompiledSystem(System([eqns_ab_s; cut], variables = vbls))
|
||||
##cut_result = HomotopyContinuation.solve(cut_system)
|
||||
##println("non-singular solutions:")
|
||||
##for soln in solutions(cut_result)
|
||||
## display(soln)
|
||||
##end
|
||||
##println("singular solutions:")
|
||||
##for sing in singular(cut_result)
|
||||
## display(sing.solution)
|
||||
##end
|
||||
|
||||
# test corresponding witness set
|
||||
cut_matrix = [1 1 1 1 0 1 1 0 1 1 0; 1 2 1 2 0 1 1 0 1 1 0; 1 1 0 1 0 1 2 0 2 0 0]
|
||||
cut_subspace = LinearSubspace(cut_matrix, [1, 1, 2])
|
||||
witness = witness_set(System(eqns_ab_s, variables = vbls), cut_subspace)
|
||||
println("witness solutions:")
|
||||
for wtns in solutions(witness)
|
||||
display(wtns)
|
||||
# test a random witness set
|
||||
system = CompiledSystem(System(eqns_ab_s, variables = vbls))
|
||||
max_slope = 2
|
||||
binom = Binomial(2max_slope, 1/2)
|
||||
Random.seed!(6071)
|
||||
samples = []
|
||||
for _ in 1:3
|
||||
cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope
|
||||
##cut_matrix = [
|
||||
## 1 1 1 1 0 1 1 0 1 1 0;
|
||||
## 1 2 1 2 0 1 1 0 1 1 0;
|
||||
## 1 1 0 1 0 1 2 0 2 0 0
|
||||
##]
|
||||
sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
|
||||
cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
|
||||
println("sphere z variables: ", vbls[sph_z_ind])
|
||||
display(cut_matrix)
|
||||
display(cut_offset)
|
||||
cut_subspace = LinearSubspace(cut_matrix, cut_offset)
|
||||
wtns = witness_set(system, cut_subspace)
|
||||
append!(samples, solution.(filter(isreal, results(wtns))))
|
||||
end
|
||||
println("$(length(samples)) sample solutions:")
|
||||
for soln in samples
|
||||
display([vbls round.(soln, digits = 6)])
|
||||
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
|
Loading…
Reference in New Issue
Block a user