using ArgParse

using Groups
using GroupAlgebras
using PropertyT

import SCS.SCSSolver

#=
Note that the element
    α(i,j,k) = ϱ(i,j)*ϱ(i,k)*inv(ϱ(i,j))*inv(ϱ(i,k)),
which surely belongs to ball of radius 4 in Aut(F₄) becomes trivial under the representation
    Aut(F₄) → GL₄(ℤ)⋉ℤ⁴ → GL₅(ℂ).
Moreover, due to work of Potapchik and Rapinchuk [1] every real representation of Aut(Fₙ) into GLₘ(ℂ) (for m ≤ 2n-2) factors through GLₙ(ℤ)⋉ℤⁿ, so will have the same problem.

We need a different approach: Here we actually compute in Aut(𝔽₄)
=#

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

function generating_set_of_AutF(N::Int)

    indexing = [[i,j] for i in 1:N for j in 1:N if i≠j]
    σs = [symmetric_AutSymbol(perm) for perm in SymmetricGroup(N)[2:end]];
    ϱs = [rmul_AutSymbol(i,j) for (i,j) in indexing]
    λs = [lmul_AutSymbol(i,j) for (i,j) in indexing]
    ɛs = [flip_AutSymbol(i) for i in 1:N];

    S = vcat(ϱs,λs)
    S = vcat(S..., σs..., ɛs)
    S = vcat(S..., [inv(g) for g in S])
    return Vector{AutWord}(unique(S)), one(AutWord)
end

function generating_set_of_OutF(N::Int)
    indexing = [[i,j] for i in 1:N for j in 1:N if i≠j]
    ϱs = [rmul_AutSymbol(i,j) for (i,j) in indexing]
    λs = [lmul_AutSymbol(i,j) for (i,j) in indexing]
    ɛs = [flip_AutSymbol(i) for i in 1:N];

    S = ϱs
    push!(S, λs..., ɛs...)
    push!(S,[inv(g) for g in S]...)

    return Vector{AutWord}(unique(S)), one(AutWord)
end

function generating_set_of_Sym(N::Int)
    σs = [symmetric_AutSymbol(perm) for perm in SymmetricGroup(N)[2:end]];

    S = σs
    push!(S, [inv(s) for s in S]...)
    return Vector{AutWord}(unique(S)), one(AutWord)
end


function products(S1::Vector{AutWord}, S2::Vector{AutWord})
    result = Vector{AutWord}()
    seen = Set{Vector{FGWord}}()
    n = length(S1)
    for (i,x) in enumerate(S1)
        for y in S2
            z::AutWord = x*y
            v::Vector{FGWord} = z(domain)
            if !in(v, seen)
                push!(seen, v)
                push!(result, z)
            end
        end
    end
    return result
end

function products_images(S1::Vector{AutWord}, S2::Vector{AutWord})
    result = Vector{Vector{FGWord}}()
    seen = Set{Vector{FGWord}}()
    n = length(S1)
    for (i,x) in enumerate(S1)
        z = x(domain)
        for y in S2
            v = y(z)
            if !in(v, seen)
                push!(seen, v)
                push!(result, v)
            end
        end
    end
    return result
end

function hashed_product{T}(image::T, B, images_dict::Dict{T, Int})
    n = size(B,1)
    column = zeros(Int,n)
    Threads.@threads for j in 1:n
        w = (B[j])(image)
        k = images_dict[w]
        k ≠ 0 || throw(ArgumentError(
            "($i,$j): $(x^-1)*$y don't seem to be supported on basis!"))
        column[j] = k
    end
    return column
end

function create_product_matrix(basis::Vector{AutWord}, images)
    n = length(basis)
    product_matrix = zeros(Int, (n, n));
    print("Creating hashtable of images...")
    @time images_dict = Dict{Vector{FGWord}, Int}(x => i
        for (i,x) in enumerate(images))

    for i in 1:n
        z = (inv(basis[i]))(domain)
        product_matrix[i,:] = hashed_product(z, basis, images_dict)
    end
    return product_matrix
