add sqadjop.jl and unit tests

src/PropertyT.jl

@@ -22,6 +22,9 @@ include("RGprojections.jl")
 include("orbitdata.jl")
 include("sos_sdps.jl")
 include("checksolution.jl")
+include("sqadjop.jl")
+
 include("1712.07167.jl")
 
+
 end # module Property(T)

src/sqadjop.jl

@@ -0,0 +1,125 @@
+indexing(n) = [(i,j) for i in 1:n for j in 1:n if i≠j]
+
+function generating_set(G::AutGroup{N}, n=N) where N
+
+    rmuls = [Groups.rmul_autsymbol(i,j) for (i,j) in indexing(n)]
+    lmuls = [Groups.lmul_autsymbol(i,j) for (i,j) in indexing(n)]
+    gen_set = G.([rmuls; lmuls])
+
+    return [gen_set; inv.(gen_set)]
+end
+
+function E(M::MatSpace, i::Integer, j::Integer, val=1)
+    @assert i ≠ j
+    @assert 1 ≤ i ≤ nrows(M)
+    @assert 1 ≤ j ≤ ncols(M)
+    m = one(M)
+    m[i,j] = val
+    return m
+end
+
+function generating_set(M::MatSpace, n=nrows(M))
+    @assert nrows(M) == ncols(M)
+
+    elts = [E(M, i,j) for (i,j) in indexing(n)]
+    return elem_type(M)[elts; inv.(elts)]
+end
+
+isopposite(σ::perm, τ::perm, i=1, j=2) =
+    σ[i] ≠ τ[i] && σ[i] ≠ τ[j] &&
+    σ[j] ≠ τ[i] && σ[j] ≠ τ[j]
+
+isadjacent(σ::perm, τ::perm, i=1, j=2) =
+    (σ[i] == τ[i] && σ[j] ≠ τ[j]) || # first equal, second differ
+    (σ[j] == τ[j] && σ[i] ≠ τ[i]) || # sedond equal, first differ
+    (σ[i] == τ[j] && σ[j] ≠ τ[i]) || # first σ equal to second τ
+    (σ[j] == τ[i] && σ[i] ≠ τ[j])    # second σ equal to first τ
+
+Base.div(X::GroupRingElem, x::Number) = parent(X)(X.coeffs.÷x)
+
+function Sq(RG::GroupRing, N::Integer)
+    S₂ = generating_set(RG.group, 2)
+    ℤ = Int64
+    Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
+
+    Alt_N = [g for g in PermutationGroup(N) if parity(g) == 0]
+
+    sq = RG()
+    for σ in Alt_N
+        GroupRings.addeq!(sq, *(σ(Δ₂), σ(Δ₂), false))
+    end
+    return sq÷factorial(N-2)
+end
+
+function Adj(RG::GroupRing, N::Integer)
+    S₂ = generating_set(RG.group, 2)
+    ℤ = Int64
+    Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
+
+    Alt_N = [g for g in PermutationGroup(N) if parity(g) == 0]
+    Δ₂s = Dict(σ=>σ(Δ₂) for σ in Alt_N)
+    adj = RG()
+
+    for σ in Alt_N
+        for τ in Alt_N
+            if isadjacent(σ, τ)
+                GroupRings.addeq!(adj, *(Δ₂s[σ], Δ₂s[τ], false))
+            end
+        end
+    end
+    return adj÷factorial(N-2)^2
+end
+
+function Op(RG::GroupRing, N::Integer)
+    if N < 4
+        return RG()
+    end
+    S₂ = generating_set(RG.group, 2)
+    ℤ = Int64
+    Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
+
+    Alt_N = [g for g in PermutationGroup(N) if parity(g) == 0]
+    Δ₂s = Dict(σ=>σ(Δ₂) for σ in Alt_N)
+    op = RG()
+
+    for σ in Alt_N
+        for τ in Alt_N
+            if isopposite(σ, τ)
+                GroupRings.addeq!(op, *(Δ₂s[σ], Δ₂s[τ], false))
+            end
+        end
+    end
+    return op÷factorial(N-2)^2
+end
+
+for Elt in [:Sq, :Adj, :Op]
+    @eval begin
+        $Elt(RG::GroupRing{AutGroup{N}}) where N = $Elt(RG, N)
+        $Elt(RG::GroupRing{<:MatSpace}) = $Elt(RG, nrows(RG.group))
+    end
+end
+
+function SqAdjOp(RG::GroupRing, N::Integer)
+    S₂ = generating_set(RG.group, 2)
+    ℤ = Int64
+    Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ);
+
+    Alt_N = [σ for σ in PermutationGroup(N) if parity(σ) == 0]
+    sq, adj, op = RG(), RG(), RG()
+
+    Δ₂s = Dict(σ=>σ(Δ₂) for σ in Alt_N)
+
+    for σ in Alt_N
+        GroupRings.addeq!(sq, *(Δ₂s[σ], Δ₂s[σ], false))
+        for τ in Alt_N
+            if isopposite(σ, τ)
+                GroupRings.addeq!(op, *(Δ₂s[σ], Δ₂s[τ], false))
+            elseif isadjacent(σ, τ)
+                GroupRings.addeq!(adj, *(Δ₂s[σ], Δ₂s[τ], false))
+            end
+        end
+    end
+
+    k = factorial(N-2)
+    return sq÷k, adj÷k^2, op÷k^2
+end

