Stow algebraic engine prototype

We're using the Gram matrix engine for the next stage of development,
so the algebraic engine shouldn't be at the top level anymore.

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)

engine-proto/ganja-test/ganja-test.html

@@ -0,0 +1,96 @@
+<!DOCTYPE html>
+<html>
+<head>
+  <style>
+    body {
+      background-color: #ffe0f0;
+    }
+    
+    /* needed to keep Ganja canvas from blowing up */
+    canvas {
+      min-width: 600px;
+      max-width: 600px;
+      min-height: 600px;
+      max-height: 600px;
+    }
+  </style>
+  <script src="https://unpkg.com/ganja.js"></script>
+</head>
+<body>
+  <p><button onclick="flip()">Flip</button></p>
+  <script>
+  // in the default view, e4 + e5 is the point at infinity
+  let CGA3 = Algebra(4, 1);
+  let elements = [
+    CGA3.inline(() => Math.sqrt(0.5)*( 1e1 + 1e2 + 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*( 1e1 - 1e2 - 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*(-1e1 + 1e2 - 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*(-1e1 - 1e2 + 1e3 + 1e5))(),
+    CGA3.inline(() => -Math.sqrt(3)*1e4 + Math.sqrt(2)*1e5)()
+  ];
+  /*
+    these blocks of commented-out code can be used to confirm that a spacelike
+    vector and its Hodge dual represent the same generalized sphere
+  */
+  /*let elements = [
+    CGA3.inline(() => Math.sqrt(0.5)*!( 1e1 + 1e2 + 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*!( 1e1 - 1e2 - 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*!(-1e1 + 1e2 - 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*!(-1e1 - 1e2 + 1e3 + 1e5))(),
+    CGA3.inline(() => !(-Math.sqrt(3)*1e4 + Math.sqrt(2)*1e5))()
+  ];*/
+  /*let elements = [
+    CGA3.inline(() =>  1e1 + 1e5)(),
+    CGA3.inline(() =>  1e2 + 1e5)(),
+    CGA3.inline(() =>  1e3 + 1e5)(),
+    CGA3.inline(() => -1e4 + 1e5)(),
+    CGA3.inline(() => Math.sqrt(0.5)*(1e1 + 1e2 + 1e3 + 1e5))(),
+    CGA3.inline(() => Math.sqrt(0.5)*!(1e1 + 1e2 + 1e3 - 0.01e4 + 1e5))()
+  ];*/
+  
+  // set up palette
+  var colorIndex;
+  var palette = [0xff00b0, 0x00ffb0, 0x00b0ff, 0x8040ff, 0xc0c0c0];
+  function nextColor() {
+    colorIndex = (colorIndex + 1) % palette.length;
+    return palette[colorIndex];
+  }
+  function resetColorCycle() {
+    colorIndex = palette.length - 1;
+  }
+  resetColorCycle();
+  
+  // create scene function
+  function scene() {
+    commands = [];
+    resetColorCycle();
+    elements.forEach((elt) => commands.push(nextColor(), elt));
+    return commands;
+  }
+  
+  // initialize graph
+  let graph = CGA3.graph(
+    scene,
+    {
+      conformal: true, gl: true, grid: true
+    }
+  )
+  document.body.appendChild(graph);
+  
+  function flip() {
+    let last = elements.length - 1;
+    for (let n = 0; n < last; ++n) {
+      // reflect
+      elements[n] = CGA3.Mul(CGA3.Mul(elements[last], elements[n]), elements[last]);
+      
+      // de-noise
+      for (let k = 6; k < elements[n].length; ++k) {
+      /*for (let k = 0; k < 26; ++k) {*/
+        elements[n][k] = 0;
+      }
+    }
+    requestAnimationFrame(graph.update.bind(graph, scene));
+  }
+  </script>
+</body>
+</html>

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, """
+  <p><button id="flip-button" onclick="flip()">Flip</button></p>
+""", async = false)
+
+# === set up visualization ===
+
+# list elements. in the default view, e4 + e5 is the point at infinity
+elements = sqrt(0.5) * BigFloat[
+  1  1 -1 -1  0;
+  1 -1  1 -1  0;
+  1 -1 -1  1  0;
+  0  0  0  0 -sqrt(6);
+  1  1  1  1  2
+]
+
+# load elements
+for vec in eachcol(elements)
+  add_element!(vec)
+end
+
+# initialize visualization
+@js win begin
+  graph = CGA3.graph(
+    scene,
+    Dict(
+      "conformal" => true,
+      "gl" => true,
+      "grid" => true
+    )
+  )
+  document.body.appendChild(graph)
+end

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

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

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")

engine-proto/gram-test/gram-test.jl

@@ -0,0 +1,85 @@
+using LinearAlgebra
+using AbstractAlgebra
+
+function printgood(msg)
+  printstyled("✓", color = :green)
+  println(" ", msg)
+end
+
+function printbad(msg)
+  printstyled("✗", color = :red)
+  println(" ", msg)
+end
+
+F, gens = rational_function_field(AbstractAlgebra.Rationals{BigInt}(), ["a₁", "a₂", "b₁", "b₂", "c₁", "c₂"])
+a = gens[1:2]
+b = gens[3:4]
+c = gens[5:6]
+
+# three mutually tangent spheres which are all perpendicular to the x, y plane
+gram = [
+  -1  1  1;
+   1 -1  1;
+   1  1 -1
+]
+
+eig = eigen(gram)
+n_pos = count(eig.values .> 0.5)
+n_neg = count(eig.values .< -0.5)
+if n_pos + n_neg == size(gram, 1)
+  printgood("Non-degenerate subspace")
+else
+  printbad("Degenerate subspace")
+end
+sig_rem = Int64[ones(1-n_pos); -ones(4-n_neg)]
+unk = hcat(a, b, c)
+M = matrix_space(F, 5, 5)
+big_gram = M(F.([
+  diagm(sig_rem) unk;
+  transpose(unk) gram
+]))
+
+r, p, L, U = lu(big_gram)
+if isone(p)
+  printgood("Found a solution")
+else
+  printbad("Didn't find a solution")
+end
+solution = transpose(L)
+mform = U * inv(solution)
+
+vals = [0, 0, 0, 1, 0, -3//4]
+solution_ex = [evaluate(entry, vals) for entry in solution]
+mform_ex = [evaluate(entry, vals) for entry in mform]
+
+std_basis = [
+  0 0 0 1  1;
+  0 0 0 1 -1;
+  1 0 0 0  0;
+  0 1 0 0  0;
+  0 0 1 0  0
+]
+std_solution = M(F.(std_basis)) * solution
+std_solution_ex = std_basis * solution_ex
+
+println("Minkowski form:")
+display(mform_ex)
+
+big_gram_recovered = transpose(solution_ex) * mform_ex * solution_ex
+valid = all(iszero.(
+  [evaluate(entry, vals) for entry in big_gram] - big_gram_recovered
+))
+if valid
+  printgood("Recovered Gram matrix:")
+else
+  printbad("Didn't recover Gram matrix. Instead, got:")
+end
+display(big_gram_recovered)
+
+# this should be a solution
+hand_solution = [0 0 1 0 0; 0 0 -1 2 2; 0 0 0 1 -1; 1 0 0 0 0; 0 1 0 0 0]
+unmix = Rational{Int64}[[1//2 1//2; 1//2 -1//2] zeros(Int64, 2, 3); zeros(Int64, 3, 2) Matrix{Int64}(I, 3, 3)]
+hand_solution_diag = unmix * hand_solution
+big_gram_hand_recovered = transpose(hand_solution_diag) * diagm([1; -ones(Int64, 4)]) * hand_solution_diag
+println("Gram matrix from hand-written solution:")
+display(big_gram_hand_recovered)

engine-proto/gram-test/gram-test.sage

@@ -0,0 +1,27 @@
+F = QQ['a', 'b', 'c'].fraction_field()
+a, b, c = F.gens()
+
+# three mutually tangent spheres which are all perpendicular to the x, y plane
+gram = matrix([
+  [-1, 0,  0,  0,  0],
+  [0, -1,  a,  b,  c],
+  [0,  a, -1,  1,  1],
+  [0,  b,  1, -1,  1],
+  [0,  c,  1,  1, -1]
+])
+
+P, L, U = gram.LU()
+solution = (P * L).transpose()
+mform = U * L.transpose().inverse()
+
+concrete = solution.subs({a: 0, b: 1, c: -3/4})
+
+std_basis = matrix([
+  [0, 0, 0, 1,  1],
+  [0, 0, 0, 1, -1],
+  [1, 0, 0, 0,  0],
+  [0, 1, 0, 0,  0],
+  [0, 0, 1, 0,  0]
+])
+std_solution = std_basis * solution
+std_concrete = std_basis * concrete

engine-proto/gram-test/irisawa-hexlet.jl

@@ -0,0 +1,77 @@
+include("Engine.jl")
+
+using SparseArrays
+
+# this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article below
+# includes a nice translation of the problem statement, which was recorded in
+# Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and Present_)
+#
+#   "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki
+#   https://www.nippon.com/en/japan-topics/c12801/
+#
+
+# initialize the partial gram matrix
+J = Int64[]
+K = Int64[]
+values = BigFloat[]
+for s in 1:9
+  # each sphere is represented by a spacelike vector
+  push!(J, s)
+  push!(K, s)
+  push!(values, 1)
+  
+  # the circumscribing sphere is internally tangent to all of the other spheres
+  if s > 1
+    append!(J, [1, s])
+    append!(K, [s, 1])
+    append!(values, [1, 1])
+  end
+  
+  if s > 3
+    # each chain sphere is externally tangent to the "sun" and "moon" spheres
+    for n in 2:3
+      append!(J, [s, n])
+      append!(K, [n, s])
+      append!(values, [-1, -1])
+    end
+    
+    # each chain sphere is externally tangent to the next chain sphere
+    s_next = 4 + mod(s-3, 6)
+    append!(J, [s, s_next])
+    append!(K, [s_next, s])
+    append!(values, [-1, -1])
+  end
+end
+gram = sparse(J, K, values)
+
+# make an initial guess
+guess = hcat(
+  Engine.sphere(BigFloat[0, 0, 0], BigFloat(15)),
+  Engine.sphere(BigFloat[0, 0, -9], BigFloat(5)),
+  Engine.sphere(BigFloat[0, 0, 11], BigFloat(3)),
+  (
+    Engine.sphere(9*BigFloat[cos(k*π/3), sin(k*π/3), 0], BigFloat(2.5))
+    for k in 1:6
+  )...
+)
+frozen = [CartesianIndex(4, k) for k in 1:4]
+
+# 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")
+if success
+  println("Chain diameters:")
+  println("  ", 1 / L[4,4], " sun (given)")
+  for k in 5:9
+    println("  ", 1 / L[4,k], " sun")
+  end
+end

engine-proto/gram-test/low-rank-test.jl

@@ -0,0 +1,49 @@
+using LowRankModels
+using LinearAlgebra
+using SparseArrays
+
+# testing Gram matrix recovery using the LowRankModels package
+
+# initialize the partial gram matrix for an arrangement of seven spheres in
+# which spheres 1 through 5 are mutually tangent, and spheres 3 through 7 are
+# also mutually tangent
+I = Int64[]
+J = Int64[]
+values = Float64[]
+for i in 1:7
+  for j in 1:7
+    if (i <= 5 && j <= 5) || (i >= 3 && j >= 3)
+      push!(I, i)
+      push!(J, j)
+      push!(values, i == j ? 1 : -1)
+    end
+  end
+end
+gram = sparse(I, J, values)
+
+# in this initial guess, the mutual tangency condition is satisfied for spheres
+# 1 through 5
+X₀ = sqrt(0.5) * [
+  1 0  1  1  1;
+  1 0  1 -1 -1;
+  1 0 -1  1 -1;
+  1 0 -1 -1  1;
+  2 -sqrt(6) 0 0 0;
+  0.2  0.3 -0.1 -0.2  0.1;
+  0.1 -0.2  0.3  0.4 -0.1
+]'
+Y₀ = diagm([-1, 1, 1, 1, 1]) * X₀
+
+# search parameters
+search_params = ProxGradParams(
+  1.0;
+  max_iter = 100,
+  inner_iter = 1,
+  abs_tol = 1e-16,
+  rel_tol = 1e-9,
+  min_stepsize = 0.01
+)
+
+# complete gram matrix
+model = GLRM(gram, QuadLoss(), ZeroReg(), ZeroReg(), 5, X = X₀, Y = Y₀)
+X, Y, history = fit!(model, search_params)

engine-proto/gram-test/overlap-test.jl

@@ -0,0 +1,37 @@
+using LinearAlgebra
+using AbstractAlgebra
+
+function printgood(msg)
+  printstyled("✓", color = :green)
+  println(" ", msg)
+end
+
+function printbad(msg)
+  printstyled("✗", color = :red)
+  println(" ", msg)
+end
+
+F, gens = rational_function_field(AbstractAlgebra.Rationals{BigInt}(), ["x", "t₁", "t₂", "t₃"])
+x = gens[1]
+t = gens[2:4]
+
+# three mutually tangent spheres which are all perpendicular to the x, y plane
+M = matrix_space(F, 7, 7)
+gram = M(F[
+   1    -1    -1 -1 -1  t[1] t[2];
+  -1     1    -1 -1 -1  x    t[3]
+  -1    -1     1 -1 -1 -1   -1;
+  -1    -1    -1  1 -1 -1   -1;
+  -1    -1    -1 -1  1 -1   -1;
+   t[1]  x    -1 -1 -1  1   -1;
+   t[2]  t[3] -1 -1 -1 -1    1
+])
+
+r, p, L, U = lu(gram)
+if isone(p)
+  printgood("Found a solution")
+else
+  printbad("Didn't find a solution")
+end
+solution = transpose(L)
+mform = U * inv(solution)

engine-proto/gram-test/overlapping-pyramids.jl

@@ -0,0 +1,90 @@
+include("Engine.jl")
+
+using SparseArrays
+using AbstractAlgebra
+using PolynomialRoots
+using Random
+
+# initialize the partial gram matrix for an arrangement of seven spheres in
+# which spheres 1 through 5 are mutually tangent, and spheres 3 through 7 are
+# also mutually tangent
+J = Int64[]
+K = Int64[]
+values = BigFloat[]
+for j in 1:7
+  for k in 1:7
+    if (j <= 5 && k <= 5) || (j >= 3 && k >= 3)
+      push!(J, j)
+      push!(K, k)
+      push!(values, j == k ? 1 : -1)
+    end
+  end
+end
+gram = sparse(J, K, values)
+
+# set the independent variable
+indep_val = -9//5
+gram[6, 1] = BigFloat(indep_val)
+gram[1, 6] = gram[6, 1]
+
+# in this initial guess, the mutual tangency condition is satisfied for spheres
+# 1 through 5
+Random.seed!(50793)
+guess = let
+  a = sqrt(BigFloat(3)/2)
+  hcat(
+    sqrt(1/BigFloat(2)) * 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
+    ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
+    Engine.rand_on_shell(fill(BigFloat(-1), 2))
+  )
+end
+
+# complete the gram matrix using Newton's method with backtracking
+L, success, history = Engine.realize_gram(gram, guess)
+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")
+
+# === algebraic check ===
+
+#=
+R, gens = polynomial_ring(AbstractAlgebra.Rationals{BigInt}(), ["x", "t₁", "t₂", "t₃"])
+x = gens[1]
+t = gens[2:4]
+
+S, u = polynomial_ring(AbstractAlgebra.Rationals{BigInt}(), "u")
+
+M = matrix_space(R, 7, 7)
+gram_symb = M(R[
+   1    -1    -1 -1 -1  t[1] t[2];
+  -1     1    -1 -1 -1  x    t[3]
+  -1    -1     1 -1 -1 -1   -1;
+  -1    -1    -1  1 -1 -1   -1;
+  -1    -1    -1 -1  1 -1   -1;
+   t[1]  x    -1 -1 -1  1   -1;
+   t[2]  t[3] -1 -1 -1 -1    1
+])
+rank_constraints = det.([
+  gram_symb[1:6, 1:6],
+  gram_symb[2:7, 2:7],
+  gram_symb[[1, 3, 4, 5, 6, 7], [1, 3, 4, 5, 6, 7]]
+])
+
+# solve for x and t
+x_constraint = 25//16 * to_univariate(S, evaluate(rank_constraints[1], [2], [indep_val]))
+t₂_constraint = 25//16 * to_univariate(S, evaluate(rank_constraints[3], [2], [indep_val]))
+x_vals = PolynomialRoots.roots(x_constraint.coeffs)
+t₂_vals = PolynomialRoots.roots(t₂_constraint.coeffs)
+=#

engine-proto/gram-test/sphere-in-tetrahedron.jl

@@ -0,0 +1,67 @@
+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:6
+  for k in 1:6
+    filled = false
+    if j == 6
+      if k <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif k == 6
+      if j <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == k
+      push!(values, 1)
+      filled = true
+    elseif j <= 4 && k <= 4
+      push!(values, -1/BigFloat(3))
+      filled = true
+    else
+      push!(values, -1)
+      filled = true
+    end
+    if filled
+      push!(J, j)
+      push!(K, k)
+    end
+  end
+end
+gram = sparse(J, K, values)
+
+# set initial guess
+Random.seed!(99230)
+guess = hcat(
+  sqrt(1/BigFloat(3)) * BigFloat[
+    1  1 -1 -1  0
+    1 -1  1 -1  0
+    1 -1 -1  1  0
+    0  0  0  0  1.5
+    1  1  1  1 -0.5
+  ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
+  BigFloat[0, 0, 0, 0, 1]
+)
+frozen = [CartesianIndex(j, 6) 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")

engine-proto/gram-test/tetrahedron-radius-ratio.jl

@@ -0,0 +1,96 @@
+include("Engine.jl")
+
+using LinearAlgebra
+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:11
+  for k in 1:11
+    filled = false
+    if j == 11
+      if k <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif k == 11
+      if j <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == k
+      push!(values, j <= 6 ? 1 : 0)
+      filled = true
+    elseif j <= 4
+      if k <= 4
+        push!(values, -1/BigFloat(3))
+        filled = true
+      elseif k == 5
+        push!(values, -1)
+        filled = true
+      elseif 7 <= k <= 10 && k - j != 6
+        push!(values, 0)
+        filled = true
+      end
+    elseif k <= 4
+      if j == 5
+        push!(values, -1)
+        filled = true
+      elseif 7 <= j <= 10 && j - k != 6
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == 6 && 7 <= k <= 10 || k == 6 && 7 <= j <= 10
+      push!(values, 0)
+      filled = true
+    end
+    if filled
+      push!(J, j)
+      push!(K, k)
+    end
+  end
+end
+gram = sparse(J, K, values)
+
+# set initial guess
+Random.seed!(99230)
+guess = hcat(
+  sqrt(1/BigFloat(3)) * BigFloat[
+    1  1 -1 -1  0    0
+    1 -1  1 -1  0    0
+    1 -1 -1  1  0    0
+    0  0  0  0  1.5  0.5
+    1  1  1  1 -0.5 -1.5
+  ] + 0.0*Engine.rand_on_shell(fill(BigFloat(-1), 6)),
+  Engine.point([-0.5, -0.5, -0.5] + 0.3*randn(3)),
+  Engine.point([-0.5,  0.5,  0.5] + 0.3*randn(3)),
+  Engine.point([ 0.5, -0.5,  0.5] + 0.3*randn(3)),
+  Engine.point([ 0.5,  0.5, -0.5] + 0.3*randn(3)),
+  BigFloat[0, 0, 0, 0, 1]
+)
+frozen = vcat(
+  [CartesianIndex(4, k) for k in 7:10],
+  [CartesianIndex(j, 11) 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])
+if success
+  infty = BigFloat[0, 0, 0, 0, 1]
+  radius_ratio = dot(infty, Engine.Q * L[:,5]) / dot(infty, Engine.Q * L[:,6])
+  println("\nCircumradius / inradius: ", radius_ratio)
+end

notes/inversive.md

@@ -6,22 +6,22 @@ These coordinates are of form $I=(c, r, x, y, z)$ where we think of $c$ as the c
 
 | Entity or Relationship                                                         | Representation                                                                                                                                                                                                                                                                                                                                                                                                                       | Comments/questions                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
 | ------------------------------------------------------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
-| Sphere s with radius r>0 centered on P = (x,y,z)                               | $I_s = (1/c, 1/r, x/r, y/r, z/r)$ satisfying $Q(I_s,I_s) = -1$, i.e., $c = r/(\|P\|^2 - r^2)$.                                                                                                                                                                                                                                                                                                                                       | Can also write $I_s = (\|P\|^2/r - r, 1/r, x/r. y/r, z/r)$ -- so there is no trouble if $\|E_{I_s}\| = r$, just get first coordinate to be 0.                                                                                                                                                                                                                                                                                                                                                                                                                                            |
-| Plane p with unit normal (x,y,z), a distance s from origin                     | $I_p = (2s, 0, x, y, z)$                                                                                                                                                                                                                                                                                                                                                                                                             | Note $Q(I_p, I_p)$ is still -1. Also, there are two representations for each plane through the origin, namely $(0,0,x,y,z)$ and $(0,0,-x,-y,-z)$                                                                                                                                                                                                                                                                                                                                                                                                                                         |
+| Sphere s with radius r>0 centered on P = (x,y,z)                               | $I_s = (1/c, 1/r, x/r, y/r, z/r)$ satisfying $Q(I_s,I_s) = -1$, i.e., $c = r/(\|P\|^2 - r^2)$.                                                                                                                                                                                                                                                                                                                                       | Can also write $I_s = (\|P\|^2/r - r, 1/r, x/r. y/r, z/r)$—so there is no trouble if $\|E_{I_s}\| = r$, just get first coordinate to be 0.                                                                                                                                                                                                                                                                                                                                                                                                                                               |
+| Plane p with unit normal (x,y,z) through the point s(x,y,z)                    | $I_p = (-2s, 0, -x, -y, -z)$                                                                                                                                                                                                                                                                                                                                                                                                         | Note $Q(I_p, I_p)$ is still −1.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
 | Point P with Euclidean coordinates (x,y,z)                                     | $I_P = (\|P\|^2, 1, x, y, z)$                                                                                                                                                                                                                                                                                                                                                                                                        | Note $Q(I_P,I_P) = 0$.  Because of this we might choose  some other scaling of the inversive coordinates, say $(\||P\||,1/\||P\||,x/\||P\||,y/\||P\||,z/\||P\||)$ instead, but that fails at the origin, and likely won't have some of the other nice properties listed below.  Note that scaling just the co-radius by $s$ and the radius by $1/s$ (which still preserves $Q=0$) dilates by a factor of $s$ about the origin, so that $(\|P\|, \|P\|, x, y, z)$, which might look symmetric, would actually have to represent the Euclidean point $(x/\||P\||, y/\||P\||, z/\||P\||)$ . |
 | ∞, the "point at infinity"                                                     | $I_\infty = (1,0,0,0,0)$                                                                                                                                                                                                                                                                                                                                                                                                             | The only solution to $Q(I,I) = 0$ not covered by the above case.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
 | P lies on sphere or plane given by I                                           | $Q(I_P, I) = 0$                                                                                                                                                                                                                                                                                                                                                                                                                      |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
-| Sphere/planes represented by I and J are tangent                               | $Q(I,J) = 1$ (??, see note at right)                                                                                                                                                                                                                                                                                                                                                                                                 | Seems as though this must be $Q(I,J) = \pm1$  ? For example, the $xy$ plane represented by (0,0,0,0,1)  is tangent to the unit circle centered at (0,0,1) rep'd by (0,1,0,0,1), but their Q-product is -1. And in general you can reflect any sphere tangent to any plane through the plane and it should flip the sign of $Q(I,J)$, if I am not mistaken.                                                                                                                                                                                                                               |
+| Sphere/planes represented by I and J are tangent                               | If $I$ and $J$ have the same orientation where they touch, $Q(I,J) = -1$. If they have opposing orientations, $Q(I,J) = 1$.                                                                                                                                                                                                                                                                                                          | For example, the $xy$ plane with normal $-e_z$, represented by $(0,0,0,0,1)$, is tangent with matching orientation to the unit sphere centered at $(0,0,1)$ with outward normals, represented by $(0,1,0,0,1)$. Accordingly, their $Q$-product is −1.                                                                                                                                                                                                                                                                                                                                    |
 | Sphere/planes represented by I and J intersect (respectively, don't intersect) | $\|Q(I,J)\| < (\text{resp. }>)\; 1$                                                                                                                                                                                                                                                                                                                                                                                                  | Follows from the angle formula, at least conceptually.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   |
 | P is center of sphere represented by I                                         | Well, $Q(I_P, I)$ comes out to be $(\|P\|^2/r - r + \|P\|^2/r)/2 - \|P\|^2/r$ or just $-r/2$ .                                                                                                                                                                                                                                                                                                                                       | Is it if and only if ?   No this probably doesn't work because center is not conformal quantity.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
 | Distance between P and R is d                                                  | $Q(I_P, I_R) = d^2/2$                                                                                                                                                                                                                                                                                                                                                                                                                |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
 | Distance between P and sphere/plane rep by I                                   |                                                                                                                                                                                                                                                                                                                                                                                                                                      | In the very simple case of a plane $I$ rep'd by $(2s, 0, x, y, z)$ and a point $P$ that lies on its perpendicular through the origin, rep'd by $(r^2, 1, rx, ry, rz)$ we get $Q(I, I_p) = s-r$, which is indeed the signed distance between $I$ and $P$. Not sure if this generalizes to other combinations?                                                                                                                                                                                                                                                                             |
-| Distance between sphere/planes rep by I and J                                  | Note that for any two Euclidean-concentric spheres rep by $I$ and $J$ with radius $r$ and $s,$ $Q(I,J) = -\frac12\left(\frac rs  + \frac sr\right)$ depends only on the ratio of $r$ and $s$. So this can't give something that determines the Euclidean distance between the two spheres, which presumably grows as the two spheres are blown up proportionally. For another example, for any two parallel planes, $Q(I,J) = \pm1$. | Alex had said: Q(I,J)=cosh^2 (d/2) maybe where d is distance in usual hyperbolic metric. Or maybe cosh d. That may be right depending on what's meant by the hyperbolic metric there, but it seems like it won't determine a reasonable Euclidean distance between planes, which should differ between different pairs of parallel planes.                                                                                                                                                                                                                                               |
+| Distance between sphere/planes rep by I and J                                  | Note that for any two Euclidean-concentric spheres rep by $I$ and $J$ with radius $r$ and $s,$ $Q(I,J) = -\frac12\left(\frac rs  + \frac sr\right)$ depends only on the ratio of $r$ and $s$. So this can't give something that determines the Euclidean distance between the two spheres, which presumably grows as the two spheres are blown up proportionally. For another example, for any two parallel planes, $Q(I,J) = \pm1$. | Alex had said: $Q(I,J)=\cosh(d/2)^2$ maybe where d is distance in usual hyperbolic metric. Or maybe $\cosh(d)$. That may be right depending on what's meant by the hyperbolic metric there, but it seems like it won't determine a reasonable Euclidean distance between planes, which should differ between different pairs of parallel planes.                                                                                                                                                                                                                                         |
 | Sphere centered on P through R                                                 |                                                                                                                                                                                                                                                                                                                                                                                                                                      | Probably just calculate distance etc.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
-| Plane rep'd by I goes through center of sphere rep'd by J                      | I think this is equivalent to the plane being perpendicular to the sphere, i.e.$Q(I,J) = 0$.                                                                                                                                                                                                                                                                                                                                         |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
-| Dihedral angle between planes (or spheres?) rep by I and J                     | $\theta = \arccos(Q(I,J))$                                                                                                                                                                                                                                                                                                                                                                                                           | Aaron Fenyes points out: The angle between spheres in $S^3$ matches the angle between the planes they bound in $R^{(1,4)}$, which matches the angle between the spacelike vectors perpendicular to those planes. So we should have $Q(I,J) = \cos\theta$. Note that when the spheres do not intersect, we can interpret this as the "imaginary angle" between them, via $\cosh t = \cos it$.                                                                                                                                                                                             |
-| R, P, S are collinear                                                          | Maybe just cross product of two differences is 0. Or, $R,P,S,\infty$ lie on a circle, or equivalently, $I_R,I_P,I_S,I_\infty$ span a plane (rather than a three-space).                                                                                                                                                                                                                                                              | Not a conformal property, but $R,P,S,\infty$ lying on a circle _is_.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
-| Plane through noncollinear R, P, S                                             | Should be, just solve Q(I, I_R) = 0 etc.                                                                                                                                                                                                                                                                                                                                                                                             |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
+| Plane rep'd by I goes through center of sphere rep'd by J                      | I think this is equivalent to the plane being perpendicular to the sphere, i.e. $Q(I,J) = 0$.                                                                                                                                                                                                                                                                                                                                        |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
+| Dihedral angle between planes (or spheres?) rep by I and J                     | $\theta = \arccos(Q(I,J))$                                                                                                                                                                                                                                                                                                                                                                                                           | Aaron Fenyes points out: The angle between spheres in $S^3$ matches the angle between the planes they bound in $R^{(1,4)}$, which matches the angle between the spacelike vectors perpendicular to those planes. So we should have $Q(I,J) = \cos(\theta)$. Note that when the spheres do not intersect, we can interpret this as the "imaginary angle" between them, via $\cosh(t) = \cos(it)$.                                                                                                                                                                                         |
+| R, P, S are collinear                                                          | Maybe just cross product of two differences is 0. Or, $R,P,S,\infty$ lie on a circle, or equivalently, $I_R,I_P,I_S,I_\infty$ span a plane (rather than a three-space).                                                                                                                                                                                                                                                              | $R,P,S$ lying on a line isn't a conformal property, but $R,P,S,\infty$ lying on a circle is.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             |
+| Plane through noncollinear R, P, S                                             | Should be, just solve $Q(I, I_R) = 0$ etc.                                                                                                                                                                                                                                                                                                                                                                                           |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
 | circle                                                                         | Maybe concentric sphere and the containing plane? Note it is easy to constrain the relationship between those two: they must be perpendicular.                                                                                                                                                                                                                                                                                       | Defn: circle is intersection of two spheres. That does cover lines. But you lose the canonicalness                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
 | line                                                                           | Maybe two perpendicular containing planes? Maybe the plane perpendicular to the line and through origin, together with the point of the line on that plane? Or maybe just as a bag of collinear points?                                                                                                                                                                                                                              | The first is the limiting case of the possible circle rep, but it is not canonical. The second appears to be canonical, but I don't see a circle rep that corresponds to it.                                                                                                                                                                                                                                                                                                                                                                                                             |