Engine prototype (#13)
This PR adds code for a Julia-language prototype of a configuration solver, in the `engine-proto` folder. It uses Julia version 1.10.0. ### Approaches Development of this PR tried two broad approaches to the constraint geometry problem. Each one suggested various solution techniques. The Gram matrix approach, with the low-rank factorization technique, seems the most promising. - **Algebraic** *(In the `alg-test` subfolder).* Write the constraints as polynomials in the inversive coordinates of the elements, and use computational algebraic geometry techniques to solve the resulting system. We tried the following techniques. - **Gröbner bases** *(`Engine.Algebraic.jl`).* Symbolic. Find a Gröbner basis for the ideal generated by the constraint equations. Information about the solution variety, like its codimension, is then relatively easy to extract. - **Homotopy continuation** *(`Engine.Numerical.jl`).* Numerical. Cut the solution set along a random hyperplane to get a generic zero-dimensional slice, and then use a fancy homotopy technique to approximate the points in that slice. A few notes about our experiences can be found on the [engine prototype](wiki/Engine-prototype) wiki page. - **Gram matrix** *(in the `gram-test` subfolder).* A construction is described completely, up to conformal transformations, by the Gram matrix of the vectors representing its elements. Express the constraints as fixed entries of the Gram matrix, and use numerical linear algebra techniques to find a list of vectors whose Gram matrix fits the bill. We tried the following techniques. - **LDL decomposition** *(`gram-test.sage`, `gram-test.jl`, `overlap-test.jl`).* Find a cluster of up to five elements whose Gram matrix is completely filled in by the constraints. Use LDL decomposition to find a list of vectors with that Gram matrix. This technique can be made algebraic, as seen in `overlap-test.jl`. - **Low-rank factorization** *(source files listed in findings section).* Write down a quadratic loss function that says how far a set of vectors is from meeting the Gram matrix constraints. Use a smooth optimization technique like Newton's method or gradient descent to find a zero of the loss function. In addition to the polished prototype described in the results section, we have an early prototype using an off-the-shelf factorization package (`low-rank-test.jl`) and an visualization of the loss function landscape near global minima (`basin-shapes.jl`). The [Gram matrix parameterization](wiki/Gram-matrix-parameterization) wiki page contains detailed notes on this approach. ### Findings With the algebraic approach, we hit a performance wall pretty quickly as our constructions grew. It was often hard to find real solutions of the polynomial system, since the techniques we use work most naturally in the complex world. With the Gram matrix approach, on the other hand, we could solve interesting problems in acceptably short times using the low-rank factorization technique. We put the optimization routine in its own module (`Engine.jl`) and used it to solve five example problems: - `overlapping-pyramids.jl` - `circles-in-triangle.jl` - `sphere-in-tetrahedron.jl` - `tetrahedron-radius-ratio.jl` - `irisawa-hexlet.jl` We plan to use low-rank factorization of the Gram matrix in our first app prototype. ### Visualizations We used the visualizer in the `ganja-test` folder to visually check our low-rank factorization results. The visualizer runs [Ganja.js](https://enkimute.github.io/ganja.js/) in an Electron app, made with [Blink](https://github.com/JuliaGizmos/Blink.jl). Although Ganja.js makes beautiful pictures under most circumstances, we found two obstacles to using it in production. - It seems to have precision problems with low-curvature spheres. - We couldn't figure out how to customize its clipping and transparency settings, and the default settings often obscure construction details. Co-authored-by: Aaron Fenyes <aaron.fenyes@fareycircles.ooo> Co-authored-by: Glen Whitney <glen@studioinfinity.org> Reviewed-on: #13 Co-authored-by: Vectornaut <vectornaut@nobody@nowhere.net> Co-committed-by: Vectornaut <vectornaut@nobody@nowhere.net>
This commit is contained in:
parent
c48d685ad6
commit
b92be312e8
223
engine-proto/alg-test/ConstructionViewer.jl
Normal file
223
engine-proto/alg-test/ConstructionViewer.jl
Normal 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
203
engine-proto/alg-test/Engine.Algebraic.jl
Normal 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
53
engine-proto/alg-test/Engine.Numerical.jl
Normal 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
76
engine-proto/alg-test/Engine.jl
Normal 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
111
engine-proto/alg-test/HittingSet.jl
Normal 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
|
||||
96
engine-proto/ganja-test/ganja-test.html
Normal file
96
engine-proto/ganja-test/ganja-test.html
Normal file
@ -0,0 +1,96 @@
|
||||
<!DOCTYPE html>
|
||||
<html>
|
||||
<head>
|
||||
<style>
|
||||
body {
|
||||
background-color: #ffe0f0;
|
||||
}
|
||||
|
||||
/* needed to keep Ganja canvas from blowing up */
|
||||
canvas {
|
||||
min-width: 600px;
|
||||
max-width: 600px;
|
||||
min-height: 600px;
|
||||
max-height: 600px;
|
||||
}
|
||||
</style>
|
||||
<script src="https://unpkg.com/ganja.js"></script>
|
||||
</head>
|
||||
<body>
|
||||
<p><button onclick="flip()">Flip</button></p>
|
||||
<script>
|
||||
// in the default view, e4 + e5 is the point at infinity
|
||||
let CGA3 = Algebra(4, 1);
|
||||
let elements = [
|
||||
CGA3.inline(() => Math.sqrt(0.5)*( 1e1 + 1e2 + 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*( 1e1 - 1e2 - 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*(-1e1 + 1e2 - 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*(-1e1 - 1e2 + 1e3 + 1e5))(),
|
||||
CGA3.inline(() => -Math.sqrt(3)*1e4 + Math.sqrt(2)*1e5)()
|
||||
];
|
||||
/*
|
||||
these blocks of commented-out code can be used to confirm that a spacelike
|
||||
vector and its Hodge dual represent the same generalized sphere
|
||||
*/
|
||||
/*let elements = [
|
||||
CGA3.inline(() => Math.sqrt(0.5)*!( 1e1 + 1e2 + 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*!( 1e1 - 1e2 - 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*!(-1e1 + 1e2 - 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*!(-1e1 - 1e2 + 1e3 + 1e5))(),
|
||||
CGA3.inline(() => !(-Math.sqrt(3)*1e4 + Math.sqrt(2)*1e5))()
|
||||
];*/
|
||||
/*let elements = [
|
||||
CGA3.inline(() => 1e1 + 1e5)(),
|
||||
CGA3.inline(() => 1e2 + 1e5)(),
|
||||
CGA3.inline(() => 1e3 + 1e5)(),
|
||||
CGA3.inline(() => -1e4 + 1e5)(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*(1e1 + 1e2 + 1e3 + 1e5))(),
|
||||
CGA3.inline(() => Math.sqrt(0.5)*!(1e1 + 1e2 + 1e3 - 0.01e4 + 1e5))()
|
||||
];*/
|
||||
|
||||
// set up palette
|
||||
var colorIndex;
|
||||
var palette = [0xff00b0, 0x00ffb0, 0x00b0ff, 0x8040ff, 0xc0c0c0];
|
||||
function nextColor() {
|
||||
colorIndex = (colorIndex + 1) % palette.length;
|
||||
return palette[colorIndex];
|
||||
}
|
||||
function resetColorCycle() {
|
||||
colorIndex = palette.length - 1;
|
||||
}
|
||||
resetColorCycle();
|
||||
|
||||
// create scene function
|
||||
function scene() {
|
||||
commands = [];
|
||||
resetColorCycle();
|
||||
elements.forEach((elt) => commands.push(nextColor(), elt));
|
||||
return commands;
|
||||
}
|
||||
|
||||
// initialize graph
|
||||
let graph = CGA3.graph(
|
||||
scene,
|
||||
{
|
||||
conformal: true, gl: true, grid: true
|
||||
}
|
||||
)
|
||||
document.body.appendChild(graph);
|
||||
|
||||
function flip() {
|
||||
let last = elements.length - 1;
|
||||
for (let n = 0; n < last; ++n) {
|
||||
// reflect
|
||||
elements[n] = CGA3.Mul(CGA3.Mul(elements[last], elements[n]), elements[last]);
|
||||
|
||||
// de-noise
|
||||
for (let k = 6; k < elements[n].length; ++k) {
|
||||
/*for (let k = 0; k < 26; ++k) {*/
|
||||
elements[n][k] = 0;
|
||||
}
|
||||
}
|
||||
requestAnimationFrame(graph.update.bind(graph, scene));
|
||||
}
|
||||
</script>
|
||||
</body>
|
||||
</html>
|
||||
127
engine-proto/ganja-test/ganja-test.jl
Normal file
127
engine-proto/ganja-test/ganja-test.jl
Normal file
@ -0,0 +1,127 @@
|
||||
using Blink
|
||||
using Colors
|
||||
|
||||
# === 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)
|
||||
|
||||
function add_element!(vec)
|
||||
# add element
|
||||
full_vec = [0; vec; fill(0, 26)]
|
||||
n = @js win elements.push(@new CGA3($full_vec))
|
||||
|
||||
# generate palette. this is Gadfly's `default_discrete_colors` palette,
|
||||
# available under the MIT license
|
||||
palette = distinguishable_colors(
|
||||
n,
|
||||
[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 win palette = $palette_packed
|
||||
end
|
||||
|
||||
# === build page ===
|
||||
|
||||
# create window and open developer console
|
||||
win = Window()
|
||||
opentools(win)
|
||||
|
||||
# set stylesheet
|
||||
style!(win, """
|
||||
body {
|
||||
background-color: #ffe0f0;
|
||||
}
|
||||
|
||||
/* needed to keep Ganja canvas from blowing up */
|
||||
canvas {
|
||||
min-width: 600px;
|
||||
max-width: 600px;
|
||||
min-height: 600px;
|
||||
max-height: 600px;
|
||||
}
|
||||
""")
|
||||
|
||||
# load Ganja.js
|
||||
loadjs!(win, "https://unpkg.com/ganja.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 visualization handle
|
||||
var graph;
|
||||
|
||||
// create scene function
|
||||
function scene() {
|
||||
commands = [];
|
||||
for (let n = 0; n < elements.length; ++n) {
|
||||
commands.push(palette[n], elements[n]);
|
||||
}
|
||||
return commands;
|
||||
}
|
||||
|
||||
function flip() {
|
||||
let last = elements.length - 1;
|
||||
for (let n = 0; n < last; ++n) {
|
||||
// reflect
|
||||
elements[n] = CGA3.Mul(CGA3.Mul(elements[last], elements[n]), elements[last]);
|
||||
|
||||
// de-noise
|
||||
for (let k = 6; k < elements[n].length; ++k) {
|
||||
elements[n][k] = 0;
|
||||
}
|
||||
}
|
||||
requestAnimationFrame(graph.update.bind(graph, scene));
|
||||
}
|
||||
""")
|
||||
|
||||
# set up controls
|
||||
body!(win, """
|
||||
<p><button id="flip-button" onclick="flip()">Flip</button></p>
|
||||
""", async = false)
|
||||
|
||||
# === set up visualization ===
|
||||
|
||||
# list elements. in the default view, e4 + e5 is the point at infinity
|
||||
elements = sqrt(0.5) * BigFloat[
|
||||
1 1 -1 -1 0;
|
||||
1 -1 1 -1 0;
|
||||
1 -1 -1 1 0;
|
||||
0 0 0 0 -sqrt(6);
|
||||
1 1 1 1 2
|
||||
]
|
||||
|
||||
# load elements
|
||||
for vec in eachcol(elements)
|
||||
add_element!(vec)
|
||||
end
|
||||
|
||||
# initialize visualization
|
||||
@js win begin
|
||||
graph = CGA3.graph(
|
||||
scene,
|
||||
Dict(
|
||||
"conformal" => true,
|
||||
"gl" => true,
|
||||
"grid" => true
|
||||
)
|
||||
)
|
||||
document.body.appendChild(graph)
|
||||
end
|
||||
450
engine-proto/gram-test/Engine.jl
Normal file
450
engine-proto/gram-test/Engine.jl
Normal file
@ -0,0 +1,450 @@
|
||||
module Engine
|
||||
|
||||
using LinearAlgebra
|
||||
using GenericLinearAlgebra
|
||||
using SparseArrays
|
||||
using Random
|
||||
using Optim
|
||||
|
||||
export
|
||||
rand_on_shell, Q, DescentHistory,
|
||||
realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram
|
||||
|
||||
# === guessing ===
|
||||
|
||||
sconh(t, u) = 0.5*(exp(t) + u*exp(-t))
|
||||
|
||||
function rand_on_sphere(rng::AbstractRNG, ::Type{T}, n) where T
|
||||
out = randn(rng, T, n)
|
||||
tries_left = 2
|
||||
while dot(out, out) < 1e-6 && tries_left > 0
|
||||
out = randn(rng, T, n)
|
||||
tries_left -= 1
|
||||
end
|
||||
normalize(out)
|
||||
end
|
||||
|
||||
##[TO DO] write a test to confirm that the outputs are on the correct shells
|
||||
function rand_on_shell(rng::AbstractRNG, shell::T) where T <: Number
|
||||
space_part = rand_on_sphere(rng, T, 4)
|
||||
rapidity = randn(rng, T)
|
||||
sig = sign(shell)
|
||||
nullmix * [sconh(rapidity, sig)*space_part; sconh(rapidity, -sig)]
|
||||
end
|
||||
|
||||
rand_on_shell(rng::AbstractRNG, shells::Array{T}) where T <: Number =
|
||||
hcat([rand_on_shell(rng, sh) for sh in shells]...)
|
||||
|
||||
rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), shells)
|
||||
|
||||
# === elements ===
|
||||
|
||||
point(pos) = [pos; 0.5; 0.5 * dot(pos, pos)]
|
||||
|
||||
plane(normal, offset) = [-normal; 0; -offset]
|
||||
|
||||
function sphere(center, radius)
|
||||
dist_sq = dot(center, center)
|
||||
[
|
||||
center / radius;
|
||||
0.5 / radius;
|
||||
0.5 * (dist_sq / radius - radius)
|
||||
]
|
||||
end
|
||||
|
||||
# === Gram matrix realization ===
|
||||
|
||||
# basis changes
|
||||
nullmix = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [-1 1; 1 1]//2]
|
||||
unmix = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [-1 1; 1 1]]
|
||||
|
||||
# the Lorentz form
|
||||
Q = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [0 -2; -2 0]]
|
||||
|
||||
# project a matrix onto the subspace of matrices whose entries vanish away from
|
||||
# the given indices
|
||||
function proj_to_entries(mat, indices)
|
||||
result = zeros(size(mat))
|
||||
for (j, k) in indices
|
||||
result[j, k] = mat[j, k]
|
||||
end
|
||||
result
|
||||
end
|
||||
|
||||
# the difference between the matrices `target` and `attempt`, projected onto the
|
||||
# subspace of matrices whose entries vanish at each empty index of `target`
|
||||
function proj_diff(target::SparseMatrixCSC{T, <:Any}, attempt::Matrix{T}) where T
|
||||
J, K, values = findnz(target)
|
||||
result = zeros(size(target))
|
||||
for (j, k, val) in zip(J, K, values)
|
||||
result[j, k] = val - attempt[j, k]
|
||||
end
|
||||
result
|
||||
end
|
||||
|
||||
# a type for keeping track of gradient descent history
|
||||
struct DescentHistory{T}
|
||||
scaled_loss::Array{T}
|
||||
neg_grad::Array{Matrix{T}}
|
||||
base_step::Array{Matrix{T}}
|
||||
hess::Array{Hermitian{T, Matrix{T}}}
|
||||
slope::Array{T}
|
||||
stepsize::Array{T}
|
||||
positive::Array{Bool}
|
||||
backoff_steps::Array{Int64}
|
||||
last_line_L::Array{Matrix{T}}
|
||||
last_line_loss::Array{T}
|
||||
|
||||
function DescentHistory{T}(
|
||||
scaled_loss = Array{T}(undef, 0),
|
||||
neg_grad = Array{Matrix{T}}(undef, 0),
|
||||
hess = Array{Hermitian{T, Matrix{T}}}(undef, 0),
|
||||
base_step = Array{Matrix{T}}(undef, 0),
|
||||
slope = Array{T}(undef, 0),
|
||||
stepsize = Array{T}(undef, 0),
|
||||
positive = Bool[],
|
||||
backoff_steps = Int64[],
|
||||
last_line_L = Array{Matrix{T}}(undef, 0),
|
||||
last_line_loss = Array{T}(undef, 0)
|
||||
) where T
|
||||
new(scaled_loss, neg_grad, hess, base_step, slope, stepsize, positive, backoff_steps, last_line_L, last_line_loss)
|
||||
end
|
||||
end
|
||||
|
||||
# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
|
||||
# explicit entry of `gram`. use gradient descent starting from `guess`
|
||||
function realize_gram_gradient(
|
||||
gram::SparseMatrixCSC{T, <:Any},
|
||||
guess::Matrix{T};
|
||||
scaled_tol = 1e-30,
|
||||
min_efficiency = 0.5,
|
||||
init_stepsize = 1.0,
|
||||
backoff = 0.9,
|
||||
max_descent_steps = 600,
|
||||
max_backoff_steps = 110
|
||||
) where T <: Number
|
||||
# start history
|
||||
history = DescentHistory{T}()
|
||||
|
||||
# scale tolerance
|
||||
scale_adjustment = sqrt(T(nnz(gram)))
|
||||
tol = scale_adjustment * scaled_tol
|
||||
|
||||
# initialize variables
|
||||
stepsize = init_stepsize
|
||||
L = copy(guess)
|
||||
|
||||
# do gradient descent
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
loss = dot(Δ_proj, Δ_proj)
|
||||
for _ in 1:max_descent_steps
|
||||
# stop if the loss is tolerably low
|
||||
if loss < tol
|
||||
break
|
||||
end
|
||||
|
||||
# find negative gradient of loss function
|
||||
neg_grad = 4*Q*L*Δ_proj
|
||||
slope = norm(neg_grad)
|
||||
dir = neg_grad / slope
|
||||
|
||||
# store current position, loss, and slope
|
||||
L_last = L
|
||||
loss_last = loss
|
||||
push!(history.scaled_loss, loss / scale_adjustment)
|
||||
push!(history.neg_grad, neg_grad)
|
||||
push!(history.slope, slope)
|
||||
|
||||
# find a good step size using backtracking line search
|
||||
push!(history.stepsize, 0)
|
||||
push!(history.backoff_steps, max_backoff_steps)
|
||||
empty!(history.last_line_L)
|
||||
empty!(history.last_line_loss)
|
||||
for backoff_steps in 0:max_backoff_steps
|
||||
history.stepsize[end] = stepsize
|
||||
L = L_last + stepsize * dir
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
loss = dot(Δ_proj, Δ_proj)
|
||||
improvement = loss_last - loss
|
||||
push!(history.last_line_L, L)
|
||||
push!(history.last_line_loss, loss / scale_adjustment)
|
||||
if improvement >= min_efficiency * stepsize * slope
|
||||
history.backoff_steps[end] = backoff_steps
|
||||
break
|
||||
end
|
||||
stepsize *= backoff
|
||||
end
|
||||
|
||||
# [DEBUG] if we've hit a wall, quit
|
||||
if history.backoff_steps[end] == max_backoff_steps
|
||||
break
|
||||
end
|
||||
end
|
||||
|
||||
# return the factorization and its history
|
||||
push!(history.scaled_loss, loss / scale_adjustment)
|
||||
L, history
|
||||
end
|
||||
|
||||
function basis_matrix(::Type{T}, j, k, dims) where T
|
||||
result = zeros(T, dims)
|
||||
result[j, k] = one(T)
|
||||
result
|
||||
end
|
||||
|
||||
# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
|
||||
# explicit entry of `gram`. use Newton's method starting from `guess`
|
||||
function realize_gram_newton(
|
||||
gram::SparseMatrixCSC{T, <:Any},
|
||||
guess::Matrix{T};
|
||||
scaled_tol = 1e-30,
|
||||
rate = 1,
|
||||
max_steps = 100
|
||||
) where T <: Number
|
||||
# start history
|
||||
history = DescentHistory{T}()
|
||||
|
||||
# find the dimension of the search space
|
||||
dims = size(guess)
|
||||
element_dim, construction_dim = dims
|
||||
total_dim = element_dim * construction_dim
|
||||
|
||||
# list the constrained entries of the gram matrix
|
||||
J, K, _ = findnz(gram)
|
||||
constrained = zip(J, K)
|
||||
|
||||
# scale the tolerance
|
||||
scale_adjustment = sqrt(T(length(constrained)))
|
||||
tol = scale_adjustment * scaled_tol
|
||||
|
||||
# use Newton's method
|
||||
L = copy(guess)
|
||||
for step in 0:max_steps
|
||||
# evaluate the loss function
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
loss = dot(Δ_proj, Δ_proj)
|
||||
|
||||
# store the current loss
|
||||
push!(history.scaled_loss, loss / scale_adjustment)
|
||||
|
||||
# stop if the loss is tolerably low
|
||||
if loss < tol || step > max_steps
|
||||
break
|
||||
end
|
||||
|
||||
# find the negative gradient of loss function
|
||||
neg_grad = 4*Q*L*Δ_proj
|
||||
|
||||
# find the negative Hessian of the loss function
|
||||
hess = Matrix{T}(undef, total_dim, total_dim)
|
||||
indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
|
||||
for (j, k) in indices
|
||||
basis_mat = basis_matrix(T, j, k, dims)
|
||||
neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
|
||||
neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
|
||||
deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
|
||||
hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
|
||||
end
|
||||
hess = Hermitian(hess)
|
||||
push!(history.hess, hess)
|
||||
|
||||
# compute the Newton step
|
||||
step = hess \ reshape(neg_grad, total_dim)
|
||||
L += rate * reshape(step, dims)
|
||||
end
|
||||
|
||||
# return the factorization and its history
|
||||
L, history
|
||||
end
|
||||
|
||||
LinearAlgebra.eigen!(A::Symmetric{BigFloat, Matrix{BigFloat}}; sortby::Nothing) =
|
||||
eigen!(Hermitian(A))
|
||||
|
||||
function convertnz(type, mat)
|
||||
J, K, values = findnz(mat)
|
||||
sparse(J, K, type.(values))
|
||||
end
|
||||
|
||||
function realize_gram_optim(
|
||||
gram::SparseMatrixCSC{T, <:Any},
|
||||
guess::Matrix{T}
|
||||
) where T <: Number
|
||||
# find the dimension of the search space
|
||||
dims = size(guess)
|
||||
element_dim, construction_dim = dims
|
||||
total_dim = element_dim * construction_dim
|
||||
|
||||
# list the constrained entries of the gram matrix
|
||||
J, K, _ = findnz(gram)
|
||||
constrained = zip(J, K)
|
||||
|
||||
# scale the loss function
|
||||
scale_adjustment = length(constrained)
|
||||
|
||||
function loss(L_vec)
|
||||
L = reshape(L_vec, dims)
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
dot(Δ_proj, Δ_proj) / scale_adjustment
|
||||
end
|
||||
|
||||
function loss_grad!(storage, L_vec)
|
||||
L = reshape(L_vec, dims)
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
storage .= reshape(-4*Q*L*Δ_proj, total_dim) / scale_adjustment
|
||||
end
|
||||
|
||||
function loss_hess!(storage, L_vec)
|
||||
L = reshape(L_vec, dims)
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
|
||||
for (j, k) in indices
|
||||
basis_mat = basis_matrix(T, j, k, dims)
|
||||
neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
|
||||
neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
|
||||
deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj) / scale_adjustment
|
||||
storage[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
|
||||
end
|
||||
end
|
||||
|
||||
optimize(
|
||||
loss, loss_grad!, loss_hess!,
|
||||
reshape(guess, total_dim),
|
||||
Newton()
|
||||
)
|
||||
end
|
||||
|
||||
# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
|
||||
# explicit entry of `gram`. use gradient descent starting from `guess`
|
||||
function realize_gram(
|
||||
gram::SparseMatrixCSC{T, <:Any},
|
||||
guess::Matrix{T},
|
||||
frozen = nothing;
|
||||
scaled_tol = 1e-30,
|
||||
min_efficiency = 0.5,
|
||||
init_rate = 1.0,
|
||||
backoff = 0.9,
|
||||
reg_scale = 1.1,
|
||||
max_descent_steps = 200,
|
||||
max_backoff_steps = 110
|
||||
) where T <: Number
|
||||
# start history
|
||||
history = DescentHistory{T}()
|
||||
|
||||
# find the dimension of the search space
|
||||
dims = size(guess)
|
||||
element_dim, construction_dim = dims
|
||||
total_dim = element_dim * construction_dim
|
||||
|
||||
# list the constrained entries of the gram matrix
|
||||
J, K, _ = findnz(gram)
|
||||
constrained = zip(J, K)
|
||||
|
||||
# scale the tolerance
|
||||
scale_adjustment = sqrt(T(length(constrained)))
|
||||
tol = scale_adjustment * scaled_tol
|
||||
|
||||
# list the un-frozen indices
|
||||
has_frozen = !isnothing(frozen)
|
||||
if has_frozen
|
||||
is_unfrozen = fill(true, size(guess))
|
||||
is_unfrozen[frozen] .= false
|
||||
unfrozen = findall(is_unfrozen)
|
||||
unfrozen_stacked = reshape(is_unfrozen, total_dim)
|
||||
end
|
||||
|
||||
# initialize variables
|
||||
grad_rate = init_rate
|
||||
L = copy(guess)
|
||||
|
||||
# use Newton's method with backtracking and gradient descent backup
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
loss = dot(Δ_proj, Δ_proj)
|
||||
for step in 1:max_descent_steps
|
||||
# stop if the loss is tolerably low
|
||||
if loss < tol
|
||||
break
|
||||
end
|
||||
|
||||
# find the negative gradient of loss function
|
||||
neg_grad = 4*Q*L*Δ_proj
|
||||
|
||||
# find the negative Hessian of the loss function
|
||||
hess = Matrix{T}(undef, total_dim, total_dim)
|
||||
indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
|
||||
for (j, k) in indices
|
||||
basis_mat = basis_matrix(T, j, k, dims)
|
||||
neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
|
||||
neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
|
||||
deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
|
||||
hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
|
||||
end
|
||||
hess = Hermitian(hess)
|
||||
push!(history.hess, hess)
|
||||
|
||||
# regularize the Hessian
|
||||
min_eigval = minimum(eigvals(hess))
|
||||
push!(history.positive, min_eigval > 0)
|
||||
if min_eigval <= 0
|
||||
hess -= reg_scale * min_eigval * I
|
||||
end
|
||||
|
||||
# compute the Newton step
|
||||
neg_grad_stacked = reshape(neg_grad, total_dim)
|
||||
if has_frozen
|
||||
hess = hess[unfrozen_stacked, unfrozen_stacked]
|
||||
neg_grad_compressed = neg_grad_stacked[unfrozen_stacked]
|
||||
else
|
||||
neg_grad_compressed = neg_grad_stacked
|
||||
end
|
||||
base_step_compressed = hess \ neg_grad_compressed
|
||||
if has_frozen
|
||||
base_step_stacked = zeros(total_dim)
|
||||
base_step_stacked[unfrozen_stacked] .= base_step_compressed
|
||||
else
|
||||
base_step_stacked = base_step_compressed
|
||||
end
|
||||
base_step = reshape(base_step_stacked, dims)
|
||||
push!(history.base_step, base_step)
|
||||
|
||||
# store the current position, loss, and slope
|
||||
L_last = L
|
||||
loss_last = loss
|
||||
push!(history.scaled_loss, loss / scale_adjustment)
|
||||
push!(history.neg_grad, neg_grad)
|
||||
push!(history.slope, norm(neg_grad))
|
||||
|
||||
# find a good step size using backtracking line search
|
||||
push!(history.stepsize, 0)
|
||||
push!(history.backoff_steps, max_backoff_steps)
|
||||
empty!(history.last_line_L)
|
||||
empty!(history.last_line_loss)
|
||||
rate = one(T)
|
||||
step_success = false
|
||||
for backoff_steps in 0:max_backoff_steps
|
||||
history.stepsize[end] = rate
|
||||
L = L_last + rate * base_step
|
||||
Δ_proj = proj_diff(gram, L'*Q*L)
|
||||
loss = dot(Δ_proj, Δ_proj)
|
||||
improvement = loss_last - loss
|
||||
push!(history.last_line_L, L)
|
||||
push!(history.last_line_loss, loss / scale_adjustment)
|
||||
if improvement >= min_efficiency * rate * dot(neg_grad, base_step)
|
||||
history.backoff_steps[end] = backoff_steps
|
||||
step_success = true
|
||||
break
|
||||
end
|
||||
rate *= backoff
|
||||
end
|
||||
|
||||
# if we've hit a wall, quit
|
||||
if !step_success
|
||||
return L_last, false, history
|
||||
end
|
||||
end
|
||||
|
||||
# return the factorization and its history
|
||||
push!(history.scaled_loss, loss / scale_adjustment)
|
||||
L, loss < tol, history
|
||||
end
|
||||
|
||||
end
|
||||
99
engine-proto/gram-test/basin-shapes.jl
Normal file
99
engine-proto/gram-test/basin-shapes.jl
Normal file
@ -0,0 +1,99 @@
|
||||
include("Engine.jl")
|
||||
|
||||
using LinearAlgebra
|
||||
using SparseArrays
|
||||
|
||||
function sphere_in_tetrahedron_shape()
|
||||
# initialize the partial gram matrix for a sphere inscribed in a regular
|
||||
# tetrahedron
|
||||
J = Int64[]
|
||||
K = Int64[]
|
||||
values = BigFloat[]
|
||||
for j in 1:5
|
||||
for k in 1:5
|
||||
push!(J, j)
|
||||
push!(K, k)
|
||||
if j == k
|
||||
push!(values, 1)
|
||||
elseif (j <= 4 && k <= 4)
|
||||
push!(values, -1/BigFloat(3))
|
||||
else
|
||||
push!(values, -1)
|
||||
end
|
||||
end
|
||||
end
|
||||
gram = sparse(J, K, values)
|
||||
|
||||
# plot loss along a slice
|
||||
loss_lin = []
|
||||
loss_sq = []
|
||||
mesh = range(0.9, 1.1, 101)
|
||||
for t in mesh
|
||||
L = hcat(
|
||||
Engine.plane(normalize(BigFloat[ 1, 1, 1]), BigFloat(1)),
|
||||
Engine.plane(normalize(BigFloat[ 1, -1, -1]), BigFloat(1)),
|
||||
Engine.plane(normalize(BigFloat[-1, 1, -1]), BigFloat(1)),
|
||||
Engine.plane(normalize(BigFloat[-1, -1, 1]), BigFloat(1)),
|
||||
Engine.sphere(BigFloat[0, 0, 0], BigFloat(t))
|
||||
)
|
||||
Δ_proj = Engine.proj_diff(gram, L'*Engine.Q*L)
|
||||
push!(loss_lin, norm(Δ_proj))
|
||||
push!(loss_sq, dot(Δ_proj, Δ_proj))
|
||||
end
|
||||
mesh, loss_lin, loss_sq
|
||||
end
|
||||
|
||||
function circles_in_triangle_shape()
|
||||
# initialize the partial gram matrix for a sphere inscribed in a regular
|
||||
# tetrahedron
|
||||
J = Int64[]
|
||||
K = Int64[]
|
||||
values = BigFloat[]
|
||||
for j in 1:8
|
||||
for k in 1:8
|
||||
filled = false
|
||||
if j == k
|
||||
push!(values, 1)
|
||||
filled = true
|
||||
elseif (j == 1 || k == 1)
|
||||
push!(values, 0)
|
||||
filled = true
|
||||
elseif (j == 2 || k == 2)
|
||||
push!(values, -1)
|
||||
filled = true
|
||||
end
|
||||
#=elseif (j <= 5 && j != 2 && k == 9 || k == 9 && k <= 5 && k != 2)
|
||||
push!(values, 0)
|
||||
filled = true
|
||||
end=#
|
||||
if filled
|
||||
push!(J, j)
|
||||
push!(K, k)
|
||||
end
|
||||
end
|
||||
end
|
||||
append!(J, [6, 4, 6, 5, 7, 5, 7, 3, 8, 3, 8, 4])
|
||||
append!(K, [4, 6, 5, 6, 5, 7, 3, 7, 3, 8, 4, 8])
|
||||
append!(values, fill(-1, 12))
|
||||
|
||||
# plot loss along a slice
|
||||
loss_lin = []
|
||||
loss_sq = []
|
||||
mesh = range(0.99, 1.01, 101)
|
||||
for t in mesh
|
||||
L = hcat(
|
||||
Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
|
||||
Engine.sphere(BigFloat[0, 0, 0], BigFloat(t)),
|
||||
Engine.plane(BigFloat[1, 0, 0], BigFloat(1)),
|
||||
Engine.plane(BigFloat[cos(2pi/3), sin(2pi/3), 0], BigFloat(1)),
|
||||
Engine.plane(BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(1)),
|
||||
Engine.sphere(4//3*BigFloat[-1, 0, 0], BigFloat(1//3)),
|
||||
Engine.sphere(4//3*BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//3)),
|
||||
Engine.sphere(4//3*BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//3))
|
||||
)
|
||||
Δ_proj = Engine.proj_diff(gram, L'*Engine.Q*L)
|
||||
push!(loss_lin, norm(Δ_proj))
|
||||
push!(loss_sq, dot(Δ_proj, Δ_proj))
|
||||
end
|
||||
mesh, loss_lin, loss_sq
|
||||
end
|
||||
76
engine-proto/gram-test/circles-in-triangle.jl
Normal file
76
engine-proto/gram-test/circles-in-triangle.jl
Normal file
@ -0,0 +1,76 @@
|
||||
include("Engine.jl")
|
||||
|
||||
using SparseArrays
|
||||
using Random
|
||||
|
||||
# initialize the partial gram matrix for a sphere inscribed in a regular
|
||||
# tetrahedron
|
||||
J = Int64[]
|
||||
K = Int64[]
|
||||
values = BigFloat[]
|
||||
for j in 1:9
|
||||
for k in 1:9
|
||||
filled = false
|
||||
if j == 9
|
||||
if k <= 5 && k != 2
|
||||
push!(values, 0)
|
||||
filled = true
|
||||
end
|
||||
elseif k == 9
|
||||
if j <= 5 && j != 2
|
||||
push!(values, 0)
|
||||
filled = true
|
||||
end
|
||||
elseif j == k
|
||||
push!(values, 1)
|
||||
filled = true
|
||||
elseif j == 1 || k == 1
|
||||
push!(values, 0)
|
||||
filled = true
|
||||
elseif j == 2 || k == 2
|
||||