2017-06-22 14:12:35 +02:00
|
|
|
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)
|
|
|
|
|
|
|
|
###############################################################################
|
|
|
|
#
|
|
|
|
# Orbit stuff
|
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
|
|
|
function orbit_decomposition(G::Nemo.Group, E::Vector, rdict=GroupRings.reverse_dict(E))
|
|
|
|
|
|
|
|
elts = collect(elements(G))
|
|
|
|
|
|
|
|
tovisit = trues(E);
|
|
|
|
orbits = Vector{Vector{Int}}()
|
|
|
|
|
|
|
|
for i in 1:endof(E)
|
|
|
|
if tovisit[i]
|
2017-08-27 18:32:19 +02:00
|
|
|
orbit = zeros(Int, length(elts))
|
2017-06-22 14:12:35 +02:00
|
|
|
a = E[i]
|
2017-08-27 18:32:19 +02:00
|
|
|
Threads.@threads for i in 1:length(elts)
|
|
|
|
orbit[i] = rdict[elts[i](a)]
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
2017-08-27 18:32:19 +02:00
|
|
|
tovisit[orbit] = false
|
2017-06-22 14:12:35 +02:00
|
|
|
push!(orbits, unique(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
|
|
|
|
|
2017-11-08 09:25:40 +01:00
|
|
|
function orbit_constraint(constraints::Vector{Vector{Tuple{Int,Int}}}, n)
|
2017-06-22 14:12:35 +02:00
|
|
|
result = spzeros(n,n)
|
|
|
|
for cnstr in constraints
|
|
|
|
for p in cnstr
|
2017-11-08 09:25:40 +01:00
|
|
|
result[p[2], p[1]] += 1.0/length(constraints)
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
end
|
2017-11-08 09:25:40 +01:00
|
|
|
return result
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
|
|
|
###############################################################################
|
|
|
|
#
|
2017-11-08 09:28:04 +01:00
|
|
|
# Matrix-, Permutation- and C*-representations
|
2017-06-22 14:12:35 +02:00
|
|
|
#
|
|
|
|
###############################################################################
|
|
|
|
|
2017-11-08 09:28:04 +01:00
|
|
|
function matrix_repr(p::perm)
|
|
|
|
N = parent(p).n
|
|
|
|
return sparse(1:N, p.d, [1.0 for _ in 1:N])
|
|
|
|
end
|
|
|
|
|
|
|
|
function matrix_reps{T<:GroupElem}(preps::Dict{T,perm})
|
|
|
|
kk = collect(keys(preps))
|
|
|
|
mreps = Vector{SparseMatrixCSC{Float64, Int}}(length(kk))
|
|
|
|
Threads.@threads for i in 1:length(kk)
|
|
|
|
mreps[i] = matrix_repr(preps[kk[i]])
|
|
|
|
end
|
|
|
|
return Dict(kk[i] => mreps[i] for i in 1:length(kk))
|
|
|
|
end
|
2017-08-27 18:35:35 +02:00
|
|
|
|
2017-11-08 09:28:04 +01:00
|
|
|
function perm_repr(g::GroupElem, E::Vector, E_dict)
|
|
|
|
p = Vector{Int}(length(E))
|
|
|
|
for (i,elt) in enumerate(E)
|
|
|
|
p[i] = E_dict[g(elt)]
|
2017-08-27 18:35:35 +02:00
|
|
|
end
|
2017-11-08 09:28:04 +01:00
|
|
|
return p
|
2017-08-27 18:35:35 +02:00
|
|
|
end
|
|
|
|
|
2017-11-08 09:28:04 +01:00
|
|
|
function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E))
|
|
|
|
elts = collect(elements(G))
|
2017-08-27 18:35:35 +02:00
|
|
|
l = length(elts)
|
|
|
|
preps = Vector{Nemo.perm}(l)
|
|
|
|
|
2017-11-08 09:28:04 +01:00
|
|
|
permG = Nemo.PermutationGroup(length(E))
|
2017-08-27 18:35:35 +02:00
|
|
|
|
|
|
|
Threads.@threads for i in 1:l
|
2017-11-08 09:28:04 +01:00
|
|
|
preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict))
|
2017-08-27 18:35:35 +02:00
|
|
|
end
|
|
|
|
|
2017-11-08 09:28:04 +01:00
|
|
|
return Dict(elts[i]=>preps[i] for i in 1:l)
|
2017-08-27 18:35:35 +02:00
|
|
|
end
|
|
|
|
|
2017-11-08 09:56:28 +01:00
|
|
|
function perm_reps(S::Vector, autS::Group, radius::Int)
|
2017-11-08 09:28:04 +01:00
|
|
|
E, _ = Groups.generate_balls(S, radius=radius)
|
2017-11-08 09:56:28 +01:00
|
|
|
return perm_reps(autS, E)
|
2017-08-27 18:35:35 +02:00
|
|
|
end
|
2017-08-03 11:35:34 +02:00
|
|
|
|
2017-08-27 18:35:35 +02:00
|
|
|
function reconstruct_sol{T<:GroupElem, S<:Nemo.perm}(preps::Dict{T, S},
|
2017-11-13 10:57:57 +01:00
|
|
|
Us::Vector, Ps::Vector, dims::Vector)
|
2017-08-27 18:37:09 +02:00
|
|
|
|
|
|
|
l = length(Us)
|
|
|
|
transfP = [dims[π].*Us[π]*Ps[π]*Us[π]' for π in 1:l]
|
|
|
|
tmp = [zeros(Float64, size(first(transfP))) for _ in 1:l]
|
|
|
|
perms = collect(keys(preps))
|
|
|
|
|
|
|
|
@inbounds Threads.@threads for π in 1:l
|
|
|
|
for p in perms
|
|
|
|
BLAS.axpy!(1.0, view(transfP[π], preps[p].d, preps[p].d), tmp[π])
|
2017-08-03 11:35:34 +02:00
|
|
|
end
|
|
|
|
end
|
2017-08-27 18:37:09 +02:00
|
|
|
|
|
|
|
recP = 1/length(perms) .* sum(tmp)
|
2017-08-03 11:35:34 +02:00
|
|
|
recP[abs.(recP) .< eps(eltype(recP))] = zero(eltype(recP))
|
|
|
|
return recP
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2017-07-27 22:02:07 +02:00
|
|
|
function Cstar_repr{T}(x::GroupRingElem{T}, mreps::Dict)
|
2017-11-08 09:29:12 +01:00
|
|
|
return sum(x[g].*mreps[g] for g in parent(x).basis if x[g] != zero(T))
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2017-11-08 09:50:58 +01:00
|
|
|
function orthSVD{T}(M::AbstractMatrix{T})
|
|
|
|
M = full(M)
|
|
|
|
fact = svdfact(M)
|
|
|
|
M_rank = sum(fact[:S] .> maximum(size(M))*eps(T))
|
|
|
|
return fact[:U][:,1:M_rank]
|
2017-06-22 14:12:35 +02:00
|
|
|
end
|
|
|
|
|
2017-11-08 09:56:28 +01:00
|
|
|
function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S::Vector{T}, autS::Nemo.Group; radius=2)
|
2017-06-22 14:12:35 +02:00
|
|
|
isdir(name) || mkdir(name)
|
|
|
|
|
|
|
|
info(logger, "Generating ball of radius $(2*radius)")
|
|
|
|
|
|
|
|
# TODO: Fix that by multiple dispatch?
|
|
|
|
Id = (isa(G, Nemo.Ring) ? one(G) : G())
|
|
|
|
|
2017-11-08 09:38:56 +01:00
|
|
|
@logtime logger E_2R, sizes = Groups.generate_balls(S, Id, radius=2*radius);
|
2017-06-22 14:12:35 +02:00
|
|
|
info(logger, "Balls of sizes $sizes.")
|
|
|
|
info(logger, "Reverse dict")
|
2017-11-08 09:38:56 +01:00
|
|
|
@logtime logger E_rdict = GroupRings.reverse_dict(E_2R)
|
2017-06-22 14:12:35 +02:00
|
|
|
|
|
|
|
info(logger, "Product matrix")
|
2017-11-08 09:38:56 +01:00
|
|
|
@logtime logger pm = GroupRings.create_pm(E_2R, E_rdict, sizes[radius], twisted=true)
|
|
|
|
RG = GroupRing(G, E_2R, E_rdict, pm)
|
2017-06-22 14:12:35 +02:00
|
|
|
Δ = PropertyT.splaplacian(RG, S)
|
|
|
|
@assert GroupRings.augmentation(Δ) == 0
|
|
|
|
save(joinpath(name, "delta.jld"), "Δ", Δ.coeffs)
|
|
|
|
save(joinpath(name, "pm.jld"), "pm", pm)
|
|
|
|
|
2017-11-08 09:56:28 +01:00
|
|
|
info(logger, "Decomposing E into orbits of $(autS)")
|
|
|
|
@logtime logger orbs = orbit_decomposition(autS, E_2R, E_rdict)
|
2017-11-08 09:38:56 +01:00
|
|
|
@assert sum(length(o) for o in orbs) == length(E_2R)
|
2017-11-08 15:40:38 +01:00
|
|
|
info(logger, "E consists of $(length(orbs)) orbits!")
|
2017-06-22 14:12:35 +02:00
|
|
|
save(joinpath(name, "orbits.jld"), "orbits", orbs)
|
|
|
|
|
|
|
|
info(logger, "Action matrices")
|
2017-11-08 09:56:28 +01:00
|
|
|
@logtime logger reps = perm_reps(autS, E_2R[1:sizes[radius]], E_rdict)
|
2017-11-08 09:33:00 +01:00
|
|
|
save_preps(joinpath(name, "preps.jld"), reps)
|
|
|
|
reps = matrix_reps(reps)
|
2017-06-22 14:12:35 +02:00
|
|
|
|
|
|
|
info(logger, "Projections")
|
2017-11-08 09:56:28 +01:00
|
|
|
@logtime logger autS_mps = rankOne_projections(autS);
|
2017-06-22 14:12:35 +02:00
|
|
|
|
2017-11-08 09:56:28 +01:00
|
|
|
@logtime logger π_E_projections = [Cstar_repr(p, reps) for p in autS_mps]
|
2017-06-22 14:12:35 +02:00
|
|
|
|
|
|
|
info(logger, "Uπs...")
|
2017-11-05 16:38:51 +01:00
|
|
|
@logtime logger Uπs = orthSVD.(π_E_projections)
|
2017-06-22 14:12:35 +02:00
|
|
|
|
|
|
|
multiplicities = size.(Uπs,2)
|
|
|
|
info(logger, "multiplicities = $multiplicities")
|
2017-11-08 09:56:28 +01:00
|
|
|
dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps];
|
2017-06-22 14:12:35 +02:00
|
|
|
info(logger, "dimensions = $dimensions")
|
|
|
|
@assert dot(multiplicities, dimensions) == sizes[radius]
|
|
|
|
|
2017-08-27 21:25:12 +02:00
|
|
|
save(joinpath(name, "U_pis.jld"),
|
|
|
|
"Uπs", Uπs,
|
|
|
|
"dims", dimensions)
|
2017-06-22 14:12:35 +02:00
|
|
|
return 0
|
|
|
|
end
|