using Pkg
Pkg.activate("..")
using Dates
now()

using LinearAlgebra
using AbstractAlgebra
using Groups
using GroupRings
using PropertyT

N = 4
G = MatrixAlgebra(zz, N)

S = PropertyT.generating_set(G)

halfradius = 2
E_R, sizes = Groups.generate_balls(S, radius=2*halfradius);
E_rdict = GroupRings.reverse_dict(E_R)
pm = GroupRings.create_pm(E_R, E_rdict, sizes[halfradius]; twisted=true);
RG = GroupRing(G, E_R, E_rdict, pm)
@show sizes;
Δ = length(S)*one(RG) - sum(RG(s) for s in S)

od = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
orbit_data = PropertyT.decimate(od);

@time AdjN = PropertyT.Adj(RG, N)
@time OpN = PropertyT.Op(RG, N);

if N == 3
    k = 0
elseif N == 4
    k = 1
elseif N == 5
    k = 1.5
end
elt = AdjN + k*OpN;
elt.coeffs

# @time SDP_problem, varλ, varP = PropertyT.SOS_problem(elt, Δ, orbit_data)
if N == 3
    UB = 0.158
elseif N == 4
    UB = 0.82005
elseif N == 5
    UB = 1.5005
end
SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=UB)

using JuMP
using SCS
λ = Ps = warm = nothing

with_SCS = with_optimizer(SCS.Optimizer, 
    linear_solver=SCS.Direct, 
    eps=3e-13,
    max_iters=10000,
    alpha=1.5,
    acceleration_lookback=10,
    warm_start=true)

status, warm = PropertyT.solve(SDP_problem, with_SCS, warm);

λ = value(SDP_problem[:λ])
Ps = [value.(P) for P in varP]
@show(status, λ);

Qs = real.(sqrt.(Ps));
Q = PropertyT.reconstruct(Qs, orbit_data);

function SOS_residual(x::GroupRingElem, Q::Matrix)
    RG = parent(x)
    @time sos = PropertyT.compute_SOS(RG, Q);
    return x - sos
end

residual = SOS_residual(elt - λ*Δ, Q)
@show norm(residual, 1);

using IntervalArithmetic
IntervalArithmetic.setrounding(Interval, :tight)
IntervalArithmetic.setformat(sigfigs=12);

Q_aug, check_columns_augmentation = PropertyT.augIdproj(Interval, Q);
@assert check_columns_augmentation
elt_int = elt - @interval(λ)*Δ;
residual_int = SOS_residual(elt_int, Q_aug)
@show norm(residual_int, 1);

certified_λ = @interval(λ) - 2^2*norm(residual_int,1)

certified_λ.lo

using Dates
now()

versioninfo()