GroupsWithPropertyT/OrbitDecomposition.jl

303 lines
8.0 KiB
Julia
Raw Normal View History

push!(LOAD_PATH, "./")
using Nemo
using Groups
using WreathProducts
using GroupRings
2017-06-06 16:20:44 +02:00
using PropertyT
import Nemo.elements
using JLD
using ProgressMeter
2017-06-06 16:19:17 +02:00
import Base.convert
convert(::Type{Int}, x::Nemo.fmpz) = x.d
function products{T<:GroupElem}(X::AbstractVector{T}, Y::AbstractVector{T})
result = Vector{T}()
seen = Set{T}()
@showprogress for x in X
for y in Y
z = x*y
if !in(z, seen)
push!(seen, z)
push!(result, z)
end
end
end
return result
end
function generate_balls{T<:GroupElem}(S::Vector{T}, Id::T; radius=2)
sizes = Vector{Int}()
S = unshift!(S, Id)
B = [Id]
for i in 1:radius
B = products(B, S);
push!(sizes, length(B))
end
return B, sizes
end
function elements(F::Nemo.FqNmodFiniteField)
deg = Int(degree(F))
char = Int(characteristic(F))
z = (gen(F)^i for i in 0:deg-1)
crtsn_prd = Base.product([0:char-1 for i in 1:deg]...)
function _it()
for crt in crtsn_prd
g = sum([b*a for (a,b) in zip(z,crt)])
produce(g)
end
end
return Task(_it)
end
function AutFG_emb(A::AutGroup, g::WreathProductElem)
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(FN)")
parent(g).P.n == length(A.objectGroup.gens) || throw("No natural action of $(parent(g)) on $A")
powers = [(elt == parent(elt)() ? 0: 1) for elt in g.n.elts]
elt = reduce(*, [A(Groups.flip_autsymbol(i))^pow for (i,pow) in enumerate(powers)])
Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])
return elt
end
function (g::WreathProductElem)(a::AutGroupElem)
g = AutFG_emb(parent(a),g)
return g*a*g^-1
end
function orbit_decomposition(G::Nemo.Group, E::Vector, rdict=GroupRings.reverse_dict(E))
elts = collect(elements(G))
tovisit = trues(E);
orbits = Vector{Set{Int}}()
@showprogress "Orbit decomposition... " for i in 1:endof(E)
if tovisit[i]
orbit = Set{Int}()
a = E[i]
for g in elts
idx = rdict[g(a)]
tovisit[idx] = false
push!(orbit,idx)
end
push!(orbits, orbit)
end
end
return orbits
end
function matrix_repr(g::WreathProductElem, E, E_dict)
rep_matrix = spzeros(Int, length(E), length(E))
for (i,e) in enumerate(E)
j = E_dict[g(e)]
rep_matrix[i,j] = 1
end
return rep_matrix
end
function action_mreps(G::Nemo.Group, E, E_dict)
result = @showprogress [matrix_repr(g, E, E_dict) for g in elements(G)]
return result
end
function chars(G::PermutationGroup)
permtype_unsorted(σ::Nemo.perm) = [length(c) for c in cycles(σ)]
permtype(σ::Nemo.perm) = sort(permtype_unsorted(σ))
χ_id(σ::Nemo.perm) = 1
χ_sgn(σ::Nemo.perm) = (-1)^parity(σ)
function χ_reg(σ::Nemo.perm)
fixed_points = countnz([(x == y? 1 : 0) for (x,y) in enumerate(σ.d)])
return fixed_points - 1
end
χ_regsgn(σ::Nemo.perm) = (-1)^parity(σ)*χ_reg(σ)
function χ_regviaS3(σ::Nemo.perm)
@assert parent(σ).n == 4
t = permtype(σ)
if t == [1,1,1,1]
result = 2
elseif t == [2,2]
result = 2
elseif t == [1,3]
result = -1
else
result = 0
end
return result
end
chars = [χ_id, χ_sgn, χ_regviaS3, χ_reg, χ_regsgn]
if G.n == 1
return chars[1:1]
elseif G.n == 2
return chars[1:2]
elseif G.n == 3
return [chars[1:2]..., chars[4]]
elseif G.n == 4
return chars[1:5]
else
throw("Characters for $G unknown!")
end
end
function epsilon(i, g::DirectProducts.DirectProductGroupElem)
return reduce(*, 1, ((-1)^isone(g.elts[j]) for j in 1:i))
end
function central_projection(RG::GroupRing, char::Function, T::Type=Rational{Int})
result = RG(T)
for g in RG.basis
result[g] = char(g)
end
return convert(T, char(RG.group())//Int(order(RG.group))*result)
end
function rankOne_projections(G::PermutationGroup, T::Type=Rational{Int})
RG = GroupRing(G)
projections = [central_projection(RG, χ, T) for χ in chars(G)]
if G.n == 1 || G.n == 2
return projections
elseif G.n == 3
rankone_projs = [
projections[1],
projections[2],
1//2*(one(RG) - RG(RG.group([2,1,3])))*projections[3]
]
return rankone_projs
elseif G.n == 4
rankone_projs = [
projections[1],
projections[2],
1//2*(one(RG) - RG(RG.group([2,1,3,4])))*projections[3],
1//2*(one(RG) - RG(RG.group([2,1,3,4])))*projections[4],
1//2*(one(RG) + RG(RG.group([2,1,3,4])))*projections[5]]
return rankone_projs
else
throw("Rank-one projections for $G unknown!")
end
end
function rankOne_projections(BN::WreathProducts.WreathProduct, T::Type=Rational{Int})
N = BN.P.n
# projections as elements of the group rings RSₙ
SNprojs_nc = [rankOne_projections(PermutationGroup(i), T) for i in 1:N]
# embedding into group ring of BN
RBN = GroupRing(BN)
RFFFF_projs = [central_projection(GroupRing(BN.N), g->epsilon(i,g), T)
for i in 0:BN.P.n]
Qs = [RBN(q, g -> BN(g)) for q in RFFFF_projs]
function incl(k::Int, g::perm, WP::WreathProduct=BN)
@assert length(g.d) + k <= WP.P.n
arr = [1:k; g.d .+ k; (length(g.d)+k+1):WP.P.n]
return WP(WP.P(arr))
end
all_projs=[Qs[1]*RBN(p, g-> incl(0,g)) for p in SNprojs_nc[N]]
for i in 1:N-1
Sk_first = [RBN(p, g->incl(0,g)) for p in SNprojs_nc[i]]
Sk_last = [RBN(p, g->incl(i,g)) for p in SNprojs_nc[N-i]]
append!(all_projs, [Qs[i+1]*p1*p2
for (p1,p2) in Base.product(Sk_first,Sk_last)])
end
append!(all_projs, [Qs[N+1]*RBN(p, g-> incl(0,g)) for p in SNprojs_nc[N]])
return all_projs
end
2017-06-06 16:19:46 +02:00
function Cstar_repr(x::GroupRingElem, matrix_reps)
res = zeros(matrix_reps[1])
for i in 1:length(parent(x).basis)
res += x.coeffs[i]*matrix_reps[i]
end
return res
end
function orthSVD(M::AbstractMatrix)
M = full(M)
# matrixRank = rank(M)
fact = svdfact(M)
sings = fact[:S]
# @show sings[sings.>1e-14]
M_rank = sum(fact[:S] .> maximum(size(M))*eps(eltype(fact[:S])))
Ufactor = fact[:U]
return Ufactor[:,1:M_rank]
end
function Uπ_matrices(P_matrices; orth=orthSVD)
U_p_matrices = Vector{Array{Float64,2}}()
@showprogress "Computing Uπ mats..." for (i,p_mat) in enumerate(P_matrices)
U_p = orth(p_mat)
push!(U_p_matrices, U_p)
end
return U_p_matrices
end
2017-06-06 16:20:20 +02:00
function compute_orbit_data(name::String, N::Int)
isdir(name) || mkdir(name)
SOutFN = AutGroup(FreeGroup(N), special=true, outer=true)
S = generators(SOutFN);
S = [S; [inv(s) for s in S]]
@show S
info("Generating ball of radius 4")
@time E4, sizes = generate_balls(S, SOutFN(), radius=4);
info("Reverse dict")
@time E_dict = GroupRings.reverse_dict(E4)
BN = WreathProduct(Nemo.FiniteField(2,1, "a")[1], PermutationGroup(N))
if !isfile(joinpath(name, "orbits.jld"))
info("Decomposing E into orbits of B$(N)")
@time orbs = orbit_decomposition(BN, E4, E_dict)
@assert sum(length(o) for o in orbs) == length(E4)
save(joinpath(name, "orbits.jld"), "orbits", orbs)
end
info("Action matrices")
pm = GroupRings.create_pm(E4, E_dict, sizes[2], twisted=true)
RSOutFN = GroupRing(SOutFN, E4, E_dict, pm)
PropertyT.splaplacian(RSOutFN, S)
E2 = E4[1:sizes[2]]
@time BNactionE_mreps = action_mreps(BN, E2, E_dict)
save(joinpath(name, "matrix_reps.jld"), "matrix_reps", BNactionE_mreps)
info("Projections")
@time BN_mps = rankOne_projections(BN);
@time π_E_projections = @showprogress [Cstar_repr(p, BNactionE_mreps) for p in BN_mps]
info("Uπs...")
@time Uπs = Uπ_matrices(π_E_projections);
@show multiplicities = [size(U,2) for U in Uπs];
@show dimensions = [Int(p[BN()]*Int(order(BN))) for p in BN_mps];
@assert dot(multiplicities, dimensions) == sizes[2]
save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "dims", dimensions)
end