GroupsWithPropertyT/OrbitDecomposition.jl

push!(LOAD_PATH, "./")

using Nemo
using Groups
using WreathProducts

using GroupRings
using PropertyT

import Nemo.elements

using JLD

include("Projections.jl")

###############################################################################
#
#  Iterator protocol for Nemo.FinField
#
###############################################################################

type FFEltsIter{T<:Nemo.FinField}
    all::Int
    field::T

    function FFEltsIter(F::T)
        return new(Int(characteristic(F)^degree(F)), F)
    end
end
FFEltsIter{T<:Nemo.FinField}(F::T) = FFEltsIter{T}(F)

import Base: start, next, done, eltype, length

Base.start(A::FFEltsIter) = (zero(A.field), 0)
Base.next(A::FFEltsIter, state) = next_ffelem(state...)
Base.done(A::FFEltsIter, state) = state[2] >= A.all
Base.eltype(::Type{FFEltsIter}) = elem_type(A.field)
Base.length(A::FFEltsIter) = A.all

function next_ffelem(f::Nemo.FinFieldElem, c::Int)
    if c == 0
        return (f, (f, 1))
    elseif c == 1
        f = one(parent(f))
        return (f, (f, 2))
    else
        f = gen(parent(f))*f
        return (f, (f, c+1))
    end
end

import Nemo.elements
elements(F::Nemo.FinField) = FFEltsIter(F)

###############################################################################
#
#  Action of Premutations on Nemo.MatElem
#
###############################################################################

function (p::Nemo.perm)(A::Nemo.MatElem)
    length(p.d) == A.r == A.c || throw("Can't act via $p on matrix of size ($(A.r), $(A.c))")
    R = parent(A)
    inv_p = inv(p)
    return R(Nemo.matrix_repr(p))*A*R(Nemo.matrix_repr(inv_p))
end

###############################################################################
#
#  Orbit stuff
#
###############################################################################

function orbit_decomposition(G::Nemo.Group, E::Vector, rdict=GroupRings.reverse_dict(E))

    elts = collect(elements(G))

    tovisit = trues(E);
    orbits = Vector{Set{Int}}()

    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 orbit_spvector(vect::AbstractVector, orbits)
    orb_vector = spzeros(length(orbits))

    for (i,o) in enumerate(orbits)
        k = vect[collect(o)]
        val = k[1]
        @assert all(k .== val)
        orb_vector[i] = val
    end

    return orb_vector
end

function orbit_constraint(cnstrs::Vector{Vector{Vector{Int64}}}, n)
    result = spzeros(n,n)
    for cnstr in cnstrs
        for p in cnstr
            result[p[1],p[2]] += 1.0
        end
    end
    return 1/length(cnstrs)*result
end

###############################################################################
#
#  Matrix- and C*-representations
#
###############################################################################

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 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)
    fact = svdfact(M)
    sings = fact[:S]
    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}}()
    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}, AutS; radius=2)
   isdir(name) || mkdir(name)

   info(logger, "Generating ball of radius $(2*radius)")
   if isa(G, Nemo.Ring)
      Id = one(G)
   else
      Id = G()
   end
   @time E4, sizes = Groups.generate_balls(S, Id, radius=2*radius);
   info(logger, "Balls of sizes $sizes.")
   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"), "Δ", Δ.coeffs)
   save(joinpath(name, "pm.jld"), "pm", pm)

   info(logger, "Decomposing E into orbits of $(AutS)")
   @time orbs = orbit_decomposition(AutS, 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 AutS_matrixreps = [matrix_repr(g, E2, E_dict) for g in elements(AutS)]

   info(logger, "Projections")
   @time AutS_mps = rankOne_projections(AutS);

   @time π_E_projections = [Cstar_repr(p, AutS_matrixreps) for p in AutS_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[AutS()]*Int(order(AutS))) for p in AutS_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