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)