PropertyT.jl/AutF4.jl

using Combinatorics

using JuMP
import SCS: SCSSolver
import Mosek: MosekSolver

push!(LOAD_PATH, "./")
using SemiDirectProduct
using GroupAlgebras
include("property(T).jl")

const N = 4

const VERBOSE = true


function permutation_matrix(p::Vector{Int})
    n = length(p)
    sort(p) == collect(1:n) || throw(ArgumentError("Input array must be a permutation of 1:n"))
    A = eye(n)
    return A[p,:]
end

SymmetricGroup(n::Int) = [nthperm(collect(1:n), k) for k in 1:factorial(n)]

# const SymmetricGroup = [permutation_matrix(x) for x in SymmetricGroup_perms]

function E(i::Int, j::Int; dim::Int=N)
    @assert i≠j
    k = eye(dim)
    k[i,j] = 1
    return k
end

function eltary_basis_vector(i::Int; dim::Int=N)
    result = zeros(dim)
    if 0 < i ≤ dim
        result[i] = 1
    end
    return result
end

v(i::Int; dim=N) = eltary_basis_vector(i,dim=dim)

ϱ(i::Int,j::Int,n=N) = SemiDirectProductElement(E(i,j,dim=n), v(j,dim=n))
λ(i::Int,j::Int,n=N) = SemiDirectProductElement(E(i,j,dim=n), -v(j,dim=n))

function ɛ(i::Int, n::Int=N)
    result = eye(n)
    result[i,i] = -1
    return SemiDirectProductElement(result)
end

σ(permutation::Vector{Int}) =
     SemiDirectProductElement(permutation_matrix(permutation))

function AutF_generating_set(n::Int=N)
    indexing = [[i,j] for i in 1:n for j in 1:n if i≠j]
    ϱs = [ϱ(ij...) for ij in indexing]
    λs = [λ(ij...) for ij in indexing]
    ɛs = [ɛ(i) for i in 1:N]
    σs = [σ(perm) for perm in SymmetricGroup(n)]
    S = vcat(ϱs, λs, ɛs, σs);
    S = unique(vcat(S, [inv(x) for x in S]));
    return S
end

const ID = eye(N+1)

const S₁ = AutF_generating_set(N)

matrix_S₁ = [matrix_repr(x) for x in S₁]

const TOL=10.0^-7

matrix_S₁[1:10,:][:,1]

Δ, cm = prepare_Laplacian_and_constraints(matrix_S₁)

#solver = SCSSolver(eps=TOL, max_iters=ITERATIONS, verbose=true);
solver = MosekSolver(MSK_DPAR_INTPNT_CO_TOL_REL_GAP=TOL,
#                      MSK_DPAR_INTPNT_CO_TOL_PFEAS=1e-15,
#                      MSK_DPAR_INTPNT_CO_TOL_DFEAS=1e-15,
#                      MSK_IPAR_PRESOLVE_USE=0,
                  QUIET=!VERBOSE)

# κ, A = solve_for_property_T(S₁, solver, verbose=VERBOSE)

product_matrix = readdlm("SL₃Z.product_matrix", Int)
L = readdlm("SL₃Z.Δ.coefficients")[:, 1]
Δ = GroupAlgebraElement(L, product_matrix)

A = readdlm("matrix.A.Mosek")
κ = readdlm("kappa.Mosek")[1]

# @show eigvals(A)
@assert isapprox(eigvals(A), abs(eigvals(A)), atol=TOL)
@assert A == Symmetric(A)


const A_sqrt = real(sqrtm(A))

SOS_EOI_fp_L₁, Ω_fp_dist = check_solution(κ, A_sqrt, Δ)

κ_rational = rationalize(BigInt, κ;)
A_sqrt_rational = rationalize(BigInt, A_sqrt)
Δ_rational = rationalize(BigInt, Δ)

SOS_EOI_rat_L₁, Ω_rat_dist = check_solution(κ_rational, A_sqrt_rational, Δ_rational)