end

function ΔandSDPconstraints(identity::AutWord, S::Vector{AutWord})

    println("Generating Balls of increasing radius...")
    @time B₁ = vcat([identity], S)
    @time B₂ = products(B₁,B₁);
    @show length(B₂)
    if length(B₂) != length(B₁)
        @time B₃ = products(B₁, B₂)
        @show length(B₃)
        if length(B₃) != length(B₂)
            @time B₄_images = products_images(B₁, B₃)
        else
            B₄_images = unique([f(domain) for f in B₃])
        end
    else
        B₃ = B₂
        B₄ = B₂
        B₄_images = unique([f(domain) for f in B₃])
    end

    @show length(B₄_images)
    # @assert length(B₄_images) == 3425657

    println("Creating product matrix...")
    @time pm = create_product_matrix(B₂, B₄_images)
    println("Creating sdp_constratints...")
    @time sdp_constraints = PropertyT.constraints_from_pm(pm)

    L_coeff = PropertyT.splaplacian_coeff(S, B₂, length(B₄_images))
    Δ = PropertyT.GroupAlgebraElement(L_coeff, Array{Int,2}(pm))

    return Δ, sdp_constraints
end

ΔandSDPconstraints(identity::AutWord, S::Vector{AutWord}, r::Int) = ΔandSDPconstraints(identity, S)

const symbols =  [FGSymbol("x₁",1), FGSymbol("x₂",1), FGSymbol("x₃",1), FGSymbol("x₄",1), FGSymbol("x₅",1), FGSymbol("x₆",1)]

const TOL=1e-8
const N = 4
const domain = Vector{FGWord}(symbols[1:N])

function cpuinfo_physicalcores()
    maxcore = -1
    for line in eachline("/proc/cpuinfo")
        if startswith(line, "core id")
            maxcore = max(maxcore, parse(Int, split(line, ':')[2]))
        end
    end
    maxcore < 0 && error("failure to read core ids from /proc/cpuinfo")
    return maxcore + 1
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--tol"
            help = "set numerical tolerance for the SDP solver (default: 1e-5)"
            arg_type = Float64
            default = 1e-5
        "--iterations"
            help = "set maximal number of iterations for the SDP solver (default: 20000)"
            arg_type = Int
            default = 20000
        "--upper-bound"
            help = "Set an upper bound for the spectral gap (default: Inf)"
            arg_type = Float64
            default = Inf
        "--cpus"
            help = "Set number of cpus used by solver (default: auto)"
            arg_type = Int
            required = false
        "-N"
            help = "Consider automorphisms of free group on N generators (default: N=3)"
            arg_type = Int
            default = 3
    end

    return parse_args(s)
end


# const name = "SYM$N"
# const upper_bound=factorial(N)-TOL^(1/5)
# S() = generating_set_of_Sym(N)

# name = "AutF$N"
# S() = generating_set_of_AutF(N)

function main()
    parsed_args = parse_commandline()

    tol = parsed_args["tol"]
    iterations = parsed_args["iterations"]

    solver = SCSSolver(eps=tol, max_iters=iterations, verbose=true, linearsolver=SCS.Indirect)

    N = parsed_args["N"]
    upper_bound = parsed_args["upper-bound"]

    name = "OutF$N"
    name = name*"-$(string(upper_bound))"
    S() = generating_set_of_OutF(N)

    if parsed_args["cpus"] ≠ nothing
        if parsed_args["cpus"] > cpuinfo_physicalcores()
            warn("Number of specified cores exceeds the physical core cound. Performance will suffer.")
        end
        Blas.set_num_threads(parsed_args["cpus"])
    end
    @time PropertyT.check_property_T(name, S, solver, upper_bound, tol, 2)
    return 0
end

main()