kacmoody.jl

# using Pkg
# Pkg.activate(joinpath(@__DIR__, ".."))

using PropertyT
using PropertyT.GroupRings
using PropertyT.Groups
using PropertyT.AbstractAlgebra
using PropertyT.JuMP

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

import LinearAlgebra: norm

using SCS
with_SCS(;iters=30_000, acceleration=50, eps=1e-10) = PropertyT.JuMP.with_optimizer(SCS.Optimizer,
    linear_solver=SCS.Direct,
    max_iters=iters,
    eps=eps,
    alpha=(acceleration == 0 ? 1.95 : 1.5),
    acceleration_lookback=acceleration,
    warm_start=true)

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

function check_obtained_solution(model, elt, orderunit)
    λ = value(model[:λ])
    @info "floating_λ = $λ"

    P = value.(model[:P])
    Q = real.(sqrt((P .+ P')./2))
    residual = SOS_residual(elt - λ*orderunit, Q)
    @show norm(residual, 1);

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

    certified_λ = @interval(λ) - 2^2*norm(residual_int, 1)
    @info "certified λ = $certified_λ"
    return certified_λ
end

function create_and_solve(elt, u; upper_bound = Inf, warmstart=nothing)

    sdp_problem = PropertyT.SOS_problem(elt, u, upper_bound = 0.01);
    @show sdp_problem

    status, warmstart = PropertyT.solve(
        sdp_problem,
        with_SCS(iters=5_000, acceleration=25, eps=1e-7),
        warmstart)

    for i in 1:3
        status, ws_new  = PropertyT.solve(sdp_problem,
            with_SCS(acceleration=0, eps=4e-10), warmstart)
        if status != JuMP.MOI.INTERRUPTED
            warmstart = ws_new
        end
        status == JuMP.MOI.OPTIMAL && break
    end

    status == JuMP.MOI.INTERRUPTED && throw("Optimisation step has been interrupted!")

    sdp_problem, warmstart
end

status, ws = nothing, nothing

begin
    const N = 4
    SL, RG = let halfradius = 2, RR = GF(2)
        G = MatrixAlgebra(RR, N)
        S = PropertyT.generating_set(G)
        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)
        G, RG
    end
    @show RG

    Δs = let RSL=RG
        Δs = fill(RSL(0), (N-1,))
        for shift in 1:N-1
            S = [
            PropertyT.EltaryMat(RSL.group, i+shift-1,j+shift-1) for (i,j) in PropertyT.indexing(2)]
            S = unique([S; inv.(S)])
            Δs[shift] = RSL(length(S)) - sum(RSL(s) for s in S)
        end

        let elt = Δs[1]*Δs[2] + Δs[2]*Δs[1]
            Alt_N = [g for g in PermutationGroup(N) if parity(g) == 0]
            @info "" sum(σ(elt) for σ in Alt_N) - PropertyT.Adj(RSL) == RSL(0)
        end

        Δs
    end
end

function Adj(deltas)
    n = length(deltas)
    res = zero(parent(first(deltas)))
    for (i,j) in zip(1:n-1, 2:n)
        res += deltas[i]*deltas[j] + deltas[j]*deltas[i]
    end
    return res
end

function Op(deltas)
    n = length(deltas)
    res = zero(parent(first(deltas)))
    for i in 1:n-2
        for j in i+2:n
            x = deltas[i]*deltas[j]
            @assert x == deltas[j]*deltas[i]
            res += 2x
        end
    end
    return res
end

Sq(deltas) = sum(d^2 for d in deltas)

if N == 4
    @assert Op(Δs) == Δs[1]*Δs[3] + Δs[3]*Δs[1]
    @assert Adj(Δs) == Δs[1]*Δs[2] + Δs[2]*Δs[1] + Δs[2]*Δs[3] + Δs[3]*Δs[2]
end

Δ = sum(Δs)
ws = nothing

Δ123 = let S = PropertyT.generating_set(SL, 3)
    RG(length(S)) - sum(RG(s) for s in S)
end

Δ234 = let S = PropertyT.generating_set(SL, 3)
    S = perm"(1,2,3,4)".(S)
    RG(length(S)) - sum(RG(s) for s in S)
end

begin
    elt = 0.75Sq(Δs) + Adj(Δs) + Op(Δs)
    # elt = Δ^2
    # elt = Sq(Δs) + 0.5*(Δ123*Δ234 + Δ234*Δ123)
    u = Δ

    m, ws = create_and_solve(elt, u, upper_bound=0.01, warmstart=ws)
    check_obtained_solution(m, elt, u)
end