diff --git a/OrbitDecomposition.jl b/OrbitDecomposition.jl
new file mode 100644
index 0000000..603f85b
--- /dev/null
+++ b/OrbitDecomposition.jl
@@ -0,0 +1,306 @@
+push!(LOAD_PATH, "./")
+
+using Nemo
+using Groups
+using WreathProducts
+
+using GroupRings
+
+import Nemo.elements
+
+using JLD
+using ProgressMeter
+
+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
+
+import Base.convert
+convert(::Type{Int}, x::Nemo.fmpz) = x.d
+
+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 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 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 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 main(N::Int)
+   name = "test"
+   # SOutF$(N)_E4
+
+   isdir(name) || mkdir(name)
+
+   SOutFN = AutGroup(FreeGroup(N), special=true, outer=true)
+   S = generators(SOutFN);
+   S = [S; [inv(s) for s in S]]
+   @show S
+
+   info("Generating ball of radius 4")
+   @time E4, sizes = generate_balls(S, SOutFN(), radius=4);
+   info("Reverse dict")
+   @time E_dict = GroupRings.reverse_dict(E4)
+
+
+   BN = WreathProduct(Nemo.FiniteField(2,1, "a")[1], PermutationGroup(N))
+
+   if !isfile(joinpath(name, "orbits.jld"))
+      info("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)
+   end
+
+   info("Action matrices")
+
+   pm = GroupRings.create_pm(E4, E_dict, sizes[2], twisted=true)
+
+   RSOutFN = GroupRing(SOutFN, E4, E_dict, pm)
+
+   PropertyT.splaplacian(RSOutFN, S)
+
+   E2 = E4[1:sizes[2]]
+   @time BNactionE_mreps = action_mreps(BN, E2, E_dict)
+
+   save(joinpath(name, "matrix_reps.jld"), "matrix_reps", BNactionE_mreps)
+
+   info("Projections")
+   @time BN_mps = rankOne_projections(BN);
+
+   @time π_E_projections = @showprogress [Cstar_repr(p, BNactionE_mreps) for p in BN_mps]
+
+   info("Uπs...")
+   @time Uπs = Uπ_matrices(π_E_projections);
+
+   @show multiplicities = [size(U,2) for U in Uπs];
+   @show dimensions = [Int(p[BN()]*Int(order(BN))) for p in BN_mps];
+   @assert dot(multiplicities, dimensions) == sizes[2]
+
+   save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "dims", dimensions)
+end
+
+@time main(4)