Initial Commit

AutFN.jl

@@ -0,0 +1,183 @@
+using JLD
+using JuMP
+import SCS: SCSSolver
+import Mosek: MosekSolver
+
+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 = create_product_matrix(B₂, B₄_images)
+    println("Creating sdp_constratints...")
+    @time sdp_constraints = constraints_from_pm(pm)
+
+    L_coeff = splaplacian_coeff(S, B₂, length(B₄_images))
+    Δ = GroupAlgebraElement(L_coeff, Array{Int,2}(pm))
+
+    return Δ, sdp_constraints
+end
+
+
+@everywhere push!(LOAD_PATH, "./")
+using GroupAlgebras
+include("property(T).jl")
+
+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)

GroupAlgebras.jl

@@ -0,0 +1,133 @@
+module GroupAlgebras
+
+import Base: convert, show, isequal, ==
+import Base: +, -, *, //
+import Base: size, length, norm, rationalize
+
+export GroupAlgebraElement
+
+
+immutable GroupAlgebraElement{T<:Number}
+    coefficients::AbstractVector{T}
+    product_matrix::Array{Int,2}
+    # basis::Array{Any,1}
+
+    function GroupAlgebraElement(coefficients::AbstractVector,
+        product_matrix::Array{Int,2})
+
+        size(product_matrix, 1) == size(product_matrix, 2) ||
+            throw(ArgumentError("Product matrix has to be square"))
+        new(coefficients, product_matrix)
+    end
+end
+
+# GroupAlgebraElement(c,pm,b) = GroupAlgebraElement(c,pm)
+GroupAlgebraElement{T}(c::AbstractVector{T},pm) = GroupAlgebraElement{T}(c,pm)
+
+convert{T<:Number}(::Type{T}, X::GroupAlgebraElement) =
+    GroupAlgebraElement(convert(AbstractVector{T}, X.coefficients), X.product_matrix)
+
+show{T}(io::IO, X::GroupAlgebraElement{T}) = print(io,
+    "Element of Group Algebra over $T of length $(length(X)):\n $(X.coefficients)")
+
+
+function isequal{T, S}(X::GroupAlgebraElement{T}, Y::GroupAlgebraElement{S})
+    if T != S
+        warn("Comparing elements with different coefficients Rings!")
+    end
+    X.product_matrix == Y.product_matrix || return false
+    X.coefficients == Y.coefficients || return false
+    return true
+end
+
+(==)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = isequal(X,Y)
+
+function add{T<:Number}(X::GroupAlgebraElement{T}, Y::GroupAlgebraElement{T})
+    X.product_matrix == Y.product_matrix || throw(ArgumentError(
+    "Elements don't seem to belong to the same Group Algebra!"))
+    return GroupAlgebraElement(X.coefficients+Y.coefficients, X.product_matrix)
+end
+
+function add{T<:Number, S<:Number}(X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{S})
+    warn("Adding elements with different base rings!")
+    return GroupAlgebraElement(+(promote(X.coefficients, Y.coefficients)...),
+    X.product_matrix)
+end
+
+(+)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = add(X,Y)
+(-)(X::GroupAlgebraElement) = GroupAlgebraElement(-X.coefficients, X.product_matrix)
+(-)(X::GroupAlgebraElement, Y::GroupAlgebraElement) = add(X,-Y)
+
+function algebra_multiplication{T<:Number}(X::AbstractVector{T}, Y::AbstractVector{T}, pm::Array{Int,2})
+    result = zeros(X)
+    for (j,y) in enumerate(Y)
+        if y != zero(T)
+            for (i, index) in enumerate(pm[:,j])
+                if X[i] != zero(T)
+                    index == 0 && throw(ArgumentError("The product don't seem to belong to the span of basis!"))
+                    result[index] += X[i]*y
+                end
+            end
+        end
+    end
+    return result
+end
+
+function group_star_multiplication{T<:Number}(X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{T})
+    X.product_matrix == Y.product_matrix || ArgumentError(
+    "Elements don't seem to belong to the same Group Algebra!")
+    result = algebra_multiplication(X.coefficients, Y.coefficients, X.product_matrix)
+    return GroupAlgebraElement(result, X.product_matrix)
+end
+
+function group_star_multiplication{T<:Number, S<:Number}(
+    X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{S})
+    S == T || warn("Multiplying elements with different base rings!")
+    return group_star_multiplication(promote(X,Y)...)
+end
+
+(*){T<:Number, S<:Number}(X::GroupAlgebraElement{T},
+    Y::GroupAlgebraElement{S}) = group_star_multiplication(X,Y);
+
+(*){T<:Number}(a::T, X::GroupAlgebraElement{T}) = GroupAlgebraElement(
+    a*X.coefficients, X.product_matrix)
+
+function scalar_multiplication{T<:Number, S<:Number}(a::T,
+    X::GroupAlgebraElement{S})
+    promote_type(T,S) == S || warn("Scalar and coefficients are in different rings! Promoting result to $(promote_type(T,S))")
+    return GroupAlgebraElement(a*X.coefficients, X.product_matrix)
+end
+
+(*){T<:Number}(a::T,X::GroupAlgebraElement) = scalar_multiplication(a, X)
+
+//{T<:Rational, S<:Rational}(X::GroupAlgebraElement{T}, a::S) =
+    GroupAlgebraElement(X.coefficients//a, X.product_matrix)
+
+//{T<:Rational, S<:Integer}(X::GroupAlgebraElement{T}, a::S) =
+    X//convert(T,a)
+
+length(X::GroupAlgebraElement) = length(X.coefficients)
+size(X::GroupAlgebraElement) = size(X.coefficients)
+
+function norm(X::GroupAlgebraElement, p=2)
+    if p == 1
+        return sum(abs(X.coefficients))
+    elseif p == Inf
+        return max(abs(X.coefficients))
+    else
+        return norm(X.coefficients, p)
+    end
+end
+
+ɛ(X::GroupAlgebraElement) = sum(X.coefficients)
+
+function rationalize{T<:Integer, S<:Number}(
+    ::Type{T}, X::GroupAlgebraElement{S}; tol=eps(S))
+    v = rationalize(T, X.coefficients, tol=tol)
+    return GroupAlgebraElement(v, X.product_matrix)
+end
+
+end

SL3Z.jl

@@ -0,0 +1,142 @@
+using JLD
+using JuMP
+import Primes: isprime
+import SCS: SCSSolver
+import Mosek: MosekSolver
+
+using Mods
+
+using Groups
+
+function SL_generatingset(n::Int)
+
+    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
+
+    S = [E(i,j,N=n) for (i,j) in indexing];
+    S = vcat(S, [convert(Array{Int,2},x') for x in S]);
+    S = vcat(S, [convert(Array{Int,2},inv(x)) for x in S]);
+    return unique(S)
+end
+
+function E(i::Int, j::Int; val=1, N::Int=3, mod=Inf)
+   @assert i≠j
+   m = eye(Int, N)
+   m[i,j] = val
+   if mod == Inf
+       return m
+   else
+       return [Mod(x,mod) for x in m]
+   end
+end
+
+function cofactor(i,j,M)
+    z1 = ones(Bool,size(M,1))
+    z1[i] = false
+
+    z2 = ones(Bool,size(M,2))
+    z2[j] = false
+
+    return M[z1,z2]
+end
+
+import Base.LinAlg.det
+
+function det(M::Array{Mod,2})
+    if size(M,1) ≠ size(M,2)
+        d = Mod(0,M[1,1].mod)
+    elseif size(M,1) == 2
+        d =  M[1,1]*M[2,2] - M[1,2]*M[2,1]
+    else
+        d = zero(eltype(M))
+        for i in 1:size(M,1)
+            d += (-1)^(i+1)*M[i,1]*det(cofactor(i,1,M))
+        end
+    end
+#     @show (M, d)
+    return d
+end
+
+function adjugate(M)
+    K = similar(M)
+    for i in 1:size(M,1), j in 1:size(M,2)
+        K[j,i] = (-1)^(i+j)*det(cofactor(i,j,M))
+    end
+    return K
+end
+
+import Base: inv, one, zero, *
+
+one(::Type{Mod}) = 1
+zero(::Type{Mod}) = 0
+zero(x::Mod) = Mod(x.mod)
+
+function inv(M::Array{Mod,2})
+    d = det(M)
+    d ≠ 0*d || thow(ArgumentError("Matrix is not invertible!"))
+    return inv(det(M))*adjugate(M)
+    return adjugate(M)
+end
+
+function SL_generatingset(n::Int, p::Int)
+    (p > 1 && n > 1) || throw(ArgumentError("Both n and p should be integers!"))
+    isprime(p) || throw(ArgumentError("p should be a prime number!"))
+
+    indexing = [(i,j) for i in 1:n for j in 1:n if i≠j]
+    S = [E(i,j, N=n, mod=p) for (i,j) in indexing]
+    S = vcat(S, [inv(s) for s in S])
+    S = vcat(S, [permutedims(x, [2,1]) for x in S]);
+
+    return unique(S)
+end
+
+function products{T}(U::AbstractVector{T}, V::AbstractVector{T})
+    result = Vector{T}()
+    for u in U
+        for v in V
+            push!(result, u*v)
+        end
+    end
+    return unique(result)
+end
+
+function ΔandSDPconstraints(identity, S)
+    B₁ = vcat([identity], S)
+    B₂ = products(B₁, B₁);
+    B₃ = products(B₁, B₂);
+    B₄ = products(B₁, B₃);
+    @assert B₄[1:length(B₂)] == B₂
+
+    product_matrix = create_product_matrix(B₄,length(B₂));
+    sdp_constraints = constraints_from_pm(product_matrix, length(B₄))
+    L_coeff = splaplacian_coeff(S, B₂, length(B₄));
+    Δ = GroupAlgebraElement(L_coeff, product_matrix)
+
+    return Δ, sdp_constraints
+end
+
+
+
+
+
+
+@everywhere push!(LOAD_PATH, "./")
+using GroupAlgebras
+include("property(T).jl")
+
+const N = 3
+
+const name = "SL$(N)Z"
+const ID = eye(Int, N)
+S() = SL_generatingset(N)
+const upper_bound=0.27
+
+
+# const p = 7
+# const upper_bound=0.738 # (N,p) = (3,7)
+
+# const name = "SL($N,$p)"
+# const ID = [Mod(x,p) for x in eye(Int,N)]
+# S() = SL_generatingset(N, p)
+
+BLAS.set_num_threads(4)
+@time check_property_T(name, ID, S; verbose=true, tol=1e-10, upper_bound=upper_bound)