engine-proto/ConstructionViewer.jl

@@ -65,12 +65,8 @@ mutable struct ConstructionViewer
       }
     """)
     
-    # load Ganja.js. for an automatically updated web-hosted version, load from
+    # load Ganja.js
-    #
+    loadjs!(win, "https://unpkg.com/ganja.js")
-    #   https://unpkg.com/ganja.js
-    #
-    # instead
-    loadjs!(win, "http://localhost:8000/ganja-1.0.204.js")
     
     # create global functions and variables
     script!(win, """
@@ -127,24 +123,12 @@ mutable struct ConstructionViewer
   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 = []
+  elements_full = [
+    [0; elt; fill(0, 26)]
     for elt in eachcol(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,

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

@@ -28,25 +28,6 @@
     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;
@@ -85,7 +66,6 @@
       
       // de-noise
       for (let k = 6; k < elements[n].length; ++k) {
-      /*for (let k = 0; k < 26; ++k) {*/
         elements[n][k] = 0;
       }
     }

engine-proto/gram-test/Engine.jl

@@ -1,146 +0,0 @@
-module Engine
-
-using LinearAlgebra
-using SparseArrays
-using Random
-
-export rand_on_shell, Q, DescentHistory, 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)
-  [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 ===
-
-plane(normal, offset) = [normal; offset; offset]
-
-function sphere(center, radius)
-  dist_sq = dot(center, center)
-  return [
-    center / radius;
-    0.5 * ((dist_sq - 1) / radius - radius);
-    0.5 * ((dist_sq + 1) / radius - radius)
-  ]
-end
-
-# === Gram matrix realization ===
-
-# the Lorentz form
-Q = diagm([1, 1, 1, 1, -1])
-
-# 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}
-  slope::Array{T}
-  stepsize::Array{T}
-  backoff_steps::Array{Int64}
-  
-  function DescentHistory{T}(
-    scaled_loss = Array{T}(undef, 0),
-    slope = Array{T}(undef, 0),
-    stepsize = Array{T}(undef, 0),
-    backoff_steps = Int64[]
-  ) where T
-    new(scaled_loss, slope, stepsize, backoff_steps)
-  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(
-  gram::SparseMatrixCSC{T, <:Any},
-  guess::Matrix{T};
-  scaled_tol = 1e-30,
-  target_improvement = 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 = norm(Δ_proj)
-  for step 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)
-    
-    # store current position, loss, and slope
-    L_last = L
-    loss_last = loss
-    push!(history.scaled_loss, loss / scale_adjustment)
-    push!(history.slope, slope)
-    
-    # find a good step size using backtracking line search
-    push!(history.stepsize, 0)
-    push!(history.backoff_steps, max_backoff_steps)
-    for backoff_steps in 0:max_backoff_steps
-      history.stepsize[end] = stepsize
-      L = L_last + stepsize * neg_grad
-      Δ_proj = proj_diff(gram, L'*Q*L)
-      loss = norm(Δ_proj)
-      improvement = loss_last - loss
-      if improvement >= target_improvement * stepsize * slope
-        history.backoff_steps[end] = backoff_steps
-        break
-      end
-      stepsize *= backoff
-    end
-  end
-  
-  # return the factorization and its history
-  push!(history.scaled_loss, loss / scale_adjustment)
-  L, history
-end
-
-end

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

@@ -1,79 +0,0 @@
-include("Engine.jl")
-
-using SparseArrays
-using AbstractAlgebra
-using PolynomialRoots
-
-# 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
-    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)
-## Random.seed!(58271) # stuck; step size collapses on step 48
-## Random.seed!(58272) # good convergence
-## Random.seed!(58273) # stuck; step size collapses on step 18
-## Random.seed!(58274) # stuck
-## Random.seed!(58275) #
-## guess = Engine.rand_on_shell(fill(BigFloat(-1), 8))
-
-# set initial guess
-#=
-guess = hcat(
-  Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
-  Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)),
-  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(BigFloat[-1, 0, 0], BigFloat(1//5)),
-  Engine.sphere(BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//5)),
-  Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5))
-)
-=#
-guess = hcat(
-  Engine.plane(BigFloat[0, 0, 1], BigFloat(0)),
-  Engine.sphere(BigFloat[0, 0, 0], BigFloat(0.9)),
-  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))
-)
-
-# complete the gram matrix using gradient descent
-L, history = Engine.realize_gram(gram, guess, max_descent_steps = 200)
-completed_gram = L'*Engine.Q*L
-println("Completed Gram matrix:\n")
-display(completed_gram)
-println("\nSteps: ", size(history.stepsize, 1))
-println("Loss: ", history.scaled_loss[end], "\n")

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

