92 lines
2.5 KiB
Julia
92 lines
2.5 KiB
Julia
using Pkg
|
||
Pkg.activate(".")
|
||
using Test
|
||
|
||
using PropertyT
|
||
using Nemo
|
||
|
||
using PropertyT.LinearAlgebra
|
||
using PropertyT.SparseArrays
|
||
using PropertyT.JuMP
|
||
|
||
using PropertyT.AbstractAlgebra
|
||
using PropertyT.Groups
|
||
using PropertyT.GroupRings
|
||
|
||
using PropertyT.JLD
|
||
|
||
|
||
# include("../1712.07167/Roots.jl")
|
||
include(joinpath("..", "src", "SpNs.jl"))
|
||
|
||
using .SpNs
|
||
using .SpNs.Roots
|
||
|
||
BLAS.set_num_threads(Threads.nthreads())
|
||
ENV["OMP_NUM_THREADS"] = Threads.nthreads()
|
||
|
||
const N = 2
|
||
const RADIUS = 3
|
||
const UPPER_BOUND = 0.801
|
||
|
||
root_system = SpNs.gens_roots(N)
|
||
S = [g.matrix for g in root_system]
|
||
G = parent(S[1])
|
||
autS = WreathProduct(PermGroup(2), PermGroup(N))
|
||
|
||
using SCS
|
||
with_SCS = with_optimizer(SCS.Optimizer, linear_solver=SCS.Direct, max_iters=400000, eps=1e-11, alpha=1.95, acceleration_lookback=1, warm_start=true)
|
||
|
||
sett = PropertyT.Settings("Sp($(2N),Z)_r$RADIUS", G, S, autS, with_SCS;
|
||
radius=RADIUS, upper_bound=UPPER_BOUND, warmstart=true)
|
||
|
||
PropertyT.print_summary(sett)
|
||
|
||
fp = PropertyT.fullpath(sett)
|
||
isdir(fp) || mkpath(fp)
|
||
|
||
Δ = PropertyT.Laplacian(S, RADIUS)
|
||
RG = parent(Δ)
|
||
Δs = SpNs.laplacians(RG, root_system)
|
||
|
||
Sq = sum( Δα^2 for (α, Δα) in Δs );
|
||
# Op = sum( Δα * sum(Δβ for (β, Δβ) in Δs if isorthogonal(α, β)) for (α, Δα) in Δs );
|
||
Adj = sum( Δα * sum(Δβ for (β, Δβ) in Δs if !isproportional(α, β)) for (α, Δα) in Δs );
|
||
|
||
const elt = Adj
|
||
|
||
solver_logfile(sett) = joinpath(PropertyT.fullpath(sett), "Adj_solver_$(PropertyT.Dates.now()).log")
|
||
|
||
λ, P = PropertyT.approximate_by_SOS(sett, elt, Δ,
|
||
solverlog=solver_logfile(sett))
|
||
|
||
save(PropertyT.filename(sett, :solution), "λ", λ, "P", P)
|
||
λ < 0 && @warn "Solver did not produce a valid solution!"
|
||
|
||
Q = real(sqrt(P));
|
||
certified_λ = PropertyT.certify_SOS_decomposition(elt, Δ, λ, Q, R=RADIUS)
|
||
|
||
@info "The obtained SOS can be use to certify λ ≥" certified_λ
|
||
|
||
|
||
|
||
|
||
#
|
||
# od = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
|
||
# orbit_data = PropertyT.decimate(od);
|
||
#
|
||
# using SCS
|
||
# using JuMP
|
||
# warmstart = nothing
|
||
#
|
||
# UB = 0.88 # for elt = Δ²;
|
||
# SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data, upper_bound=UB);
|
||
# # SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data);
|
||
#
|
||
# sett = PropertyT.Settings("Sp($(2N),Z)", G, S, autS, solver(2000, accel=20);
|
||
# upper_bound=1.3, warmstart=false)
|
||
#
|
||
# with_SCS = with_optimizer(SCS.Optimizer, linear_solver=SCS.Direct, max_iters=100000, eps=1e-12, alpha=1.95, acceleration_lookback=1, warm_start=true)
|
||
#
|
||
# status, warmstart = PropertyT.solve(SDP_problem, with_SCS, warmstart)
|