doc: expand on representation of collinearity

In the process, clarify the signed distance from a point to a sphere and
add inversion across a sphere.

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/alg-test/Engine.Algebraic.jl

@@ -0,0 +1,203 @@
+module Algebraic
+
+export
+  codimension, dimension,
+  Construction, realize,
+  Element, Point, Sphere,
+  Relation, LiesOn, AlignsWithBy, mprod
+
+import Subscripts
+using LinearAlgebra
+using AbstractAlgebra
+using Groebner
+using ...HittingSet
+
+# --- commutative algebra ---
+
+# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate
+# polynomial rings when the coefficients are integers. we use Groebner to extend
+# support to rationals and to finite fields of prime order
+Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} =
+  Generic.Ideal{U}(base_ring(I), groebner(gens(I)))
+
+function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}}
+  leading = [exponent_vector(f, 1) for f in gens(I)]
+  targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading]
+  length(HittingSet.solve(HittingSetProblem(targets), maxdepth))
+end
+
+dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} =
+  length(gens(base_ring(I))) - codimension(I, maxdepth)
+
+# --- primitve elements ---
+
+abstract type Element{T} end
+
+mutable struct Point{T} <: Element{T}
+  coords::Vector{MPolyRingElem{T}}
+  vec::Union{Vector{MPolyRingElem{T}}, Nothing}
+  rel::Nothing
+  
+  ## [to do] constructor argument never needed?
+  Point{T}(
+    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
+    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
+  ) where T = new(coords, vec, nothing)
+end
+
+function buildvec!(pt::Point)
+  coordring = parent(pt.coords[1])
+  pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
+end
+
+mutable struct Sphere{T} <: Element{T}
+  coords::Vector{MPolyRingElem{T}}
+  vec::Union{Vector{MPolyRingElem{T}}, Nothing}
+  rel::Union{MPolyRingElem{T}, Nothing}
+  
+  ## [to do] constructor argument never needed?
+  Sphere{T}(
+    coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
+    vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
+    rel::Union{MPolyRingElem{T}, Nothing} = nothing
+  ) where T = new(coords, vec, rel)
+end
+
+function buildvec!(sph::Sphere)
+  coordring = parent(sph.coords[1])
+  sph.vec = sph.coords
+  sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
+end
+
+const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}(
+  nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ],
+  nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
+)
+
+coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index]
+
+function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex)
+  eltindex, elt = indexed_elt
+  name = coordname(elt, coordindex)
+  if !isnothing(name)
+    subscript = Subscripts.sub(string(eltindex))
+    push!(coordnamelist, Symbol(name, subscript))
+  end
+end
+
+function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex)
+  elt = indexed_elt[2]
+  if !isnothing(coordname(elt, coordindex))
+    push!(elt.coords, popfirst!(coordlist))
+  end
+end
+
+# --- primitive relations ---
+
+abstract type Relation{T} end
+
+mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end])
+
+# elements: point, sphere
+struct LiesOn{T} <: Relation{T}
+  elements::Vector{Element{T}}
+  
+  LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph])
+end
+
+equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec)
+
+# elements: sphere, sphere
+struct AlignsWithBy{T} <: Relation{T}
+  elements::Vector{Element{T}}
+  cos_angle::T
+  
+  AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle)
+end
+
+equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle
+
+# --- constructions ---
+
+mutable struct Construction{T}
+  points::Vector{Point{T}}
+  spheres::Vector{Sphere{T}}
+  relations::Vector{Relation{T}}
+  
+  function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{Relation{T}}()) where T
+    allelements = union(elements, (rel.elements for rel in relations)...)
+    new{T}(
+      filter(elt -> isa(elt, Point), allelements),
+      filter(elt -> isa(elt, Sphere), allelements),
+      relations
+    )
+  end
+end
+
+function Base.push!(ctx::Construction{T}, elt::Point{T}) where T
+  push!(ctx.points, elt)
+end
+
+function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T
+  push!(ctx.spheres, elt)
+end
+
+function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
+  push!(ctx.relations, rel)
+  for elt in rel.elements
+    push!(ctx, elt)
+  end
+end
+
+function realize(ctx::Construction{T}) where T
+  # collect coordinate names
+  coordnamelist = Symbol[]
+  eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points)))
+  for coordindex in 1:5
+    for indexed_elt in eltenum
+      pushcoordname!(coordnamelist, indexed_elt, coordindex)
+    end
+  end
+  
+  # construct coordinate ring
+  coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
+  
+  # retrieve coordinates
+  for (_, elt) in eltenum
+    empty!(elt.coords)
+  end
+  for coordindex in 1:5
+    for indexed_elt in eltenum
+      takecoord!(coordqueue, indexed_elt, coordindex)
+    end
+  end
+  
+  # construct coordinate vectors
+  for (_, elt) in eltenum
+    buildvec!(elt)
+  end
+  
+  # turn relations into equations
+  eqns = vcat(
+    equation.(ctx.relations),
+    [elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)]
+  )
+  
+  # add relations to center, orient, and scale the construction
+  # [to do] the scaling constraint, as written, can be impossible to satisfy
+  #         when all of the spheres have to go through the origin
+  if !isempty(ctx.points)
+    append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3])
+  end
+  if !isempty(ctx.spheres)
+    append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
+  end
+  n_elts = length(ctx.points) + length(ctx.spheres)
+  if n_elts > 0
+    push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
+  end
+  
+  (Generic.Ideal(coordring, eqns), eqns)
+end
+
+end