@@ -0,0 +1,137 @@
+using LinearAlgebra
+using SparseArrays
+using AbstractAlgebra
+using PolynomialRoots
+
+# testing Gram matrix recovery using a homemade gradient descent routine
+
+# === gradient descent ===
+
+# 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, attempt)
+  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
+
+# === example ===
+
+# the Lorentz form
+Q = diagm([1, 1, 1, 1, -1])
+
+# 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
+#
+# using gram[6, 2] or gram[7, 1] as the independent variable seems to stall
+# convergence, even if its value comes from a known solution, like
+#
+#   gram[6, 2] = 0.9936131705272925
+#
+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
+guess = sqrt(0.5) * BigFloat[
+  1  1 -1 -1  0       -0.1  0.3;
+  1 -1  1 -1  0       -0.5  0.4;
+  1 -1 -1  1  0        0.1 -0.2;
+  0  0  0  0 -sqrt(6)  0.3 -0.2;
+  1  1  1  1  2        0.2  0.1;
+]
+
+# search parameters
+steps = 600
+line_search_max_steps = 100
+init_stepsize = BigFloat(1)
+step_shrink_factor = BigFloat(0.5)
+target_improvement_factor = BigFloat(0.5)
+
+# complete the gram matrix using gradient descent
+loss_history = Array{BigFloat}(undef, steps + 1)
+stepsize_history = Array{BigFloat}(undef, steps)
+line_search_depth_history = fill(line_search_max_steps, steps)
+stepsize = init_stepsize
+L = copy(guess)
+Δ_proj = proj_diff(gram, L'*Q*L)
+loss = norm(Δ_proj)
+for step in 1:steps
+  # find negative gradient of loss function
+  neg_grad = 4*Q*L*Δ_proj
+  slope = norm(neg_grad)
+  
+  # store current position and loss
+  L_last = L
+  loss_last = loss
+  loss_history[step] = loss
+  
+  # find a good step size using backtracking line search
+  for line_search_depth in 1:line_search_max_steps
+    stepsize_history[step] = stepsize
+    global L = L_last + stepsize * neg_grad
+    global Δ_proj = proj_diff(gram, L'*Q*L)
+    global loss = norm(Δ_proj)
+    improvement = loss_last - loss
+    if improvement >= target_improvement_factor * stepsize * slope
+      line_search_depth_history[step] = line_search_depth
+      break
+    end
+    global stepsize *= step_shrink_factor
+  end
+end
+completed_gram = L'*Q*L
+loss_history[steps + 1] = loss
+println("Completed Gram matrix:\n")
+display(completed_gram)
+println("\nLoss: ", loss, "\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/overlapping-pyramids.jl

@@ -1,86 +0,0 @@
-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
-#
-# using gram[6, 2] or gram[7, 1] as the independent variable seems to stall
-# convergence, even if its value comes from a known solution, like
-#
-#   gram[6, 2] = 0.9936131705272925
-#
-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 = hcat(
-  sqrt(1/BigFloat(2)) * BigFloat[
-    1  1 -1 -1  0;
-    1 -1  1 -1  0;
-    1 -1 -1  1  0;
-    0  0  0  0 -sqrt(BigFloat(6));
-    1  1  1  1  2;
-  ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5)),
-  Engine.rand_on_shell(fill(BigFloat(-1), 2))
-)
-
-# complete the gram matrix using gradient descent
-L, history = Engine.realize_gram(gram, guess)
-completed_gram = L'*Engine.Q*L
-println("Completed Gram matrix:\n")
-display(completed_gram)
-println("\nSteps: ", size(history.stepsize, 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

@@ -1,44 +0,0 @@
-include("Engine.jl")
-
-using SparseArrays
-using AbstractAlgebra
-using PolynomialRoots
-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: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)
-
-# set initial guess
-Random.seed!(99230)
-guess = sqrt(1/BigFloat(3)) * BigFloat[
-  1  1 -1 -1  0
-  1 -1  1 -1  0
-  1 -1 -1  1  0
-  1  1  1  1 -2
-  1  1  1  1  1
-] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5))
-
-# complete the gram matrix using gradient descent
-L, history = Engine.realize_gram(gram, guess)
-completed_gram = L'*Engine.Q*L
-println("Completed Gram matrix:\n")
-display(completed_gram)
-println("\nSteps: ", size(history.stepsize, 1))
-println("Loss: ", history.scaled_loss[end], "\n")