add orbit-related code

src/OrbitDecomposition.jl

@@ -0,0 +1,206 @@
+push!(LOAD_PATH, "./")
+
+using Nemo
+using Groups
+
+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)
+
+###############################################################################
+#
+#  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]
+            orbit = Vector{Int}()
+            a = E[i]
+            for g in elts
+                idx = rdict[g(a)]
+                tovisit[idx] = false
+                push!(orbit,idx)
+            end
+            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
+
+function orbit_constraint(constraints::Vector{Vector{Vector{Int64}}}, n)
+    result = spzeros(n,n)
+    for cnstr in constraints
+        for p in cnstr
+            result[p[2], p[1]] += 1.0
+        end
+    end
+    return 1/length(constraints)*result
+end
+
+###############################################################################
+#
+#  Matrix- and C*-representations
+#
+###############################################################################
+
+function matrix_repr(g::GroupElem, 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 matrix_reps{T<:GroupElem}(G::Nemo.Group, S::Vector{T}, AutS::Nemo.Group, radius::Int)
+   Id = (isa(G, Nemo.Ring) ? one(G) : G())
+   E2, _ = Groups.generate_balls(S, Id, radius=radius)
+   Edict = GroupRings.reverse_dict(E2)
+
+   mreps = Dict(g=>matrix_repr(g, E2, Edict) for g in elements(AutS))
+   return mreps
+end
+
+function reconstruct_sol(mreps::Dict, Us::Vector, Ps::Vector, dims::Vector)
+   recP = zeros(size(Us[1],1), size(Us[1],1))
+
+   for g in keys(mreps)
+      for π in 1:endof(Us)
+         recP .+= dims[π] .* mreps[g]*transform(Us[π]', Ps[π])*mreps[inv(g)]
+      end
+   end
+   recP .*= 1/length(collect(keys(mreps)))
+
+   return recP
+end
+
+function Cstar_repr(x::GroupRingElem, mreps)
+    k = collect(keys(mreps))[1]
+    res = zeros(size(mreps[k])...)
+
+    for g in parent(x).basis
+        res .+= x[g]*mreps[g]
+    end
+    return res
+end
+
+function orthSVD(M::AbstractMatrix)
+    M = full(M)
+    fact = svdfact(M)
+    singv = fact[:S]
+    M_rank = sum(singv .> maximum(size(M))*eps(eltype(singv)))
+    return fact[:U][:,1:M_rank]
+end
+
+function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S::Vector{T}, AutS; radius=2)
+   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())
+
+   @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_mreps = Dict(g=>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_mreps) for p in AutS_mps]
+
+   info(logger, "Uπs...")
+   @time Uπs = orthSVD.(π_E_projections)
+
+   multiplicities = size.(Uπs,2)
+   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

src/Projections.jl

