From b864cf7866e3399da3fcf861d752849a2c8dc5b5 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Wed, 24 Jan 2024 11:16:24 -0500 Subject: [PATCH 01/32] Start drafting engine prototype --- engine-proto/engine.jl | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 engine-proto/engine.jl diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl new file mode 100644 index 0000000..8b82b47 --- /dev/null +++ b/engine-proto/engine.jl @@ -0,0 +1,34 @@ +module Engine + +export Construction, Sphere, mprod, point + +using LinearAlgebra +using Groebner + +mutable struct Construction + nextid::Int64 + + Construction(; nextid = 0) = new(nextid) +end + +struct Sphere{T<:Number} + vec::Vector{T} + id + + function Sphere(vec::Vector{T}, ctx::Construction) where T <: Number + id = ctx.nextid + ctx.nextid += 1 + new{T}(vec, id) + end +end + +function mprod(sv::Sphere, sw::Sphere) + v = sv.vec + w = sw.vec + v[1]*w[2] + v[2]*w[1] - dot(v[3:end], w[3:end]) +end + +point(pt::Vector{<:Number}, ctx::Construction) = + Sphere([one(eltype(pt)), dot(pt, pt), pt...], ctx) + +end \ No newline at end of file From 4d5aa3b327f680cb7d1e2478bd499faefdab1004 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 26 Jan 2024 11:14:32 -0500 Subject: [PATCH 02/32] Realize geometric elements as symbolic vectors --- engine-proto/engine.jl | 107 ++++++++++++++++++++++++++++++++++------- 1 file changed, 89 insertions(+), 18 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index 8b82b47..df75fbe 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -1,34 +1,105 @@ module Engine -export Construction, Sphere, mprod, point +export Construction, mprod +import Subscripts using LinearAlgebra +using AbstractAlgebra using Groebner -mutable struct Construction - nextid::Int64 +# --- primitve elements --- + +mutable struct Point{T} + coords::Union{Vector{MPolyRingElem{T}}, Nothing} + vec::Union{Vector{MPolyRingElem{T}}, Nothing} - Construction(; nextid = 0) = new(nextid) + ## [to do] constructor argument never needed? + Point{T}(vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing) where T = new(vec) end -struct Sphere{T<:Number} - vec::Vector{T} - id +coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] + +function buildvec(pt::Point, coordqueue) + pt.coords = splice!(coordqueue, 1:3) + coordring = parent(coordqueue[1]) + pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] +end + +mutable struct Sphere{T} + coords::Union{Vector{MPolyRingElem{T}}, Nothing} + vec::Union{Vector{MPolyRingElem{T}}, Nothing} - function Sphere(vec::Vector{T}, ctx::Construction) where T <: Number - id = ctx.nextid - ctx.nextid += 1 - new{T}(vec, id) + Sphere{T}(vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing) where T = new(vec) +end + +coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] + +function buildvec(sph::Sphere, coordqueue) + sph.coords = splice!(coordqueue, 1:5) + sph.vec = sph.coords +end + +# --- primitive relations --- + +abstract type Relation{T} end + +mprod(v, w) = v[1]*w[2] + w[1]*v[2] - dot(v[3:end], w[3:end]) + +struct LiesOn{T} <: Relation{T} + pt::Point{T} + sph::Sphere{T} +end + +struct AlignsWithBy{T} <: Relation{T} + sph_v::Sphere{T} + sph_w::Sphere{T} + cos_angle::T +end + +# --- constructions --- + +mutable struct Construction{T} + points::Vector{Point{T}} + spheres::Vector{Sphere{T}} + + Construction{T}(; points = Point{T}[], spheres = Sphere{T}[]) where T = new{T}(points, spheres) +end + +function Base.push!(ctx::Construction{T}, elem::Point{T}) where T + push!(ctx.points, elem) +end + +function Base.push!(ctx::Construction{T}, elem::Sphere{T}) where T + push!(ctx.spheres, elem) +end + +function realize(ctx::Construction{T}) where T + # collect variable names + allcoordnames = Symbol[] + elements = vcat(ctx.points, ctx.spheres) + for (index, elem) in enumerate(elements) + subscript = Subscripts.sub(string(index)) + append!(allcoordnames, + [Symbol(name, subscript) for name in coordnames(elem)] + ) + end + + # construct coordinate ring + coordring, coordqueue = polynomial_ring(parent_type(T)(), allcoordnames) + + # construct coordinate vectors + for elem in elements + buildvec(elem, coordqueue) end end -function mprod(sv::Sphere, sw::Sphere) - v = sv.vec - w = sw.vec - v[1]*w[2] + v[2]*w[1] - dot(v[3:end], w[3:end]) end -point(pt::Vector{<:Number}, ctx::Construction) = - Sphere([one(eltype(pt)), dot(pt, pt), pt...], ctx) +# ~~~ sandbox setup ~~~ -end \ No newline at end of file +a = Engine.Point{Rational{Int64}}() +b = Engine.Point{Rational{Int64}}() +s = Engine.Sphere{Rational{Int64}}() +ctx = Engine.Construction{Rational{Int64}}(points = [a]) +Engine.push!(ctx, b) +Engine.push!(ctx, s) \ No newline at end of file From 463a3b21e1fa1c1a531645b4f9cfafbb2fbd2034 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 27 Jan 2024 12:28:29 -0500 Subject: [PATCH 03/32] Realize relations as equations --- engine-proto/engine.jl | 87 ++++++++++++++++++++++++++++++------------ 1 file changed, 62 insertions(+), 25 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index df75fbe..f672745 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -9,34 +9,47 @@ using Groebner # --- primitve elements --- -mutable struct Point{T} +abstract type Element{T} end + +mutable struct Point{T} <: Element{T} coords::Union{Vector{MPolyRingElem{T}}, Nothing} vec::Union{Vector{MPolyRingElem{T}}, Nothing} + rel::Nothing ## [to do] constructor argument never needed? - Point{T}(vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing) where T = new(vec) + Point{T}( + coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing + ) where T = new(coords, vec, nothing) end coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] function buildvec(pt::Point, coordqueue) - pt.coords = splice!(coordqueue, 1:3) coordring = parent(coordqueue[1]) + pt.coords = splice!(coordqueue, 1:3) pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] end -mutable struct Sphere{T} +mutable struct Sphere{T} <: Element{T} coords::Union{Vector{MPolyRingElem{T}}, Nothing} vec::Union{Vector{MPolyRingElem{T}}, Nothing} + rel::Union{MPolyRingElem{T}, Nothing} - Sphere{T}(vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing) where T = new(vec) + Sphere{T}( + coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + rel::Union{MPolyRingElem{T}, Nothing} = nothing + ) where T = new(coords, vec, rel) end coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] function buildvec(sph::Sphere, coordqueue) + coordring = parent(coordqueue[1]) sph.coords = splice!(coordqueue, 1:5) sph.vec = sph.coords + sph.rel = mprod(sph.coords, sph.coords) + one(coordring) end # --- primitive relations --- @@ -45,52 +58,70 @@ abstract type Relation{T} end mprod(v, w) = v[1]*w[2] + w[1]*v[2] - dot(v[3:end], w[3:end]) +# elements: point, sphere struct LiesOn{T} <: Relation{T} - pt::Point{T} - sph::Sphere{T} + elements::Vector{Element{T}} + + LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph]) end +equation(rel::LiesOn) = dot(rel.elements[1].vec, rel.elements[2].vec) + +# elements: sphere, sphere struct AlignsWithBy{T} <: Relation{T} - sph_v::Sphere{T} - sph_w::Sphere{T} + elements::Vector{Element{T}} cos_angle::T + + LiesOn{T}(sph1::Point{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle) end +equation(rel::AlignsWithBy) = dot(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle + # --- constructions --- mutable struct Construction{T} - points::Vector{Point{T}} - spheres::Vector{Sphere{T}} + elements::Set{Element{T}} + relations::Set{Relation{T}} - Construction{T}(; points = Point{T}[], spheres = Sphere{T}[]) where T = new{T}(points, spheres) + function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T + allelements = union(elements, (rel.elements for rel in relations)...) + new{T}(allelements, relations) + end end -function Base.push!(ctx::Construction{T}, elem::Point{T}) where T - push!(ctx.points, elem) +function Base.push!(ctx::Construction{T}, elem::Element{T}) where T + push!(ctx.elements, elem) end -function Base.push!(ctx::Construction{T}, elem::Sphere{T}) where T - push!(ctx.spheres, elem) +function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T + push!(ctx.relations, rel) + union!(ctx.elements, rel.elements) end function realize(ctx::Construction{T}) where T # collect variable names - allcoordnames = Symbol[] - elements = vcat(ctx.points, ctx.spheres) - for (index, elem) in enumerate(elements) + coordnamelist = Symbol[] + elemenum = enumerate(ctx.elements) + for (index, elem) in elemenum subscript = Subscripts.sub(string(index)) - append!(allcoordnames, + append!(coordnamelist, [Symbol(name, subscript) for name in coordnames(elem)] ) end # construct coordinate ring - coordring, coordqueue = polynomial_ring(parent_type(T)(), allcoordnames) + coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) # construct coordinate vectors - for elem in elements + for (_, elem) in elemenum buildvec(elem, coordqueue) end + + # turn relations into equations + vcat( + equation.(ctx.relations), + [elem.rel for elem in ctx.elements if !isnothing(elem.rel)] + ) end end @@ -98,8 +129,14 @@ end # ~~~ sandbox setup ~~~ a = Engine.Point{Rational{Int64}}() -b = Engine.Point{Rational{Int64}}() s = Engine.Sphere{Rational{Int64}}() -ctx = Engine.Construction{Rational{Int64}}(points = [a]) +a_on_s = Engine.LiesOn{Rational{Int64}}(a, s) +ctx = Engine.Construction{Rational{Int64}}(elements = Set([a]), relations= Set([a_on_s])) +eqns_a_s = Engine.realize(ctx) + +b = Engine.Point{Rational{Int64}}() +b_on_s = Engine.LiesOn{Rational{Int64}}(b, s) Engine.push!(ctx, b) -Engine.push!(ctx, s) \ No newline at end of file +Engine.push!(ctx, s) +Engine.push!(ctx, b_on_s) +eqns_ab_s = Engine.realize(ctx) \ No newline at end of file From 86dbd9ea45ce9908c6895355c936a29ffe00629c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 27 Jan 2024 14:21:03 -0500 Subject: [PATCH 04/32] Order variables by coordinate and then element MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In other words, order coordinates like (rₛ₁, rₛ₂, sₛ₁, sₛ₂, xₛ₁, xₛ₂, xₚ₃, yₛ₁, yₛ₂, yₚ₃, zₛ₁, zₛ₂, zₚ₃) instead of like (rₛ₁, sₛ₁, xₛ₁, yₛ₁, zₛ₁, rₛ₂, sₛ₂, xₛ₂, yₛ₂, zₛ₂, xₚ₃, yₚ₃, zₚ₃). In the test cases, this really cuts down the size of the Gröbner basis. --- engine-proto/engine.jl | 76 +++++++++++++++++++++++++++++++----------- 1 file changed, 57 insertions(+), 19 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index f672745..7e68fe4 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -12,46 +12,68 @@ using Groebner abstract type Element{T} end mutable struct Point{T} <: Element{T} - coords::Union{Vector{MPolyRingElem{T}}, Nothing} + coords::Vector{MPolyRingElem{T}} vec::Union{Vector{MPolyRingElem{T}}, Nothing} rel::Nothing ## [to do] constructor argument never needed? Point{T}( - coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing ) where T = new(coords, vec, nothing) end -coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] +##coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] -function buildvec(pt::Point, coordqueue) - coordring = parent(coordqueue[1]) - pt.coords = splice!(coordqueue, 1:3) +function buildvec!(pt::Point) + coordring = parent(pt.coords[1]) pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] end mutable struct Sphere{T} <: Element{T} - coords::Union{Vector{MPolyRingElem{T}}, Nothing} + coords::Vector{MPolyRingElem{T}} vec::Union{Vector{MPolyRingElem{T}}, Nothing} rel::Union{MPolyRingElem{T}, Nothing} + ## [to do] constructor argument never needed? Sphere{T}( - coords::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, rel::Union{MPolyRingElem{T}, Nothing} = nothing ) where T = new(coords, vec, rel) end -coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] +##coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] -function buildvec(sph::Sphere, coordqueue) - coordring = parent(coordqueue[1]) - sph.coords = splice!(coordqueue, 1:5) +function buildvec!(sph::Sphere) + coordring = parent(sph.coords[1]) sph.vec = sph.coords sph.rel = mprod(sph.coords, sph.coords) + one(coordring) end +const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}( + nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ], + nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] +) + +coordname(elem::Element, index) = coordnames[nameof(typeof(elem))][index] + +function pushcoordname!(coordnamelist, indexed_elem::Tuple{Any, Element}, coordindex) + elemindex, elem = indexed_elem + name = coordname(elem, coordindex) + if !isnothing(name) + subscript = Subscripts.sub(string(elemindex)) + push!(coordnamelist, Symbol(name, subscript)) + end +end + +function takecoord!(coordlist, indexed_elem::Tuple{Any, Element}, coordindex) + elem = indexed_elem[2] + if !isnothing(coordname(elem, coordindex)) + push!(elem.coords, popfirst!(coordlist)) + end +end + # --- primitive relations --- abstract type Relation{T} end @@ -99,22 +121,38 @@ function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T end function realize(ctx::Construction{T}) where T - # collect variable names + # collect coordinate names coordnamelist = Symbol[] elemenum = enumerate(ctx.elements) - for (index, elem) in elemenum - subscript = Subscripts.sub(string(index)) - append!(coordnamelist, - [Symbol(name, subscript) for name in coordnames(elem)] - ) + for coordindex in 1:5 + for indexed_elem in elemenum + pushcoordname!(coordnamelist, indexed_elem, coordindex) + end end + display(collect(elemenum)) + display(coordnamelist) + println() + # construct coordinate ring coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) + # retrieve coordinates + for (_, elem) in elemenum + empty!(elem.coords) + end + for coordindex in 1:5 + for indexed_elem in elemenum + takecoord!(coordqueue, indexed_elem, coordindex) + end + end + # construct coordinate vectors for (_, elem) in elemenum - buildvec(elem, coordqueue) + buildvec!(elem) + display(elem.coords) + display(elem.vec) + println() end # turn relations into equations From c29000d912a15896f98c33af5aefd691899721ad Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 28 Jan 2024 01:34:13 -0500 Subject: [PATCH 05/32] Write a simple solver for the hitting set problem I think we need this to find the dimension of the solution variety. --- engine-proto/hitting-set.jl | 110 ++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 engine-proto/hitting-set.jl diff --git a/engine-proto/hitting-set.jl b/engine-proto/hitting-set.jl new file mode 100644 index 0000000..e9aacf4 --- /dev/null +++ b/engine-proto/hitting-set.jl @@ -0,0 +1,110 @@ +module HittingSet + +HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}} + +# `subsets` should be a collection of Set objects +function HittingSetProblem(subsets, chosen = Set()) + wholeset = union(subsets...) + T = eltype(wholeset) + unsorted_moves = [ + elt => Set(filter(s -> elt ∉ s, subsets)) + for elt in wholeset + ] + moves = sort(unsorted_moves, by = pair -> length(pair.second)) + Set{T}(chosen) => moves +end + +function Base.display(problem::HittingSetProblem{T}) where T + println("HittingSetProblem{$T}") + + chosen = problem.first + println(" {", join(string.(chosen), ", "), "}") + + moves = problem.second + for (choice, missed) in moves + println(" | ", choice) + for s in missed + println(" | | {", join(string.(s), ", "), "}") + end + end + println() +end + +function solve(pblm::HittingSetProblem{T}, maxdepth = Inf) where T + problems = Dict(pblm) + println(typeof(problems)) + while length(first(problems).first) < maxdepth + subproblems = typeof(problems)() + for (chosen, moves) in problems + if isempty(moves) + return chosen + else + for (choice, missed) in moves + to_be_chosen = union(chosen, Set([choice])) + if isempty(missed) + return to_be_chosen + elseif !haskey(subproblems, to_be_chosen) + push!(subproblems, HittingSetProblem(missed, to_be_chosen)) + end + end + end + end + problems = subproblems + end + problems +end + +function test(n = 1) + T = [Int64, Int64, Symbol, Symbol][n] + subsets = Set{T}.([ + [ + [1, 3, 5], + [2, 3, 4], + [1, 4], + [2, 3, 4, 5], + [4, 5] + ], + # example from Amit Chakrabarti's graduate-level algorithms class (CS 105) + # notes by Valika K. Wan and Khanh Do Ba, Winter 2005 + # https://www.cs.dartmouth.edu/~ac/Teach/CS105-Winter05/ + [ + [1, 3], [1, 4], [1, 5], + [1, 3], [1, 2, 4], [1, 2, 5], + [4, 3], [ 2, 4], [ 2, 5], + [6, 3], [6, 4], [ 5] + ], + [ + [:w, :x, :y], + [:x, :y, :z], + [:w, :z], + [:x, :y] + ], + # Wikipedia showcases this as an example of a problem where the greedy + # algorithm performs especially poorly + [ + [:a, :x, :t1], + [:a, :y, :t2], + [:a, :y, :t3], + [:a, :z, :t4], + [:a, :z, :t5], + [:a, :z, :t6], + [:a, :z, :t7], + [:b, :x, :t8], + [:b, :y, :t9], + [:b, :y, :t10], + [:b, :z, :t11], + [:b, :z, :t12], + [:b, :z, :t13], + [:b, :z, :t14] + ] + ][n]) + problem = HittingSetProblem(subsets) + if isa(problem, HittingSetProblem{T}) + println("Correct type") + else + println("Wrong type: ", typeof(problem)) + end + problem +end + +end \ No newline at end of file From 59a527af43b869da239a4c3789a52dea52a0da09 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 29 Jan 2024 12:28:45 -0500 Subject: [PATCH 06/32] Correct Minkowski product; build chain of three spheres --- engine-proto/engine.jl | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index 7e68fe4..f2dc981 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -78,7 +78,7 @@ end abstract type Relation{T} end -mprod(v, w) = v[1]*w[2] + w[1]*v[2] - dot(v[3:end], w[3:end]) +mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end]) # elements: point, sphere struct LiesOn{T} <: Relation{T} @@ -94,7 +94,7 @@ struct AlignsWithBy{T} <: Relation{T} elements::Vector{Element{T}} cos_angle::T - LiesOn{T}(sph1::Point{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle) + AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle) end equation(rel::AlignsWithBy) = dot(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle @@ -166,15 +166,28 @@ end # ~~~ sandbox setup ~~~ -a = Engine.Point{Rational{Int64}}() -s = Engine.Sphere{Rational{Int64}}() -a_on_s = Engine.LiesOn{Rational{Int64}}(a, s) -ctx = Engine.Construction{Rational{Int64}}(elements = Set([a]), relations= Set([a_on_s])) +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])) eqns_a_s = Engine.realize(ctx) -b = Engine.Point{Rational{Int64}}() -b_on_s = Engine.LiesOn{Rational{Int64}}(b, s) +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) -eqns_ab_s = Engine.realize(ctx) \ No newline at end of file +eqns_ab_s = Engine.realize(ctx) + +spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] +tangencies = [ + Engine.AlignsWithBy{CoeffType}( + spheres[n], + spheres[mod1(n+1, length(spheres))], + -1//1 + ) + for n in 1:3 +] +ctx_chain = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) \ No newline at end of file From 0731c7aac18fee58daa0675e84eb065c029984d7 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 29 Jan 2024 12:41:07 -0500 Subject: [PATCH 07/32] Correct relation equations --- engine-proto/engine.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index f2dc981..2f7294a 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -87,7 +87,7 @@ struct LiesOn{T} <: Relation{T} LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph]) end -equation(rel::LiesOn) = dot(rel.elements[1].vec, rel.elements[2].vec) +equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec) # elements: sphere, sphere struct AlignsWithBy{T} <: Relation{T} @@ -97,7 +97,7 @@ struct AlignsWithBy{T} <: Relation{T} AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle) end -equation(rel::AlignsWithBy) = dot(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle +equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle # --- constructions --- From 6349f298ae723c21361f7450207b64f8e7f5e19c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 29 Jan 2024 19:11:21 -0500 Subject: [PATCH 08/32] Extend AbstractAlgebra ideals to rational coefficients The extension should also let us work over finite fields of prime order, although we don't need to do that. --- engine-proto/engine.jl | 43 ++++++++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index 2f7294a..b282fdd 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -7,6 +7,26 @@ using LinearAlgebra using AbstractAlgebra using Groebner +# --- commutative algebra --- + +# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate +# polynomial rings when coefficients are integers. in `reduce_gens`, the +# `lmnode` constructor requires < to be defined on the coefficients, and the +# `reducer_size` heuristic requires `ndigits` to be defined on the coefficients. +# this patch for `reducer_size` removes the `ndigits` dependency +##function Generic.reducer_size(f::T) where {U <: MPolyRingElem{<:FieldElement}, V, N, T <: Generic.lmnode{U, V, N}} +## if f.size != 0.0 +## return f.size +## end +## return 0.0 + sum(j^2 for j in 1:length(f.poly)) +##end + +# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate +# polynomial rings when the coefficients are integers. we use Groebner to extend +# support to rationals and to finite fields of prime order +Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} = + Generic.Ideal{U}(base_ring(I), groebner(gens(I))) + # --- primitve elements --- abstract type Element{T} end @@ -23,8 +43,6 @@ mutable struct Point{T} <: Element{T} ) where T = new(coords, vec, nothing) end -##coordnames(_::Point) = [:xₚ, :yₚ, :zₚ] - function buildvec!(pt::Point) coordring = parent(pt.coords[1]) pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] @@ -43,8 +61,6 @@ mutable struct Sphere{T} <: Element{T} ) where T = new(coords, vec, rel) end -##coordnames(_::Sphere) = [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] - function buildvec!(sph::Sphere) coordring = parent(sph.coords[1]) sph.vec = sph.coords @@ -130,10 +146,6 @@ function realize(ctx::Construction{T}) where T end end - display(collect(elemenum)) - display(coordnamelist) - println() - # construct coordinate ring coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) @@ -150,16 +162,14 @@ function realize(ctx::Construction{T}) where T # construct coordinate vectors for (_, elem) in elemenum buildvec!(elem) - display(elem.coords) - display(elem.vec) - println() end # turn relations into equations - vcat( + eqns = vcat( equation.(ctx.relations), [elem.rel for elem in ctx.elements if !isnothing(elem.rel)] ) + Generic.Ideal(coordring, eqns) end end @@ -172,22 +182,23 @@ 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])) -eqns_a_s = Engine.realize(ctx) +ideal_a_s = Engine.realize(ctx) 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) -eqns_ab_s = Engine.realize(ctx) +ideal_ab_s = Engine.realize(ctx) spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] tangencies = [ Engine.AlignsWithBy{CoeffType}( spheres[n], spheres[mod1(n+1, length(spheres))], - -1//1 + CoeffType(-1//1) ) for n in 1:3 ] -ctx_chain = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) \ No newline at end of file +ctx_chain = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) +ideal_chain = Engine.realize(ctx_chain) \ No newline at end of file From 4e02ee16fc32519d0d08773c58e59c44afa79ec8 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 30 Jan 2024 02:45:14 -0500 Subject: [PATCH 09/32] Find dimension of solution variety --- engine-proto/engine.jl | 31 +++++++++++++++++-------------- engine-proto/hitting-set.jl | 15 ++++++++------- 2 files changed, 25 insertions(+), 21 deletions(-) diff --git a/engine-proto/engine.jl b/engine-proto/engine.jl index b282fdd..6d3636d 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/engine.jl @@ -1,3 +1,5 @@ +include("hitting-set.jl") + module Engine export Construction, mprod @@ -6,27 +8,25 @@ import Subscripts using LinearAlgebra using AbstractAlgebra using Groebner +using ..HittingSet # --- commutative algebra --- -# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate -# polynomial rings when coefficients are integers. in `reduce_gens`, the -# `lmnode` constructor requires < to be defined on the coefficients, and the -# `reducer_size` heuristic requires `ndigits` to be defined on the coefficients. -# this patch for `reducer_size` removes the `ndigits` dependency -##function Generic.reducer_size(f::T) where {U <: MPolyRingElem{<:FieldElement}, V, N, T <: Generic.lmnode{U, V, N}} -## if f.size != 0.0 -## return f.size -## end -## return 0.0 + sum(j^2 for j in 1:length(f.poly)) -##end - # as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate # polynomial rings when the coefficients are integers. we use Groebner to extend # support to rationals and to finite fields of prime order Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} = Generic.Ideal{U}(base_ring(I), groebner(gens(I))) +function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} + leading = [exponent_vector(f, 1) for f in gens(I)] + targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading] + length(HittingSet.solve(HittingSetProblem(targets), maxdepth)) +end + +dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} = + length(gens(base_ring(I))) - codimension(I, maxdepth) + # --- primitve elements --- abstract type Element{T} end @@ -183,6 +183,7 @@ 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 freeom") b = Engine.Point{CoeffType}() b_on_s = Engine.LiesOn{CoeffType}(b, s) @@ -190,6 +191,7 @@ Engine.push!(ctx, b) Engine.push!(ctx, s) Engine.push!(ctx, b_on_s) ideal_ab_s = Engine.realize(ctx) +println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freeom") spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] tangencies = [ @@ -200,5 +202,6 @@ tangencies = [ ) for n in 1:3 ] -ctx_chain = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) -ideal_chain = Engine.realize(ctx_chain) \ No newline at end of file +ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) +ideal_tan_sph = Engine.realize(ctx_tan_sph) +println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freeom") \ No newline at end of file diff --git a/engine-proto/hitting-set.jl b/engine-proto/hitting-set.jl index e9aacf4..347c4d2 100644 --- a/engine-proto/hitting-set.jl +++ b/engine-proto/hitting-set.jl @@ -1,13 +1,15 @@ module HittingSet +export HittingSetProblem, solve + HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}} -# `subsets` should be a collection of Set objects -function HittingSetProblem(subsets, chosen = Set()) - wholeset = union(subsets...) +# `targets` should be a collection of Set objects +function HittingSetProblem(targets, chosen = Set()) + wholeset = union(targets...) T = eltype(wholeset) unsorted_moves = [ - elt => Set(filter(s -> elt ∉ s, subsets)) + elt => Set(filter(s -> elt ∉ s, targets)) for elt in wholeset ] moves = sort(unsorted_moves, by = pair -> length(pair.second)) @@ -32,7 +34,6 @@ end function solve(pblm::HittingSetProblem{T}, maxdepth = Inf) where T problems = Dict(pblm) - println(typeof(problems)) while length(first(problems).first) < maxdepth subproblems = typeof(problems)() for (chosen, moves) in problems @@ -56,7 +57,7 @@ end function test(n = 1) T = [Int64, Int64, Symbol, Symbol][n] - subsets = Set{T}.([ + targets = Set{T}.([ [ [1, 3, 5], [2, 3, 4], @@ -98,7 +99,7 @@ function test(n = 1) [:b, :z, :t14] ] ][n]) - problem = HittingSetProblem(subsets) + problem = HittingSetProblem(targets) if isa(problem, HittingSetProblem{T}) println("Correct type") else From 65d23fb6676c73aad5817941c94d75573b6e4d10 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 30 Jan 2024 02:49:33 -0500 Subject: [PATCH 10/32] Use module names as filenames You're right: this naming convention seems to be standard for Julia modules now. --- engine-proto/{engine.jl => Engine.jl} | 2 +- engine-proto/{hitting-set.jl => HittingSet.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename engine-proto/{engine.jl => Engine.jl} (99%) rename engine-proto/{hitting-set.jl => HittingSet.jl} (100%) diff --git a/engine-proto/engine.jl b/engine-proto/Engine.jl similarity index 99% rename from engine-proto/engine.jl rename to engine-proto/Engine.jl index 6d3636d..a632581 100644 --- a/engine-proto/engine.jl +++ b/engine-proto/Engine.jl @@ -1,4 +1,4 @@ -include("hitting-set.jl") +include("HittingSet.jl") module Engine diff --git a/engine-proto/hitting-set.jl b/engine-proto/HittingSet.jl similarity index 100% rename from engine-proto/hitting-set.jl rename to engine-proto/HittingSet.jl From a3f3f6a31bde0f28c01cb35c9ad18f7719232d90 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 1 Feb 2024 16:13:22 -0500 Subject: [PATCH 11/32] Order spheres before points within each coordinate block MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In the cases I've tried so far, this leads to substantially smaller Gröbner bases. --- engine-proto/Engine.jl | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index a632581..1977e99 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -118,28 +118,39 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - # --- constructions --- mutable struct Construction{T} - elements::Set{Element{T}} + points::Set{Point{T}} + spheres::Set{Sphere{T}} relations::Set{Relation{T}} function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T allelements = union(elements, (rel.elements for rel in relations)...) - new{T}(allelements, relations) + new{T}( + filter(elt -> isa(elt, Point), allelements), + filter(elt -> isa(elt, Sphere), allelements), + relations + ) end end -function Base.push!(ctx::Construction{T}, elem::Element{T}) where T - push!(ctx.elements, elem) +function Base.push!(ctx::Construction{T}, elem::Point{T}) where T + push!(ctx.points, elem) +end + +function Base.push!(ctx::Construction{T}, elem::Sphere{T}) where T + push!(ctx.spheres, elem) end function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T push!(ctx.relations, rel) - union!(ctx.elements, rel.elements) + for elt in rel.elements + push!(ctx, elt) + end end function realize(ctx::Construction{T}) where T # collect coordinate names coordnamelist = Symbol[] - elemenum = enumerate(ctx.elements) + elemenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points))) for coordindex in 1:5 for indexed_elem in elemenum pushcoordname!(coordnamelist, indexed_elem, coordindex) @@ -167,7 +178,7 @@ function realize(ctx::Construction{T}) where T # turn relations into equations eqns = vcat( equation.(ctx.relations), - [elem.rel for elem in ctx.elements if !isnothing(elem.rel)] + [elem.rel for (_, elem) in elemenum if !isnothing(elem.rel)] ) Generic.Ideal(coordring, eqns) end From 21f09c4a4dfbea36ed28088030480bd7b6be22cb Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sun, 4 Feb 2024 16:08:13 -0500 Subject: [PATCH 12/32] Switch element abbreviation from "elem" to "elt" --- engine-proto/Engine.jl | 46 +++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 1977e99..72b5923 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -72,21 +72,21 @@ const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}( nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] ) -coordname(elem::Element, index) = coordnames[nameof(typeof(elem))][index] +coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index] -function pushcoordname!(coordnamelist, indexed_elem::Tuple{Any, Element}, coordindex) - elemindex, elem = indexed_elem - name = coordname(elem, coordindex) +function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex) + eltindex, elt = indexed_elt + name = coordname(elt, coordindex) if !isnothing(name) - subscript = Subscripts.sub(string(elemindex)) + subscript = Subscripts.sub(string(eltindex)) push!(coordnamelist, Symbol(name, subscript)) end end -function takecoord!(coordlist, indexed_elem::Tuple{Any, Element}, coordindex) - elem = indexed_elem[2] - if !isnothing(coordname(elem, coordindex)) - push!(elem.coords, popfirst!(coordlist)) +function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex) + elt = indexed_elt[2] + if !isnothing(coordname(elt, coordindex)) + push!(elt.coords, popfirst!(coordlist)) end end @@ -132,12 +132,12 @@ mutable struct Construction{T} end end -function Base.push!(ctx::Construction{T}, elem::Point{T}) where T - push!(ctx.points, elem) +function Base.push!(ctx::Construction{T}, elt::Point{T}) where T + push!(ctx.points, elt) end -function Base.push!(ctx::Construction{T}, elem::Sphere{T}) where T - push!(ctx.spheres, elem) +function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T + push!(ctx.spheres, elt) end function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T @@ -150,10 +150,10 @@ end function realize(ctx::Construction{T}) where T # collect coordinate names coordnamelist = Symbol[] - elemenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points))) + eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points))) for coordindex in 1:5 - for indexed_elem in elemenum - pushcoordname!(coordnamelist, indexed_elem, coordindex) + for indexed_elt in eltenum + pushcoordname!(coordnamelist, indexed_elt, coordindex) end end @@ -161,24 +161,24 @@ function realize(ctx::Construction{T}) where T coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) # retrieve coordinates - for (_, elem) in elemenum - empty!(elem.coords) + for (_, elt) in eltenum + empty!(elt.coords) end for coordindex in 1:5 - for indexed_elem in elemenum - takecoord!(coordqueue, indexed_elem, coordindex) + for indexed_elt in eltenum + takecoord!(coordqueue, indexed_elt, coordindex) end end # construct coordinate vectors - for (_, elem) in elemenum - buildvec!(elem) + for (_, elt) in eltenum + buildvec!(elt) end # turn relations into equations eqns = vcat( equation.(ctx.relations), - [elem.rel for (_, elem) in elemenum if !isnothing(elem.rel)] + [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)], ) Generic.Ideal(coordring, eqns) end From 43cbf8a3a0c3e2152306876e8481fbcab797f7f8 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 5 Feb 2024 00:10:13 -0500 Subject: [PATCH 13/32] Add relations to center and orient the construction --- engine-proto/Engine.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 72b5923..ac9ed35 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -178,8 +178,17 @@ function realize(ctx::Construction{T}) where T # turn relations into equations eqns = vcat( equation.(ctx.relations), - [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)], + [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)] ) + + # add relations to center and orient the construction + if !isempty(ctx.points) + append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) + end + if !isempty(ctx.spheres) + append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) + end + Generic.Ideal(coordring, eqns) end From 45aaaafc8f44a33cf6fb74c1cbc3a09771785015 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 8 Feb 2024 01:53:55 -0500 Subject: [PATCH 14/32] Seek sample solutions by cutting with a hyperplane The example hyperplane yields a single solution, with multiplicity six. You can find it analytically by hand, and homotopy continuation finds it numerically. --- engine-proto/Engine.jl | 99 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 82 insertions(+), 17 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index ac9ed35..b5eee96 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -2,12 +2,13 @@ include("HittingSet.jl") module Engine -export Construction, mprod +export Construction, mprod, codimension, dimension import Subscripts using LinearAlgebra using AbstractAlgebra using Groebner +using HomotopyContinuation: Variable, Expression, System using ..HittingSet # --- commutative algebra --- @@ -27,6 +28,34 @@ end dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} = length(gens(base_ring(I))) - codimension(I, maxdepth) +# 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 + +## [to do] not needed right now +# create a ModelKit.System from a list of elements of a multivariate polynomial +# ring. the variable ordering is taken from the polynomial ring +##function System(eqns::AbstractVector{MPolyRingElem}) +## if isempty(eqns) +## return System([]) +## else +## variables = Variable.(symbols(parent(f))) +## return System(Expression.(eqns), variables = variables) +## end +##end + # --- primitve elements --- abstract type Element{T} end @@ -189,39 +218,75 @@ function realize(ctx::Construction{T}) where T append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) end - Generic.Ideal(coordring, eqns) + (Generic.Ideal(coordring, eqns), eqns) end end # ~~~ sandbox setup ~~~ +using AbstractAlgebra +using HomotopyContinuation + 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 freeom") +##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 = Engine.realize(ctx) -println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of freeom") +ideal_ab_s, eqns_ab_s = Engine.realize(ctx) +println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " 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//1) - ) - for n in 1:3 +##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] +##tangencies = [ +## Engine.AlignsWithBy{CoeffType}( +## spheres[n], +## spheres[mod1(n+1, length(spheres))], +## CoeffType(-1//1) +## ) +## for n in 1:3 +##] +##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) +##ideal_tan_sph = Engine.realize(ctx_tan_sph) +##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom") + +# --- test rational cut --- + +cut = [ + sum(vcat(a.coords, (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])) ] -ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) -ideal_tan_sph = Engine.realize(ctx_tan_sph) -println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freeom") \ No newline at end of file +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 + + # 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, 1]) + witness = witness_set(System(eqns_ab_s, variables = vbls), cut_subspace) + println("witness solutions:") + for wtns in solutions(witness) + display(wtns) + end +end \ No newline at end of file From f97090c9974cb7d3b33b06396f5aeb84dcf620de Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 8 Feb 2024 01:58:12 -0500 Subject: [PATCH 15/32] Try a cut that goes through the trivial solution The previous cut was supposed to do this, but I was missing some parentheses. --- engine-proto/Engine.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index b5eee96..41d3ed7 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -262,8 +262,8 @@ println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of f cut = [ sum(vcat(a.coords, (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])) + 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) @@ -283,7 +283,7 @@ if cut_dim == 0 # 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, 1]) + 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) From 95c0ff14b249f4033864859fde61336f2ced5962 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 9 Feb 2024 17:09:43 -0500 Subject: [PATCH 16/32] Show explicitly that all coefficients are 1 in first cut equation --- engine-proto/Engine.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 41d3ed7..38ed672 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -261,7 +261,7 @@ println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of f # --- test rational cut --- cut = [ - sum(vcat(a.coords, (s.coords - [0, 0, 0, 0, 1]))) + 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]))) ] From 34358a872800810975e23f8499efd201982d4641 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 9 Feb 2024 23:44:10 -0500 Subject: [PATCH 17/32] Find witnesses on random rational hyperplanes Choose hyperplanes that go through the trivial solution. --- engine-proto/Engine.jl | 83 ++++++++++++++++++++++++++++++------------ 1 file changed, 60 insertions(+), 23 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 38ed672..546bf21 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -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,67 @@ println("Two points on a sphere: ", Engine.dimension(ideal_ab_s), " degrees of f # --- test rational cut --- +cut_coeffs = [ + 1 1 1 0 0 0 1 1 1 1 1; + 2 1 1 0 0 0 1 2 1 1 1; + 1 2 0 0 0 0 1 1 0 1 2 +] 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]))) + sum(vcat(cf[1:3] .* a.coords, cf[4:6] .* b.coords, cf[7:end] .* (s.coords - [0, 0, 0, 0, 1]))) + for cf in eachrow(cut_coeffs) ] 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_freedom = Engine.dimension(cut_ideal_ab_s) +println("Two points on a sphere, after cut: ", cut_freedom, " degrees of freedom") +if cut_freedom == 0 + coordring = base_ring(ideal_ab_s) + vbls = Variable.(symbols(coordring)) 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 + ##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) + # test a random witness set + 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(eqns_ab_s, variables = vbls), cut_subspace) + append!(samples, solution.(filter(isreal, results(wtns)))) + end println("witness solutions:") - for wtns in solutions(witness) - display(wtns) + 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 end \ No newline at end of file From becefe0c47253f0d7017f8e045a8dbadb761a52f Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 00:59:50 -0500 Subject: [PATCH 18/32] Try switching to compiled system --- engine-proto/Engine.jl | 107 ++++++++++++++++++----------------------- 1 file changed, 47 insertions(+), 60 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 546bf21..4bce0d7 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -263,67 +263,54 @@ println("Two points on a sphere: ", freedom, " degrees of freedom") # --- test rational cut --- -cut_coeffs = [ - 1 1 1 0 0 0 1 1 1 1 1; - 2 1 1 0 0 0 1 2 1 1 1; - 1 2 0 0 0 0 1 1 0 1 2 -] -cut = [ - sum(vcat(cf[1:3] .* a.coords, cf[4:6] .* b.coords, cf[7:end] .* (s.coords - [0, 0, 0, 0, 1]))) - for cf in eachrow(cut_coeffs) -] -cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut]) -cut_freedom = Engine.dimension(cut_ideal_ab_s) -println("Two points on a sphere, after cut: ", cut_freedom, " degrees of freedom") -if cut_freedom == 0 - coordring = base_ring(ideal_ab_s) - vbls = Variable.(symbols(coordring)) - 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 - - # test a random witness set - 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(eqns_ab_s, variables = vbls), cut_subspace) - append!(samples, solution.(filter(isreal, results(wtns)))) - end - println("witness 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 +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 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("witness 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 = 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))") + 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 \ No newline at end of file From 06872a04afb2b211550b471012d2ade5b92bd0d6 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 01:06:06 -0500 Subject: [PATCH 19/32] Say how many sample solutions we found --- engine-proto/Engine.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 4bce0d7..9841223 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -298,7 +298,7 @@ for _ in 1:3 wtns = witness_set(system, cut_subspace) append!(samples, solution.(filter(isreal, results(wtns)))) end -println("witness solutions:") +println("$(length(samples)) sample solutions:") for soln in samples display([vbls round.(soln, digits = 6)]) k_sq = abs2(soln[1]) From 8e33987f596ae9572b53885f5b4642105d45165c Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 13:46:01 -0500 Subject: [PATCH 20/32] Systematically try out different cut planes --- engine-proto/Engine.jl | 80 +++++++++++++++++++++++++++--------------- 1 file changed, 51 insertions(+), 29 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 9841223..ba4300c 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -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 - 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))}") +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 + ##] + ## [verbose] display(cut_matrix) + if through_trivial + cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)] + ## [verbose] display(cut_offset) + cut_subspace = LinearSubspace(cut_matrix, cut_offset) 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))") + cut_subspace = LinearSubspace(cut_matrix, fill(0, 3)) end + wtns = witness_set(system, cut_subspace) + 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 - 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("--- 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))") + 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 end \ No newline at end of file From af1d31f6e61552b0320c7127ef3c720f75e40dea Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 14:21:52 -0500 Subject: [PATCH 21/32] Test a scale constraint MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In all but a few cases (for example, a single point on a plane), we should be able to us the radius-coradius boost symmetry to make the average co-radius—representing the "overall scale"—roughly one. --- engine-proto/Engine.jl | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index ba4300c..0b6162e 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -210,13 +210,17 @@ function realize(ctx::Construction{T}) where T [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)] ) - # add relations to center and orient the construction + # add relations to center, orient, and scale the construction if !isempty(ctx.points) append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) end if !isempty(ctx.spheres) append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) end + n_elts = length(ctx.points) + length(ctx.spheres) + if n_elts > 0 + push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) + end (Generic.Ideal(coordring, eqns), eqns) end @@ -305,10 +309,14 @@ for through_trivial in [false, true] ## [verbose] display(cut_offset) cut_subspace = LinearSubspace(cut_matrix, cut_offset) else - cut_subspace = LinearSubspace(cut_matrix, fill(0, 3)) + cut_subspace = LinearSubspace(cut_matrix, fill(0, freedom)) end wtns = witness_set(system, cut_subspace) - for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns)))) + real_solns = solution.(filter(isreal, results(wtns))) + nontrivial_solns = filter(is_nontrivial, real_solns) + println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") + for soln in nontrivial_solns + ##[test] 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 From b3b7c2026daa828bdc119cdb37e8ca7d0a948a61 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 14:50:50 -0500 Subject: [PATCH 22/32] Separate the algebraic and numerical parts of the engine --- engine-proto/Engine.Algebraic.jl | 201 +++++++++++++++++++++++++++ engine-proto/Engine.Numerical.jl | 25 ++++ engine-proto/Engine.jl | 229 +------------------------------ 3 files changed, 233 insertions(+), 222 deletions(-) create mode 100644 engine-proto/Engine.Algebraic.jl create mode 100644 engine-proto/Engine.Numerical.jl diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl new file mode 100644 index 0000000..b9b790a --- /dev/null +++ b/engine-proto/Engine.Algebraic.jl @@ -0,0 +1,201 @@ +module Algebraic + +export + codimension, dimension, + Construction, realize, + Element, Point, Sphere, + Relation, LiesOn, AlignsWithBy, mprod + +import Subscripts +using LinearAlgebra +using AbstractAlgebra +using Groebner +using ...HittingSet + +# --- commutative algebra --- + +# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate +# polynomial rings when the coefficients are integers. we use Groebner to extend +# support to rationals and to finite fields of prime order +Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} = + Generic.Ideal{U}(base_ring(I), groebner(gens(I))) + +function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} + leading = [exponent_vector(f, 1) for f in gens(I)] + targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading] + length(HittingSet.solve(HittingSetProblem(targets), maxdepth)) +end + +dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} = + length(gens(base_ring(I))) - codimension(I, maxdepth) + +# --- primitve elements --- + +abstract type Element{T} end + +mutable struct Point{T} <: Element{T} + coords::Vector{MPolyRingElem{T}} + vec::Union{Vector{MPolyRingElem{T}}, Nothing} + rel::Nothing + + ## [to do] constructor argument never needed? + Point{T}( + coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], + vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing + ) where T = new(coords, vec, nothing) +end + +function buildvec!(pt::Point) + coordring = parent(pt.coords[1]) + pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] +end + +mutable struct Sphere{T} <: Element{T} + coords::Vector{MPolyRingElem{T}} + vec::Union{Vector{MPolyRingElem{T}}, Nothing} + rel::Union{MPolyRingElem{T}, Nothing} + + ## [to do] constructor argument never needed? + Sphere{T}( + coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], + vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, + rel::Union{MPolyRingElem{T}, Nothing} = nothing + ) where T = new(coords, vec, rel) +end + +function buildvec!(sph::Sphere) + coordring = parent(sph.coords[1]) + sph.vec = sph.coords + sph.rel = mprod(sph.coords, sph.coords) + one(coordring) +end + +const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}( + nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ], + nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] +) + +coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index] + +function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex) + eltindex, elt = indexed_elt + name = coordname(elt, coordindex) + if !isnothing(name) + subscript = Subscripts.sub(string(eltindex)) + push!(coordnamelist, Symbol(name, subscript)) + end +end + +function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex) + elt = indexed_elt[2] + if !isnothing(coordname(elt, coordindex)) + push!(elt.coords, popfirst!(coordlist)) + end +end + +# --- primitive relations --- + +abstract type Relation{T} end + +mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end]) + +# elements: point, sphere +struct LiesOn{T} <: Relation{T} + elements::Vector{Element{T}} + + LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph]) +end + +equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec) + +# elements: sphere, sphere +struct AlignsWithBy{T} <: Relation{T} + elements::Vector{Element{T}} + cos_angle::T + + AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle) +end + +equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle + +# --- constructions --- + +mutable struct Construction{T} + points::Set{Point{T}} + spheres::Set{Sphere{T}} + relations::Set{Relation{T}} + + function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T + allelements = union(elements, (rel.elements for rel in relations)...) + new{T}( + filter(elt -> isa(elt, Point), allelements), + filter(elt -> isa(elt, Sphere), allelements), + relations + ) + end +end + +function Base.push!(ctx::Construction{T}, elt::Point{T}) where T + push!(ctx.points, elt) +end + +function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T + push!(ctx.spheres, elt) +end + +function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T + push!(ctx.relations, rel) + for elt in rel.elements + push!(ctx, elt) + end +end + +function realize(ctx::Construction{T}) where T + # collect coordinate names + coordnamelist = Symbol[] + eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points))) + for coordindex in 1:5 + for indexed_elt in eltenum + pushcoordname!(coordnamelist, indexed_elt, coordindex) + end + end + + # construct coordinate ring + coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) + + # retrieve coordinates + for (_, elt) in eltenum + empty!(elt.coords) + end + for coordindex in 1:5 + for indexed_elt in eltenum + takecoord!(coordqueue, indexed_elt, coordindex) + end + end + + # construct coordinate vectors + for (_, elt) in eltenum + buildvec!(elt) + end + + # turn relations into equations + eqns = vcat( + equation.(ctx.relations), + [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)] + ) + + # add relations to center, orient, and scale the construction + if !isempty(ctx.points) + append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) + end + if !isempty(ctx.spheres) + append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) + end + n_elts = length(ctx.points) + length(ctx.spheres) + if n_elts > 0 + push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) + end + + (Generic.Ideal(coordring, eqns), eqns) +end + +end \ No newline at end of file diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl new file mode 100644 index 0000000..669bbda --- /dev/null +++ b/engine-proto/Engine.Numerical.jl @@ -0,0 +1,25 @@ +module Numerical + +using AbstractAlgebra +using HomotopyContinuation: Variable, Expression, System +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 + +end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 0b6162e..b1b0b30 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -2,229 +2,14 @@ include("HittingSet.jl") module Engine +include("Engine.Algebraic.jl") +include("Engine.Numerical.jl") + +using .Algebraic +using .Numerical + export Construction, mprod, codimension, dimension -import Subscripts -using LinearAlgebra -using AbstractAlgebra -using Groebner -using HomotopyContinuation: Variable, Expression, System -using ..HittingSet - -# --- commutative algebra --- - -# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate -# polynomial rings when the coefficients are integers. we use Groebner to extend -# support to rationals and to finite fields of prime order -Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} = - Generic.Ideal{U}(base_ring(I), groebner(gens(I))) - -function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} - leading = [exponent_vector(f, 1) for f in gens(I)] - targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading] - length(HittingSet.solve(HittingSetProblem(targets), maxdepth)) -end - -dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} = - length(gens(base_ring(I))) - codimension(I, maxdepth) - -# 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 - -## [to do] not needed right now -# create a ModelKit.System from a list of elements of a multivariate polynomial -# ring. the variable ordering is taken from the polynomial ring -##function System(eqns::AbstractVector{MPolyRingElem}) -## if isempty(eqns) -## return System([]) -## else -## variables = Variable.(symbols(parent(f))) -## return System(Expression.(eqns), variables = variables) -## end -##end - -# --- primitve elements --- - -abstract type Element{T} end - -mutable struct Point{T} <: Element{T} - coords::Vector{MPolyRingElem{T}} - vec::Union{Vector{MPolyRingElem{T}}, Nothing} - rel::Nothing - - ## [to do] constructor argument never needed? - Point{T}( - coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], - vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing - ) where T = new(coords, vec, nothing) -end - -function buildvec!(pt::Point) - coordring = parent(pt.coords[1]) - pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...] -end - -mutable struct Sphere{T} <: Element{T} - coords::Vector{MPolyRingElem{T}} - vec::Union{Vector{MPolyRingElem{T}}, Nothing} - rel::Union{MPolyRingElem{T}, Nothing} - - ## [to do] constructor argument never needed? - Sphere{T}( - coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[], - vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing, - rel::Union{MPolyRingElem{T}, Nothing} = nothing - ) where T = new(coords, vec, rel) -end - -function buildvec!(sph::Sphere) - coordring = parent(sph.coords[1]) - sph.vec = sph.coords - sph.rel = mprod(sph.coords, sph.coords) + one(coordring) -end - -const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}( - nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ], - nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ] -) - -coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index] - -function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex) - eltindex, elt = indexed_elt - name = coordname(elt, coordindex) - if !isnothing(name) - subscript = Subscripts.sub(string(eltindex)) - push!(coordnamelist, Symbol(name, subscript)) - end -end - -function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex) - elt = indexed_elt[2] - if !isnothing(coordname(elt, coordindex)) - push!(elt.coords, popfirst!(coordlist)) - end -end - -# --- primitive relations --- - -abstract type Relation{T} end - -mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end]) - -# elements: point, sphere -struct LiesOn{T} <: Relation{T} - elements::Vector{Element{T}} - - LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph]) -end - -equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec) - -# elements: sphere, sphere -struct AlignsWithBy{T} <: Relation{T} - elements::Vector{Element{T}} - cos_angle::T - - AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle) -end - -equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle - -# --- constructions --- - -mutable struct Construction{T} - points::Set{Point{T}} - spheres::Set{Sphere{T}} - relations::Set{Relation{T}} - - function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T - allelements = union(elements, (rel.elements for rel in relations)...) - new{T}( - filter(elt -> isa(elt, Point), allelements), - filter(elt -> isa(elt, Sphere), allelements), - relations - ) - end -end - -function Base.push!(ctx::Construction{T}, elt::Point{T}) where T - push!(ctx.points, elt) -end - -function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T - push!(ctx.spheres, elt) -end - -function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T - push!(ctx.relations, rel) - for elt in rel.elements - push!(ctx, elt) - end -end - -function realize(ctx::Construction{T}) where T - # collect coordinate names - coordnamelist = Symbol[] - eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points))) - for coordindex in 1:5 - for indexed_elt in eltenum - pushcoordname!(coordnamelist, indexed_elt, coordindex) - end - end - - # construct coordinate ring - coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex) - - # retrieve coordinates - for (_, elt) in eltenum - empty!(elt.coords) - end - for coordindex in 1:5 - for indexed_elt in eltenum - takecoord!(coordqueue, indexed_elt, coordindex) - end - end - - # construct coordinate vectors - for (_, elt) in eltenum - buildvec!(elt) - end - - # turn relations into equations - eqns = vcat( - equation.(ctx.relations), - [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)] - ) - - # add relations to center, orient, and scale the construction - if !isempty(ctx.points) - append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) - end - if !isempty(ctx.spheres) - append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) - end - n_elts = length(ctx.points) + length(ctx.spheres) - if n_elts > 0 - push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) - end - - (Generic.Ideal(coordring, eqns), eqns) -end - end # ~~~ sandbox setup ~~~ @@ -293,7 +78,7 @@ is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(ge max_slope = 5 binom = Binomial(2max_slope, 1/2) Random.seed!(6071) -n_planes = 36 +n_planes = 3 for through_trivial in [false, true] samples = [] for _ in 1:n_planes From 621c4c577630b3dc7e351bdd0e521c7e2a31d8ba Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 15:02:26 -0500 Subject: [PATCH 23/32] Try uniformly distributed hyperplane orientations Unit normals are uniformly distributed over the sphere. --- engine-proto/Engine.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index b1b0b30..e6b326d 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -75,26 +75,29 @@ 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) +##max_slope = 5 +##binom = Binomial(2max_slope, 1/2) Random.seed!(6071) n_planes = 3 for through_trivial in [false, true] samples = [] for _ in 1:n_planes - cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope + cut_matrix = transpose(hcat( + (normalize(randn(length(gens(coordring)))) for _ in 1:freedom)... + )) + ##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 ##] - ## [verbose] display(cut_matrix) + display(cut_matrix) ## [verbose] if through_trivial cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)] - ## [verbose] display(cut_offset) + display(cut_offset) ## [verbose] cut_subspace = LinearSubspace(cut_matrix, cut_offset) else - cut_subspace = LinearSubspace(cut_matrix, fill(0, freedom)) + cut_subspace = LinearSubspace(cut_matrix, fill(0., freedom)) end wtns = witness_set(system, cut_subspace) real_solns = solution.(filter(isreal, results(wtns))) From 6f18d4efcc560290366e4083fb7c81a349fce5a8 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 15:10:48 -0500 Subject: [PATCH 24/32] Test lots of uniformly distributed hyperplanes --- engine-proto/Engine.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index e6b326d..d8d4b52 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -78,7 +78,7 @@ is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(ge ##max_slope = 5 ##binom = Binomial(2max_slope, 1/2) Random.seed!(6071) -n_planes = 3 +n_planes = 36 for through_trivial in [false, true] samples = [] for _ in 1:n_planes @@ -91,10 +91,10 @@ for through_trivial in [false, true] ## 1 2 1 2 0 1 1 0 1 1 0; ## 1 1 0 1 0 1 2 0 2 0 0 ##] - display(cut_matrix) ## [verbose] + ## display(cut_matrix) ## [verbose] if through_trivial cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)] - display(cut_offset) ## [verbose] + ## display(cut_offset) ## [verbose] cut_subspace = LinearSubspace(cut_matrix, cut_offset) else cut_subspace = LinearSubspace(cut_matrix, fill(0., freedom)) @@ -104,7 +104,7 @@ for through_trivial in [false, true] nontrivial_solns = filter(is_nontrivial, real_solns) println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") for soln in nontrivial_solns - ##[test] for soln in filter(is_nontrivial, solution.(filter(isreal, results(wtns)))) + ## [test] 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 @@ -117,7 +117,7 @@ for through_trivial in [false, true] end println("$(length(samples)) sample solutions, not including the trivial one:") for soln in samples - ## [verbose] display([vbls round.(soln, digits = 6)]) + ## display([vbls round.(soln, digits = 6)]) ## [verbose] k_sq = abs2(soln[1]) if abs2(soln[end-2]) > 1e-12 if k_sq < 1e-12 From 1f173708eb57fc2305b9a17ab88bdf1375cc26fd Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Sat, 10 Feb 2024 17:39:26 -0500 Subject: [PATCH 25/32] Move random cut routine into engine --- engine-proto/Engine.Numerical.jl | 19 +++++++++++++++++-- engine-proto/Engine.jl | 20 +------------------- 2 files changed, 18 insertions(+), 21 deletions(-) diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl index 669bbda..1ae7b97 100644 --- a/engine-proto/Engine.Numerical.jl +++ b/engine-proto/Engine.Numerical.jl @@ -1,7 +1,8 @@ module Numerical +using LinearAlgebra using AbstractAlgebra -using HomotopyContinuation: Variable, Expression, System +using HomotopyContinuation using ..Algebraic # --- polynomial conversion --- @@ -10,7 +11,7 @@ using ..Algebraic # 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)) + f_data = zip(AbstractAlgebra.coefficients(f), exponent_vectors(f)) sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data) end @@ -22,4 +23,18 @@ function System(I::Generic.Ideal) System(eqns, variables = variables) end +# --- sampling --- + +function real_samples(F::AbstractSystem, dim) + # 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(nvariables(F))) for _ in 1:dim)... + )) + cut = LinearSubspace(normals, fill(0., dim)) + filter(isreal, results(witness_set(F, cut))) +end + end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index d8d4b52..f058085 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -82,25 +82,7 @@ n_planes = 36 for through_trivial in [false, true] samples = [] for _ in 1:n_planes - cut_matrix = transpose(hcat( - (normalize(randn(length(gens(coordring)))) for _ in 1:freedom)... - )) - ##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 - ##] - ## display(cut_matrix) ## [verbose] - if through_trivial - cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)] - ## display(cut_offset) ## [verbose] - cut_subspace = LinearSubspace(cut_matrix, cut_offset) - else - cut_subspace = LinearSubspace(cut_matrix, fill(0., freedom)) - end - wtns = witness_set(system, cut_subspace) - real_solns = solution.(filter(isreal, results(wtns))) + real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) nontrivial_solns = filter(is_nontrivial, real_solns) println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") for soln in nontrivial_solns From 6cf07dc6a10766026cf30365c3cd431d87f17967 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 12 Feb 2024 20:34:12 -0500 Subject: [PATCH 26/32] Evaluate and display elements --- engine-proto/Engine.Numerical.jl | 16 +++++- engine-proto/Engine.jl | 90 ++++++++++++++++---------------- 2 files changed, 58 insertions(+), 48 deletions(-) diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl index 1ae7b97..48fb682 100644 --- a/engine-proto/Engine.Numerical.jl +++ b/engine-proto/Engine.Numerical.jl @@ -2,7 +2,10 @@ module Numerical using LinearAlgebra using AbstractAlgebra -using HomotopyContinuation +using HomotopyContinuation: + Variable, Expression, AbstractSystem, System, LinearSubspace, + nvariables, isreal, witness_set, results +import GLMakie using ..Algebraic # --- polynomial conversion --- @@ -11,7 +14,7 @@ using ..Algebraic # 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(AbstractAlgebra.coefficients(f), exponent_vectors(f)) + f_data = zip(coefficients(f), exponent_vectors(f)) sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data) end @@ -37,4 +40,13 @@ function real_samples(F::AbstractSystem, dim) filter(isreal, results(witness_set(F, cut))) 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 \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index f058085..099a5f0 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -19,6 +19,7 @@ using Distributions using LinearAlgebra using AbstractAlgebra using HomotopyContinuation +using GLMakie CoeffType = Rational{Int64} @@ -55,62 +56,59 @@ println("Two points on a sphere: ", freedom, " degrees of freedom") 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 a random witness set system = CompiledSystem(System(eqns_ab_s, variables = vbls)) 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") +## [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)) -is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring)) -##max_slope = 5 -##binom = Binomial(2max_slope, 1/2) +## [old] is_nontrivial = soln -> norm2(abs.(real.(soln)) - trivial_soln) > 1e-4*length(gens(coordring)) Random.seed!(6071) -n_planes = 36 -for through_trivial in [false, true] - samples = [] - for _ in 1:n_planes - real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) - nontrivial_solns = filter(is_nontrivial, real_solns) - println("$(length(real_solns) - length(nontrivial_solns)) trivial solutions found") - for soln in nontrivial_solns - ## [test] 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 +n_planes = 3 +samples = [] +for _ in 1:n_planes + 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 + if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) + push!(samples, soln) 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 - ## 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 +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 = 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))") + 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(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 \ No newline at end of file From a450f701fbfa1818c7795d654e6f2a65d2b75b15 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 12 Feb 2024 21:14:07 -0500 Subject: [PATCH 27/32] Try displaying a chain of spheres For three mutually tangent spheres, I couldn't find real solutions. --- engine-proto/Engine.jl | 56 ++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 29 deletions(-) diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 099a5f0..113c6a3 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -23,23 +23,23 @@ 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])) +##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") +##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") +##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] +spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] ##tangencies = [ ## Engine.AlignsWithBy{CoeffType}( ## spheres[n], @@ -48,31 +48,29 @@ println("Two points on a sphere: ", freedom, " degrees of freedom") ## ) ## for n in 1:3 ##] -##ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) -##ideal_tan_sph = Engine.realize(ctx_tan_sph) -##println("Three mutually tangent spheres: ", Engine.dimension(ideal_tan_sph), " degrees of freedom") +tangencies = [ + Engine.AlignsWithBy{CoeffType}(spheres[1], spheres[2], CoeffType(-1)), + 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 --- -coordring = base_ring(ideal_ab_s) +coordring = base_ring(ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_ab_s, 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") +system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) 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) -n_planes = 3 +n_planes = 16 samples = [] for _ in 1:n_planes 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 if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples) push!(samples, soln) @@ -97,7 +95,7 @@ for soln in samples end # show a sample solution -function show_solution(vals) +function show_solution(ctx, vals) # evaluate elements real_vals = real.(vals) disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points] From 31d5e7e864704895dbf0273c46185476e53e40f1 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Mon, 12 Feb 2024 22:48:16 -0500 Subject: [PATCH 28/32] Play with two points on two spheres Guess conditions that make the scaling constraint impossible to satisfy. --- engine-proto/Engine.Algebraic.jl | 2 ++ engine-proto/Engine.jl | 43 +++++++++++++++++++++++--------- 2 files changed, 33 insertions(+), 12 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index b9b790a..ca39967 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -184,6 +184,8 @@ function realize(ctx::Construction{T}) where T ) # add relations to center, orient, and scale the construction + # [to do] the scaling constraint, as written, can be impossible to satisfy + # when all of the spheres have to go through the origin if !isempty(ctx.points) append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) end diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 113c6a3..7ddf72d 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -39,7 +39,7 @@ CoeffType = Rational{Int64} ##freedom = Engine.dimension(ideal_ab_s) ##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 = [ ## Engine.AlignsWithBy{CoeffType}( ## spheres[n], @@ -48,26 +48,45 @@ spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] ## ) ## for n in 1:3 ##] -tangencies = [ - Engine.AlignsWithBy{CoeffType}(spheres[1], spheres[2], CoeffType(-1)), - 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) +##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 = 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") + +p = Engine.Point{CoeffType}() +q = Engine.Point{CoeffType}() +a = Engine.Sphere{CoeffType}() +b = Engine.Sphere{CoeffType}() +p_on_a = Engine.LiesOn{CoeffType}(p, a) +p_on_b = Engine.LiesOn{CoeffType}(p, b) +q_on_a = Engine.LiesOn{CoeffType}(q, a) +q_on_b = Engine.LiesOn{CoeffType}(q, b) +ctx_joined = Engine.Construction{CoeffType}( + elements = Set([p, q, a, b]), + relations= Set([p_on_a, p_on_b, q_on_a, q_on_b]) +) +ideal_joined, eqns_joined = Engine.realize(ctx_joined) +freedom = Engine.dimension(ideal_joined) +println("Two points on two spheres: $freedom degrees of freedom") # --- test rational cut --- -coordring = base_ring(ideal_tan_sph) +coordring = base_ring(ideal_joined) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) +system = CompiledSystem(System(eqns_joined, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) Random.seed!(6071) -n_planes = 16 +n_planes = 3 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) From e41bcc7e13d29217d51ed5e7810b2c8f4e904fed Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Tue, 13 Feb 2024 04:02:14 -0500 Subject: [PATCH 29/32] Explore the performance wall Three points on two spheres is too much. --- engine-proto/Engine.Algebraic.jl | 3 ++- engine-proto/Engine.jl | 19 ++++++++----------- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index ca39967..380cee1 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -197,7 +197,8 @@ function realize(ctx::Construction{T}) where T push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) end - (Generic.Ideal(coordring, eqns), eqns) + ## [test] (Generic.Ideal(coordring, eqns), eqns) + (nothing, eqns) end end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 7ddf72d..49011c6 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -61,21 +61,18 @@ CoeffType = Rational{Int64} ##freedom = Engine.dimension(ideal_tan_sph) ##println("Three mutually tangent spheres: $freedom degrees of freedom") -p = Engine.Point{CoeffType}() -q = Engine.Point{CoeffType}() -a = Engine.Sphere{CoeffType}() -b = Engine.Sphere{CoeffType}() -p_on_a = Engine.LiesOn{CoeffType}(p, a) -p_on_b = Engine.LiesOn{CoeffType}(p, b) -q_on_a = Engine.LiesOn{CoeffType}(q, a) -q_on_b = Engine.LiesOn{CoeffType}(q, b) +points = [Engine.Point{CoeffType}() for _ in 1:3] +spheres = [Engine.Sphere{CoeffType}() for _ in 1:2] ctx_joined = Engine.Construction{CoeffType}( - elements = Set([p, q, a, b]), - relations= Set([p_on_a, p_on_b, q_on_a, q_on_b]) + 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("Two points on two spheres: $freedom degrees of freedom") +println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom") # --- test rational cut --- From 291d5c8ff685b642fd3535b6c46179c9f134c469 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 13:28:01 -0800 Subject: [PATCH 30/32] Study mutually tangent spheres with two fixed --- engine-proto/Engine.Algebraic.jl | 24 ++++----- engine-proto/Engine.jl | 93 +++++++++++++++++--------------- 2 files changed, 62 insertions(+), 55 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index 380cee1..b6fc7a7 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -186,19 +186,19 @@ function realize(ctx::Construction{T}) where T # add relations to center, orient, and scale the construction # [to do] the scaling constraint, as written, can be impossible to satisfy # when all of the spheres have to go through the origin - if !isempty(ctx.points) - append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) - end - if !isempty(ctx.spheres) - append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) - end - n_elts = length(ctx.points) + length(ctx.spheres) - if n_elts > 0 - push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) - end + ##if !isempty(ctx.points) + ## append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3]) + ##end + ##if !isempty(ctx.spheres) + ## append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4]) + ##end + ##n_elts = length(ctx.points) + length(ctx.spheres) + ##if n_elts > 0 + ## push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts) + ##end - ## [test] (Generic.Ideal(coordring, eqns), eqns) - (nothing, eqns) + (Generic.Ideal(coordring, eqns), eqns) + ## [test] (nothing, eqns) end end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 49011c6..0410f6d 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -39,15 +39,15 @@ CoeffType = Rational{Int64} ##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//1) -## ) -## for n in 1:3 -##] +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]), @@ -56,34 +56,41 @@ CoeffType = Rational{Int64} ##Engine.LiesOn{CoeffType}(points[3], spheres[1]), ##Engine.LiesOn{CoeffType}(points[3], spheres[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") +ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(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(small_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") +##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_joined) +coordring = base_ring(small_ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_joined, variables = vbls)) +system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) Random.seed!(6071) -n_planes = 3 +n_planes = 36 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) @@ -94,21 +101,21 @@ for _ in 1:n_planes 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 +##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) From a6da6f9925590cf9382fab0e6d20be470f638abe Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 14:17:03 -0800 Subject: [PATCH 31/32] Investigate why witness sets aren't reproducible All the random number generators seems to be seeded, so why aren't the results reproducible? --- engine-proto/Engine.Numerical.jl | 15 +++++++++------ engine-proto/Engine.jl | 24 ++++++++++++------------ 2 files changed, 21 insertions(+), 18 deletions(-) diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl index 48fb682..c54e71c 100644 --- a/engine-proto/Engine.Numerical.jl +++ b/engine-proto/Engine.Numerical.jl @@ -1,5 +1,6 @@ module Numerical +using Random: default_rng using LinearAlgebra using AbstractAlgebra using HomotopyContinuation: @@ -28,16 +29,18 @@ end # --- sampling --- -function real_samples(F::AbstractSystem, dim) +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(nvariables(F))) for _ in 1:dim)... - )) - cut = LinearSubspace(normals, fill(0., dim)) - filter(isreal, results(witness_set(F, cut))) + ##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}) = diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 0410f6d..be201cf 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -59,13 +59,13 @@ tangencies = [ ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(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(small_ideal_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] @@ -83,17 +83,17 @@ println("Three mutually tangent spheres, with two fixed: $freedom degrees of fre # --- test rational cut --- -coordring = base_ring(small_ideal_tan_sph) +coordring = base_ring(ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls)) +system = CompiledSystem(System(eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) -Random.seed!(6071) -n_planes = 36 +rng = MersenneTwister(6701) +n_planes = 6 samples = [] for _ in 1:n_planes - real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) + 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) From 2a505c1f5994c96bb7ed36912467ca4ef0dae4e6 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 14:27:41 -0800 Subject: [PATCH 32/32] Store elements in arrays to keep order stable This seems to restore reproducibility. --- engine-proto/Engine.Algebraic.jl | 8 ++++---- engine-proto/Engine.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index b6fc7a7..f14f58c 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -120,11 +120,11 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - # --- constructions --- mutable struct Construction{T} - points::Set{Point{T}} - spheres::Set{Sphere{T}} - relations::Set{Relation{T}} + points::Vector{Point{T}} + spheres::Vector{Sphere{T}} + relations::Vector{Relation{T}} - function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T + function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{Relation{T}}()) where T allelements = union(elements, (rel.elements for rel in relations)...) new{T}( filter(elt -> isa(elt, Point), allelements), diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index be201cf..9f7d150 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -56,7 +56,7 @@ tangencies = [ ##Engine.LiesOn{CoeffType}(points[3], spheres[1]), ##Engine.LiesOn{CoeffType}(points[3], spheres[2]) ##] -ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies)) +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 = [