engine-proto/alg-test/Engine.Numerical.jl

@@ -0,0 +1,53 @@
+module Numerical
+
+using Random: default_rng
+using LinearAlgebra
+using AbstractAlgebra
+using HomotopyContinuation:
+  Variable, Expression, AbstractSystem, System, LinearSubspace,
+  nvariables, isreal, witness_set, results
+import GLMakie
+using ..Algebraic
+
+# --- polynomial conversion ---
+
+# hat tip Sascha Timme
+#   https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/issues/520#issuecomment-1317681521
+function Base.convert(::Type{Expression}, f::MPolyRingElem)
+  variables = Variable.(symbols(parent(f)))
+  f_data = zip(coefficients(f), exponent_vectors(f))
+  sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data)
+end
+
+# create a ModelKit.System from an ideal in a multivariate polynomial ring. the
+# variable ordering is taken from the polynomial ring
+function System(I::Generic.Ideal)
+  eqns = Expression.(gens(I))
+  variables = Variable.(symbols(base_ring(I)))
+  System(eqns, variables = variables)
+end
+
+# --- sampling ---
+
+function real_samples(F::AbstractSystem, dim; rng = default_rng())
+  # choose a random real hyperplane of codimension `dim` by intersecting
+  # hyperplanes whose normal vectors are uniformly distributed over the unit
+  # sphere
+  # [to do] guard against the unlikely event that one of the normals is zero
+  normals = transpose(hcat(
+    (normalize(randn(rng, nvariables(F))) for _ in 1:dim)...
+  ))
+  cut = LinearSubspace(normals, fill(0., dim))
+  filter(isreal, results(witness_set(F, cut, seed = 0x1974abba)))
+end
+
+AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
+  GLMakie.Point3f([evaluate(u, vals) for u in pt.coords])
+
+function AbstractAlgebra.evaluate(sph::Sphere, vals::Vector{<:RingElement})
+  radius = 1 / evaluate(sph.coords[1], vals)
+  center = radius * [evaluate(u, vals) for u in sph.coords[3:end]]
+  GLMakie.Sphere(GLMakie.Point3f(center), radius)
+end
+
+end

engine-proto/alg-test/Engine.jl

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

engine-proto/alg-test/HittingSet.jl