@@ -0,0 +1,135 @@
+using DirectProducts
+using WreathProducts
+
+###############################################################################
+#
+#  Characters of PermutationGroup
+#
+###############################################################################
+
+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
+
+###############################################################################
+#
+#  Character of DirectProducts
+#
+###############################################################################
+
+function epsilon(i, g::DirectProducts.DirectProductGroupElem)
+    return reduce(*, 1, ((-1)^isone(g.elts[j]) for j in 1:i))
+end
+
+###############################################################################
+#
+#  Projections
+#
+###############################################################################
+
+function central_projection(RG::GroupRing, char::Function, T::Type=Rational{Int})
+    result = RG(T)
+    for g in RG.basis
+        result[g] = char(inv(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

src/oPropertyT.jl

@@ -0,0 +1,234 @@
+using JLD
+using JuMP
+using SCS
+
+using GroupRings
+using PropertyT
+
+using ValidatedNumerics
+using ArgParse
+
+import Nemo: Group, GroupElem
+
+immutable Settings
+   name::String
+   N::Int
+   G::Group
+   S::Vector
+   AutS::Group
+   radius::Int
+   solver::SCSSolver
+   upper_bound::Float64
+   tol::Float64
+end
+
+immutable OrbitData
+   name::String
+   Us::Vector
+   Ps::Vector{Array{JuMP.Variable,2}}
+   cnstr::Vector
+   laplacian::Vector
+   laplacianSq::Vector
+   dims::Vector{Int}
+end
+
+include("OrbitDecomposition.jl")
+
+function sparsify{T}(U::Array{T}, eps=eps(T))
+    n = rank(U)
+    W = deepcopy(U)
+    W[abs.(W) .< eps] = zero(T)
+    if rank(W) != n
+        warn("Sparsification would decrease the rank!")
+        W = U
+    end
+    W = sparse(W)
+    dropzeros!(W)
+    return W
+end
+
+function sparsify!{T}(U::SparseMatrixCSC{T}, eps=eps(T))
+    U[abs.(U) .< eps] = zero(T)
+    dropzeros!(U)
+    return U
+end
+
+sparsify{T}(U::SparseMatrixCSC{T}, eps=eps(T)) = sparsify!(deepcopy(U), eps)
+
+function init_model(Uπs)
+    m = JuMP.Model();
+    l = size(Uπs,1)
+    P = Vector{Array{JuMP.Variable,2}}(l)
+
+    for k in 1:l
+        s = size(Uπs[k],2)
+        P[k] = JuMP.@variable(m, [i=1:s, j=1:s])
+        JuMP.@SDconstraint(m, P[k] >= 0.0)
+    end
+
+    JuMP.@variable(m, λ >= 0.0)
+    JuMP.@objective(m, Max, λ)
+    return m, P
+end
+
+function init_OrbitData(name::String)
+   splap = load(joinpath(name, "delta.jld"), "Δ");
+   pm = load(joinpath(name, "pm.jld"), "pm");
+   cnstr = PropertyT.constraints_from_pm(pm);
+   splap² = GroupRings.mul(splap, splap, pm);
+
+   Uπs = load(joinpath(name, "U_pis.jld"), "Uπs");
+   # Uπs = sparsify.(Uπs);
+   #dimensions of the corresponding πs:
+   dims = load(joinpath(name, "U_pis.jld"), "dims")
+
+   m, P = init_model(Uπs);
+
+   orbits = load(joinpath(name, "orbits.jld"), "orbits");
+   n = size(Uπs[1],1)
+   orb_spcnstrm = [orbit_constraint(cnstr[collect(orb)], n) for orb in orbits]
+   orb_splap = orbit_spvector(splap, orbits)
+   orb_splap² = orbit_spvector(splap², orbits)
+
+   orbData = OrbitData(name, Uπs, P, orb_spcnstrm, orb_splap, orb_splap², dims);
+
+   # orbData = OrbitData(name, Uπs, P, orb_spcnstrm, splap, splap², dims);
+
+   return m, orbData
+end
+
+function transform(U::AbstractArray, V::AbstractArray; sparse=false)
+    w = U'*V*U
+    sparse && sparsify!(w)
+    return w
+end
+
+A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t])
+
+function constrLHS(m::JuMP.Model, data::OrbitData, t)
+    l = endof(data.Us)
+    lhs = @expression(m, sum(vecdot(A(data, π, t), data.Ps[π]) for π in 1:l))
+    return lhs
+end
+
+function addconstraints!(m::JuMP.Model, data::OrbitData, l::Int=length(data.laplacian); var::Symbol = :λ)
+    λ = m[var]
+   #  orbits = load(joinpath(data.name, "orbits.jld"), "orbits");
+   #  locate(t, orb=orbits) = findfirst(x->t in x, orb)
+    for t in 1:l
+      #   lhs = constrLHS(m, data, locate(t))
+        lhs = constrLHS(m, data, t)
+
+        d, d² = data.laplacian[t], data.laplacianSq[t]
+        if lhs == zero(lhs)
+            if d == 0 && d² == 0
+                info("Detected empty constraint")
+                continue
+            else
+                warn("Adding unsatisfiable constraint!")
+            end
+        end
+        JuMP.@constraint(m, lhs == d² - λ*d)
+    end
+end
+
+function create_SDP_problem(name::String; upper_bound=Inf)
+   info(PropertyT.logger, "Loading orbit data....")
+   t = @timed SDP_problem, orb_data = init_OrbitData(name);
+   info(PropertyT.logger, PropertyT.timed_msg(t))
+
+   if upper_bound < Inf
+      λ = JuMP.getvariable(SDP_problem, :λ)
+      JuMP.@constraint(SDP_problem, λ <= upper_bound)
+   end
+
+   info(PropertyT.logger, "Adding constraints... ")
+   t = @timed addconstraints!(SDP_problem, orb_data)
+   info(PropertyT.logger, PropertyT.timed_msg(t))
+
+   return SDP_problem, orb_data
+end
+
+function λandP(m::JuMP.Model, data::OrbitData)
+   varλ = m[:λ]
+   varP = data.Ps
+   λ, Ps = PropertyT.λandP(data.name, m, varλ, varP)
+   return λ, Ps
+end
+
+function λandP(m::JuMP.Model, data::OrbitData, sett::Settings)
+   info(PropertyT.logger, "Solving SDP problem...")
+   λ, Ps = λandP(m, data)
+
+   info(PropertyT.logger, "Reconstructing P...")
+   mreps = matrix_reps(sett.G, sett.S, sett.AutS, sett.radius)
+
+   recP = reconstruct_sol(mreps, data.Us, Ps, data.dims)
+
+   fname = PropertyT.λSDPfilenames(data.name)[2]
+   save(fname, "origP", Ps, "P", recP)
+   return λ, recP
+end
+
+function init_orbit_data(logger, sett::Settings; radius=2)
+
+   ex(fname) = isfile(joinpath(sett.name, fname))
+
+   files_exists = ex.(["delta.jld", "pm.jld", "U_pis.jld", "orbits.jld"])
+
+   if !all(files_exists)
+      compute_orbit_data(logger, sett.name, sett.G, sett.S, sett.AutS, radius=radius)
+   end
+
+   return 0
+end
+
+function orbit_check_propertyT(logger, sett::Settings)
+
+   init_orbit_data(logger, sett, radius=sett.radius)
+
+   Δ = PropertyT.ΔandSDPconstraints(sett.name, sett.G)[1]
+
+   fnames = PropertyT.λSDPfilenames(sett.name)
+
+   if all(isfile.(fnames))
+      λ, P = PropertyT.λandP(sett.name)
+   else
+      info(logger, "Creating SDP problem...")
+      SDP_problem, orb_data = create_SDP_problem(sett.name, upper_bound=sett.upper_bound)
+      JuMP.setsolver(SDP_problem, sett.solver)
+
+      λ, P = λandP(SDP_problem, orb_data, sett)
+   end
+
+   info(logger, "λ = $λ")
+   info(logger, "sum(P) = $(sum(P))")
+   info(logger, "maximum(P) = $(maximum(P))")
+   info(logger, "minimum(P) = $(minimum(P))")
+
+   if λ > 0
+
+      isapprox(eigvals(P), abs(eigvals(P)), atol=sett.tol) ||
+          warn("The solution matrix doesn't seem to be positive definite!")
+     #  @assert P == Symmetric(P)
+      Q = real(sqrtm(Symmetric(P)))
+
+      sgap = PropertyT.check_distance_to_positive_cone(Δ, λ, Q, 2*sett.radius, tol=sett.tol, rational=false)
+      if isa(sgap, Interval)
+           sgap = sgap.lo
+      end
+      if sgap > 0
+           info(logger, "λ ≥ $(Float64(trunc(sgap,12)))")
+            Kazhdan_κ = PropertyT.Kazhdan_from_sgap(sgap, length(sett.S))
+            Kazhdan_κ = Float64(trunc(Kazhdan_κ, 12))
+            info(logger, "κ($(sett.name), S) ≥ $Kazhdan_κ: Group HAS property (T)!")
+            return true
+      else
+           sgap = Float64(trunc(sgap, 12))
+           info(logger, "λ($(sett.name), S) ≥ $sgap: Group may NOT HAVE property (T)!")
+           return false
+      end
+   end
+   info(logger, "κ($(sett.name), S) ≥ $λ < 0: Tells us nothing about property (T)")
+   return false
+end