diff --git a/engine-proto/alg-test/ConstructionViewer.jl b/engine-proto/alg-test/ConstructionViewer.jl
new file mode 100644
index 0000000..b9c8ffb
--- /dev/null
+++ b/engine-proto/alg-test/ConstructionViewer.jl
@@ -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)
\ No newline at end of file
diff --git a/engine-proto/alg-test/Engine.Algebraic.jl b/engine-proto/alg-test/Engine.Algebraic.jl
new file mode 100644
index 0000000..a9b6667
--- /dev/null
+++ b/engine-proto/alg-test/Engine.Algebraic.jl
@@ -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
\ No newline at end of file
diff --git a/engine-proto/alg-test/Engine.Numerical.jl b/engine-proto/alg-test/Engine.Numerical.jl
new file mode 100644
index 0000000..d1e14bd
--- /dev/null
+++ b/engine-proto/alg-test/Engine.Numerical.jl
@@ -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
\ No newline at end of file
diff --git a/engine-proto/alg-test/Engine.jl b/engine-proto/alg-test/Engine.jl
new file mode 100644
index 0000000..f6f92c5
--- /dev/null
+++ b/engine-proto/alg-test/Engine.jl
@@ -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
\ No newline at end of file
diff --git a/engine-proto/alg-test/HittingSet.jl b/engine-proto/alg-test/HittingSet.jl
new file mode 100644
index 0000000..347c4d2
--- /dev/null
+++ b/engine-proto/alg-test/HittingSet.jl
@@ -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
\ No newline at end of file
diff --git a/engine-proto/ganja-test/ganja-test.html b/engine-proto/ganja-test/ganja-test.html
new file mode 100644
index 0000000..0207dcc
--- /dev/null
+++ b/engine-proto/ganja-test/ganja-test.html
@@ -0,0 +1,96 @@
+
+
+
+
+
+
+
+
+
+
+
diff --git a/engine-proto/ganja-test/ganja-test.jl b/engine-proto/ganja-test/ganja-test.jl
new file mode 100644
index 0000000..6c55061
--- /dev/null
+++ b/engine-proto/ganja-test/ganja-test.jl
@@ -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, """
+
+""", 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
\ No newline at end of file
diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl
new file mode 100644
index 0000000..ac5fe54
--- /dev/null
+++ b/engine-proto/gram-test/Engine.jl
@@ -0,0 +1,451 @@
+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
+## [old] Q = diagm([1, 1, 1, 1, -1])
+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 at 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
\ No newline at end of file
diff --git a/engine-proto/gram-test/basin-shapes.jl b/engine-proto/gram-test/basin-shapes.jl
new file mode 100644
index 0000000..5c03c01
--- /dev/null
+++ b/engine-proto/gram-test/basin-shapes.jl
@@ -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
\ No newline at end of file
diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl
new file mode 100644
index 0000000..1bd22a7
--- /dev/null
+++ b/engine-proto/gram-test/circles-in-triangle.jl
@@ -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
+ push!(values, -1)
+ 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))
+#= make construction rigid
+append!(J, [3, 4, 4, 5])
+append!(K, [4, 3, 5, 4])
+append!(values, fill(-0.5, 4))
+=#
+gram = sparse(J, K, values)
+
+# set initial guess
+Random.seed!(58271)
+guess = hcat(
+ Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
+ Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ Engine.plane(-BigFloat[1, 0, 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ Engine.plane(-BigFloat[cos(2pi/3), sin(2pi/3), 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ Engine.plane(-BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ Engine.sphere(BigFloat[-1, 0, 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ Engine.sphere(BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+ BigFloat[0, 0, 0, 0, 1]
+)
+frozen = [CartesianIndex(j, 9) for j in 1:5]
+
+# complete the gram matrix using Newton's method with backtracking
+L, success, history = Engine.realize_gram(gram, guess, frozen)
+completed_gram = L'*Engine.Q*L
+println("Completed Gram matrix:\n")
+display(completed_gram)
+if success
+ println("\nTarget accuracy achieved!")
+else
+ println("\nFailed to reach target accuracy")
+end
+println("Steps: ", size(history.scaled_loss, 1))
+println("Loss: ", history.scaled_loss[end], "\n")
\ No newline at end of file
diff --git a/engine-proto/gram-test/ganja-1.0.204.js b/engine-proto/gram-test/ganja-1.0.204.js
new file mode 100644
index 0000000..1e95e42
--- /dev/null
+++ b/engine-proto/gram-test/ganja-1.0.204.js
@@ -0,0 +1,1913 @@
+/** Ganja.js - Geometric Algebra - Not Just Algebra.
+ * @author Enki
+ * @link https://github.com/enkimute/ganja.js
+ */
+
+/*********************************************************************************************************************/
+//
+// Ganja.js is an Algebra generator for javascript. It generates a wide variety of Algebra's and supports operator
+// overloading, algebraic literals and a variety of graphing options.
+//
+// Ganja.js is designed with prototyping and educational purposes in mind. Clean mathematical syntax is the primary
+// target.
+//
+// Ganja.js exports only one function called *Algebra*. This function is used to generate Algebra classes. (say complex
+// numbers, minkowski or 3D CGA). The returned class can be used to create, add, multiply etc, but also to upgrade
+// javascript functions with algebraic literals, operator overloading, vectors, matrices and much more.
+//
+// As a simple example, multiplying two complex numbers 3+2i and 1+4i could be done like this :
+//
+// var complex = Algebra(0,1);
+// var a = new complex([3,2]);
+// var b = new complex([1,3]);
+// var result = a.Mul(b);
+//
+// But the same can be written using operator overloading and algebraic literals. (where scientific notation with
+// lowercase e is overloaded to directly specify generators (e1, e2, e12, ...))
+//
+// var result = Algebra(0,1,()=>(3+2e1)*(1+4e1));
+//
+// Please see github for user documentation and examples.
+//
+/*********************************************************************************************************************/
+
+// Documentation below is for implementors. I'll assume you know about Clifford Algebra's, grades, its products, etc ..
+// I'll also assume you are familiar with ES6. My style may feel a bith mathematical, advise is to read slow.
+
+(function (name, context, definition) {
+ if (typeof module != 'undefined' && module.exports) module.exports = definition();
+ else if (typeof define == 'function' && define.amd) define(name, definition);
+ else context[name] = definition();
+}('Algebra', this, function () {
+
+/** Some helpers for eigenvalues for bivector split in high-d spaces **/
+ function QR(M) {
+ // helpers
+ const {abs,sqrt} = Math;
+ const hyp = (a,b)=>abs(a)>abs(b)?abs(a)*sqrt(1+(b/a)**2):b==0?0:abs(b)*sqrt(1+(a/b)**2);
+ const [m,n] = [M.length, M[0].length];
+ var qr = M.map(r=>r.map(c=>c)), Q = M.map(r=>r.map(c=>0)), R = M.map(r=>r.map(c=>0)), d = [], k, i, j, nrm;
+ // helper matrix
+ for (k=0; k=0; --k) {
+ for (i=0; i{
+ var res = A.map(r=>r.map(c=>0));
+ for(let i=0;iA[i][i]);
+ }
+
+/** The Algebra class generator. Possible calling signatures :
+ * Algebra([func]) => algebra with no dimensions, i.e. R. Optional function for the translator.
+ * Algebra(p,[func]) => 'p' positive dimensions and an optional function to pass to the translator.
+ * Algebra(p,q,[func]) => 'p' positive and 'q' negative dimensions and optional function.
+ * Algebra(p,q,r,[func]) => 'p' positive, 'q' negative and 'r' zero dimensions and optional function.
+ * Algebra({ => for custom basis, cayley, mixing, etc pass in an object as first parameter.
+ * [p:p], => optional 'p' for # of positive dimensions
+ * [q:q], => optional 'q' for # of negative dimensions
+ * [r:r], => optional 'r' for # of zero dimensions
+ * [metric:array], => alternative for p,q,r. e.g. ([1,1,1,-1] for spacetime)
+ * [basis:array], => array of strings with basis names. (e.g. ['1','e1','e2','e12'])
+ * [Cayley:Cayley], => optional custom Cayley table (strings). (e.g. [['1','e1'],['e1','-1']])
+ * [mix:boolean], => Allows mixing of various algebras. (for space efficiency).
+ * [graded:boolean], => Use a graded algebra implementation. (automatic for +6D)
+ * [baseType:Float32Array] => optional basetype to use. (only for flat generator)
+ * },[func]) => optional function for the translator.
+ **/
+ return function Algebra(p,q,r) {
+ // Resolve possible calling signatures so we know the numbers for p,q,r. Last argument can always be a function.
+ var fu=arguments[arguments.length-1],options=p; if (options instanceof Object) {
+ q = (p.q || (p.metric && p.metric.filter(x=>x==-1).length))| 0;
+ r = (p.r || (p.metric && p.metric.filter(x=>x==0).length)) | 0;
+ p = p.p === undefined ? (p.metric && p.metric.filter(x=>x==1).length) || 0 : p.p || 0;
+ } else { options={}; p=p|0; r=r|0; q=q|0; };
+
+ // Support for multi-dual-algebras
+ if (options.dual || (p==0 && q==0 && r<0)) { r=options.dual=options.dual||-r; // Create a dual number algebra if r<0 (old) or options.dual set(new)
+ options.basis = [...Array(r+1)].map((a,i)=>i?'e0'+i:'1'); options.metric = [1,...Array(r)]; options.tot=r+1;
+ options.Cayley = [...Array(r+1)].map((a,i)=>[...Array(r+1)].map((y,j)=>i*j==0?((i+j)?'e0'+(i+j):'1'):'0'));
+ }
+ if (options.over) options.baseType = Array;
+
+ // Calculate the total number of dimensions.
+ var tot = options.tot = (options.tot||(p||0)+(q||0)+(r||0)||(options.basis&&options.basis.length))|0;
+
+ // Unless specified, generate a full set of Clifford basis names. We generate them as an array of strings by starting
+ // from numbers in binary representation and changing the set bits into their relative position.
+ // Basis names are ordered first per grade, then lexically (not cyclic!).
+ // For 10 or more dimensions all names will be double digits ! 1e01 instead of 1e1 ..
+ var basis=(options.basis&&(options.basis.length==2**tot||r<0||options.Cayley)&&options.basis)||[...Array(2**tot)] // => [undefined, undefined, undefined, undefined, undefined, undefined, undefined, undefined]
+ .map((x,xi)=>(((1<<30)+xi).toString(2)).slice(-tot||-1) // => ["000", "001", "010", "011", "100", "101", "110", "111"] (index of array in base 2)
+ .replace(/./g,(a,ai)=>a=='0'?'':String.fromCharCode(66+ai-(r!=0)))) // => ["", "3", "2", "23", "1", "13", "12", "123"] (1 bits replaced with their positions, 0's removed)
+ .sort((a,b)=>(a.toString().length==b.toString().length)?(a>b?1:b>a?-1:0):a.toString().length-b.toString().length) // => ["", "1", "2", "3", "12", "13", "23", "123"] (sorted numerically)
+ .map(x=>x&&'e'+(x.replace(/./g,x=>('0'+(x.charCodeAt(0)-65)).slice(tot>9?-2:-1) ))||'1') // => ["1", "e1", "e2", "e3", "e12", "e13", "e23", "e123"] (converted to commonly used basis names)
+
+ // See if the basis names start from 0 or 1, store grade per component and lowest component per grade.
+ var low=basis.length==1?1:basis[1].match(/\d+/g)[0]*1,
+ grades=options.grades||(options.dual&&basis.map((x,i)=>i?1:0))||basis.map(x=>tot>9?(x.length-1)/2:x.length-1),
+ grade_start=grades.map((a,b,c)=>c[b-1]!=a?b:-1).filter(x=>x+1).concat([basis.length]);
+
+ // String-simplify a concatenation of two basis blades. (and supports custom basis names e.g. e21 instead of e12)
+ // This is the function that implements e1e1 = +1/-1/0 and e1e2=-e2e1. The brm function creates the remap dictionary.
+ var simplify = (s,p,q,r)=>{
+ var sign=1,c,l,t=[],f=true,ss=s.match(tot>9?/(\d\d)/g:/(\d)/g);if (!ss) return s; s=ss; l=s.length;
+ while (f) { f=false;
+ // implement Ex*Ex = metric.
+ for (var i=0; i=(p+r)) sign*=-1; else if ((s[i]-low)t[i+1]) { c=t[i];t[i]=t[i+1];t[i+1]=c;sign*=-1;f=true; break;} if (f) { s=t;t=[];l=s.length; }
+ }
+ var ret=(sign==0)?'0':((sign==1)?'':'-')+(t.length?'e'+t.join(''):'1'); return (brm&&brm[ret])||(brm&&brm['-'+ret]&&'-'+brm['-'+ret])||ret;
+ },
+ brm=(x=>{ var ret={}; for (var i in basis) ret[basis[i]=='1'?'1':simplify(basis[i],p,q,r)] = basis[i]; return ret; })(basis);
+
+ // As an alternative to the string fiddling, one can also bit-fiddle. In this case the basisvectors are represented by integers with 1 bit per generator set.
+ var simplify_bits = (A,B,p2)=>{ var n=p2||(p+q+r),t=0,ab=A&B,res=A^B; if (ab&((1<>1); t&=B; t^=ab>>(p+r); t^=t>>16; t^=t>>8; t^=t>>4; return [1-2*(27030>>(t&15)&1),res]; },
+ bc = (v)=>{ v=v-((v>>1)& 0x55555555); v=(v&0x33333333)+((v>>2)&0x33333333); var c=((v+(v>>4)&0xF0F0F0F)*0x1010101)>>24; return c };
+
+ if (!options.graded && tot <= 6 || options.graded===false || options.Cayley) {
+ // Faster and degenerate-metric-resistant dualization. (a remapping table that maps items into their duals).
+ var drm=basis.map((a,i)=>{ return {a:a,i:i} })
+ .sort((a,b)=>a.a.length>b.a.length?1:a.a.lengthx.i).reverse(),
+ drms=drm.map((x,i)=>(x==0||i==0)?1:simplify(basis[x]+basis[i])[0]=='-'?-1:1);
+
+ /// Store the full metric (also for bivectors etc ..)
+ var metric = options.Cayley&&options.Cayley.map((x,i)=>x[i]) || basis.map((x,xi)=>simplify(x+x,p,q,r)|0); metric[0]=1;
+
+ /// Generate multiplication tables for the outer and geometric products.
+ var mulTable = options.Cayley||basis.map(x=>basis.map(y=>(x==1)?y:(y==1)?x:simplify(x+y,p,q,r)));
+
+ // subalgebra support. (must be bit-order basis blades, does no error checking.)
+ if (options.even) options.basis = basis.filter(x=>x.length%2==1);
+ if (options.basis && !options.Cayley && r>=0 && options.basis.length != 2**tot) {
+ metric = metric.filter((x,i)=>options.basis.indexOf(basis[i])!=-1);
+ mulTable = mulTable.filter((x,i)=>options.basis.indexOf(basis[i])!=-1).map(x=>x.filter((x,i)=>options.basis.indexOf(basis[i])!=-1));
+ basis = options.basis;
+ }
+
+ /// Convert Cayley table to product matrices. The outer product selects the strict sum of the GP (but without metric), the inner product
+ /// is the left contraction.
+ var gp=basis.map(x=>basis.map(x=>'0')), cp=gp.map(x=>gp.map(x=>'0')), cps=gp.map(x=>gp.map(x=>'0')), op=gp.map(x=>gp.map(x=>'0')), gpo={}; // Storage for our product tables.
+ basis.forEach((x,xi)=>basis.forEach((y,yi)=>{ var n = mulTable[xi][yi].replace(/^-/,''); if (!gpo[n]) gpo[n]=[]; gpo[n].push([xi,yi]); }));
+ basis.forEach((o,oi)=>{
+ gpo[o].forEach(([xi,yi])=>op[oi][xi]=(grades[oi]==grades[xi]+grades[yi])?((mulTable[xi][yi]=='0')?'0':((mulTable[xi][yi][0]!='-')?'':'-')+'b['+yi+']*this['+xi+']'):'0');
+ gpo[o].forEach(([xi,yi])=>{
+ gp[oi][xi] =((gp[oi][xi]=='0')?'':gp[oi][xi]+'+') + ((mulTable[xi][yi]=='0')?'0':((mulTable[xi][yi][0]!='-')?'':'-')+'b['+yi+']*this['+xi+']');
+ cp[oi][xi] =((cp[oi][xi]=='0')?'':cp[oi][xi]+'+') + ((grades[oi]==grades[yi]-grades[xi])?gp[oi][xi]:'0');
+ cps[oi][xi]=((cps[oi][xi]=='0')?'':cps[oi][xi]+'+') + ((grades[oi]==Math.abs(grades[yi]-grades[xi]))?gp[oi][xi]:'0');
+ });
+ });
+
+ /// Flat Algebra Multivector Base Class.
+ var generator = class MultiVector extends (options.baseType||Float32Array) {
+ /// constructor - create a floating point array with the correct number of coefficients.
+ constructor(a) { super(a||basis.length); return this; }
+
+ /// grade selection - return a only the part of the input with the specified grade.
+ Grade(grade,res) { res=res||new this.constructor(); for (var i=0,l=res.length; i1e-10) res.push(((this[i]==1)&&i?'':((this[i]==-1)&&i)?'-':(this[i].toFixed(10)*1))+(i==0?'':tot==1&&q==1?'i':basis[i].replace('e','e_'))); return res.join('+').replace(/\+-/g,'-')||'0'; }
+
+ /// Reversion, Involutions, Conjugation for any number of grades, component acces shortcuts.
+ get Negative (){ var res = new this.constructor(); for (var i=0; ia[drm[i]]*drms[i]); var res = new this.constructor(); res[res.length-1]=1; return this.Mul(res); };
+ get UnDual (){ if (r) return this.map((x,i,a)=>a[drm[i]]*drms[a.length-i-1]); var res = new this.constructor(); res[res.length-1]=1; return this.Div(res); };
+ get Length (){ return options.over?Math.sqrt(Math.abs(this.Mul(this.Conjugate).s.s)):Math.sqrt(Math.abs(this.Mul(this.Conjugate).s)); };
+ get VLength (){ var res = 0; for (var i=0; i'res['+xi+']=b['+xi+']+this['+xi+']').join(';\n').replace(/(b|this)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';\nreturn res')
+ generator.prototype.Scale = new Function('b,res','res=res||new this.constructor();\n'+basis.map((x,xi)=>'res['+xi+']=b*this['+xi+']').join(';\n').replace(/(b|this)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';\nreturn res')
+ generator.prototype.Sub = new Function('b,res','res=res||new this.constructor();\n'+basis.map((x,xi)=>'res['+xi+']=this['+xi+']-b['+xi+']').join(';\n').replace(/(b|this)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';\nreturn res')
+ generator.prototype.Mul = new Function('b,res','res=res||new this.constructor();\n'+gp.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\+\-/g,'-').replace(/(\w*?)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a).replace(/\+0/g,'')+';').join('\n')+'\nreturn res;');
+ generator.prototype.LDot = new Function('b,res','res=res||new this.constructor();\n'+cp.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\+\-/g,'-').replace(/\+0/g,'').replace(/(\w*?)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\n')+'\nreturn res;');
+ generator.prototype.Dot = new Function('b,res','res=res||new this.constructor();\n'+cps.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\+\-/g,'-').replace(/\+0/g,'').replace(/(\w*?)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\n')+'\nreturn res;');
+ generator.prototype.Wedge = new Function('b,res','res=res||new this.constructor();\n'+op.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\+\-/g,'-').replace(/\+0/g,'').replace(/(\w*?)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\n')+'\nreturn res;');
+// generator.prototype.Vee = new Function('b,res','res=res||new this.constructor();\n'+op.map((r,ri)=>'res['+drm[ri]+']='+r.map(x=>x.replace(/\[(.*?)\]/g,function(a,b){return '['+(drm[b|0])+']'})).join('+').replace(/\+\-/g,'-').replace(/\+0/g,'').replace(/(\w*?)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\n')+'\nreturn res;');
+ /// Conforms to the new Chapter 11 now.
+ generator.prototype.Vee = new Function('b,res',('res=res||new this.constructor();\n'+op.map((r,ri)=>'res['+drm[ri]+']='+drms[ri]+'*('+r.map(x=>x.replace(/\[(.*?)\]/g,function(a,b){return '['+(drm[b|0])+']'+(drms[b|0]>0?"":"*-1")})).join('+').replace(/\+\-/g,'-').replace(/\+0/g,'').replace(/(\w*?)\[(.*?)\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+');').join('\n')+'\nreturn res;').replace(/(b\[)|(this\[)/g,a=>a=='b['?'this[':'b['));
+ generator.prototype.eigenValues = eigenValues;
+
+ /// Add getter and setters for the basis vectors/bivectors etc ..
+ basis.forEach((b,i)=>Object.defineProperty(generator.prototype, i?b:'s', {
+ configurable: true, get(){ return this[i] }, set(x){ this[i]=x; }
+ }));
+
+ /// Graded generator for high-dimensional algebras.
+ } else {
+
+ /// extra graded lookups.
+ var basisg = grade_start.slice(0,grade_start.length-1).map((x,i)=>basis.slice(x,grade_start[i+1]));
+ var counts = grade_start.map((x,i,a)=>i==a.length-1?0:a[i+1]-x).slice(0,tot+1);
+ var basis_bits = basis.map(x=>x=='1'?0:x.slice(1).match(tot>9?/\d\d/g:/\d/g).reduce((a,b)=>a+(1<<(b-low)),0)),
+ bits_basis = []; basis_bits.forEach((b,i)=>bits_basis[b]=i);
+ var metric = basisg.map((x,xi)=>x.map((y,yi)=>simplify_bits(basis_bits[grade_start[xi]+yi],basis_bits[grade_start[xi]+yi])[0]));
+ var drms = basisg.map((x,xi)=>x.map((y,yi)=>simplify_bits(basis_bits[grade_start[xi]+yi],(~basis_bits[grade_start[xi]+yi])&((1<(typeof x=="string")?"-"+x:-x):undefined):this[i];
+ else { if (r[i]==undefined) r[i]=[]; for(var j=0,m=Math.max(this[i].length,b[i].length);jx&&x.map(y=>typeof y=="string"?y+"*"+s:y*s)); }
+
+ // geometric product.
+ Mul(b,r) {
+ r=r||new this.constructor(); var gotstring=false;
+ for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],ig.map(e=>e&&(!(e+'').match(/-{0,1}\w+/))?'('+e+')':e))
+ return r;
+ }
+ // outer product.
+ Wedge(b,r) {
+ r=r||new this.constructor();
+ for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],ix).sort((a,b)=>((a.match(/\d+/)[0]|0)-(b.match(/\d+/)[0]|0))||((a.match(/\d+$/)[0]|0)-(b.match(/\d+$/)[0]|0))).map(x=>x.replace(/\/\/\d+$/,''));
+ var r2 = 'float sum=0.0; float res=0.0;\n', g=0;
+ r.forEach(x=>{
+ var cg = x.match(/\d+/)[0]|0;
+ if (cg != g) r2 += "sum += res*res;\nres = 0.0;\n";
+ r2 += x.replace(/\[\d+\]/,'') + '\n';
+ g=cg;
+ });
+ r2+= "sum += res*res;\n";
+ return r2;
+ }
+ // Inner product glsl output.
+ IPNS_GLSL(b,point_source) {
+ var r='',count=0,curg;
+ for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],ix).sort((a,b)=>((a.match(/\d+/)[0]|0)-(b.match(/\d+/)[0]|0))||((a.match(/\d+$/)[0]|0)-(b.match(/\d+$/)[0]|0))).map(x=>x.replace(/\/\/\d+$/,''));
+ var r2 = 'float sum=0.0; float res=0.0;\n', g=0;
+ r.forEach(x=>{
+ var cg = x.match(/\d+/)[0]|0;
+ if (cg != g) r2 += "sum += res*res;\nres = 0.0;\n";
+ r2 += x.replace(/\[\d+\]/,'') + '\n';
+ g=cg;
+ });
+ r2+= "sum += res*res;\n";
+ return r2;
+ }
+ // Left contraction.
+ LDot(b,r) {
+ r=r||new this.constructor();
+ for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],ig&&g.map((c,ci)=>!c?undefined:((c+'').match(/[\+\-\*]/)?'('+c+')':c)+(gi==0?"":basisg[gi][ci])).filter(x=>x).join('+')).filter(x=>x).join('+').replace(/\+\-/g,'-')||"0"; }
+ get s () { if (this[0]) return this[0][0]||0; return 0; }
+ get Length () { var res=0; this.forEach((g,gi)=>g&&g.forEach((e,ei)=>res+=(e||0)**2*metric[gi][ei])); return Math.abs(res)**.5; }
+ get VLength () { var res=0; this.forEach((g,gi)=>g&&g.forEach((e,ei)=>res+=(e||0)**2)); return Math.abs(res)**.5; }
+ get Reverse () { var r=new this.constructor(); this.forEach((x,gi)=>x&&x.forEach((e,ei)=>{if(!r[gi])r[gi]=[]; r[gi][ei] = this[gi][ei]*[1,1,-1,-1][gi%4]; })); return r; }
+ get Involute () { var r=new this.constructor(); this.forEach((x,gi)=>x&&x.forEach((e,ei)=>{if(!r[gi])r[gi]=[]; r[gi][ei] = this[gi][ei]*[1,-1,1,-1][gi%4]; })); return r; }
+ get Conjugate () { var r=new this.constructor(); this.forEach((x,gi)=>x&&x.forEach((e,ei)=>{if(!r[gi])r[gi]=[]; r[gi][ei] = this[gi][ei]*[1,-1,-1,1][gi%4]; })); return r; }
+ get Dual() { var r=new this.constructor(); this.forEach((g,gi)=>{ if (!g) return; r[tot-gi]=[]; g.forEach((e,ei)=>r[tot-gi][counts[gi]-1-ei]=drms[gi][ei]*e); }); return r; }
+ get Normalized () { return this.Scale(1/this.Length); }
+ }
+
+
+ // This generator is UNDER DEVELOPMENT - I'm publishing it so I can test on observable.
+ }
+
+ // Generate a new class for our algebra. It extends the javascript typed arrays (default float32 but can be specified in options).
+ var res = class Element extends generator {
+
+ // constructor - create a floating point array with the correct number of coefficients.
+ constructor(a) { super(a); if (this.upgrade) this.upgrade(); return this; }
+
+ // Grade selection. (implemented by parent class).
+ Grade(grade,res) { res=res||new Element(); return super.Grade(grade,res); }
+
+ // Right and Left divide - Defined on the elements, shortcuts to multiplying with the inverse.
+ Div (b,res) { return this.Mul(b.Inverse,res); }
+ LDiv (b,res) { return b.Inverse.Mul(this,res); }
+
+
+ // Bivector split - we handle all real cases, still have to add the complex cases for those exception scenarios.
+ Split (iter=50) {
+ var k = Math.floor((p+q+r)/2), OB = this.map(x=>x), B = this.map(x=>x), m = 1;
+ var Wi = [...Array(k)].map((r,i)=>{ m = m*(i+1); var Wi = B.Scale(1/m); B = B.Wedge(OB); return Wi; });
+ if (k<3) { // The quadratic case is easy to solve. (for spaces <6D)
+ var TDT = this.Dot(this).s, TWT = this.Wedge(this);
+ if (TWT.VLength < 1E-5) return [this]; // bivector was simple.
+ var D = 0.5*Math.sqrt( TDT**2 - TWT.Mul(TWT).s );
+ var eigen = [0.5*TDT + D, 0.5*TDT - D].sort((a,b)=>Math.abs(a)6D, closed form solutions of the characteristic polyn. are impossible, use eigenvalues of companion matrix.
+ var Wis = Wi.map((W,i)=>W.Mul(W).s*(-1)**(k-i+(k%2)) );
+ var matrix = [...Array(k)].map((r,i)=>[...Array(k)].map((c,j)=>(j == k-1)?Wis[k-i-1]:(i-1==j)?1:0));
+ var eigen = eigenValues(matrix,iter).sort((a,b)=>Math.abs(a){
+ var r = Math.floor(k/2), N = Element.Scalar(0), DN = Element.Scalar(0);
+ for (var i=0; i<=r; ++i) { N.Add( Wi[2*i+1].Scale(v**(r-i)), N); DN.Add( Wi[2*i].Scale(v**(r-i)), DN); }
+ if (DN.VLength == 0) return Element.Scalar(0);
+ var ret = N.Div(DN); sum.Add(ret, sum); return ret;
+ });
+ return [this.Sub(sum),...res]; // Smallest eigvalue becomes B-rest
+ }
+
+ // Factorize a motor
+ Factorize (iter=50) {
+ var S = this.Grade(2).Split(iter);
+ var P = this.Scale(1);
+ // if (P.s) {
+ var R = S.slice(0,S.length-1).map((Si,i)=>{
+ var Mi = Element.Scalar(P.s).Add(Si);
+ var scale = Math.sqrt(Mi.Reverse.Mul(Mi).s);
+ return Mi.Scale(1/scale);
+ });
+ R.push( R.reduce((tot,fact)=>tot.Mul(fact.Reverse), Element.Scalar(1)).Mul(P) );
+ // }
+ return R;
+ }
+
+ // exp - closed form exp.
+ Exp (taylor = false) {
+ if (r==1 && tot<=4 && Math.abs(this[0])<1E-9 && !options.over && !taylor) {
+ // https://www.researchgate.net/publication/360528787_Normalization_Square_Roots_and_the_Exponential_and_Logarithmic_Maps_in_Geometric_Algebras_of_Less_than_6D
+ // 0 1 2 3 4 5
+ // 5 6 7 8 9 10
+ var l = (this[8]*this[8] + this[9]*this[9] + this[10]*this[10]);
+ if (l==0) return new Element([1, 0,0,0,0, this[5], this[6], this[7], 0, 0, 0, 0,0,0,0, 0]);
+ var m = (this[5]*this[10] + this[6]*this[9] + this[7]*this[8]), a = Math.sqrt(l), c = Math.cos(a), s = Math.sin(a)/a, t = m/l*(c-s);
+ var test = Element.Element(c, 0,0,0,0, s*this[5] + t*this[10], s*this[6] + t*this[9], s*this[7] + t*this[8], s*this[8], s*this[9], s*this[10], 0,0,0,0, m*s);
+ //return test; // tbc .. investigate pss coeff??
+
+ var u = Math.sqrt(Math.abs(this.Dot(this).s)); if (Math.abs(u)<1E-5) return this.Add(Element.Scalar(1));
+ var v = this.Wedge(this).Scale(-1/(2*u));
+ var res2 = Element.Add(Element.Sub(Math.cos(u),v.Scale(Math.sin(u))),Element.Div(Element.Mul((Element.Add(Math.sin(u),v.Scale(Math.cos(u)))),this),(Element.Add(u,v))));
+ //if ([...test].map(x=>x.toFixed(1))+'' != [...res2].map(x=>x.toFixed(1))+'') { console.log(test, res2); debugger }
+
+ return res2;
+ }
+ if (!taylor && Math.abs(this[0])<1E-9 && !options.over) {
+ return this.Grade(2).Split().reduce((total,simple)=>{
+ var square = simple.Mul(simple).s, len = Math.sqrt(Math.abs(square));
+ if (len <= 1E-5) return total.Mul(Element.Scalar(1).Add(simple));
+ if (square < 0) return total.Mul(Element.Scalar(Math.cos(len)).Add(simple.Scale(Math.sin(len)/len)) );
+ return total.Mul(Element.Scalar(Math.cosh(len)).Add(simple.Scale(Math.sinh(len)/len)) );
+ },Element.Scalar(1));
+ }
+ if (options.dual) { var f=Math.exp(this.s); return this.map((x,i)=>i?x*f:f); }
+ var res = Element.Scalar(1), y=1, M= this.Scale(1), N=this.Scale(1); for (var x=1; x<15; x++) { res=res.Add(M.Scale(1/y)); M=M.Mul(N); y=y*(x+1); };
+ return res;
+ }
+
+ // Log - only for up to 3D PGA for now
+ Log (compat = false) {
+ if (options.over) return;
+ if (!compat) {
+ if (tot==4 && q==0 && r==1 && !options.over) { // https://www.researchgate.net/publication/360528787_Normalization_Square_Roots_and_the_Exponential_and_Logarithmic_Maps_in_Geometric_Algebras_of_Less_than_6D
+ if (Math.abs(this.s)>=.99999) return Element.Bivector(this[5],this[6],this[7],0,0,0).Scale(Math.sign(this.s));
+ var a = 1/(1 - this[0]*this[0]), b = Math.acos(this[0])*Math.sqrt(a), c = a*this[15]*(1 - this[0]*b);
+ return Element.Bivector( c*this[10] + b*this[5], -c*this[9] + b*this[6], c*this[8] + b*this[7], b*this[8], b*this[9], b*this[10] );
+ }
+ return this.Factorize().reduce((sum,bi)=>{
+ var [ci,si] = [bi.s, bi.Grade(2)];
+ var square = si.Mul(si).s;
+ var len = Math.sqrt(Math.abs(square));
+ if (Math.abs(square) < 1E-5) return sum.Add(si);
+ if (square < 0) return sum.Add(si.Scale(Math.acos(ci)/len));
+ return sum.Add(si.Scale(Math.acosh(ci)/len));
+ },Element.Scalar(0));
+ }
+ var b = this.Grade(2), bdb = Element.Dot(b,b).s;
+ if (Math.abs(bdb)<=1E-5) return this.s<0?b.Scale(-1):b;
+ var s = Math.sqrt(-bdb), bwb = Element.Wedge(b,b);
+ if (Math.abs(bwb[bwb.length-1])<=1E-5 || Math.abs(this.s)<=1E-5) return b.Scale(Math.atan2(s,this.s)/s);
+ var p = bwb.Scale(-1/(2*s));
+ return Element.Mul(Element.Mul((Element.Add(Math.atan2(s,this.s),p.Scale(1/this.s))),b),Element.Sub(s,p)).Scale(1/(s*s));
+ }
+
+ // Helper for efficient inverses. (custom involutions - negates grades in arguments).
+ Map () { var res=new Element(); return super.Map(res,...arguments); }
+
+ // Factories - Make it easy to generate vectors, bivectors, etc when using the functional API. None of the examples use this but
+ // users that have used other GA libraries will expect these calls. The Coeff() is used internally when translating algebraic literals.
+ static Element() { return new Element([...arguments]); };
+ static Coeff() { return (new Element()).Coeff(...arguments); }
+ static Scalar(x) { return (new Element()).Coeff(0,x); }
+ static Vector() { return (new Element()).nVector(1,...arguments); }
+ static Bivector() { return (new Element()).nVector(2,...arguments); }
+ static Trivector() { return (new Element()).nVector(3,...arguments); }
+ static nVector(n) { return (new Element()).nVector(...arguments); }
+
+ // Static operators. The parser will always translate operators to these static calls so that scalars, vectors, matrices and other non-multivectors can also be handled.
+ // The static operators typically handle functions and matrices, calling through to element methods for multivectors. They are intended to be flexible and allow as many
+ // types of arguments as possible. If performance is a consideration, one should use the generated element methods instead. (which only accept multivector arguments)
+ static toEl(x) { if (x instanceof Function) x=x(); if (!(x instanceof Element)) x=Element.Scalar(x); return x; }
+
+ // Addition and subtraction. Subtraction with only one parameter is negation.
+ static Add(a,b,res) {
+ // Resolve expressions passed in.
+ while(a.call)a=a(); while(b.call)b=b(); if (a.Add && b.Add) return a.Add(b,res);
+ // If either is a string, the result is a string.
+ if ((typeof a=='string')||(typeof b=='string')) return a.toString()+b.toString();
+ // If only one is an array, add the other element to each of the elements.
+ if ((a instanceof Array && !a.Add)^(b instanceof Array && !b.Add)) return (a instanceof Array)?a.map(x=>Element.Add(x,b)):b.map(x=>Element.Add(a,x));
+ // If both are equal length arrays, add elements one-by-one
+ if ((a instanceof Array)&&(b instanceof Array)&&a.length==b.length) return a.map((x,xi)=>Element.Add(x,b[xi]));
+ // If they're both not elements let javascript resolve it.
+ if (!(a instanceof Element || b instanceof Element)) return a+b;
+ // Here we're left with scalars and multivectors, call through to generated code.
+ a=Element.toEl(a); b=Element.toEl(b); return a.Add(b,res);
+ }
+
+ static Sub(a,b,res) {
+ // Resolve expressions passed in.
+ while(a.call)a=a(); while(b&&b.call) b=b(); if (a.Sub && b && b.Sub) return a.Sub(b,res);
+ // If only one is an array, add the other element to each of the elements.
+ if (b&&((a instanceof Array)^(b instanceof Array))) return (a instanceof Array)?a.map(x=>Element.Sub(x,b)):b.map(x=>Element.Sub(a,x));
+ // If both are equal length arrays, add elements one-by-one
+ if (b&&(a instanceof Array)&&(b instanceof Array)&&a.length==b.length) return a.map((x,xi)=>Element.Sub(x,b[xi]));
+ // Negation
+ if (arguments.length==1) return Element.Mul(a,-1);
+ // If none are elements here, let js do it.
+ if (!(a instanceof Element || b instanceof Element)) return a-b;
+ // Here we're left with scalars and multivectors, call through to generated code.
+ a=Element.toEl(a); b=Element.toEl(b); return a.Sub(b,res);
+ }
+
+ // The geometric product. (or matrix*matrix, matrix*vector, vector*vector product if called with 1D and 2D arrays)
+ static Mul(a,b,res) {
+ // Resolve expressions
+ while(a.call&&!a.length)a=a(); while(b.call&&!b.length)b=b(); if (a.Mul && b.Mul) return a.Mul(b,res);
+ // still functions -> experimental curry style (dont use this.)
+ if (a.call && b.call) return (ai,bi)=>Element.Mul(a(ai),b(bi));
+ // scalar mul.
+ if (Number.isFinite(a) && b.Scale) return b.Scale(a); else if (Number.isFinite(b) && a.Scale) return a.Scale(b);
+ // Handle matrices and vectors.
+ if ((a instanceof Array)&&(b instanceof Array)) {
+ // vector times vector performs a dot product. (which internally uses the GP on each component)
+ if((!(a[0] instanceof Array) || (a[0] instanceof Element)) &&(!(b[0] instanceof Array) || (b[0] instanceof Element))) { var r=tot?Element.Scalar(0):0; a.forEach((x,i)=>r=Element.Add(r,Element.Mul(x,b[i]),r)); return r; }
+ // Array times vector
+ if(!(b[0] instanceof Array)) return a.map((x,i)=>Element.Mul(a[i],b));
+ // Array times Array
+ var r=a.map((x,i)=>b[0].map((y,j)=>{ var r=tot?Element.Scalar(0):0; x.forEach((xa,k)=>r=Element.Add(r,Element.Mul(xa,b[k][j]))); return r; }));
+ // Return resulting array or scalar if 1 by 1.
+ if (r.length==1 && r[0].length==1) return r[0][0]; else return r;
+ }
+ // Only one is an array multiply each of its elements with the other.
+ if ((a instanceof Array)^(b instanceof Array)) return (a instanceof Array)?a.map(x=>Element.Mul(x,b)):b.map(x=>Element.Mul(a,x));
+ // Try js multiplication, else call through to geometric product.
+ var r=a*b; if (!isNaN(r)) return r;
+ a=Element.toEl(a); b=Element.toEl(b); return a.Mul(b,res);
+ }
+
+ // The inner product. (default is left contraction).
+ static LDot(a,b,res) {
+ // Expressions
+ while(a.call)a=a(); while(b.call)b=b(); //if (a.LDot) return a.LDot(b,res);
+ // Map elements in array
+ if (b instanceof Array && !(a instanceof Array)) return b.map(x=>Element.LDot(a,x));
+ if (a instanceof Array && !(b instanceof Array)) return a.map(x=>Element.LDot(x,b));
+ // js if numbers, else contraction product.
+ if (!(a instanceof Element || b instanceof Element)) return a*b;
+ a=Element.toEl(a);b=Element.toEl(b); return a.LDot(b,res);
+ }
+
+ // The symmetric inner product. (default is left contraction).
+ static Dot(a,b,res) {
+ // Expressions
+ while(a.call)a=a(); while(b.call)b=b(); //if (a.LDot) return a.LDot(b,res);
+ // js if numbers, else contraction product.
+ if (!(a instanceof Element || b instanceof Element)) return a|b;
+ a=Element.toEl(a);b=Element.toEl(b); return a.Dot(b,res);
+ }
+
+ // The outer product. (Grassman product - no use of metric)
+ static Wedge(a,b,res) {
+ // normal behavior for booleans/numbers
+ if (typeof a in {boolean:1,number:1} && typeof b in {boolean:1,number:1}) return a^b;
+ // Expressions
+ while(a.call)a=a(); while(b.call)b=b(); if (a.Wedge) return a.Wedge(Element.toEl(b),res);
+ // The outer product of two vectors is a matrix .. internally Mul not Wedge !
+ if (a instanceof Array && b instanceof Array) return a.map(xa=>b.map(xb=>Element.Mul(xa,xb)));
+ // js, else generated wedge product.
+ if (!(a instanceof Element || b instanceof Element)) return a*b;
+ a=Element.toEl(a);b=Element.toEl(b); return a.Wedge(b,res);
+ }
+
+ // The regressive product. (Dual of the outer product of the duals).
+ static Vee(a,b,res) {
+ // normal behavior for booleans/numbers
+ if (typeof a in {boolean:1,number:1} && typeof b in {boolean:1,number:1}) return a&b;
+ // Expressions
+ while(a.call)a=a(); while(b.call)b=b(); if (a.Vee) return a.Vee(Element.toEl(b),res);
+ // js, else generated vee product. (shortcut for dual of wedge of duals)
+ if (!(a instanceof Element || b instanceof Element)) return 0;
+ a=Element.toEl(a);b=Element.toEl(b); return a.Vee(b,res);
+ }
+
+ // The sandwich product. Provided for convenience (>>> operator)
+ static sw(a,b) {
+ // Skip strings/colors
+ if (typeof b == "string" || typeof b =="number") return b;
+ // Expressions
+ while(a.call)a=a(); while(b.call)b=b(); if (a.sw) return a.sw(b);
+ // Map elements in array
+ if (b instanceof Array && !b.Add) return b.map(x=>Element.sw(a,x));
+ // Call through. no specific generated code for it so just perform the muls.
+ a=Element.toEl(a); b=Element.toEl(b); return a.Mul(b).Mul(a.Reverse);
+ }
+
+ // Division - scalars or cal through to element method.
+ static Div(a,b,res) {
+ // Expressions
+ while(a.call)a=a(); while(b.call)b=b();
+ // For DDG experiments, I'll include a quick cholesky on matrices here. (vector/matrix)
+ if ((a instanceof Array) && (b instanceof Array) && (b[0] instanceof Array)) {
+ // factor
+ var R = b.flat(), i, j, k, sum, i_n, j_n, n=b[0].length, s=new Array(n), x=new Array(n), yi;
+ for (i=0;i=0; i--) for (x[i] /= R[i*n+i], j=i+1; ja.map((r,ri)=>Element.Conjugate(a[ri][ci]))); return Element.toEl(a).Conjugate; }
+ static Normalize(a) { return Element.toEl(a).Normalized; };
+ static Length(a) { return Element.toEl(a).Length };
+
+ // Comparison operators always use length. Handle expressions, then js or length comparison
+ static eq(a,b) { if (!(a instanceof Element)||!(b instanceof Element)) return a==b; while(a.call)a=a(); while(b.call)b=b(); for (var i=0; i(b instanceof Element?b.Length:b); }
+ static lte(a,b) { while(a.call)a=a(); while(b.call)b=b(); return (a instanceof Element?a.Length:a)<=(b instanceof Element?b.Length:b); }
+ static gte(a,b) { while(a.call)a=a(); while(b.call)b=b(); return (a instanceof Element?a.Length:a)>=(b instanceof Element?b.Length:b); }
+
+ // Debug output and printing multivectors.
+ static describe(x) { if (x===true) console.log(`Basis\n${basis}\nMetric\n${metric.slice(1,1+tot)}\nCayley\n${mulTable.map(x=>(x.map(x=>(' '+x).slice(-2-tot)))).join('\n')}\nMatrix Form:\n`+gp.map(x=>x.map(x=>x.match(/(-*b\[\d+\])/)).map(x=>x&&((x[1].match(/-/)||' ')+String.fromCharCode(65+1*x[1].match(/\d+/)))||' 0')).join('\n')); return {basis:basisg||basis,metric,mulTable,matrix:gp.map(x=>x.map(x=>x.replace(/\*this\[.+\]/,'').replace(/b\[(\d+)\]/,(a,x)=>(metric[x]==-1||metric[x]==0&&grades[x]>1&&(-1)**grades[x]==(metric[basis.indexOf(basis[x].replace('0',''))]||(-1)**grades[x])?'-':'')+basis[x]).replace('--','')))} }
+
+ // Direct sum of algebras - experimental
+ static sum(B){
+ var A = Element;
+ // Get the multiplication tabe and basis.
+ var T1 = A.describe().mulTable, T2 = B.describe().mulTable;
+ var B1 = A.describe().basis, B2 = B.describe().basis;
+ // Get the maximum index of T1, minimum of T2 and rename T2 if needed.
+ var max_T1 = B1.filter(x=>x.match(/e/)).map(x=>x.match(/\d/g)).flat().map(x=>x|0).sort((a,b)=>b-a)[0];
+ var max_T2 = B2.filter(x=>x.match(/e/)).map(x=>x.match(/\d/g)).flat().map(x=>x|0).sort((a,b)=>b-a)[0];
+ var min_T2 = B2.filter(x=>x.match(/e/)).map(x=>x.match(/\d/g)).flat().map(x=>x|0).sort((a,b)=>a-b)[0];
+ // remapping ..
+ T2 = T2.map(x=>x.map(y=>y.match(/e/)?y.replace(/(\d)/g,(x)=>(x|0)+max_T1):y.replace("1","e"+(1+max_T2+max_T1))));
+ B2 = B2.map((y,i)=>i==0?y.replace("1","e"+(1+max_T2+max_T1)):y.replace(/(\d)/g,(x)=>(x|0)+max_T1));
+ // Build the new basis and multable..
+ var basis = [...B1,...B2];
+ var Cayley = T1.map((x,i)=>[...x,...T2[0].map(x=>"0")]).concat(T2.map((x,i)=>[...T1[0].map(x=>"0"),...x]))
+ // Build the new algebra.
+ var grades = [...B1.map(x=>x=="1"?0:x.length-1),...B2.map((x,i)=>i?x.length-1:0)];
+ var a = Algebra({basis,Cayley,grades,tot:Math.log2(B1.length)+Math.log2(B2.length)})
+ // And patch up ..
+ a.Scalar = function(x) {
+ var res = new a();
+ for (var i=0; i function of 1 parameter will be called with that parameter from -1 to 1 and graphed on a canvas. Returned values should also be in the [-1 1] range
+ // graph(function(x,y)) => functions of 2 parameters will be called from -1 to 1 on both arguments. Returned values can be 0-1 for greyscale or an array of three RGB values.
+ // graph(array) => array of algebraic elements (points, lines, circles, segments, texts, colors, ..) is graphed.
+ // graph(function=>array) => same as above, for animation scenario's this function is called each frame.
+ // An optional second parameter is an options object { width, height, animate, camera, scale, grid, canvas }
+ static graph(f,options) {
+ // Store the original input
+ if (!f) return; var origf=f;
+ // generate default options.
+ options=options||{}; options.scale=options.scale||1; options.camera=options.camera||(tot!=4?Element.Scalar(1): ( Element.Bivector(0,0,0,0,0,options.p||0).Exp() ).Mul( Element.Bivector(0,0,0,0,options.h||0,0).Exp()) );
+ if (options.conformal && tot==4) var ni = options.ni||this.Coeff(4,1,3,1), no = options.no||this.Coeff(4,0.5,3,-0.5), minus_no = no.Scale(-1);
+ var ww=options.width, hh=options.height, cvs=options.canvas, tpcam=new Element([0,0,0,0,0,0,0,0,0,0,0,-5,0,0,1,0]),tpy=this.Coeff(4,1),tp=new Element(),
+ // project 3D to 2D. This allows to render 3D and 2D PGA with the same code.
+ project=(o)=>{ if (!o) return o; while (o.call) o=o();
+// if (o instanceof Element && o.length == 32) o = new Element([o[0],o[1],o[2],o[3],o[4],o[6],o[7],o[8],o[10],o[11],o[13],o[16],o[17],o[19],o[22],o[26]]);
+ // Clip 3D lines so they don't go past infinity.
+ if (o instanceof Element && o.length == 16 && o[8]**2+o[9]**2+o[10]**2>0.0001) {
+ o = [[options.clip||2,1,0,0],[-(options.clip||2),1,0,0],[options.clip||2,0,1,0],[-(options.clip||2),0,1,0],[options.clip||2,0,0,1],[-(options.clip||2),0,0,1]].map(v=>{
+ var r = Element.Vector(...v).Wedge(o); return r[14]?r.Scale(1/r[14], r):undefined;
+ }).filter(x=>x && Math.abs(x[13])<= (options.clip||2)+0.001 && Math.abs(x[12]) <= (options.clip||2)+0.001 && Math.abs(x[11]) <= (options.clip||2) + 0.001).slice(0,2);
+ return o.map(o=>(tpcam).Vee(options.camera.Mul(o).Mul(options.camera.Conjugate)).Wedge(tpy));
+ }
+ // Convert 3D planes to polies.
+ if (o instanceof Element && o.length == 16 && o.Grade(1).Length>0.01) {
+ var m = Element.Add(1, Element.Mul(o.Normalized, Element.Coeff(3,1))).Normalized, e0 = 0;
+ o=Element.sw(m,[[-1,-1],[-1,1],[1,1],[-1,-1],[1,1],[1,-1]].map(([x,z])=>Element.Trivector(x*o.Length,e0,z*o.Length,1)));
+ return o.map(o=>(tpcam).Vee(options.camera.Mul(o).Mul(options.camera.Conjugate)).Wedge(tpy));
+ }
+ return (tot==4 && o instanceof Element && o.length==16)?(tpcam).Vee(options.camera.Mul(o).Mul(options.camera.Conjugate)).Wedge(tpy):(o.length==2**tot)?Element.sw(options.camera,o):o;
+ };
+ // gl escape.
+ if (options.gl && !(tot==4 && options.conformal)) return Element.graphGL(f,options); if (options.up) return Element.graphGL2(f,options);
+ // if we get an array or function without parameters, we render c2d or p2d SVG points/lines/circles/etc
+ if (!(f instanceof Function) || f.length===0) {
+ // Our current cursor, color, animation state and 2D mapping.
+ var lx,ly,lr,color,res,anim=false,to2d=(tot==5)?[0, 8, 11, 13, 19, 17, 22, 26]:(tot==3)?[0,1,2,3,4,5,6,7]:[0,7,9,10,13,12,14,15];
+ // Make sure we have an array of elements. (if its an object, convert to array with elements and names.)
+ if (f instanceof Function) f=f(); if (!(f instanceof Array)) f=[].concat.apply([],Object.keys(f).map((k)=>typeof f[k]=='number'?[f[k]]:[f[k],k]));
+ // The build function generates the actual SVG. It will be called everytime the user interacts or the anim flag is set.
+ function build(f,or) {
+ // Make sure we have an aray.
+ if (or && f && f instanceof Function) f=f();
+ // Reset position and color for cursor.
+ lx=-2;ly=options.conformal?-1.85:1.85;lr=0;color='#444';
+ // Create the svg element. (master template string till end of function)
+ var svg=new DOMParser().parseFromString(`