Vectornaut/dyna3
1
0
Fork 0
forked from StudioInfinity/dyna3
Code Pull requests Activity

Stow algebraic engine prototype

Browse source 
We're using the Gram matrix engine for the next stage of development,
so the algebraic engine shouldn't be at the top level anymore.
This commit is contained in:
Aaron Fenyes 2024-07-28 20:50:04 -07:00
parent 9d69a900e2
commit d7dbee4c05
5 changed files with 0 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

223
engine-proto/alg-test/ConstructionViewer.jl Normal file
View file

@ -0,0 +1,223 @@
module Viewer
using Blink
using Colors
using Printf
using Main.Engine
export ConstructionViewer, display!, opentools!, closetools!
# === Blink utilities ===
append_to_head!(w, type, content) = @js w begin
@var element = document.createElement($type)
element.appendChild(document.createTextNode($content))
document.head.appendChild(element)
end
style!(w, stylesheet) = append_to_head!(w, "style", stylesheet)
script!(w, code) = append_to_head!(w, "script", code)
# === construction viewer ===
mutable struct ConstructionViewer
win::Window
function ConstructionViewer()
# create window and open developer console
win = Window(Blink.Dict(:width => 620, :height => 830))
# set stylesheet
style!(win, """
body {
background-color: #ccc;
}
/* the maximum dimensions keep Ganja from blowing up the canvas */
#view {
display: block;
width: 600px;
height: 600px;
margin-top: 10px;
margin-left: 10px;
border-radius: 10px;
background-color: #f0f0f0;
}
#control-panel {
width: 600px;
height: 200px;
box-sizing: border-box;
padding: 5px 10px 5px 10px;
margin-top: 10px;
margin-left: 10px;
overflow-y: scroll;
border-radius: 10px;
background-color: #f0f0f0;
}
#control-panel > div {
margin-top: 5px;
padding: 4px;
border-radius: 5px;
border: solid;
font-family: monospace;
}
""")
# load Ganja.js. for an automatically updated web-hosted version, load from
#
# https://unpkg.com/ganja.js
#
# instead
loadjs!(win, "http://localhost:8000/ganja-1.0.204.js")
# create global functions and variables
script!(win, """
// create algebra
var CGA3 = Algebra(4, 1);
// initialize element list and palette
var elements = [];
var palette = [];
// declare handles for the view and its options
var view;
var viewOpt;
// declare handles for the controls
var controlPanel;
var visToggles;
// create scene function
function scene() {
commands = [];
for (let n = 0; n < elements.length; ++n) {
if (visToggles[n].checked) {
commands.push(palette[n], elements[n]);
}
}
return commands;
}
function updateView() {
requestAnimationFrame(view.update.bind(view, scene));
}
""")
@js win begin
# create view
viewOpt = Dict(
:conformal => true,
:gl => true,
:devicePixelRatio => window.devicePixelRatio
)
view = CGA3.graph(scene, viewOpt)
view.setAttribute(:id, "view")
view.removeAttribute(:style)
document.body.replaceChildren(view)
# create control panel
controlPanel = document.createElement(:div)
controlPanel.setAttribute(:id, "control-panel")
document.body.appendChild(controlPanel)
end
new(win)
end
end
mprod(v, w) =
v[1]*w[1] + v[2]*w[2] + v[3]*w[3] + v[4]*w[4] - v[5]*w[5]
function display!(viewer::ConstructionViewer, elements::Matrix)
# load elements
elements_full = []
for elt in eachcol(Engine.unmix * elements)
if mprod(elt, elt) < 0.5
elt_full = [0; elt; fill(0, 26)]
else
# `elt` is a spacelike vector, representing a generalized sphere, so we
# take its Hodge dual before passing it to Ganja.js. the dual represents
# the same generalized sphere, but Ganja.js only displays planes when
# they're represented by vectors in grade 4 rather than grade 1
elt_full = [fill(0, 26); -elt[5]; -elt[4]; elt[3]; -elt[2]; elt[1]; 0]
end
push!(elements_full, elt_full)
end
@js viewer.win elements = $elements_full.map((elt) -> @new CGA3(elt))
# generate palette. this is Gadfly's `default_discrete_colors` palette,
# available under the MIT license
palette = distinguishable_colors(
length(elements_full),
[LCHab(70, 60, 240)],
transform = c -> deuteranopic(c, 0.5),
lchoices = Float64[65, 70, 75, 80],
cchoices = Float64[0, 50, 60, 70],
hchoices = range(0, stop=330, length=24)
)
palette_packed = [RGB24(c).color for c in palette]
@js viewer.win palette = $palette_packed
# create visibility toggles
@js viewer.win begin
controlPanel.replaceChildren()
visToggles = []
end
for (elt, c) in zip(eachcol(elements), palette)
vec_str = join(map(t -> @sprintf("%.3f", t), elt), ", ")
color_str = "#$(hex(c))"
style_str = "background-color: $color_str; border-color: $color_str;"
@js viewer.win begin
@var toggle = document.createElement(:div)
toggle.setAttribute(:style, $style_str)
toggle.checked = true
toggle.addEventListener(
"click",
() -> begin
toggle.checked = !toggle.checked
toggle.style.backgroundColor = toggle.checked ? $color_str : "inherit";
updateView()
end
)
toggle.appendChild(document.createTextNode($vec_str))
visToggles.push(toggle);
controlPanel.appendChild(toggle);
end
end
# update view
@js viewer.win updateView()
end
function opentools!(viewer::ConstructionViewer)
size(viewer.win, 1240, 830)
opentools(viewer.win)
end
function closetools!(viewer::ConstructionViewer)
closetools(viewer.win)
size(viewer.win, 620, 830)
end
end
# ~~~ sandbox setup ~~~
elements = let
a = sqrt(BigFloat(3)/2)
sqrt(0.5) * BigFloat[
1 1 -1 -1 0
1 -1 1 -1 0
1 -1 -1 1 0
0.5 0.5 0.5 0.5 1+a
0.5 0.5 0.5 0.5 1-a
]
end
# show construction
viewer = Viewer.ConstructionViewer()
Viewer.display!(viewer, elements)