test/1812.03456.jl

@@ -0,0 +1,93 @@
+@testset "Sq, Adj, Op" begin
+    function isconstant_on_orbit(v, orb)
+        isempty(orb) && return true
+        k = v[first(orb)]
+        return all(v[o] == k for o in orb)
+    end
+
+    @testset "unit tests" begin
+        for N in [3,4]
+            M = MatrixSpace(Nemo.ZZ, N,N)
+            A = SAut(FreeGroup(N))
+            @test length(PropertyT.generating_set(M)) == 2N*(N-1)
+            S = PropertyT.generating_set(M)
+            @test all(inv(s) ∈ S for s in S)
+            @test length(PropertyT.generating_set(A)) == 4N*(N-1)
+            S = PropertyT.generating_set(A)
+            @test all(inv(s) ∈ S for s in S)
+        end
+
+        N = 4
+        M = MatrixSpace(Nemo.ZZ, N,N)
+        S = PropertyT.generating_set(M)
+
+        @test PropertyT.E(M, 1, 2) isa MatElem
+        e12 = PropertyT.E(M, 1, 2)
+        @test e12[1,2] == 1
+        @test inv(e12)[1,2] == -1
+        @test e12 ∈ S
+
+        @test PropertyT.isopposite(perm"(1,2,3)(4)", perm"(1,4,2)")
+        @test PropertyT.isadjacent(perm"(1,2,3)", perm"(1,2)(3)")
+
+        @test !PropertyT.isopposite(perm"(1,2,3)", perm"(1,2)(3)")
+        @test !PropertyT.isadjacent(perm"(1,4)", perm"(2,3)(4)")
+
+        @test isconstant_on_orbit([1,1,1,2,2], [2,3])
+        @test !isconstant_on_orbit([1,1,1,2,2], [2,3,4])
+    end
+
+    @testset "Sq, Adj, Op" begin
+
+        N = 4
+        M = MatrixSpace(Nemo.ZZ, N,N)
+        S = PropertyT.generating_set(M)
+        Δ = PropertyT.Laplacian(S, 2)
+        RG = parent(Δ)
+
+        autS = WreathProduct(PermGroup(2), PermGroup(N))
+        orbits = PropertyT.orbit_decomposition(autS, RG.basis)
+
+        @test PropertyT.Sq(RG) isa GroupRingElem
+        sq = PropertyT.Sq(RG)
+        @test all(isconstant_on_orbit(sq, orb) for orb in orbits)
+
+        @test PropertyT.Adj(RG) isa GroupRingElem
+        adj = PropertyT.Adj(RG)
+        @test all(isconstant_on_orbit(adj, orb) for orb in orbits)
+
+        @test PropertyT.Op(RG) isa GroupRingElem
+        op = PropertyT.Op(RG)
+        @test all(isconstant_on_orbit(op, orb) for orb in orbits)
+
+        sq, adj, op = PropertyT.SqAdjOp(RG, N)
+        @test sq == PropertyT.Sq(RG)
+        @test adj == PropertyT.Adj(RG)
+        @test op == PropertyT.Op(RG)
+
+        e = one(M)
+        g = PropertyT.E(M, 1,2)
+        h = PropertyT.E(M, 1,3)
+        k = PropertyT.E(M, 3,4)
+
+        edges = N*(N-1)÷2
+        @test sq[e] == 20*edges
+        @test sq[g] == sq[h] == -8
+        @test sq[g^2] == sq[h^2] == 1
+        @test sq[g*h] == sq[h*g] == 0
+
+        #     @test adj[e] == ...
+        @test adj[g] == adj[h] # == ...
+        @test adj[g^2] == adj[h^2] == 0
+        @test adj[g*h] == adj[h*g] # == ...
+
+
+        #     @test op[e] == ...
+        @test op[g] == op[h] # == ...
+        @test op[g^2] == op[h^2] == 0
+        @test op[g*h] == op[h*g] == 0
+        @test op[g*k] == op[k*g] # == ...
+        @test op[h*k] == op[k*h] == 0
+    end
+end
+

test/runtests.jl

@@ -23,3 +23,4 @@ solver(iters; accel=1) =
 include("1703.09680.jl")
 include("1712.07167.jl")
 include("SOS_correctness.jl")
+include("1812.03456.jl")