Compare commits
No commits in common. "023759a26715e07c2dd00ad274af55ec0f369350" and "8eb1ebb8d2920617d11effd25afb668732b45419" have entirely different histories.
023759a267
...
8eb1ebb8d2
@ -65,12 +65,8 @@ mutable struct ConstructionViewer
|
||||
}
|
||||
""")
|
||||
|
||||
# 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")
|
||||
# load Ganja.js
|
||||
loadjs!(win, "https://unpkg.com/ganja.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 = []
|
||||
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
|
||||
elements_full = [
|
||||
[0; elt; fill(0, 26)]
|
||||
for elt in eachcol(elements)
|
||||
]
|
||||
@js viewer.win elements = $elements_full.map((elt) -> @new CGA3(elt))
|
||||
|
||||
# generate palette. this is Gadfly's `default_discrete_colors` palette,
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
@ -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
|
||||
@ -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")
|
||||
137
engine-proto/gram-test/descent-test.jl
Normal file
137
engine-proto/gram-test/descent-test.jl
Normal file
@ -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)
|
||||
File diff suppressed because it is too large
Load Diff
@ -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)
|
||||
@ -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")
|
||||
Loading…
Reference in New Issue
Block a user