GroupsWithPropertyT/OrbitDecomposition.jl

299 lines
8.2 KiB
Julia
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

push!(LOAD_PATH, "./")
using Nemo
using Groups
using WreathProducts
using GroupRings
using PropertyT
import Nemo.elements
using JLD
using ProgressMeter
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
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
function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Group, S::Vector{T}; radius=2)
isdir(name) || mkdir(name)
info(logger, "Generating ball of radius 4")
@time E4, sizes = generate_balls(S, G(), radius=2*radius);
info(logger, "Reverse dict")
@time E_dict = GroupRings.reverse_dict(E4)
info(logger, "Product matrix")
@time pm = GroupRings.create_pm(E4, E_dict, sizes[radius], twisted=true)
RG = GroupRing(G, E4, E_dict, pm)
Δ = PropertyT.splaplacian(RG, S)
@assert GroupRings.augmentation(Δ) == 0
save(joinpath(name, "delta.jld"), "delta", Δ.coeffs)
save(joinpath(name, "pm.jld"), "pm", pm)
N = length(G.objectGroup.gens)
info(logger, "WreathProduct")
@time BN = WreathProduct(Nemo.FiniteField(2,1, "a")[1], PermutationGroup(N))
info(logger, "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)
info(logger, "Action matrices")
E2 = E4[1:sizes[radius]]
@time BNactionE_mreps = action_mreps(BN, E2, E_dict)
info(logger, "Projections")
@time BN_mps = rankOne_projections(BN);
@time π_E_projections = @showprogress [Cstar_repr(p, BNactionE_mreps) for p in BN_mps]
info(logger, "Uπs...")
@time Uπs = Uπ_matrices(π_E_projections);
multiplicities = [size(U,2) for U in Uπs];
info(logger, "multiplicities = $multiplicities")
dimensions = [Int(p[BN()]*Int(order(BN))) for p in BN_mps];
info(logger, "dimensions = $dimensions")
@assert dot(multiplicities, dimensions) == sizes[radius]
save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "dims", dimensions)
return 0
end