Start "circles in triangle" from a very close guess

Randomly perturb the pre-solved part of the guess, and randomly choose
the unsolved part.

engine-proto/ConstructionViewer.jl

@@ -65,8 +65,12 @@ mutable struct ConstructionViewer
       }
     """)
     
-    # load Ganja.js
-    loadjs!(win, "https://unpkg.com/ganja.js")
+    # 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, """
@@ -123,12 +127,24 @@ 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 = [
-    [0; elt; fill(0, 26)]
-    for elt in eachcol(elements)
-  ]
+  elements_full = []
+  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,6 +28,25 @@
     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;
@@ -66,6 +85,7 @@
       
       // 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

@@ -0,0 +1,146 @@
+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

@@ -0,0 +1,79 @@
+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

@@ -1,137 +0,0 @@
-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

@@ -0,0 +1,86 @@
+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

@@ -0,0 +1,44 @@
+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")