using Groups
using ProgressMeter

#=
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))
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))
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))
end


function products(S1::Vector{AutWord}, S2::Vector{AutWord})
    result = Vector{AutWord}()
    seen = Set{Vector{FGWord}}()
    n = length(S1)
    p = Progress(n, 1, "Computing complete products...", 50)
    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
        next!(p)
    end
    return result
end

function products_images(S1::Vector{AutWord}, S2::Vector{AutWord})
    result = Vector{Vector{FGWord}}()
    seen = Set{Vector{FGWord}}()
    n = length(S1)

    p = Progress(n, 1, "Computing images of elts in B₄...", 50)
    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
        next!(p)
    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))

    p = Progress(n, 1, "Computing product matrix in basis...", 50)
    for i in 1:n
        z = (inv(basis[i]))(domain)
        product_matrix[i,:] = hashed_product(z, basis, images_dict)
        next!(p)
    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 = PropertyT.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


using GroupAlgebras
using PropertyT

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])
const ID = one(AutWord)

# 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)

name = "OutF$N"
S() = generating_set_of_OutF(N)
const upper_bound=0.05

BLAS.set_num_threads(4)
@time check_property_T(name, ID, S; verbose=true, tol=TOL, upper_bound=upper_bound)