add check_SpNs.jl

bin/check_SpNs.jl

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