Tidy up Rust .gitignore files

Also, drop the commented-out beginnings of a `plotters-gtk4` version. I
can't use `plotters-gtk4` on my machine because it requires GTK 4.14 or
higher, and Ubuntu 22.04 is still at GTK 4.6.

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

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

@@ -0,0 +1,99 @@
+include("Engine.jl")
+
+using LinearAlgebra
+using SparseArrays
+
+function sphere_in_tetrahedron_shape()
+  # initialize the partial gram matrix for a sphere inscribed in a regular
+  # tetrahedron
+  J = Int64[]
+  K = Int64[]
+  values = BigFloat[]
+  for j in 1:5
+    for k in 1:5
+      push!(J, j)
+      push!(K, k)
+      if j == k
+        push!(values, 1)
+      elseif (j <= 4 && k <= 4)
+        push!(values, -1/BigFloat(3))
+      else
+        push!(values, -1)
+      end
+    end
+  end
+  gram = sparse(J, K, values)
+  
+  # plot loss along a slice
+  loss_lin = []
+  loss_sq = []
+  mesh = range(0.9, 1.1, 101)
+  for t in mesh
+    L = hcat(
+      Engine.plane(normalize(BigFloat[ 1,  1,  1]), BigFloat(1)),
+      Engine.plane(normalize(BigFloat[ 1, -1, -1]), BigFloat(1)),
+      Engine.plane(normalize(BigFloat[-1,  1, -1]), BigFloat(1)),
+      Engine.plane(normalize(BigFloat[-1, -1,  1]), BigFloat(1)),
+      Engine.sphere(BigFloat[0, 0, 0], BigFloat(t))
+    )
+    Δ_proj = Engine.proj_diff(gram, L'*Engine.Q*L)
+    push!(loss_lin, norm(Δ_proj))
+    push!(loss_sq, dot(Δ_proj, Δ_proj))
+  end
+  mesh, loss_lin, loss_sq
+end
+
+function circles_in_triangle_shape()
+  # initialize the partial gram matrix for a sphere inscribed in a regular
+  # tetrahedron
+  J = Int64[]
+  K = Int64[]
+  values = BigFloat[]
+  for j in 1:8
+    for k in 1:8
+      filled = false
+      if j == k
+        push!(values, 1)
+        filled = true
+      elseif (j == 1 || k == 1)
+        push!(values, 0)
+        filled = true
+      elseif (j == 2 || k == 2)
+        push!(values, -1)
+        filled = true
+      end
+      #=elseif (j <= 5 && j != 2 && k == 9 || k == 9 && k <= 5 && k != 2)
+        push!(values, 0)
+        filled = true
+      end=#
+      if filled
+        push!(J, j)
+        push!(K, k)
+      end
+    end
+  end
+  append!(J, [6, 4, 6, 5, 7, 5, 7, 3, 8, 3, 8, 4])
+  append!(K, [4, 6, 5, 6, 5, 7, 3, 7, 3, 8, 4, 8])
+  append!(values, fill(-1, 12))
+  
+  # plot loss along a slice
+  loss_lin = []
+  loss_sq = []
+  mesh = range(0.99, 1.01, 101)
+  for t in mesh
+    L = hcat(
+      Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
+      Engine.sphere(BigFloat[0, 0, 0], BigFloat(t)),
+      Engine.plane(BigFloat[1, 0, 0], BigFloat(1)),
+      Engine.plane(BigFloat[cos(2pi/3), sin(2pi/3), 0], BigFloat(1)),
+      Engine.plane(BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(1)),
+      Engine.sphere(4//3*BigFloat[-1, 0, 0], BigFloat(1//3)),
+      Engine.sphere(4//3*BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//3)),
+      Engine.sphere(4//3*BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//3))
+    )
+    Δ_proj = Engine.proj_diff(gram, L'*Engine.Q*L)
+    push!(loss_lin, norm(Δ_proj))
+    push!(loss_sq, dot(Δ_proj, Δ_proj))
+  end
+  mesh, loss_lin, loss_sq
+end

engine-proto/gram-test/circles-in-triangle.jl

@@ -0,0 +1,76 @@
+include("Engine.jl")
+
+using SparseArrays
+using Random
+
+# initialize the partial gram matrix for a sphere inscribed in a regular
+# tetrahedron
+J = Int64[]
+K = Int64[]
+values = BigFloat[]
+for j in 1:9
+  for k in 1:9
+    filled = false
+    if j == 9
+      if k <= 5 && k != 2
+        push!(values, 0)
+        filled = true
+      end
+    elseif k == 9
+      if j <= 5 && j != 2
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == k
+      push!(values, 1)
+      filled = true
+    elseif j == 1 || k == 1
+      push!(values, 0)
+      filled = true
+    elseif j == 2 || k == 2
+      push!(values, -1)
+      filled = true
+    end
+    if filled
+      push!(J, j)
+      push!(K, k)
+    end
+  end
+end
+append!(J, [6, 4, 6, 5, 7, 5, 7, 3, 8, 3, 8, 4])
+append!(K, [4, 6, 5, 6, 5, 7, 3, 7, 3, 8, 4, 8])
+append!(values, fill(-1, 12))
+#= make construction rigid
+append!(J, [3, 4, 4, 5])
+append!(K, [4, 3, 5, 4])
+append!(values, fill(-0.5, 4))
+=#
+gram = sparse(J, K, values)
+
+# set initial guess
+Random.seed!(58271)
+guess = hcat(
+  Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
+  Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  Engine.plane(-BigFloat[1, 0, 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  Engine.plane(-BigFloat[cos(2pi/3), sin(2pi/3), 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  Engine.plane(-BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(-1)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  Engine.sphere(BigFloat[-1, 0, 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  Engine.sphere(BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5)) + 0.1*Engine.rand_on_shell([BigFloat(-1)]),
+  BigFloat[0, 0, 0, 0, 1]
+)
+frozen = [CartesianIndex(j, 9) for j in 1:5]
+
+# complete the gram matrix using Newton's method with backtracking
+L, success, history = Engine.realize_gram(gram, guess, frozen)
+completed_gram = L'*Engine.Q*L
+println("Completed Gram matrix:\n")
+display(completed_gram)
+if success
+  println("\nTarget accuracy achieved!")
+else
+  println("\nFailed to reach target accuracy")
+end
+println("Steps: ", size(history.scaled_loss, 1))
+println("Loss: ", history.scaled_loss[end], "\n")

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

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

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

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

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

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

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

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

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

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

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

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

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

@@ -0,0 +1,67 @@
+include("Engine.jl")
+
+using SparseArrays
+using Random
+
+# initialize the partial gram matrix for a sphere inscribed in a regular
+# tetrahedron
+J = Int64[]
+K = Int64[]
+values = BigFloat[]
+for j in 1:6
+  for k in 1:6
+    filled = false
+    if j == 6
+      if k <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif k == 6
+      if j <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == k
+      push!(values, 1)
+      filled = true
+    elseif j <= 4 && k <= 4
+      push!(values, -1/BigFloat(3))
+      filled = true
+    else
+      push!(values, -1)
+      filled = true
+    end
+    if filled
+      push!(J, j)
+      push!(K, k)
+    end
+  end
+end
+gram = sparse(J, K, values)
+
+# set initial guess
+Random.seed!(99230)
+guess = hcat(
+  sqrt(1/BigFloat(3)) * BigFloat[
+    1  1 -1 -1  0
+    1 -1  1 -1  0
+    1 -1 -1  1  0
+    0  0  0  0  1.5
+    1  1  1  1 -0.5
+  ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
+  BigFloat[0, 0, 0, 0, 1]
+)
+frozen = [CartesianIndex(j, 6) for j in 1:5]
+
+# complete the gram matrix using Newton's method with backtracking
+L, success, history = Engine.realize_gram(gram, guess, frozen)
+completed_gram = L'*Engine.Q*L
+println("Completed Gram matrix:\n")
+display(completed_gram)
+if success
+  println("\nTarget accuracy achieved!")
+else
+  println("\nFailed to reach target accuracy")
+end
+println("Steps: ", size(history.scaled_loss, 1))
+println("Loss: ", history.scaled_loss[end], "\n")

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

@@ -0,0 +1,96 @@
+include("Engine.jl")
+
+using LinearAlgebra
+using SparseArrays
+using Random
+
+# initialize the partial gram matrix for a sphere inscribed in a regular
+# tetrahedron
+J = Int64[]
+K = Int64[]
+values = BigFloat[]
+for j in 1:11
+  for k in 1:11
+    filled = false
+    if j == 11
+      if k <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif k == 11
+      if j <= 4
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == k
+      push!(values, j <= 6 ? 1 : 0)
+      filled = true
+    elseif j <= 4
+      if k <= 4
+        push!(values, -1/BigFloat(3))
+        filled = true
+      elseif k == 5
+        push!(values, -1)
+        filled = true
+      elseif 7 <= k <= 10 && k - j != 6
+        push!(values, 0)
+        filled = true
+      end
+    elseif k <= 4
+      if j == 5
+        push!(values, -1)
+        filled = true
+      elseif 7 <= j <= 10 && j - k != 6
+        push!(values, 0)
+        filled = true
+      end
+    elseif j == 6 && 7 <= k <= 10 || k == 6 && 7 <= j <= 10
+      push!(values, 0)
+      filled = true
+    end
+    if filled
+      push!(J, j)
+      push!(K, k)
+    end
+  end
+end
+gram = sparse(J, K, values)
+
+# set initial guess
+Random.seed!(99230)
+guess = hcat(
+  sqrt(1/BigFloat(3)) * BigFloat[
+    1  1 -1 -1  0    0
+    1 -1  1 -1  0    0
+    1 -1 -1  1  0    0
+    0  0  0  0  1.5  0.5
+    1  1  1  1 -0.5 -1.5
+  ] + 0.0*Engine.rand_on_shell(fill(BigFloat(-1), 6)),
+  Engine.point([-0.5, -0.5, -0.5] + 0.3*randn(3)),
+  Engine.point([-0.5,  0.5,  0.5] + 0.3*randn(3)),
+  Engine.point([ 0.5, -0.5,  0.5] + 0.3*randn(3)),
+  Engine.point([ 0.5,  0.5, -0.5] + 0.3*randn(3)),
+  BigFloat[0, 0, 0, 0, 1]
+)
+frozen = vcat(
+  [CartesianIndex(4, k) for k in 7:10],
+  [CartesianIndex(j, 11) for j in 1:5]
+)
+
+# complete the gram matrix using Newton's method with backtracking
+L, success, history = Engine.realize_gram(gram, guess, frozen)
+completed_gram = L'*Engine.Q*L
+println("Completed Gram matrix:\n")
+display(completed_gram)
+if success
+  println("\nTarget accuracy achieved!")
+else
+  println("\nFailed to reach target accuracy")
+end
+println("Steps: ", size(history.scaled_loss, 1))
+println("Loss: ", history.scaled_loss[end])
+if success
+  infty = BigFloat[0, 0, 0, 0, 1]
+  radius_ratio = dot(infty, Engine.Q * L[:,5]) / dot(infty, Engine.Q * L[:,6])
+  println("\nCircumradius / inradius: ", radius_ratio)
+end

lang-trials/rust-benchmark-native/.gitignore

@@ -0,0 +1,3 @@
+target
+dist
+Cargo.lock

lang-trials/rust-benchmark-native/Cargo.toml

@@ -0,0 +1,16 @@
+[package]
+name = "rust-benchmark-native"
+version = "0.1.0"
+authors = ["Aaron"]
+edition = "2021"
+
+[dependencies]
+cairo-rs = "0.20.1"
+gtk = { package = "gtk4", version = "0.9.0" }
+nalgebra = "0.33.0"
+plotters = "0.3.6"
+plotters-cairo = "0.7.0"
+
+[profile.release]
+opt-level = "s" # optimize for small code size
+debug = true # include debug symbols

lang-trials/rust-benchmark-native/src/engine.rs

@@ -0,0 +1,105 @@
+use nalgebra::{*, allocator::Allocator};
+use std::f64::consts::{PI, E};
+
+/* dynamic matrices */
+pub fn rand_eigval_series(dim: usize, time_res: usize) -> Vec<OVector<Complex<f64>, Dyn>> {
+    // initialize the random matrix
+    let mut rand_mat = DMatrix::<f64>::from_fn(dim, dim, |j, k| {
+        let n = j*dim + k;
+        E*((n*n) as f64) % 2.0 - 1.0
+    }) * (3.0 / (dim as f64)).sqrt();
+    
+    // initialize the rotation step
+    let mut rot_step = DMatrix::<f64>::identity(dim, dim);
+    let max_freq = 4;
+    for n in (0..dim).step_by(2) {
+        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
+        let ang_cos = ang.cos();
+        let ang_sin = ang.sin();
+        rot_step[(n, n)] = ang_cos;
+        rot_step[(n+1, n)] = ang_sin;
+        rot_step[(n, n+1)] = -ang_sin;
+        rot_step[(n+1, n+1)] = ang_cos;
+    }
+    
+    // find the eigenvalues
+    let mut eigval_series = Vec::<OVector<Complex<f64>, Dyn>>::with_capacity(time_res);
+    eigval_series.push(rand_mat.complex_eigenvalues());
+    for _ in 1..time_res {
+        rand_mat = &rot_step * rand_mat;
+        eigval_series.push(rand_mat.complex_eigenvalues());
+    }
+    eigval_series
+}
+
+/* dynamic single float matrices */
+/*pub fn rand_eigval_series(dim: usize, time_res: usize) -> Vec<OVector<Complex<f32>, Dyn>> {
+    // initialize the random matrix
+    let mut rand_mat = DMatrix::<f32>::from_fn(dim, dim, |j, k| {
+        let n = j*dim + k;
+        (E as f32)*((n*n) as f32) % 2.0_f32 - 1.0_f32
+    }) * (3.0_f32 / (dim as f32)).sqrt();
+    
+    // initialize the rotation step
+    let mut rot_step = DMatrix::<f32>::identity(dim, dim);
+    let max_freq = 4;
+    for n in (0..dim).step_by(2) {
+        let ang = (PI as f32) * ((n % max_freq) as f32) / (time_res as f32);
+        let ang_cos = ang.cos();
+        let ang_sin = ang.sin();
+        rot_step[(n, n)] = ang_cos;
+        rot_step[(n+1, n)] = ang_sin;
+        rot_step[(n, n+1)] = -ang_sin;
+        rot_step[(n+1, n+1)] = ang_cos;
+    }
+    
+    // find the eigenvalues
+    let mut eigval_series = Vec::<OVector<Complex<f32>, Dyn>>::with_capacity(time_res);
+    eigval_series.push(rand_mat.complex_eigenvalues());
+    for _ in 1..time_res {
+        rand_mat = &rot_step * rand_mat;
+        eigval_series.push(rand_mat.complex_eigenvalues());
+    }
+    eigval_series
+}*/
+
+/* static matrices. should only be used when the dimension is really small */
+/*pub fn rand_eigval_series<N>(time_res: usize) -> Vec<OVector<Complex<f64>, N>>
+    where
+        N: ToTypenum + DimName + DimSub<U1>,
+        DefaultAllocator:
+            Allocator<N> +
+            Allocator<N, N> +
+            Allocator<<N as DimSub<U1>>::Output> +
+            Allocator<N, <N as DimSub<U1>>::Output>
+{
+    // initialize the random matrix
+    let dim = N::try_to_usize().unwrap();
+    let mut rand_mat = OMatrix::<f64, N, N>::from_fn(|j, k| {
+        let n = j*dim + k;
+        E*((n*n) as f64) % 2.0 - 1.0
+    }) * (3.0 / (dim as f64)).sqrt();
+    /*let mut rand_mat = OMatrix::<f64, N, N>::identity();*/
+    
+    // initialize the rotation step
+    let mut rot_step = OMatrix::<f64, N, N>::identity();
+    let max_freq = 4;
+    for n in (0..dim).step_by(2) {
+        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
+        let ang_cos = ang.cos();
+        let ang_sin = ang.sin();
+        rot_step[(n, n)] = ang_cos;
+        rot_step[(n+1, n)] = ang_sin;
+        rot_step[(n, n+1)] = -ang_sin;
+        rot_step[(n+1, n+1)] = ang_cos;
+    }
+    
+    // find the eigenvalues
+    let mut eigval_series = Vec::<OVector<Complex<f64>, N>>::with_capacity(time_res);
+    eigval_series.push(rand_mat.complex_eigenvalues());
+    for _ in 1..time_res {
+        rand_mat = &rot_step * rand_mat;
+        eigval_series.push(rand_mat.complex_eigenvalues());
+    }
+    eigval_series
+}*/

lang-trials/rust-benchmark-native/src/main.rs

@@ -0,0 +1,104 @@
+// based on Olivier Pelhatre's GTK 3 example, ported to GTK 4
+//
+//   https://github.com/Ouam74/RUST_Real-time_plots_using_GTK-rs_and_Plotters-rs
+//
+// a self-contained component might draw on the example below, by StackOverflow
+// user Nicolas
+//
+//   https://stackoverflow.com/a/76548487
+//
+// here's a crash course in `plotters`
+//
+//   https://plotters-rs.github.io/book/basic/basic_data_plotting.html
+//
+
+extern crate cairo;
+use plotters::prelude::*;
+use plotters_cairo::CairoBackend;
+use gtk::{
+    glib,
+    prelude::*,
+    Adjustment,
+    Align,
+    Application,
+    ApplicationWindow,
+    Box,
+    DrawingArea,
+    Label,
+    Orientation,
+    Scale
+};
+use std::time::Instant;
+
+mod engine;
+
+fn main() -> glib::ExitCode {
+    let app = Application::builder()
+        .application_id("org.studioinfinity.rust-benchmark-native")
+        .build();
+    
+    app.connect_activate(|app| {
+        const TIME_RES: usize = 100;
+        let start_time = Instant::now();
+        let eigval_series = engine::rand_eigval_series(60, TIME_RES);
+        let run_time = start_time.elapsed().as_millis();
+        
+        // application state
+        let time_step = Adjustment::new(0.0, 0.0, TIME_RES as f64, 1.0, 0.0, 0.0);
+        
+        // create the window.
+        let window = ApplicationWindow::builder()
+            .application(app)
+            .title("The circular law")
+            .build();
+        
+        // create a vertical box
+        let container = Box::new(Orientation::Vertical, 5);
+        window.set_child(Some(&container));
+        
+        // create the run time readout
+        let run_time_readout = Label::builder()
+            .margin_top(5)
+            .margin_start(10)
+            .halign(Align::Start)
+            .label(glib::gformat!("{} ms", run_time))
+            .build();
+        container.append(&run_time_readout);
+        
+        // set up the drawing area
+        let drawing_area = DrawingArea::builder()
+            .content_width(600)
+            .content_height(600)
+            .build();
+        let time_step_for_draw = time_step.clone();
+        let draw_eigvals = move |_: &DrawingArea, context: &cairo::Context, width: i32, height: i32| {
+            let root = CairoBackend::new(&context, (width as u32, height as u32)).unwrap().into_drawing_area();
+            let _ = root.fill(&BLACK);
+            
+            const R_DISP: f64 = 1.5;
+            let mut chart = ChartBuilder::on(&root)
+                .build_cartesian_2d(-R_DISP..R_DISP, -R_DISP..R_DISP)
+                .unwrap();
+            let time_step_val = (time_step_for_draw.value() as usize).min(TIME_RES-1);
+            let eigval_iter = eigval_series[time_step_val].iter();
+            let _ = chart.draw_series(
+                eigval_iter.map(|z| Circle::new((z.re, z.im), 3, WHITE.filled()))
+            );
+            let _ = root.present();
+        };
+        DrawingAreaExtManual::set_draw_func(&drawing_area, draw_eigvals);
+        container.append(&drawing_area);
+        
+        // set up the time step slider
+        let time_step_scale = Scale::new(Orientation::Horizontal, Some(&time_step));
+        time_step_scale.connect_value_changed(move |_: &Scale| {
+            drawing_area.queue_draw();
+        });
+        container.append(&time_step_scale);
+        
+        // show the window
+        window.present();
+    });
+    
+    app.run()
+}

lang-trials/rust-benchmark/.gitignore

@@ -0,0 +1,3 @@
+target
+dist
+Cargo.lock

lang-trials/rust-benchmark/Cargo.toml

@@ -0,0 +1,34 @@
+[package]
+name = "rust-benchmark"
+version = "0.1.0"
+authors = ["Aaron"]
+edition = "2021"
+
+[features]
+default = ["console_error_panic_hook"]
+
+[dependencies]
+nalgebra = "0.33.0"
+sycamore = "0.9.0-beta.2"
+
+# The `console_error_panic_hook` crate provides better debugging of panics by
+# logging them with `console.error`. This is great for development, but requires
+# all the `std::fmt` and `std::panicking` infrastructure, so isn't great for
+# code size when deploying.
+console_error_panic_hook = { version = "0.1.7", optional = true }
+
+[dependencies.web-sys]
+version = "0.3.69"
+features = [
+  'CanvasRenderingContext2d',
+  'HtmlCanvasElement',
+  'Window',
+  'Performance'
+]
+
+[dev-dependencies]
+wasm-bindgen-test = "0.3.34"
+
+[profile.release]
+opt-level = "s" # optimize for small code size
+debug = true # include debug symbols

lang-trials/rust-benchmark/index.html

@@ -0,0 +1,9 @@
+<!DOCTYPE html>
+<html>
+  <head>
+    <meta charset="utf-8"/>
+    <title>The circular law</title>
+    <link data-trunk rel="css" href="main.css"/>
+  </head>
+  <body></body>
+</html>

lang-trials/rust-benchmark/main.css

@@ -0,0 +1,23 @@
+body {
+  margin-left: 20px;
+  margin-top: 20px;
+  color: #fcfcfc;
+  background-color: #202020;
+}
+
+#app {
+  display: flex;
+  flex-direction: column;
+  width: 600px;
+}
+
+canvas {
+  float: left;
+  background-color: #020202;
+  border-radius: 10px;
+  margin-top: 5px;
+}
+
+input {
+  margin-top: 5px;
+}

lang-trials/rust-benchmark/notes

@@ -0,0 +1,13 @@
+in profiling, most time is being spent in the `reflect` method:
+
+f64:
+  sycamore_trial-3d0aca3efee8b5fd.wasm.nalgebra::geometry::reflection::Reflection<T,D,S>::reflect::h7899977a4ba0b1d3
+  sycamore_trial-3d0aca3efee8b5fd.wasm.nalgebra::geometry::reflection::Reflection<T,D,S>::reflect::hc337c3cb6e3b4061
+  sycamore_trial-3d0aca3efee8b5fd.wasm.nalgebra::geometry::reflection::Reflection<T,D,S>::reflect_rows::h43d0f6838d0c2833
+
+f32:
+  sycamore_trial-3d0aca3efee8b5fd.wasm.nalgebra::geometry::reflection::Reflection<T,D,S>::reflect::h0e8ec322f198f847
+  sycamore_trial-3d0aca3efee8b5fd.wasm.nalgebra::geometry::reflection::Reflection<T,D,S>::reflect::h9928bdd5e72743ea
+  sycamore_trial-3d0aca3efee8b5fd.wasm.nalgebra::geometry::reflection::Reflection<T,D,S>::reflect_rows::h49f571fd8fc9b0f2
+
+in one test, we spent 4000 ms in "WASM closure", but the enveloping "VoidFunction" takes 1300 ms longer. in another test, though, there's no overhang; the 7000 ms we spent in `rand_eigval_series` accounts for basically the entire load time, and matches the clock timing

lang-trials/rust-benchmark/src/engine.rs

@@ -0,0 +1,104 @@
+use nalgebra::{*, allocator::Allocator};
+use std::f64::consts::{PI, E};
+
+/* dynamic matrices */
+pub fn rand_eigval_series(dim: usize, time_res: usize) -> Vec<OVector<Complex<f64>, Dyn>> {
+    // initialize the random matrix
+    let mut rand_mat = DMatrix::<f64>::from_fn(dim, dim, |j, k| {
+        let n = j*dim + k;
+        E*((n*n) as f64) % 2.0 - 1.0
+    }) * (3.0 / (dim as f64)).sqrt();
+    
+    // initialize the rotation step
+    let mut rot_step = DMatrix::<f64>::identity(dim, dim);
+    let max_freq = 4;
+    for n in (0..dim).step_by(2) {
+        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
+        let ang_cos = ang.cos();
+        let ang_sin = ang.sin();
+        rot_step[(n, n)] = ang_cos;
+        rot_step[(n+1, n)] = ang_sin;
+        rot_step[(n, n+1)] = -ang_sin;
+        rot_step[(n+1, n+1)] = ang_cos;
+    }
+    
+    // find the eigenvalues
+    let mut eigval_series = Vec::<OVector<Complex<f64>, Dyn>>::with_capacity(time_res);
+    eigval_series.push(rand_mat.complex_eigenvalues());
+    for _ in 1..time_res {
+        rand_mat = &rot_step * rand_mat;
+        eigval_series.push(rand_mat.complex_eigenvalues());
+    }
+    eigval_series
+}
+
+/* dynamic single float matrices */
+/*pub fn rand_eigval_series(dim: usize, time_res: usize) -> Vec<OVector<Complex<f32>, Dyn>> {
+    // initialize the random matrix
+    let mut rand_mat = DMatrix::<f32>::from_fn(dim, dim, |j, k| {
+        let n = j*dim + k;
+        (E as f32)*((n*n) as f32) % 2.0_f32 - 1.0_f32
+    }) * (3.0_f32 / (dim as f32)).sqrt();
+    
+    // initialize the rotation step
+    let mut rot_step = DMatrix::<f32>::identity(dim, dim);
+    let max_freq = 4;
+    for n in (0..dim).step_by(2) {
+        let ang = (PI as f32) * ((n % max_freq) as f32) / (time_res as f32);
+        let ang_cos = ang.cos();
+        let ang_sin = ang.sin();
+        rot_step[(n, n)] = ang_cos;
+        rot_step[(n+1, n)] = ang_sin;
+        rot_step[(n, n+1)] = -ang_sin;
+        rot_step[(n+1, n+1)] = ang_cos;
+    }
+    
+    // find the eigenvalues
+    let mut eigval_series = Vec::<OVector<Complex<f32>, Dyn>>::with_capacity(time_res);
+    eigval_series.push(rand_mat.complex_eigenvalues());
+    for _ in 1..time_res {
+        rand_mat = &rot_step * rand_mat;
+        eigval_series.push(rand_mat.complex_eigenvalues());
+    }
+    eigval_series
+}*/
+
+/* static matrices. should only be used when the dimension is really small */
+/*pub fn rand_eigval_series<N>(time_res: usize) -> Vec<OVector<Complex<f64>, N>>
+    where
+        N: ToTypenum + DimName + DimSub<U1>,
+        DefaultAllocator:
+            Allocator<N> +
+            Allocator<N, N> +
+            Allocator<<N as DimSub<U1>>::Output> +
+            Allocator<N, <N as DimSub<U1>>::Output>
+{
+    // initialize the random matrix
+    let dim = N::try_to_usize().unwrap();
+    let mut rand_mat = OMatrix::<f64, N, N>::from_fn(|j, k| {
+        let n = j*dim + k;
+        E*((n*n) as f64) % 2.0 - 1.0
+    }) * (3.0 / (dim as f64)).sqrt();
+    
+    // initialize the rotation step
+    let mut rot_step = OMatrix::<f64, N, N>::identity();
+    let max_freq = 4;
+    for n in (0..dim).step_by(2) {
+        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
+        let ang_cos = ang.cos();
+        let ang_sin = ang.sin();
+        rot_step[(n, n)] = ang_cos;
+        rot_step[(n+1, n)] = ang_sin;
+        rot_step[(n, n+1)] = -ang_sin;
+        rot_step[(n+1, n+1)] = ang_cos;
+    }
+    
+    // find the eigenvalues
+    let mut eigval_series = Vec::<OVector<Complex<f64>, N>>::with_capacity(time_res);
+    eigval_series.push(rand_mat.complex_eigenvalues());
+    for _ in 1..time_res {
+        rand_mat = &rot_step * rand_mat;
+        eigval_series.push(rand_mat.complex_eigenvalues());
+    }
+    eigval_series
+}*/

lang-trials/rust-benchmark/src/main.rs

@@ -0,0 +1,78 @@
+use std::f64::consts::PI as PI;
+use sycamore::{prelude::*, rt::{JsCast, JsValue}};
+use web_sys::window;
+
+mod engine;
+
+fn main() {
+    // set up a config option that forwards panic messages to `console.error`
+    #[cfg(feature = "console_error_panic_hook")]
+    console_error_panic_hook::set_once();
+    
+    sycamore::render(|| {
+        let time_res: usize = 100;
+        let time_step = create_signal(0.0);
+        let run_time_report = create_signal(-1.0);
+        let display = create_node_ref();
+        
+        on_mount(move || {
+            let performance = window().unwrap().performance().unwrap();
+            let start_time = performance.now();
+            /*let eigval_series = engine::rand_eigval_series::<U60>(time_res);*/
+            let eigval_series = engine::rand_eigval_series(60, time_res);
+            let run_time = performance.now() - start_time;
+            run_time_report.set(run_time);
+            
+            let canvas = display
+                .get::<DomNode>()
+                .unchecked_into::<web_sys::HtmlCanvasElement>();
+            let ctx = canvas
+                .get_context("2d")
+                .unwrap()
+                .unwrap()
+                .dyn_into::<web_sys::CanvasRenderingContext2d>()
+                .unwrap();
+            ctx.set_fill_style(&JsValue::from("white"));
+            
+            create_effect(move || {
+                // center and normalize the coordinate system
+                let width = canvas.width() as f64;
+                let height = canvas.height() as f64;
+                ctx.set_transform(1.0, 0.0, 0.0, -1.0, 0.5*width, 0.5*height).unwrap();
+                
+                // clear the previous frame
+                ctx.clear_rect(-0.5*width, -0.5*width, width, height);
+                
+                // find the resolution
+                const R_DISP: f64 = 1.5;
+                let res = width / (2.0*R_DISP);
+                
+                // draw the eigenvalues
+                let eigvals = &eigval_series[time_step.get() as usize];
+                for n in 0..eigvals.len() {
+                    ctx.begin_path();
+                    ctx.arc(
+                        /* typecast only needed for single float version */
+                        res * f64::from(eigvals[n].re),
+                        res * f64::from(eigvals[n].im),
+                        3.0,
+                        0.0, 2.0*PI
+                    ).unwrap();
+                    ctx.fill();
+                }
+            });
+        });
+        
+        view! {
+            div(id="app") {
+                div { (run_time_report.get()) " ms" }
+                canvas(ref=display, width="600", height="600")
+                input(
+                    type="range",
+                    max=(time_res - 1).to_string(),
+                    bind:valueAsNumber=time_step
+                )
+            }
+        }
+    });
+}

lang-trials/rust/.gitignore

@@ -0,0 +1,3 @@
+target
+dist
+Cargo.lock

lang-trials/rust/Cargo.toml

@@ -0,0 +1,32 @@
+[package]
+name = "sycamore-trial"
+version = "0.1.0"
+authors = ["Aaron"]
+edition = "2021"
+
+[features]
+default = ["console_error_panic_hook"]
+
+[dependencies]
+nalgebra = "0.33.0"
+sycamore = "0.9.0-beta.2"
+
+# The `console_error_panic_hook` crate provides better debugging of panics by
+# logging them with `console.error`. This is great for development, but requires
+# all the `std::fmt` and `std::panicking` infrastructure, so isn't great for
+# code size when deploying.
+console_error_panic_hook = { version = "0.1.7", optional = true }
+
+[dependencies.web-sys]
+version = "0.3.69"
+features = [
+  'CanvasRenderingContext2d',
+  'HtmlCanvasElement',
+]
+
+[dev-dependencies]
+wasm-bindgen-test = "0.3.34"
+
+[profile.release]
+# Tell `rustc` to optimize for small code size.
+opt-level = "s"

lang-trials/rust/index.html

@@ -0,0 +1,8 @@
+<!DOCTYPE html>
+<html>
+  <head>
+    <title>Lattice circle</title>
+    <link data-trunk rel="css" href="main.css"/>
+  </head>
+  <body></body>
+</html>

lang-trials/rust/main.css

@@ -0,0 +1,50 @@
+body {
+  margin-left: 20px;
+  margin-top: 20px;
+  color: #fcfcfc;
+  background-color: #202020;
+}
+
+input {
+  color: inherit;
+  background-color: #020202;
+  border: 1px solid #606060;
+  min-width: 40px;
+  border-radius: 4px;
+}
+
+input.point-1 {
+  border-color: #ba5d09;
+}
+
+input.point-2 {
+  border-color: #0e8a06;
+}
+
+input.point-3 {
+  border-color: #8951fb;
+}
+
+#data-panel {
+  float: left;
+  margin-left: 20px;
+  display: grid;
+  grid-template-columns: auto auto;
+  gap: 10px 10px;
+  width: 120px;
+}
+
+#data-panel > div {
+  text-align: center;
+}
+
+#result-display {
+  margin-top: 10px;
+  font-weight: bold;
+}
+
+canvas {
+  float: left;
+  background-color: #020202;
+  border-radius: 10px;
+}

lang-trials/rust/src/engine.rs

@@ -0,0 +1,38 @@
+use nalgebra::*;
+
+pub struct Circle {
+  pub center_x: f64,
+  pub center_y: f64,
+  pub radius: f64,
+}
+
+// construct the circle through the points given by the columns of `points`
+pub fn circ_thru(points: Matrix2x3<f64>) -> Option<Circle> {
+    // build the matrix that maps the circle's coefficient vector to the
+    // negative of the linear part of the circle's equation, evaluated at the
+    // given points
+    let neg_lin_part = stack![2.0*points.transpose(), Vector3::repeat(1.0)];
+    
+    // find the quadrdatic part of the circle's equation, evaluated at the given
+    // points
+    let quad_part = Vector3::from_iterator(
+      points.column_iter().map(|v| v.dot(&v))
+    );
+    
+    // find the circle's coefficient vector, and from there its center and
+    // radius
+    match neg_lin_part.lu().solve(&quad_part) {
+        None => None,
+        Some(coeffs) => {
+            let center_x = coeffs[0];
+            let center_y = coeffs[1];
+            Some(Circle {
+                center_x: center_x,
+                center_y: center_y,
+                radius: (
+                    coeffs[2] + center_x*center_x + center_y*center_y
+                ).sqrt(),
+            })
+        }
+    }
+}

lang-trials/rust/src/main.rs

@@ -0,0 +1,114 @@
+use nalgebra::Matrix2x3;
+use std::f64::consts::PI as PI;
+use sycamore::{prelude::*, rt::{JsCast, JsValue}};
+
+mod engine;
+
+fn main() {
+  // set up a config option that forwards panic messages to `console.error`
+  #[cfg(feature = "console_error_panic_hook")]
+  console_error_panic_hook::set_once();
+  
+  sycamore::render(|| {
+    let data = [-1.0, 0.0, 0.0, -1.0, 1.0, 0.0].map(|n| create_signal(n));
+    let display = create_node_ref();
+    
+    on_mount(move || {
+      let canvas = display
+        .get::<DomNode>()
+        .unchecked_into::<web_sys::HtmlCanvasElement>();
+      let ctx = canvas
+        .get_context("2d")
+        .unwrap()
+        .unwrap()
+        .dyn_into::<web_sys::CanvasRenderingContext2d>()
+        .unwrap();
+      
+      create_effect(move || {
+        // center and normalize the coordinate system
+        let width = canvas.width() as f64;
+        let height = canvas.height() as f64;
+        ctx.set_transform(1.0, 0.0, 0.0, -1.0, 0.5*width, 0.5*height).unwrap();
+        
+        // clear the previous frame
+        ctx.clear_rect(-0.5*width, -0.5*width, width, height);
+        
+        // find the resolution
+        const R_DISP: f64 = 5.0;
+        let res = width / (2.0*R_DISP);
+        
+        // set colors
+        let highlight_style = JsValue::from("white");
+        let grid_style = JsValue::from("#404040");
+        let point_fill_styles = ["#ba5d09", "#0e8a06", "#8951fb"];
+        let point_stroke_styles = ["#f89142", "#58c145", "#c396fc"];
+        
+        // draw the grid
+        let r_grid = (R_DISP - 0.01).floor() as i32;
+        let edge_scr = res * R_DISP;
+        ctx.set_stroke_style(&grid_style);
+        for t in -r_grid ..= r_grid {
+          let t_scr = res * (t as f64);
+          
+          // draw horizontal grid line
+          ctx.begin_path();
+          ctx.move_to(-edge_scr, t_scr);
+          ctx.line_to(edge_scr, t_scr);
+          ctx.stroke();
+          
+          // draw vertical grid line
+          ctx.begin_path();
+          ctx.move_to(t_scr, -edge_scr);
+          ctx.line_to(t_scr, edge_scr);
+          ctx.stroke();
+        }
+        
+        // find and draw the circle through the given points
+        let data_vals = data.map(|sig| sig.get()).to_vec();
+        let points = Matrix2x3::from_vec(data_vals);
+        if let Some(circ) = engine::circ_thru(points) {
+          ctx.begin_path();
+          ctx.set_stroke_style(&highlight_style);
+          ctx.arc(
+            res * circ.center_x,
+            res * circ.center_y,
+            res * circ.radius,
+            0.0, 2.0*PI
+          ).unwrap();
+          ctx.stroke();
+        }
+        
+        // draw the data points
+        for n in 0..3 {
+          ctx.begin_path();
+          ctx.set_fill_style(&JsValue::from(point_fill_styles[n]));
+          ctx.set_stroke_style(&JsValue::from(point_stroke_styles[n]));
+          let ind_x = 2*n;
+          let ind_y = ind_x + 1;
+          ctx.arc(
+            res * data[ind_x].get(),
+            res * data[ind_y].get(),
+            3.0,
+            0.0, 2.0*PI
+          ).unwrap();
+          ctx.fill();
+          ctx.stroke();
+        }
+      });
+    });
+    
+    view! {
+      canvas(ref=display, width="600", height="600")
+      div(id="data-panel") {
+        div { "x" }
+        div { "y" }
+        input(type="number", class="point-1", bind:valueAsNumber=data[0])
+        input(type="number", class="point-1", bind:valueAsNumber=data[1])
+        input(type="number", class="point-2", bind:valueAsNumber=data[2])
+        input(type="number", class="point-2", bind:valueAsNumber=data[3])
+        input(type="number", class="point-3", bind:valueAsNumber=data[4])
+        input(type="number", class="point-3", bind:valueAsNumber=data[5])
+      }
+    }
+  });
+}

lang-trials/scala-benchmark/.gitignore

@@ -0,0 +1,2 @@
+target
+sbt.json

lang-trials/scala-benchmark/build.sbt

@@ -0,0 +1,9 @@
+enablePlugins(ScalaJSPlugin)
+
+name := "Circular Law"
+scalaVersion := "3.4.2"
+scalaJSUseMainModuleInitializer := true
+
+libraryDependencies += "com.raquo" %%% "laminar" % "17.0.0"
+libraryDependencies += "ai.dragonfly" %%% "slash" % "0.3.1"
+libraryDependencies += "org.scala-js" %%% "scalajs-dom" % "2.8.0"

lang-trials/scala-benchmark/index.html

@@ -0,0 +1,10 @@
+<!DOCTYPE html>
+<html>
+  <head>
+    <meta charset="UTF-8">
+    <title>The circular law</title>
+    <script type="text/javascript" src="./target/scala-3.4.2/circular-law-opt/main.js"></script>
+    <link rel="stylesheet" href="main.css"/>
+  </head>
+  <body></body>
+</html>

lang-trials/scala-benchmark/main.css

@@ -0,0 +1,23 @@
+body {
+  margin-left: 20px;
+  margin-top: 20px;
+  color: #fcfcfc;
+  background-color: #202020;
+}
+
+#app {
+  display: flex;
+  flex-direction: column;
+  width: 600px;
+}
+
+canvas {
+  float: left;
+  background-color: #020202;
+  border-radius: 10px;
+  margin-top: 5px;
+}
+
+input {
+  margin-top: 5px;
+}

lang-trials/scala-benchmark/project/build.properties

@@ -0,0 +1 @@
+sbt.version=1.10.1

lang-trials/scala-benchmark/project/plugins.sbt

@@ -0,0 +1 @@
+addSbtPlugin("org.scala-js" % "sbt-scalajs" % "1.16.0")

lang-trials/scala-benchmark/src/main/scala/CircularLawApp.scala

@@ -0,0 +1,99 @@
+import com.raquo.laminar.api.L.{*, given}
+import narr.*
+import org.scalajs.dom
+import org.scalajs.dom.document
+import scala.collection.mutable.ArrayBuffer
+import scala.math.{cos, sin}
+import slash.matrix.Matrix
+import slash.matrix.decomposition.Eigen
+
+object CircularLawApp:
+  val canvas = canvasTag(widthAttr := 600, heightAttr := 600)
+  val ctx = canvas.ref.getContext("2d").asInstanceOf[dom.CanvasRenderingContext2D]
+  
+  val (eigvalSeries, runTimeReport) = randEigvalSeries[60]()
+  val timeStepState = Var("0")
+  
+  def draw(timeStep: String): Unit =
+    // center and normalize the coordinate system
+    val width = canvas.ref.width
+    val height = canvas.ref.height
+    ctx.setTransform(1d, 0d, 0d, -1d, 0.5*width, 0.5*height)
+    
+    // clear the previous frame
+    ctx.clearRect(-0.5*width, -0.5*width, width, height)
+    
+    // find the resolution
+    val rDisp: Double = 1.5
+    val res = width / (2*rDisp)
+    
+    // draw the eigenvalues
+    val eigvals = eigvalSeries(timeStep.toInt)
+    for n <- 0 to eigvals(0).length-1 do
+      ctx.beginPath()
+      ctx.arc(
+        res * eigvals(0)(n),
+        res * eigvals(1)(n),
+        3d,
+        0d, 2*math.Pi
+      )
+      ctx.fill()
+  
+  def complexEigenvalues[N <: Int](mat: Matrix[N, N])(using ValueOf[N]): (NArray[Double], NArray[Double]) =
+    val eigen = Eigen(mat)
+    (
+      eigen.realEigenvalues.asInstanceOf[NArray[Double]],
+      eigen.imaginaryEigenvalues.asInstanceOf[NArray[Double]]
+    )
+  
+  def randEigvalSeries[N <: Int]()(using ValueOf[N]): (ArrayBuffer[(NArray[Double], NArray[Double])], String) =
+    // start timing
+    val startTime = System.currentTimeMillis()
+    
+    // initialize the random matrix step
+    val dim: Int = valueOf[N]
+    var randMat = new Matrix[N, N](
+      NArray.tabulate(dim*dim)(k => (math.E*k*k) % 2 - 1)
+    ).times(math.sqrt(3d / dim))
+    
+    // initialize the rotation step
+    val timeRes = 100
+    val maxFreq = 4
+    val rotStep = Matrix.identity[N, N]
+    for n <- 0 to dim by 2 do
+      val ang = math.Pi * (n % maxFreq) / timeRes
+      val cos_ang = cos(ang)
+      val sin_ang = sin(ang)
+      rotStep(n, n) = cos_ang
+      rotStep(n+1, n) = sin_ang
+      rotStep(n, n+1) = -sin_ang
+      rotStep(n+1, n+1) = cos_ang
+    
+    // find the eigenvalues
+    val eigvalSeries = ArrayBuffer(complexEigenvalues(randMat))
+    for _ <- 1 to timeRes-1 do
+      randMat = rotStep * randMat
+      eigvalSeries += complexEigenvalues(randMat)
+    
+    // finish timing
+    val runTime = System.currentTimeMillis() - startTime
+    (eigvalSeries, runTime.toString() + " ms")
+  
+  def main(args: Array[String]): Unit =
+    ctx.fillStyle = "white"
+    
+    lazy val app = div(
+      idAttr := "app",
+      div(runTimeReport),
+      canvas,
+      input(
+        typ := "range",
+        maxAttr := (eigvalSeries.length-1).toString,
+        controlled(
+          value <-- timeStepState.signal,
+          onInput.mapToValue --> timeStepState.writer
+        ),
+        timeStepState.signal --> draw
+      )
+    )
+    renderOnDomContentLoaded(document.body, app)

lang-trials/scala/.gitignore

@@ -0,0 +1,2 @@
+target
+sbt.json

lang-trials/scala/build.sbt

@@ -0,0 +1,12 @@
+enablePlugins(ScalaJSPlugin)
+
+name := "Lattice Circle"
+scalaVersion := "3.4.2"
+
+// This is an application with a main method
+scalaJSUseMainModuleInitializer := true
+
+libraryDependencies += "com.raquo" %%% "laminar" % "17.0.0"
+/*libraryDependencies += "org.scalanlp" %% "breeze" % "2.1.0"*/
+libraryDependencies += "ai.dragonfly" %%% "slash" % "0.3.1"
+libraryDependencies += "org.scala-js" %%% "scalajs-dom" % "2.8.0"

lang-trials/scala/index.html

@@ -0,0 +1,10 @@
+<!DOCTYPE html>
+<html>
+  <head>
+    <meta charset="UTF-8">
+    <title>Lattice circle</title>
+    <script type="text/javascript" src="./target/scala-3.4.2/lattice-circle-fastopt/main.js"></script>
+    <link rel="stylesheet" href="main.css"/>
+  </head>
+  <body></body>
+</html>

lang-trials/scala/main.css

@@ -0,0 +1,45 @@
+body {
+  margin-left: 20px;
+  margin-top: 20px;
+  color: #fcfcfc;
+  background-color: #202020;
+}
+
+input {
+  color: inherit;
+  background-color: #020202;
+  border: 1px solid #606060;
+  min-width: 40px;
+  border-radius: 4px;
+}
+
+input.point-1 {
+  border-color: #ba5d09;
+}
+
+input.point-2 {
+  border-color: #0e8a06;
+}
+
+input.point-3 {
+  border-color: #8951fb;
+}
+
+#data-panel {
+  float: left;
+  margin-left: 20px;
+  display: grid;
+  grid-template-columns: auto auto;
+  gap: 10px 10px;
+  width: 120px;
+}
+
+#data-panel > div {
+  text-align: center;
+}
+
+canvas {
+  float: left;
+  background-color: #020202;
+  border-radius: 10px;
+}

lang-trials/scala/project/build.properties

@@ -0,0 +1 @@
+sbt.version=1.10.1

lang-trials/scala/project/plugins.sbt

@@ -0,0 +1 @@
+addSbtPlugin("org.scala-js" % "sbt-scalajs" % "1.16.0")

lang-trials/scala/src/main/scala/LatticeCircleApp.scala

@@ -0,0 +1,160 @@
+// based on the Laminar example app
+//
+//   https://github.com/raquo/laminar-examples/blob/master/src/main/scala/App.scala
+//
+// and Li Haoyi's example canvas app
+//
+//   http://www.lihaoyi.com/hands-on-scala-js/#MakingaCanvasApp
+//
+
+import com.raquo.laminar.api.L.{*, given}
+import narr.*
+import org.scalajs.dom
+import org.scalajs.dom.document
+import scala.math
+import slash.matrix.*
+
+class Circle(var centerX: Double, var centerY: Double, var radius: Double)
+
+object LatticeCircleApp:
+  val canvas = canvasTag(widthAttr := 600, heightAttr := 600)
+  val ctx = canvas.ref.getContext("2d").asInstanceOf[dom.CanvasRenderingContext2D]
+  val data = List("-1", "0", "0", "-1", "1", "0").map(Var(_))
+  
+  def circThru(points: Matrix[3, 2]): Option[Circle] =
+    // build the matrix that maps the circle's coefficient vector to the
+    // negative of the linear part of the circle's equation, evaluated at the
+    // given points
+    val negLinPart = Matrix.ones[3, 3]
+    negLinPart.setMatrix(0, 0, points * 2.0)
+    
+    // find the quadrdatic part of the circle's equation, evaluated at the given
+    // points
+    val quadPart = Matrix[3, 1](
+      NArray.tabulate[Double](3)(
+        k => points(k, 0)*points(k, 0) + points(k, 1)*points(k, 1)
+      )
+    )
+    
+    // find the circle's coefficient vector, and from there its center and
+    // radius
+    try
+      val coeffs = negLinPart.solve(quadPart)
+      val centerX = coeffs(0, 0)
+      val centerY = coeffs(1, 0)
+      Some(Circle(
+        centerX,
+        centerY,
+        math.sqrt(coeffs(2, 0) + centerX*centerX + centerY*centerY)
+      ))
+    catch
+      _ => return None
+  
+  def draw(): Unit =
+    // center and normalize the coordinate system
+    val width = canvas.ref.width
+    val height = canvas.ref.height
+    ctx.setTransform(1.0, 0.0, 0.0, -1.0, 0.5*width, 0.5*height)
+    
+    // clear the previous frame
+    ctx.clearRect(-0.5*width, -0.5*width, width, height)
+    
+    // find the resolution
+    val rDisp = 5.0
+    val res = width / (2.0*rDisp)
+    
+    // set colors
+    val highlightStyle = "white"
+    val gridStyle = "#404040"
+    val pointFillStyles = List("#ba5d09", "#0e8a06", "#8951fb")
+    val pointStrokeStyles = List("#f89142", "#58c145", "#c396fc")
+    
+    // draw the grid
+    val rGrid = (rDisp - 0.01).floor.toInt
+    val edgeScr = res * rDisp
+    ctx.strokeStyle = gridStyle
+    for t <- -rGrid to rGrid do
+      val tScr = res * t
+      
+      // draw horizontal grid line
+      ctx.beginPath();
+      ctx.moveTo(-edgeScr, tScr)
+      ctx.lineTo(edgeScr, tScr)
+      ctx.stroke()
+      
+      // draw vertical grid line
+      ctx.beginPath();
+      ctx.moveTo(tScr, -edgeScr)
+      ctx.lineTo(tScr, edgeScr)
+      ctx.stroke()
+    
+    // find and draw the circle through the given points
+    val dataNow = NArray.tabulate(6)(n =>
+      try
+        data(n).signal.now().toDouble
+      catch
+        _ => Double.NaN
+    )
+    if dataNow.forall(t => t == t.floor) then
+      // all of the coordinates are integer and non-NaN
+      val points = Matrix[3, 2](dataNow)
+      circThru(points) match
+        case Some(circ) =>
+          ctx.beginPath()
+          ctx.strokeStyle = highlightStyle
+          ctx.arc(
+            res * circ.centerX,
+            res * circ.centerY,
+            res * circ.radius,
+            0.0, 2.0*math.Pi
+          )
+          ctx.stroke()
+        case None =>
+    
+    // draw the data points
+    for n <- 0 to 2 do
+      val indX = 2*n
+      val indY = indX + 1
+      if
+        dataNow(indX) == dataNow(indX).floor &&
+        dataNow(indY) == dataNow(indY).floor
+      then
+        ctx.beginPath()
+        ctx.fillStyle = pointFillStyles(n)
+        ctx.strokeStyle = pointStrokeStyles(n)
+        ctx.arc(
+          res * dataNow(indX),
+          res * dataNow(indY),
+          3.0,
+          0.0, 2.0*math.Pi
+        )
+        ctx.fill()
+        ctx.stroke()
+  
+  def coordInput(n: Int): Input =
+    input(
+      typ := "number",
+      cls := s"point-${(1.0 + 0.5*n).floor.toInt}",
+      controlled(
+        value <-- data(n).signal,
+        onInput.mapToValue --> data(n).writer
+      ),
+      data(n).signal --> { _ => draw() }
+    )
+  
+  def main(args: Array[String]): Unit =
+    lazy val app = div(
+      canvas,
+      div(
+        idAttr := "data-panel",
+        div("x"),
+        div("y"),
+        coordInput(0),
+        coordInput(1),
+        coordInput(2),
+        coordInput(3),
+        coordInput(4),
+        coordInput(5)
+      )
+    )
+    renderOnDomContentLoaded(document.body, app)

notes/inversive.md

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