2018-09-05 09:56:04 +02:00
|
|
|
###############################################################################
|
|
|
|
#
|
2020-10-17 11:39:19 +02:00
|
|
|
# BlockDecomposition
|
2018-09-05 09:56:04 +02:00
|
|
|
#
|
|
|
|
###############################################################################
|
2017-10-09 19:45:28 +02:00
|
|
|
|
2020-10-17 11:39:19 +02:00
|
|
|
struct BlockDecomposition{T<:AbstractArray{Float64, 2}, GEl<:GroupElem, P<:Generic.Perm}
|
2018-08-20 03:54:03 +02:00
|
|
|
orbits::Vector{Vector{Int}}
|
2018-09-05 09:56:04 +02:00
|
|
|
preps::Dict{GEl, P}
|
2018-08-20 03:54:03 +02:00
|
|
|
Uπs::Vector{T}
|
2018-01-01 14:06:33 +01:00
|
|
|
dims::Vector{Int}
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2020-10-18 01:08:23 +02:00
|
|
|
function BlockDecomposition(RG::GroupRing, autS::Group; verbose=true)
|
2019-02-24 00:17:20 +01:00
|
|
|
verbose && @info "Decomposing basis of RG into orbits of" autS
|
2018-09-05 09:56:04 +02:00
|
|
|
@time orbs = orbit_decomposition(autS, RG.basis, RG.basis_dict)
|
|
|
|
@assert sum(length(o) for o in orbs) == length(RG.basis)
|
2019-02-24 00:17:20 +01:00
|
|
|
verbose && @info "The action has $(length(orbs)) orbits"
|
2017-08-27 18:58:01 +02:00
|
|
|
|
2019-02-24 00:17:20 +01:00
|
|
|
verbose && @info "Finding projections in the Group Ring of" autS
|
2019-01-11 06:37:24 +01:00
|
|
|
@time autS_mps = Projections.rankOne_projections(GroupRing(autS, collect(autS)))
|
2017-07-26 18:02:21 +02:00
|
|
|
|
2019-02-24 00:17:20 +01:00
|
|
|
verbose && @info "Finding AutS-action matrix representation"
|
2018-09-16 18:03:25 +02:00
|
|
|
@time preps = perm_reps(autS, RG.basis[1:size(RG.pm,1)], RG.basis_dict)
|
|
|
|
@time mreps = matrix_reps(preps)
|
2018-08-20 04:00:14 +02:00
|
|
|
|
2019-02-24 00:17:20 +01:00
|
|
|
verbose && @info "Computing the projection matrices Uπs"
|
2018-09-05 09:56:04 +02:00
|
|
|
@time Uπs = [orthSVD(matrix_repr(p, mreps)) for p in autS_mps]
|
2018-08-20 04:00:14 +02:00
|
|
|
|
2018-09-05 09:56:04 +02:00
|
|
|
multiplicities = size.(Uπs,2)
|
2019-07-05 18:55:25 +02:00
|
|
|
dimensions = [Int(p[one(autS)]*Int(order(autS))) for p in autS_mps]
|
2019-02-21 14:52:20 +01:00
|
|
|
if verbose
|
|
|
|
info_strs = ["",
|
|
|
|
lpad("multiplicities", 14) * " =" * join(lpad.(multiplicities, 4), ""),
|
|
|
|
lpad("dimensions", 14) * " =" * join(lpad.(dimensions, 4), "")
|
|
|
|
]
|
2019-02-24 00:17:20 +01:00
|
|
|
@info join(info_strs, "\n")
|
2019-02-21 14:52:20 +01:00
|
|
|
end
|
2018-09-05 09:56:04 +02:00
|
|
|
@assert dot(multiplicities, dimensions) == size(RG.pm,1)
|
2018-08-20 04:00:14 +02:00
|
|
|
|
2020-10-17 11:39:19 +02:00
|
|
|
return BlockDecomposition(orbs, preps, Uπs, dimensions)
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2020-10-18 01:08:23 +02:00
|
|
|
function decimate(od::BlockDecomposition; verbose=true)
|
2020-10-17 02:13:05 +02:00
|
|
|
nzros = [i for i in 1:length(od.Uπs) if !isempty(od.Uπs[i])]
|
2017-11-05 20:55:53 +01:00
|
|
|
|
2020-10-17 02:13:05 +02:00
|
|
|
Us = sparsify!.(od.Uπs, eps(Float64) * 1e4, verbose = verbose)[nzros]
|
|
|
|
#dimensions of the corresponding Uπs:
|
2018-09-05 09:56:04 +02:00
|
|
|
dims = od.dims[nzros]
|
2017-11-05 20:55:53 +01:00
|
|
|
|
2020-10-17 11:39:19 +02:00
|
|
|
return BlockDecomposition(od.orbits, od.preps, Array{Float64}.(Us), dims)
|
2017-06-22 15:11:14 +02:00
|
|
|
end
|
|
|
|
|
2018-09-05 10:40:08 +02:00
|
|
|
function orthSVD(M::AbstractMatrix{T}) where {T<:AbstractFloat}
|
2020-10-17 02:13:05 +02:00
|
|
|
fact = svd(convert(Matrix{T}, M))
|
|
|
|
M_rank = sum(fact.S .> maximum(size(M)) * eps(T))
|
|
|
|
return fact.U[:, 1:M_rank]
|
2018-09-05 10:40:08 +02:00
|
|
|
end
|
2018-08-20 03:57:33 +02:00
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
orbit_decomposition(
|
|
|
|
G::Group,
|
|
|
|
E::AbstractVector,
|
|
|
|
rdict = GroupRings.reverse_dict(E);
|
|
|
|
op = ^,
|
|
|
|
) = orbit_decomposition(collect(G), E, rdict; op=op)
|
2018-09-09 11:37:11 +02:00
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function orbit_decomposition(elts::AbstractVector{<:GroupElem}, E::AbstractVector, rdict=GroupRings.reverse_dict(E); op=^)
|
2018-09-09 11:37:11 +02:00
|
|
|
|
2019-01-11 06:32:09 +01:00
|
|
|
tovisit = trues(size(E));
|
2018-09-09 11:37:11 +02:00
|
|
|
orbits = Vector{Vector{Int}}()
|
|
|
|
|
|
|
|
orbit = zeros(Int, length(elts))
|
|
|
|
|
2019-01-11 06:32:09 +01:00
|
|
|
for i in eachindex(E)
|
2018-09-09 11:37:11 +02:00
|
|
|
if tovisit[i]
|
|
|
|
g = E[i]
|
2019-01-11 06:32:09 +01:00
|
|
|
Threads.@threads for j in eachindex(elts)
|
2020-10-10 12:55:02 +02:00
|
|
|
orbit[j] = rdict[op(g, elts[j])]
|
2018-09-09 11:37:11 +02:00
|
|
|
end
|
2019-01-11 06:32:09 +01:00
|
|
|
tovisit[orbit] .= false
|
2018-09-09 11:37:11 +02:00
|
|
|
push!(orbits, unique(orbit))
|
|
|
|
end
|
|
|
|
end
|
|
|
|
return orbits
|
|
|
|
end
|
|
|
|
|
2018-09-05 09:10:58 +02:00
|
|
|
###############################################################################
|
|
|
|
#
|
|
|
|
# Sparsification
|
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
|
|
|
dens(M::SparseMatrixCSC) = nnz(M)/length(M)
|
2019-01-11 06:32:09 +01:00
|
|
|
dens(M::AbstractArray) = count(!iszero, M)/length(M)
|
2018-01-02 03:19:57 +01:00
|
|
|
|
2020-10-17 02:13:05 +02:00
|
|
|
function sparsify!(M::SparseMatrixCSC{Tv,Ti}, tol=eps(Tv); verbose=false) where {Tv,Ti}
|
2018-09-05 09:06:00 +02:00
|
|
|
|
2018-09-05 09:10:58 +02:00
|
|
|
densM = dens(M)
|
2020-10-17 02:13:05 +02:00
|
|
|
droptol!(M, tol)
|
2018-09-05 09:06:00 +02:00
|
|
|
|
2020-10-17 02:13:05 +02:00
|
|
|
verbose && @info(
|
|
|
|
"Sparsified density:",
|
|
|
|
rpad(densM, 20),
|
|
|
|
" → ",
|
|
|
|
rpad(dens(M), 20),
|
|
|
|
" ($(nnz(M)) non-zeros)"
|
|
|
|
)
|
2018-01-01 14:06:33 +01:00
|
|
|
|
2018-09-05 09:10:58 +02:00
|
|
|
return M
|
|
|
|
end
|
2018-01-01 14:06:33 +01:00
|
|
|
|
2020-10-17 02:13:05 +02:00
|
|
|
function sparsify!(M::AbstractArray{T}, tol=eps(T); verbose=false) where T
|
2018-09-05 09:10:58 +02:00
|
|
|
densM = dens(M)
|
2020-10-17 02:13:05 +02:00
|
|
|
clamp_small!(M, tol)
|
2018-01-01 14:06:33 +01:00
|
|
|
|
2018-09-05 09:10:58 +02:00
|
|
|
if verbose
|
2019-02-21 21:10:36 +01:00
|
|
|
@info("Sparsifying $(size(M))-matrix... \n $(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(count(!iszero, M)) non-zeros)")
|
2018-09-05 09:10:58 +02:00
|
|
|
end
|
2018-01-01 14:06:33 +01:00
|
|
|
|
2018-09-05 09:10:58 +02:00
|
|
|
return sparse(M)
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
2018-09-05 09:10:58 +02:00
|
|
|
|
2020-10-17 02:13:05 +02:00
|
|
|
function clamp_small!(M::AbstractArray{T}, tol=eps(T)) where T
|
2019-01-10 04:48:30 +01:00
|
|
|
for n in eachindex(M)
|
2020-10-17 02:13:05 +02:00
|
|
|
if abs(M[n]) < tol
|
2019-01-10 04:48:30 +01:00
|
|
|
M[n] = zero(T)
|
|
|
|
end
|
|
|
|
end
|
|
|
|
return M
|
|
|
|
end
|
|
|
|
|
2019-01-11 06:32:09 +01:00
|
|
|
function sparsify(U::AbstractArray{T}, tol=eps(T); verbose=false) where T
|
|
|
|
return sparsify!(deepcopy(U), tol, verbose=verbose)
|
|
|
|
end
|
2018-09-09 11:37:11 +02:00
|
|
|
|
|
|
|
###############################################################################
|
|
|
|
#
|
|
|
|
# perm-, matrix-, representations
|
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
|
|
|
function perm_repr(g::GroupElem, E::Vector, E_dict)
|
2019-01-11 06:32:09 +01:00
|
|
|
p = Vector{Int}(undef, length(E))
|
2018-09-09 11:37:11 +02:00
|
|
|
for (i,elt) in enumerate(E)
|
2020-10-10 12:55:02 +02:00
|
|
|
p[i] = E_dict[elt^g]
|
2018-09-09 11:37:11 +02:00
|
|
|
end
|
|
|
|
return p
|
|
|
|
end
|
|
|
|
|
|
|
|
function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
|
2019-01-11 06:37:24 +01:00
|
|
|
elts = collect(G)
|
2018-09-09 11:37:11 +02:00
|
|
|
l = length(elts)
|
2019-10-30 23:58:10 +01:00
|
|
|
preps = Vector{Generic.Perm}(undef, l)
|
2018-09-09 11:37:11 +02:00
|
|
|
|
2020-06-23 16:15:01 +02:00
|
|
|
permG = SymmetricGroup(length(E))
|
2018-09-09 11:37:11 +02:00
|
|
|
|
|
|
|
Threads.@threads for i in 1:l
|
|
|
|
preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict), false)
|
|
|
|
end
|
|
|
|
|
|
|
|
return Dict(elts[i]=>preps[i] for i in 1:l)
|
|
|
|
end
|
|
|
|
|
|
|
|
function matrix_repr(x::GroupRingElem, mreps::Dict)
|
2019-01-11 06:32:09 +01:00
|
|
|
nzeros = findall(!iszero, x.coeffs)
|
2018-09-09 11:37:11 +02:00
|
|
|
return sum(x[i].*mreps[parent(x).basis[i]] for i in nzeros)
|
|
|
|
end
|
|
|
|
|
2019-10-30 23:58:10 +01:00
|
|
|
function matrix_reps(preps::Dict{T,Generic.Perm{I}}) where {T<:GroupElem, I<:Integer}
|
2018-09-09 11:37:11 +02:00
|
|
|
kk = collect(keys(preps))
|
2019-01-11 06:32:09 +01:00
|
|
|
mreps = Vector{SparseMatrixCSC{Float64, Int}}(undef, length(kk))
|
2018-09-09 11:37:11 +02:00
|
|
|
Threads.@threads for i in 1:length(kk)
|
|
|
|
mreps[i] = AbstractAlgebra.matrix_repr(preps[kk[i]])
|
|
|
|
end
|
|
|
|
return Dict(kk[i] => mreps[i] for i in 1:length(kk))
|
|
|
|
end
|
2019-01-11 07:54:36 +01:00
|
|
|
|
|
|
|
###############################################################################
|
|
|
|
#
|
|
|
|
# actions
|
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function Base.:^(y::GroupRingElem, g::GroupRingElem, op = ^)
|
2019-06-28 00:20:30 +02:00
|
|
|
res = parent(y)()
|
|
|
|
for elt in GroupRings.supp(g)
|
2020-10-10 12:55:02 +02:00
|
|
|
res += g[elt] * ^(y, elt, op)
|
2019-06-28 00:20:30 +02:00
|
|
|
end
|
|
|
|
return res
|
|
|
|
end
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function Base.:^(y::GroupRingElem, g::GroupElem, op = ^)
|
2019-06-25 17:36:37 +02:00
|
|
|
RG = parent(y)
|
|
|
|
result = zero(RG, eltype(y.coeffs))
|
2019-04-12 23:18:48 +02:00
|
|
|
|
2019-06-25 17:36:37 +02:00
|
|
|
for (idx, c) in enumerate(y.coeffs)
|
2020-06-23 16:13:42 +02:00
|
|
|
if !iszero(c)
|
2020-10-10 12:55:02 +02:00
|
|
|
result[op(RG.basis[idx], g)] = c
|
2019-01-11 07:54:36 +01:00
|
|
|
end
|
|
|
|
end
|
|
|
|
return result
|
|
|
|
end
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function Base.:^(
|
|
|
|
y::GroupRingElem{T,<:SparseVector},
|
|
|
|
g::GroupElem,
|
|
|
|
op = ^,
|
|
|
|
) where {T}
|
2019-06-25 17:36:37 +02:00
|
|
|
RG = parent(y)
|
2020-10-10 12:55:02 +02:00
|
|
|
index = [RG.basis_dict[op(RG.basis[idx], g)] for idx in y.coeffs.nzind]
|
2019-06-25 17:36:37 +02:00
|
|
|
|
|
|
|
result = GroupRingElem(sparsevec(index, y.coeffs.nzval, y.coeffs.n), RG)
|
|
|
|
|
|
|
|
return result
|
|
|
|
end
|
|
|
|
|
2020-06-23 16:13:42 +02:00
|
|
|
###############################################################################
|
|
|
|
#
|
|
|
|
# perm && WreathProductElems actions: MatAlgElem
|
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function Base.:^(A::MatAlgElem, p::Generic.Perm)
|
|
|
|
length(p.d) == size(A, 1) == size(A, 2) ||
|
|
|
|
throw("Can't act via $p on matrix of size $(size(A))")
|
2019-06-25 17:36:37 +02:00
|
|
|
result = similar(A)
|
2020-10-10 12:55:02 +02:00
|
|
|
@inbounds for i = 1:size(A, 1)
|
|
|
|
for j = 1:size(A, 2)
|
|
|
|
result[p[i], p[j]] = A[i, j] # action by permuting rows and colums/conjugation
|
2019-06-25 17:36:37 +02:00
|
|
|
end
|
2019-06-28 01:34:37 +02:00
|
|
|
end
|
2019-06-25 17:36:37 +02:00
|
|
|
return result
|
|
|
|
end
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function Base.:^(A::MatAlgElem, g::WreathProductElem{N}) where {N}
|
2019-06-25 17:37:08 +02:00
|
|
|
# @assert N == size(A,1) == size(A,2)
|
2020-10-10 12:55:02 +02:00
|
|
|
flips = ntuple(i -> (g.n[i].d[1] == 1 && g.n[i].d[2] == 2 ? 1 : -1), N)
|
2019-06-25 17:37:08 +02:00
|
|
|
result = similar(A)
|
2019-06-28 09:15:01 +02:00
|
|
|
R = base_ring(parent(A))
|
|
|
|
tmp = R(1)
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
@inbounds for i = 1:size(A, 1)
|
|
|
|
for j = 1:size(A, 2)
|
2019-06-30 11:58:32 +02:00
|
|
|
x = A[i, j]
|
2020-10-10 12:55:02 +02:00
|
|
|
if flips[i] * flips[j] == 1
|
2019-06-30 11:58:32 +02:00
|
|
|
result[g.p[i], g.p[j]] = x
|
2019-06-28 09:15:01 +02:00
|
|
|
else
|
2019-06-30 11:58:32 +02:00
|
|
|
result[g.p[i], g.p[j]] = -x
|
2019-06-28 09:15:01 +02:00
|
|
|
end
|
2019-06-25 17:37:08 +02:00
|
|
|
end
|
|
|
|
end
|
|
|
|
return result
|
2019-01-11 07:54:36 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
###############################################################################
|
|
|
|
#
|
2020-06-23 16:13:42 +02:00
|
|
|
# perm && WreathProductElems actions: Automorphism
|
2019-01-11 07:54:36 +01:00
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
2020-10-10 12:55:02 +02:00
|
|
|
function Base.:^(a::Automorphism, g::GroupElem)
|
2020-06-23 16:13:42 +02:00
|
|
|
Ag = parent(a)(g)
|
|
|
|
Ag_inv = inv(Ag)
|
|
|
|
res = append!(Ag, a, Ag_inv)
|
|
|
|
return Groups.freereduce!(res)
|
|
|
|
end
|
|
|
|
|
|
|
|
(A::AutGroup)(p::Generic.Perm) = A(Groups.AutSymbol(p))
|
|
|
|
|
|
|
|
function (A::AutGroup)(g::WreathProductElem)
|
2019-01-11 07:54:36 +01:00
|
|
|
isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
|
2020-10-10 12:55:02 +02:00
|
|
|
parent(g).P.n == length(A.objectGroup.gens) ||
|
|
|
|
throw("No natural embedding of $(parent(g)) into $A")
|
2019-07-05 18:55:25 +02:00
|
|
|
elt = one(A)
|
|
|
|
Id = one(parent(g.n.elts[1]))
|
2020-10-10 12:55:02 +02:00
|
|
|
for i = 1:length(g.p.d)
|
2020-06-23 16:13:42 +02:00
|
|
|
if g.n.elts[i] != Id
|
|
|
|
push!(elt, Groups.flip(i))
|
|
|
|
end
|
|
|
|
end
|
|
|
|
push!(elt, Groups.AutSymbol(g.p))
|
2019-01-11 07:54:36 +01:00
|
|
|
return elt
|
|
|
|
end
|
|
|
|
|
|
|
|
|
2020-06-23 16:13:42 +02:00
|
|
|
# fallback:
|
|
|
|
Base.one(p::Generic.Perm) = Perm(length(p.d))
|