@@ -0,0 +1,111 @@
+module HittingSet
+
+export HittingSetProblem, solve
+
+HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}}
+
+# `targets` should be a collection of Set objects
+function HittingSetProblem(targets, chosen = Set())
+  wholeset = union(targets...)
+  T = eltype(wholeset)
+  unsorted_moves = [
+    elt => Set(filter(s -> elt ∉ s, targets))
+    for elt in wholeset
+  ]
+  moves = sort(unsorted_moves, by = pair -> length(pair.second))
+  Set{T}(chosen) => moves
+end
+
+function Base.display(problem::HittingSetProblem{T}) where T
+  println("HittingSetProblem{$T}")
+  
+  chosen = problem.first
+  println(" {", join(string.(chosen), ", "), "}")
+  
+  moves = problem.second
+  for (choice, missed) in moves
+    println(" | ", choice)
+    for s in missed
+      println(" | | {", join(string.(s), ", "), "}")
+    end
+  end
+  println()
+end
+
+function solve(pblm::HittingSetProblem{T}, maxdepth = Inf) where T
+  problems = Dict(pblm)
+  while length(first(problems).first) < maxdepth
+    subproblems = typeof(problems)()
+    for (chosen, moves) in problems
+      if isempty(moves)
+        return chosen
+      else
+        for (choice, missed) in moves
+          to_be_chosen = union(chosen, Set([choice]))
+          if isempty(missed)
+            return to_be_chosen
+          elseif !haskey(subproblems, to_be_chosen)
+            push!(subproblems, HittingSetProblem(missed, to_be_chosen))
+          end
+        end
+      end
+    end
+    problems = subproblems
+  end
+  problems
+end
+
+function test(n = 1)
+  T = [Int64, Int64, Symbol, Symbol][n]
+  targets = Set{T}.([
+    [
+      [1, 3, 5],
+      [2, 3, 4],
+      [1, 4],
+      [2, 3, 4, 5],
+      [4, 5]
+    ],
+    # example from Amit Chakrabarti's graduate-level algorithms class (CS 105)
+    # notes by Valika K. Wan and Khanh Do Ba, Winter 2005
+    # https://www.cs.dartmouth.edu/~ac/Teach/CS105-Winter05/
+    [
+      [1, 3], [1,    4], [1,    5],
+      [1, 3], [1, 2, 4], [1, 2, 5],
+      [4, 3], [   2, 4], [   2, 5],
+      [6, 3], [6,    4], [      5]
+    ],
+    [
+      [:w, :x, :y],
+      [:x, :y, :z],
+      [:w, :z],
+      [:x, :y]
+    ],
+    # Wikipedia showcases this as an example of a problem where the greedy
+    # algorithm performs especially poorly
+    [
+      [:a, :x, :t1],
+      [:a, :y, :t2],
+      [:a, :y, :t3],
+      [:a, :z, :t4],
+      [:a, :z, :t5],
+      [:a, :z, :t6],
+      [:a, :z, :t7],
+      [:b, :x, :t8],
+      [:b, :y, :t9],
+      [:b, :y, :t10],
+      [:b, :z, :t11],
+      [:b, :z, :t12],
+      [:b, :z, :t13],
+      [:b, :z, :t14]
+    ]
+  ][n])
+  problem = HittingSetProblem(targets)
+  if isa(problem, HittingSetProblem{T})
+    println("Correct type")
+  else
+    println("Wrong type: ", typeof(problem))
+  end
+  problem
+end
+
+end

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,450 @@
+module Engine
+
+using LinearAlgebra
+using GenericLinearAlgebra
+using SparseArrays
+using Random
+using Optim
+
+export
+  rand_on_shell, Q, DescentHistory,
+  realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram
+
+# === guessing ===
+
+sconh(t, u) = 0.5*(exp(t) + u*exp(-t))
+
+function rand_on_sphere(rng::AbstractRNG, ::Type{T}, n) where T
+  out = randn(rng, T, n)
+  tries_left = 2
+  while dot(out, out) < 1e-6 && tries_left > 0
+    out = randn(rng, T, n)
+    tries_left -= 1
+  end
+  normalize(out)
+end
+
+##[TO DO] write a test to confirm that the outputs are on the correct shells
+function rand_on_shell(rng::AbstractRNG, shell::T) where T <: Number
+  space_part = rand_on_sphere(rng, T, 4)
+  rapidity = randn(rng, T)
+  sig = sign(shell)
+  nullmix * [sconh(rapidity, sig)*space_part; sconh(rapidity, -sig)]
+end
+
+rand_on_shell(rng::AbstractRNG, shells::Array{T}) where T <: Number =
+  hcat([rand_on_shell(rng, sh) for sh in shells]...)
+
+rand_on_shell(shells::Array{<:Number}) = rand_on_shell(Random.default_rng(), shells)
+
+# === elements ===
+
+point(pos) = [pos; 0.5; 0.5 * dot(pos, pos)]
+
+plane(normal, offset) = [-normal; 0; -offset]
+
+function sphere(center, radius)
+  dist_sq = dot(center, center)
+  [
+    center / radius;
+    0.5 / radius;
+    0.5 * (dist_sq / radius - radius)
+  ]
+end
+
+# === Gram matrix realization ===
+
+# basis changes
+nullmix = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [-1 1; 1 1]//2]
+unmix = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [-1 1; 1 1]]
+
+# the Lorentz form
+Q = [Matrix{Int64}(I, 3, 3) zeros(Int64, 3, 2); zeros(Int64, 2, 3) [0 -2; -2 0]]
+
+# project a matrix onto the subspace of matrices whose entries vanish away from
+# the given indices
+function proj_to_entries(mat, indices)
+  result = zeros(size(mat))
+  for (j, k) in indices
+    result[j, k] = mat[j, k]
+  end
+  result
+end
+
+# the difference between the matrices `target` and `attempt`, projected onto the
+# subspace of matrices whose entries vanish at each empty index of `target`
+function proj_diff(target::SparseMatrixCSC{T, <:Any}, attempt::Matrix{T}) where T
+  J, K, values = findnz(target)
+  result = zeros(size(target))
+  for (j, k, val) in zip(J, K, values)
+    result[j, k] = val - attempt[j, k]
+  end
+  result
+end
+
+# a type for keeping track of gradient descent history
+struct DescentHistory{T}
+  scaled_loss::Array{T}
+  neg_grad::Array{Matrix{T}}
+  base_step::Array{Matrix{T}}
+  hess::Array{Hermitian{T, Matrix{T}}}
+  slope::Array{T}
+  stepsize::Array{T}
+  positive::Array{Bool}
+  backoff_steps::Array{Int64}
+  last_line_L::Array{Matrix{T}}
+  last_line_loss::Array{T}
+  
+  function DescentHistory{T}(
+    scaled_loss = Array{T}(undef, 0),
+    neg_grad = Array{Matrix{T}}(undef, 0),
+    hess = Array{Hermitian{T, Matrix{T}}}(undef, 0),
+    base_step = Array{Matrix{T}}(undef, 0),
+    slope = Array{T}(undef, 0),
+    stepsize = Array{T}(undef, 0),
+    positive = Bool[],
+    backoff_steps = Int64[],
+    last_line_L = Array{Matrix{T}}(undef, 0),
+    last_line_loss = Array{T}(undef, 0)
+  ) where T
+    new(scaled_loss, neg_grad, hess, base_step, slope, stepsize, positive, backoff_steps, last_line_L, last_line_loss)
+  end
+end
+
+# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
+# explicit entry of `gram`. use gradient descent starting from `guess`
+function realize_gram_gradient(
+  gram::SparseMatrixCSC{T, <:Any},
+  guess::Matrix{T};
+  scaled_tol = 1e-30,
+  min_efficiency = 0.5,
+  init_stepsize = 1.0,
+  backoff = 0.9,
+  max_descent_steps = 600,
+  max_backoff_steps = 110
+) where T <: Number
+  # start history
+  history = DescentHistory{T}()
+  
+  # scale tolerance
+  scale_adjustment = sqrt(T(nnz(gram)))
+  tol = scale_adjustment * scaled_tol
+  
+  # initialize variables
+  stepsize = init_stepsize
+  L = copy(guess)
+  
+  # do gradient descent
+  Δ_proj = proj_diff(gram, L'*Q*L)
+  loss = dot(Δ_proj, Δ_proj)
+  for _ in 1:max_descent_steps
+    # stop if the loss is tolerably low
+    if loss < tol
+      break
+    end
+    
+    # find negative gradient of loss function
+    neg_grad = 4*Q*L*Δ_proj
+    slope = norm(neg_grad)
+    dir = neg_grad / slope
+    
+    # store current position, loss, and slope
+    L_last = L
+    loss_last = loss
+    push!(history.scaled_loss, loss / scale_adjustment)
+    push!(history.neg_grad, neg_grad)
+    push!(history.slope, slope)
+    
+    # find a good step size using backtracking line search
+    push!(history.stepsize, 0)
+    push!(history.backoff_steps, max_backoff_steps)
+    empty!(history.last_line_L)
+    empty!(history.last_line_loss)
+    for backoff_steps in 0:max_backoff_steps
+      history.stepsize[end] = stepsize
+      L = L_last + stepsize * dir
+      Δ_proj = proj_diff(gram, L'*Q*L)
+      loss = dot(Δ_proj, Δ_proj)
+      improvement = loss_last - loss
+      push!(history.last_line_L, L)
+      push!(history.last_line_loss, loss / scale_adjustment)
+      if improvement >= min_efficiency * stepsize * slope
+        history.backoff_steps[end] = backoff_steps
+        break
+      end
+      stepsize *= backoff
+    end
+    
+    # [DEBUG] if we've hit a wall, quit
+    if history.backoff_steps[end] == max_backoff_steps
+      break
+    end
+  end
+  
+  # return the factorization and its history
+  push!(history.scaled_loss, loss / scale_adjustment)
+  L, history
+end
+
+function basis_matrix(::Type{T}, j, k, dims) where T
+  result = zeros(T, dims)
+  result[j, k] = one(T)
+  result
+end
+
+# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
+# explicit entry of `gram`. use Newton's method starting from `guess`
+function realize_gram_newton(
+  gram::SparseMatrixCSC{T, <:Any},
+  guess::Matrix{T};
+  scaled_tol = 1e-30,
+  rate = 1,
+  max_steps = 100
+) where T <: Number
+  # start history
+  history = DescentHistory{T}()
+  
+  # find the dimension of the search space
+  dims = size(guess)
+  element_dim, construction_dim = dims
+  total_dim = element_dim * construction_dim
+  
+  # list the constrained entries of the gram matrix
+  J, K, _ = findnz(gram)
+  constrained = zip(J, K)
+  
+  # scale the tolerance
+  scale_adjustment = sqrt(T(length(constrained)))
+  tol = scale_adjustment * scaled_tol
+  
+  # use Newton's method
+  L = copy(guess)
+  for step in 0:max_steps
+    # evaluate the loss function
+    Δ_proj = proj_diff(gram, L'*Q*L)
+    loss = dot(Δ_proj, Δ_proj)
+    
+    # store the current loss
+    push!(history.scaled_loss, loss / scale_adjustment)
+    
+    # stop if the loss is tolerably low
+    if loss < tol || step > max_steps
+      break
+    end
+    
+    # find the negative gradient of loss function
+    neg_grad = 4*Q*L*Δ_proj
+    
+    # find the negative Hessian of the loss function
+    hess = Matrix{T}(undef, total_dim, total_dim)
+    indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
+    for (j, k) in indices
+      basis_mat = basis_matrix(T, j, k, dims)
+      neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
+      neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
+      deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
+      hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
+    end
+    hess = Hermitian(hess)
+    push!(history.hess, hess)
+    
+    # compute the Newton step
+    step = hess \ reshape(neg_grad, total_dim)
+    L += rate * reshape(step, dims)
+  end
+  
+  # return the factorization and its history
+  L, history
+end
+
+LinearAlgebra.eigen!(A::Symmetric{BigFloat, Matrix{BigFloat}}; sortby::Nothing) =
+  eigen!(Hermitian(A))
+
+function convertnz(type, mat)
+  J, K, values = findnz(mat)
+  sparse(J, K, type.(values))
+end
+
+function realize_gram_optim(
+  gram::SparseMatrixCSC{T, <:Any},
+  guess::Matrix{T}
+) where T <: Number
+  # find the dimension of the search space
+  dims = size(guess)
+  element_dim, construction_dim = dims
+  total_dim = element_dim * construction_dim
+  
+  # list the constrained entries of the gram matrix
+  J, K, _ = findnz(gram)
+  constrained = zip(J, K)
+  
+  # scale the loss function
+  scale_adjustment = length(constrained)
+  
+  function loss(L_vec)
+    L = reshape(L_vec, dims)
+    Δ_proj = proj_diff(gram, L'*Q*L)
+    dot(Δ_proj, Δ_proj) / scale_adjustment
+  end
+  
+  function loss_grad!(storage, L_vec)
+    L = reshape(L_vec, dims)
+    Δ_proj = proj_diff(gram, L'*Q*L)
+    storage .= reshape(-4*Q*L*Δ_proj, total_dim) / scale_adjustment
+  end
+  
+  function loss_hess!(storage, L_vec)
+    L = reshape(L_vec, dims)
+    Δ_proj = proj_diff(gram, L'*Q*L)
+    indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
+    for (j, k) in indices
+      basis_mat = basis_matrix(T, j, k, dims)
+      neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
+      neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
+      deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj) / scale_adjustment
+      storage[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
+    end
+  end
+  
+  optimize(
+    loss, loss_grad!, loss_hess!,
+    reshape(guess, total_dim),
+    Newton()
+  )
+end
+
+# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
+# explicit entry of `gram`. use gradient descent starting from `guess`
+function realize_gram(
+  gram::SparseMatrixCSC{T, <:Any},
+  guess::Matrix{T},
+  frozen = nothing;
+  scaled_tol = 1e-30,
+  min_efficiency = 0.5,
+  init_rate = 1.0,
+  backoff = 0.9,
+  reg_scale = 1.1,
+  max_descent_steps = 200,
+  max_backoff_steps = 110
+) where T <: Number
+  # start history
+  history = DescentHistory{T}()
+  
+  # find the dimension of the search space
+  dims = size(guess)
+  element_dim, construction_dim = dims
+  total_dim = element_dim * construction_dim
+  
+  # list the constrained entries of the gram matrix
+  J, K, _ = findnz(gram)
+  constrained = zip(J, K)
+  
+  # scale the tolerance
+  scale_adjustment = sqrt(T(length(constrained)))
+  tol = scale_adjustment * scaled_tol
+  
+  # list the un-frozen indices
+  has_frozen = !isnothing(frozen)
+  if has_frozen
+    is_unfrozen = fill(true, size(guess))
+    is_unfrozen[frozen] .= false
+    unfrozen = findall(is_unfrozen)
+    unfrozen_stacked = reshape(is_unfrozen, total_dim)
+  end
+  
+  # initialize variables
+  grad_rate = init_rate
+  L = copy(guess)
+  
+  # use Newton's method with backtracking and gradient descent backup
+  Δ_proj = proj_diff(gram, L'*Q*L)
+  loss = dot(Δ_proj, Δ_proj)
+  for step in 1:max_descent_steps
+    # stop if the loss is tolerably low
+    if loss < tol
+      break
+    end
+    
+    # find the negative gradient of loss function
+    neg_grad = 4*Q*L*Δ_proj
+    
+    # find the negative Hessian of the loss function
+    hess = Matrix{T}(undef, total_dim, total_dim)
+    indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
+    for (j, k) in indices
+      basis_mat = basis_matrix(T, j, k, dims)
+      neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
+      neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
+      deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
+      hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
+    end
+    hess = Hermitian(hess)
+    push!(history.hess, hess)
+    
+    # regularize the Hessian
+    min_eigval = minimum(eigvals(hess))
+    push!(history.positive, min_eigval > 0)
+    if min_eigval <= 0
+      hess -= reg_scale * min_eigval * I
+    end
+    
+    # compute the Newton step
+    neg_grad_stacked = reshape(neg_grad, total_dim)
+    if has_frozen
+      hess = hess[unfrozen_stacked, unfrozen_stacked]
+      neg_grad_compressed = neg_grad_stacked[unfrozen_stacked]
+    else
+      neg_grad_compressed = neg_grad_stacked
+    end
+    base_step_compressed = hess \ neg_grad_compressed
+    if has_frozen
+      base_step_stacked = zeros(total_dim)
+      base_step_stacked[unfrozen_stacked] .= base_step_compressed
+    else
+      base_step_stacked = base_step_compressed
+    end
+    base_step = reshape(base_step_stacked, dims)
+    push!(history.base_step, base_step)
+    
+    # store the current position, loss, and slope
+    L_last = L
+    loss_last = loss
+    push!(history.scaled_loss, loss / scale_adjustment)
+    push!(history.neg_grad, neg_grad)
+    push!(history.slope, norm(neg_grad))
+    
+    # find a good step size using backtracking line search
+    push!(history.stepsize, 0)
+    push!(history.backoff_steps, max_backoff_steps)
+    empty!(history.last_line_L)
+    empty!(history.last_line_loss)
+    rate = one(T)
+    step_success = false
+    for backoff_steps in 0:max_backoff_steps
+      history.stepsize[end] = rate
+      L = L_last + rate * base_step
+      Δ_proj = proj_diff(gram, L'*Q*L)
+      loss = dot(Δ_proj, Δ_proj)
+      improvement = loss_last - loss
+      push!(history.last_line_L, L)
+      push!(history.last_line_loss, loss / scale_adjustment)
+      if improvement >= min_efficiency * rate * dot(neg_grad, base_step)
+        history.backoff_steps[end] = backoff_steps
+        step_success = true
+        break
+      end
+      rate *= backoff
+    end
+    
+    # if we've hit a wall, quit
+    if !step_success
+      return L_last, false, history
+    end
+  end
+  
+  # return the factorization and its history
+  push!(history.scaled_loss, loss / scale_adjustment)
+  L, loss < tol, history
+end
+
+end

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

@@ -2,28 +2,29 @@
 
 (proposed by Alex Kontorovich as a practical system for doing 3D geometric calculations)
 
-These coordinates are of form $I=(c, r, x, y, z)$ where we think of $c$ as the co-radius, $r$ as the radius, and $x, y, z$ as the "Euclidean" part, which we abbreviate $E_I$.  There is an underlying basic quadratic form $Q(I_1,I_2) = (c_1r_2+c_2r_1)/2 - x_1x_2 -y_1y_2-z_1z_2$ which aids in calculation/verification of coordinates in this representation. We have:
+These coordinates are of form $I=(c, b, x, y, z)$ where we think of $c$ as the co-radius, $b$ as the "bend" (reciprocal radius), and $x, y, z$ as the "Euclidean" part, which we abbreviate $E_I$.  There is an underlying basic quadratic form $Q(I_1,I_2) = (c_1b_2+c_2b_1)/2 - x_1x_2 -y_1y_2-z_1z_2$ which aids in calculation/verification of coordinates in this representation. We have:
 
-| 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)$                                                                                                                                                                                                                                                                                                                                                                                                                                         |
-| 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 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.                                                                                                                                                                                                                                               |
-| 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.                                                                                                                                                                                                                                                                                                                                                                                             |                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
-| 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.                                                                                                                                                                                                                                                                                                                                                                                                             |
+| Entity or Relationship                                                             | Representation                                                                                                                                                                                                                                                                                                                                                                                                                       | Comments/questions                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
+| ---------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
+| Sphere $s$ with radius $r>0$ centered on $P = (x,y,z)$                             | $I_s = (\frac1{c}, \frac1{r}, \frac{x}{r}, \frac{y}{r}, \frac{z}{r})$ satisfying $Q(I_s,I_s) = -1,$ i.e., $c = r/(\|P\|^2 - r^2)$.                                                                                                                                                                                                                                                                                                   | Note that $1/c = \|P\|^2/r - r$, so there is no trouble if $\|P\| = r$; we just get first coordinate to be 0. Using the point representation $I_P$ from below, let's orient the sphere so that its normals point into the "positive side," where $Q(I_P, I_s) > 0$. The vector $I_s$ then represents a sphere with outward normals, while $-I_s$ represents one with inward normals.                                                                                                                                                                                                                                                                                                                                       |
+| Plane $p$ with unit normal $(x,y,z)$ through the (Euclidean) point $(sx,sy,sz)$    | $I_p = (-2s, 0, -x, -y, -z)$                                                                                                                                                                                                                                                                                                                                                                                                         | Note that $Q(I_p, I_p)$ is still $-1$. We orient planes using the same convention we use for spheres. For example, $(-2, 0, -1/\sqrt3, -1/\sqrt3, -1/\sqrt3)$ and $(2, 0, 1/\sqrt3, 1/\sqrt3, 1/\sqrt3)$ represent planes that coincide in space, which have their normals pointing away from and toward the origin, respectively. Note that the ray from $(sx, sy, sz) \in p$ in direction $(-x, -y, -z)$ is the ray perpendicular to the plane through the origin; since $(-x, -y, -z)$ is a unit vector, $(sx, sy, sz)$ and hence $p$ is at distance $s$ from the origin. These coordinates are essentially the limit of a sphere's coordinates as its radius goes to infinity, or equivalently, as its bend goes to 0. |
+| Point $P$ with Euclidean coordinates $(x,y,z)$                                     | $I_P = (\|P\|^2, 1, x, y, z)$                                                                                                                                                                                                                                                                                                                                                                                                        | Note $Q(I_P,I_P) = 0$. This gives us the freedom to choose a different normalization. For example, we could scale the representation shown here by $(\|P\|^2+1)^{-1}$, putting it on the sphere where the light cone intersects the plane where the first two coordinates sum to $1$.                                                                                                                                                                                                                                                                                                                                                                                                                                      |
+| ∞, the "point at infinity"                                                         | $I_\infty = (1,0,0,0,0)$                                                                                                                                                                                                                                                                                                                                                                                                             | The only solution to $Q(I,I) = 0$ not covered by (some normalization of) the above case.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   |
+| Point $P$ lies on sphere or plane given by $I$                                     | $Q(I_P, I) = 0$                                                                                                                                                                                                                                                                                                                                                                                                                      | Actually also works if $I$ is the coordinates of a point, in which case "lies on" simply means "coincides with".                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           |
+| 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) | $\lvert Q(I,J)\rvert \le (\text{resp. }>)\; 1$                                                                                                                                                                                                                                                                                                                                                                                       | Follows from the angle formula and the tangency condition, at least conceptually. One subtlety: parallel planes have $Q$ - product $\pm 1$, because they intersect at infinity (and in fact, are "tangent" there)!                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           |
+| $P$ is center of sphere rep'd by $I$                                               | $Q(I, I_P) = -r/2$, where $1/r = 2Q(I_\infty, I)$ is the signed bend of the sphere, and $I_P$ is normalized in the standard way, which is to set $Q(I_\infty, I_P) = 1/2$                                                                                                                                                                                                                                                            | This relationship is equivalent to both of the following. (1) The point $P$ has signed distance $-r$ from the sphere. (2) Inversion across the sphere maps $\infty$ to $P$.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                |
+| Distance between points $P$ and $R$ is $d$                                         | $Q(I_P, I_R) = d^2/2$                                                                                                                                                                                                                                                                                                                                                                                                                | If $P$ and $R$ are represented by non-normalized vectors $V_P$ and $V_R$, the relation becomes $Q(V_P, V_R) = 2\,Q(V_P, I_\infty)\,Q(V_R, I_\infty)\,d^2$. This version of the relation makes it easier to see why $d$ goes to infinity as $P$ or $R$ approaches the point at infinity.                                                                                                                                                                                                                                                                                                                                                                                                                                    |
+| Signed distance between point rep'd by $V$ and sphere/plane rep'd by $I$ is $d$    | In general, $\frac{Q(I, V)}{2Q(I_\infty, V)} = Q(I_\infty, I)\,d^2 + d$. When $V$ is normalized in the usual way, this simplifies to $Q(I, V) = d^2/r + d$ for a sphere of radius $r$, and to $Q(I, V) = d$ for a plane.                                                                                                                                                                                                             | We can use a Euclidean motion, represented linearly by a Lorentz transformation that fixes $I_\infty$, to put the point on the $z$ axis and put the nearest point on the sphere/plane at the origin with its normal pointing in the positive $z$ direction. Then the sphere/plane is represented by $I = (0, 1/r, 0, 0, -1)$, and the point can be represented by any multiple of $I_P = (d^2, 1, 0, 0, d)$, giving $Q(I, I_P) = d^2/2r + d.$ We turn this into a general expression by writing it in terms of Lorentz-invariant quantities and making it independent of the normalization of $I_P$.                                                                                                                       |
+| 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 point $P$ through point $R$                                     |                                                                                                                                                                                                                                                                                                                                                                                                                                      | Probably just calculate distance etc.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
+| Plane rep'd by $I$ goes through center of sphere rep'd by $J$                      | This is equivalent to the plane being perpendicular to the sphere: that is, $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)$.                                                                                                                                                                                                                                                                                                                           |
+| Points $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). Or we can add two planes constrained to be perpendicular with one constrained to contain the origin, and all three points constrained to lie on both. But that's a lot of auxiliary entities and constraints...                                                                                                                                                                                                                                                             | $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. However, there is a distinguished "standard" choice we could make: always choose one plane to contain the origin and the line, and the other to be the perpendicular plane containing the line. That choice or Plücker coordinates might be the best we can do. If we use the standardized perpendicular planes choice, then adding a line would be equivalent to adding two planes and the two constraints that one contains the origin and the other is perpendicular to it. That doesn't seem so bad. The second convention (perpendicular plane through the origin and a point on it) appears to be canonical, but there doesn't seem to be a circle representation that tends to it in the limit.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
+| Inversion of entity represented by $v$ across sphere $s$, rep'd by $I_s$           | $v \mapsto v + 2Q(I_s, v)\,I_s$                                                                                                                                                                                                                                                                                                                                                                                                      | This is just an educated guess, but its behavior is consistent with inversion in at least two ways. (1) It fixes points on $s$ and spheres perpendicular to $s$. (2) It preserves dihedral angles with $s$.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                |
 
 The unification of spheres/planes is indeed attractive for a project like Dyna3. The relationship between this representation and Geometric Algebras is a bit murky; likely it somehow fits under the Geometric Algebra umbrella.