Engine prototype #13
@ -227,6 +227,7 @@ end
|
||||
|
||||
using Random
|
||||
using Distributions
|
||||
using LinearAlgebra
|
||||
using AbstractAlgebra
|
||||
using HomotopyContinuation
|
||||
|
||||
@ -278,39 +279,60 @@ vbls = Variable.(symbols(coordring))
|
||||
|
||||
# test a random witness set
|
||||
system = CompiledSystem(System(eqns_ab_s, variables = vbls))
|
||||
max_slope = 2
|
||||
sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
|
||||
println("sphere z variables: ", vbls[sph_z_ind])
|
||||
trivial_soln = fill(0, length(gens(coordring)))
|
||||
trivial_soln[sph_z_ind] .= 1
|
||||
println("trivial solutions: $trivial_soln")
|
||||
norm2 = vec -> real(dot(conj.(vec), vec))
|
||||
is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring))
|
||||
max_slope = 5
|
||||
binom = Binomial(2max_slope, 1/2)
|
||||
Random.seed!(6071)
|
||||
samples = []
|
||||
for _ in 1:3
|
||||
n_planes = 36
|
||||
for through_trivial in [false, true]
|
||||
samples = []
|
||||
for _ in 1:n_planes
|
||||
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))
|
||||
## [verbose] display(cut_matrix)
|
||||
if through_trivial
|
||||
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)
|
||||
## [verbose] display(cut_offset)
|
||||
cut_subspace = LinearSubspace(cut_matrix, cut_offset)
|
||||
else
|
||||
cut_subspace = LinearSubspace(cut_matrix, fill(0, 3))
|
||||
end
|
||||
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)])
|
||||
for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns))))
|
||||
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
|
||||
push!(samples, soln)
|
||||
end
|
||||
end
|
||||
end
|
||||
if through_trivial
|
||||
println("--- planes through trivial solution ---")
|
||||
else
|
||||
println("--- planes through origin ---")
|
||||
end
|
||||
println("$(length(samples)) sample solutions, not including the trivial one:")
|
||||
for soln in samples
|
||||
## [verbose] 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))}")
|
||||
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))")
|
||||
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))")
|
||||
println(" center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
|
||||
end
|
||||
end
|
||||
end
|
Loading…
Reference in New Issue
Block a user