diff --git a/engine-proto/alg-test/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl
similarity index 100%
rename from engine-proto/alg-test/Engine.Algebraic.jl
rename to engine-proto/Engine.Algebraic.jl
diff --git a/engine-proto/alg-test/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl
similarity index 100%
rename from engine-proto/alg-test/Engine.Numerical.jl
rename to engine-proto/Engine.Numerical.jl
diff --git a/engine-proto/alg-test/Engine.jl b/engine-proto/Engine.jl
similarity index 100%
rename from engine-proto/alg-test/Engine.jl
rename to engine-proto/Engine.jl
diff --git a/engine-proto/alg-test/HittingSet.jl b/engine-proto/HittingSet.jl
similarity index 100%
rename from engine-proto/alg-test/HittingSet.jl
rename to engine-proto/HittingSet.jl
diff --git a/engine-proto/alg-test/ConstructionViewer.jl b/engine-proto/alg-test/ConstructionViewer.jl
deleted file mode 100644
index b9c8ffb..0000000
--- a/engine-proto/alg-test/ConstructionViewer.jl
+++ /dev/null
@@ -1,223 +0,0 @@
-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/ganja-test/ganja-test.html b/engine-proto/ganja-test/ganja-test.html
deleted file mode 100644
index 0207dcc..0000000
--- a/engine-proto/ganja-test/ganja-test.html
+++ /dev/null
@@ -1,96 +0,0 @@
-
-
-
-
-
-
-
-
-
-
-
diff --git a/engine-proto/ganja-test/ganja-test.jl b/engine-proto/ganja-test/ganja-test.jl
deleted file mode 100644
index 6c55061..0000000
--- a/engine-proto/ganja-test/ganja-test.jl
+++ /dev/null
@@ -1,127 +0,0 @@
-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
deleted file mode 100644
index ac5fe54..0000000
--- a/engine-proto/gram-test/Engine.jl
+++ /dev/null
@@ -1,451 +0,0 @@
-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
deleted file mode 100644
index 5c03c01..0000000
--- a/engine-proto/gram-test/basin-shapes.jl
+++ /dev/null
@@ -1,99 +0,0 @@
-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
deleted file mode 100644
index 1bd22a7..0000000
--- a/engine-proto/gram-test/circles-in-triangle.jl
+++ /dev/null
@@ -1,76 +0,0 @@
-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
deleted file mode 100644
index 1e95e42..0000000
--- a/engine-proto/gram-test/ganja-1.0.204.js
+++ /dev/null
@@ -1,1913 +0,0 @@
-/** 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(`