Compare commits
6 Commits
Author | SHA1 | Date | |
---|---|---|---|
2400b1d143 | |||
f5262ebeea | |||
2e2add905b | |||
f22c4c9192 | |||
6f82292b52 | |||
c4dafc635e |
@ -1,10 +1,4 @@
|
||||
# DEPRECATED!
|
||||
|
||||
This repository has not been updated for a while!
|
||||
If You are interested in replicating results for [1712.07167](https://arxiv.org/abs/1712.07167) please check [these instruction](https://kalmar.faculty.wmi.amu.edu.pl/post/1712.07176/)
|
||||
Also [this notebook](https://nbviewer.jupyter.org/gist/kalmarek/03510181bc1e7c98615e86e1ec580b2a) could be of some help. If everything else fails the [zenodo dataset](https://zenodo.org/record/1133440) should contain the last-resort instructions.
|
||||
|
||||
This repository contains some legacy code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
|
||||
This repository contains code for computations in [Certifying Numerical Estimates of Spectral Gaps](https://arxiv.org/abs/1703.09680).
|
||||
|
||||
# Installing
|
||||
To run the code You need `julia-v0.5` (should work on `v0.6`, but with warnings).
|
||||
|
157
paper_data/check_positivity.jl
Normal file
157
paper_data/check_positivity.jl
Normal file
@ -0,0 +1,157 @@
|
||||
using AbstractAlgebra
|
||||
using Groups
|
||||
using GroupRings
|
||||
using PropertyT
|
||||
using IntervalArithmetic
|
||||
|
||||
using SCS
|
||||
using JLD
|
||||
|
||||
include("sqadjop.jl")
|
||||
|
||||
invalid_use_message = """You need to call this script in the parent folder of oSAutF5_r2 folder.
|
||||
Provide also the two (numerical) parameters: `-k` and `-lambda`"""
|
||||
|
||||
if !(iseven(length(ARGS)))
|
||||
throw(invalid_use_message)
|
||||
end
|
||||
|
||||
K = LAMBDA = nothing
|
||||
|
||||
for i in 1:2:length(ARGS)
|
||||
if ARGS[i] == "-k"
|
||||
K = parse(Float64, ARGS[i+1])
|
||||
elseif ARGS[i] == "-lambda"
|
||||
LAMBDA = parse(Float64, ARGS[i+1])
|
||||
end
|
||||
end
|
||||
|
||||
if K == nothing || LAMBDA == nothing
|
||||
throw(invalid_use_message)
|
||||
end
|
||||
|
||||
info("Running checks for 36(Adj₅ + $K·Op₅) - $LAMBDA·Δ₅")
|
||||
|
||||
G = AutGroup(FreeGroup(5), special=true)
|
||||
S = generating_set(G)
|
||||
|
||||
const prefix = "oSAutF5_r2"
|
||||
isdir(prefix) || mkpath(prefix)
|
||||
|
||||
const DELTA_FILE = joinpath(prefix,"delta.jld")
|
||||
const SQADJOP_FILE = joinpath(prefix, "SqAdjOp_coeffs.jld")
|
||||
const ORBITDATA_FILE = joinpath(prefix, "OrbitData.jld")
|
||||
|
||||
const fullpath = joinpath(prefix, string(LAMBDA))
|
||||
isdir(fullpath) || mkpath(fullpath)
|
||||
const SOLUTION_FILE = joinpath(fullpath, "solution.jld")
|
||||
|
||||
info("Looking for delta.jld, SqAdjOp_coeffs.jld and OrbitData.jld in $prefix")
|
||||
|
||||
if isfile(DELTA_FILE) && isfile(SQADJOP_FILE) && isfile(ORBITDATA_FILE)
|
||||
# cached
|
||||
Δ = PropertyT.loadGRElem(DELTA_FILE, G)
|
||||
RG = parent(Δ)
|
||||
orbit_data = load(ORBITDATA_FILE, "OrbitData")
|
||||
sq_c, adj_c, op_c = load(SQADJOP_FILE, "Sq", "Adj", "Op")
|
||||
sq = GroupRingElem(sq_c, RG)
|
||||
adj = GroupRingElem(adj_c, RG)
|
||||
op = GroupRingElem(op_c, RG);
|
||||
else
|
||||
info("Computing Laplacian")
|
||||
Δ = PropertyT.Laplacian(S, 2)
|
||||
PropertyT.saveGRElem(DELTA_FILE, Δ)
|
||||
RG = parent(Δ)
|
||||
|
||||
info("Computing Sq, Adj, Op")
|
||||
@time sq, adj, op = Sq(RG), Adj(RG), Op(RG)
|
||||
|
||||
save(SQADJOP_FILE, "Sq", sq.coeffs, "Adj", adj.coeffs, "Op", op.coeffs)
|
||||
|
||||
info("Compute OrbitData")
|
||||
if !isfile(ORBITDATA_FILE)
|
||||
orbit_data = PropertyT.OrbitData(RG, sett.autS)
|
||||
save(ORBITDATA_FILE, "OrbitData", orbit_data)
|
||||
else
|
||||
orbit_data = load(ORBITDATA_FILE, "OrbitData")
|
||||
end
|
||||
end;
|
||||
|
||||
orbit_data = PropertyT.decimate(orbit_data);
|
||||
|
||||
elt = adj + K*op;
|
||||
|
||||
info("Looking for solution.jld in $fullpath")
|
||||
|
||||
if !isfile(SOLUTION_FILE)
|
||||
info("solution.jld not found, attempting to recreate one.")
|
||||
|
||||
SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=LAMBDA)
|
||||
|
||||
begin
|
||||
scs_solver = SCS.SCSSolver(linear_solver=SCS.Direct,
|
||||
eps=1e-12,
|
||||
max_iters=200_000,
|
||||
alpha=1.5,
|
||||
acceleration_lookback=1)
|
||||
|
||||
JuMP.setsolver(SDP_problem, scs_solver)
|
||||
end
|
||||
|
||||
λ = Ps = ws = nothing
|
||||
const WARMSTART_FILE = joinpath(fullpath, "warmstart.jld")
|
||||
if isfile(WARMSTART_FILE)
|
||||
ws = load(WARMSTART_FILE, "warmstart")
|
||||
end
|
||||
|
||||
i = 0
|
||||
# for i in 1:6
|
||||
status= :Unknown
|
||||
while status !=:Optimal
|
||||
i += 1
|
||||
SOLVERLOG_FILE = joinpath(fullpath, "solver_$(now()).log")
|
||||
status, (λ, Ps, ws) = PropertyT.solve(SOLVERLOG_FILE, SDP_problem, varλ, varP, ws);
|
||||
precision = abs(λ - LAMBDA)
|
||||
println("i = $i, \t precision = $precision")
|
||||
|
||||
if all((!isnan).(ws[1]))
|
||||
save(WARMSTART_FILE, "warmstart", ws, "λ", λ, "Ps", Ps)
|
||||
save(WARMSTART_FILE[1:end-4]*"_$(now()).jld", "warmstart", ws, "λ", λ, "Ps", Ps)
|
||||
else
|
||||
warn("No valid solution was saved!")
|
||||
end
|
||||
end
|
||||
|
||||
info("Reconstructing P...")
|
||||
@time P = PropertyT.reconstruct(Ps, orbit_data);
|
||||
info("Computing Q = √P")
|
||||
@time const Q = real(sqrtm(P));
|
||||
|
||||
save(SOLUTION_FILE, "λ", λ, "Q", Q)
|
||||
end
|
||||
|
||||
info("Checking the sum of squares solution for 36(Adj₅ + $K·Op₅) - $LAMBDA·Δ₅")
|
||||
Q, λ = load(SOLUTION_FILE, "Q", "λ")
|
||||
|
||||
function SOS_residual(eoi::GroupRingElem, Q::Matrix)
|
||||
RG = parent(eoi)
|
||||
@time sos = PropertyT.compute_SOS(RG, Q);
|
||||
return eoi - sos
|
||||
end
|
||||
|
||||
info("In floating point arithmetic:")
|
||||
EOI = elt - λ*Δ
|
||||
b = SOS_residual(EOI, Q)
|
||||
@show norm(b, 1);
|
||||
|
||||
info("In interval arithmetic:")
|
||||
EOI_int = elt - @interval(λ)*Δ;
|
||||
Q_int = PropertyT.augIdproj(Q);
|
||||
@assert all([zero(eltype(Q)) in sum(view(Q_int, :, i)) for i in 1:size(Q_int, 2)])
|
||||
b_int = SOS_residual(EOI_int, Q_int)
|
||||
@show norm(b_int, 1);
|
||||
|
||||
λ_cert = @interval(λ) - 2^2*norm(b_int,1)
|
||||
info("λ is certified to be > ", λ_cert.lo)
|
||||
|
||||
info("i.e Adj₅ + $K·Op₅ - $((λ_cert/36).lo)·Δ₅ ∈ Σ²₂ ISAut(F₅)")
|
51
paper_data/sqadjop.jl
Normal file
51
paper_data/sqadjop.jl
Normal file
@ -0,0 +1,51 @@
|
||||
indexing(n) = [(i,j) for i in 1:n for j in 1:n if i≠j]
|
||||
|
||||
function generating_set(G::AutGroup{N}, n=N) where N
|
||||
|
||||
rmuls = [Groups.rmul_autsymbol(i,j) for (i,j) in indexing(n)]
|
||||
lmuls = [Groups.lmul_autsymbol(i,j) for (i,j) in indexing(n)]
|
||||
gen_set = G.([rmuls; lmuls])
|
||||
|
||||
return [gen_set; inv.(gen_set)]
|
||||
end
|
||||
|
||||
function Sq(RG::GroupRing{AutGroup{N}}) where N
|
||||
S₂ = generating_set(RG.group, 2)
|
||||
ℤ = Int64
|
||||
Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
|
||||
|
||||
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
|
||||
elt = sum(σ(Δ₂)^2 for σ in Alt_N)
|
||||
# return RG(elt.coeffs÷factorial(N-2))
|
||||
return elt
|
||||
end
|
||||
|
||||
function Adj(RG::GroupRing{AutGroup{N}}) where N
|
||||
S₂ = small_generating_set(RG, 2)
|
||||
ℤ = Int64
|
||||
Δ₂ = length(T2)*one(RG, ℤ) - RG(S₂, ℤ);
|
||||
|
||||
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
|
||||
|
||||
adj(σ::perm, τ::perm, i=1, j=2) = Set([σ[i], σ[j]]) ∩ Set([τ[i], τ[j]])
|
||||
adj(σ::perm) = [τ for τ in Alt_N if length(adj(σ, τ)) == 1]
|
||||
|
||||
@time elt = sum(σ(Δ₂)*sum(τ(Δ₂) for τ in adj(σ)) for σ in Alt_N);
|
||||
# return RG(elt.coeffs÷factorial(N-2)^2)
|
||||
return elt
|
||||
end
|
||||
|
||||
function Op(RG::GroupRing{AutGroup{N}}) where N
|
||||
S₂ = small_generating_set(RG, 2)
|
||||
ℤ = Int64
|
||||
Δ₂ = length(T2)*one(RG, ℤ) - RG(S₂, ℤ);
|
||||
|
||||
Alt_N = [g for g in elements(PermutationGroup(N)) if parity(g) == 0]
|
||||
|
||||
adj(σ::perm, τ::perm, i=1, j=2) = Set([σ[i], σ[j]]) ∩ Set([τ[i], τ[j]])
|
||||
op(σ::perm) = [τ for τ in Alt_N if length(adj(σ, τ)) == 0]
|
||||
|
||||
@time elt = sum(σ(Δ₂)*sum(τ(Δ₂) for τ in op(σ)) for σ in Alt_N);
|
||||
# return RG(elt.coeffs÷factorial(N-2)^2)
|
||||
return elt
|
||||
end
|
Loading…
Reference in New Issue
Block a user