203
engine-proto/alg-test/Engine.Algebraic.jl Normal file
View file

@ -0,0 +1,203 @@
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::Vector{Point{T}}
spheres::Vector{Sphere{T}}
relations::Vector{Relation{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),
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
# [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
(Generic.Ideal(coordring, eqns), eqns)
end
end

53
engine-proto/alg-test/Engine.Numerical.jl Normal file
View file

@ -0,0 +1,53 @@
module Numerical
using Random: default_rng
using LinearAlgebra
using AbstractAlgebra
using HomotopyContinuation:
Variable, Expression, AbstractSystem, System, LinearSubspace,
nvariables, isreal, witness_set, results
import GLMakie
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
# --- sampling ---
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(rng, nvariables(F))) for _ in 1:dim)...
))
cut = LinearSubspace(normals, fill(0., dim))
filter(isreal, results(witness_set(F, cut, seed = 0x1974abba)))
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

76
engine-proto/alg-test/Engine.jl Normal file
View file

@ -0,0 +1,76 @@
include("HittingSet.jl")
module Engine
include("Engine.Algebraic.jl")
include("Engine.Numerical.jl")
using .Algebraic
using .Numerical
export Construction, mprod, codimension, dimension
end
# ~~~ sandbox setup ~~~
using Random
using Distributions
using LinearAlgebra
using AbstractAlgebra
using HomotopyContinuation
using GLMakie
CoeffType = Rational{Int64}
spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
tangencies = [
Engine.AlignsWithBy{CoeffType}(
spheres[n],
spheres[mod1(n+1, length(spheres))],
CoeffType(1)
)
for n in 1:3
]
ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = 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")
# --- test rational cut ---
coordring = base_ring(ideal_tan_sph)
vbls = Variable.(symbols(coordring))
# test a random witness set
system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
norm2 = vec -> real(dot(conj.(vec), vec))
rng = MersenneTwister(6071)
n_planes = 6
samples = []
for _ in 1:n_planes
real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng))
for soln in real_solns
if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
push!(samples, soln)
end
end
end
println("Found $(length(samples)) sample solutions")
# show a sample solution
function show_solution(ctx, vals)
# evaluate elements
real_vals = real.(vals)
disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
# create scene
scene = Scene()
cam3d!(scene)
scatter!(scene, disp_points, color = :green)
for sph in disp_spheres
mesh!(scene, sph, color = :gray)
end
scene
end

111
engine-proto/alg-test/HittingSet.jl Normal file
View file

@ -0,0 +1,111 @@
module HittingSet
export HittingSetProblem, solve
HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}}
# `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, targets))
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)
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]
targets = 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(targets)
if isa(problem, HittingSetProblem{T})
println("Correct type")
else
println("Wrong type: ", typeof(problem))
end
problem
end
end