Merge pull request #4 from kalmarek/enh/AA_MatrixAlgebra

Use Matrix Algebras

.travis.yml

@@ -6,6 +6,7 @@ os:
 julia:
   - 1.0
   - 1.1
+  - 1.2
   - nightly
 notifications:
   email: true

Manifest.toml

@@ -2,9 +2,9 @@
 
 [[AbstractAlgebra]]
 deps = ["InteractiveUtils", "LinearAlgebra", "Markdown", "Random", "SparseArrays", "Test"]
-git-tree-sha1 = "c0d57a3f0618bfbb214005860b9b5e5bceafa61c"
+git-tree-sha1 = "0d4f7283435bd7e12a703a3fd58aa11224a96019"
 uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
-version = "0.5.0"
+version = "0.5.2"
 
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
@@ -16,10 +16,10 @@ uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
 version = "0.8.10"
 
 [[BinaryProvider]]
-deps = ["Libdl", "SHA"]
+deps = ["Libdl", "Logging", "SHA"]
-git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648"
+git-tree-sha1 = "8153fd64131cd00a79544bb23788877741f627bb"
 uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.5.4"
+version = "0.5.5"
 
 [[Blosc]]
 deps = ["BinaryProvider", "CMakeWrapper", "Compat", "Libdl"]
@@ -124,19 +124,19 @@ version = "0.10.3"
 
 [[GroupRings]]
 deps = ["AbstractAlgebra", "LinearAlgebra", "Markdown", "SparseArrays"]
-git-tree-sha1 = "6d43558d0e7946d4e0e3d94a108344b2514e95b7"
+git-tree-sha1 = "9f41bad54217ce2605c13f10b79d1221c259ed32"
 repo-rev = "master"
 repo-url = "https://github.com/kalmarek/GroupRings.jl"
 uuid = "0befed6a-bd73-11e8-1e41-a1190947c2f5"
-version = "0.2.1"
+version = "0.3.0"
 
 [[Groups]]
 deps = ["AbstractAlgebra", "LinearAlgebra", "Markdown"]
-git-tree-sha1 = "64dcaa46affb4568501d35e6fae36a71a7ae5cbb"
+git-tree-sha1 = "bf267f82f4c313a6deb3db64efc7893b139d4340"
 repo-rev = "master"
 repo-url = "https://github.com/kalmarek/Groups.jl"
 uuid = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
-version = "0.2.1"
+version = "0.2.2"
 
 [[HDF5]]
 deps = ["BinDeps", "Blosc", "Homebrew", "Libdl", "Mmap", "WinRPM"]
@@ -192,9 +192,9 @@ version = "0.4.1"
 
 [[LibCURL]]
 deps = ["BinaryProvider", "Libdl"]
-git-tree-sha1 = "5ee138c679fa202ebe211b2683d1eee2a87b3dbe"
+git-tree-sha1 = "fd5fc15f2a04608fe1435a769dbbfc7959ff1daa"
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.5.1"
+version = "0.5.2"
 
 [[LibExpat]]
 deps = ["Compat"]
@@ -246,12 +246,6 @@ git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 version = "0.3.2"
 
-[[Nemo]]
-deps = ["AbstractAlgebra", "BinaryProvider", "InteractiveUtils", "Libdl", "LinearAlgebra", "Markdown", "Test"]
-git-tree-sha1 = "b359c25b708c3edcc2f17345d59da7169c207385"
-uuid = "2edaba10-b0f1-5616-af89-8c11ac63239a"
-version = "0.14.0"
-
 [[OrderedCollections]]
 deps = ["Random", "Serialization", "Test"]
 git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"

Project.toml

@@ -1,7 +1,7 @@
 name = "PropertyT"
 uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d"
 authors = ["Marek Kaluba <kalmar@amu.edu.pl>"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
@@ -14,12 +14,11 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
-Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-GroupRings = "^0.2.0"
+GroupRings = "^0.3.0"
 Groups = "^0.2.1"
 JuMP = "^0.19.0"
 

src/1712.07167.jl

@@ -10,9 +10,9 @@ abstract type Settings end
 
 struct Naive{El} <: Settings
     name::String
-    G::Group
+    G::Union{Group, NCRing}
     S::Vector{El}
-    radius::Int
+    halfradius::Int
     upper_bound::Float64
 
     solver::JuMP.OptimizerFactory
@@ -21,10 +21,10 @@ end
 
 struct Symmetrized{El} <: Settings
     name::String
-    G::Group
+    G::Union{Group, NCRing}
     S::Vector{El}
     autS::Group
-    radius::Int
+    halfradius::Int
     upper_bound::Float64
 
     solver::JuMP.OptimizerFactory
@@ -32,15 +32,15 @@ struct Symmetrized{El} <: Settings
 end
 
 function Settings(name::String,
-    G::Group, S::Vector{<:GroupElem}, solver::JuMP.OptimizerFactory;
+    G::Union{Group, NCRing}, S::AbstractVector{El}, solver::JuMP.OptimizerFactory;
-    radius::Integer=2, upper_bound::Float64=1.0, warmstart=true)
+    halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) where El <: Union{GroupElem, NCRingElem}
-    return Naive(name, G, S, radius, upper_bound, solver, warmstart)
+    return Naive(name, G, S, halfradius, upper_bound, solver, warmstart)
 end
 
 function Settings(name::String,
-    G::Group, S::Vector{<:GroupElem}, autS::Group, solver::JuMP.OptimizerFactory;
+    G::Union{Group, NCRing}, S::AbstractVector{El}, autS::Group, solver::JuMP.OptimizerFactory;
-    radius::Integer=2, upper_bound::Float64=1.0, warmstart=true)
+    halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) where El <: Union{GroupElem, NCRingElem}
-    return Symmetrized(name, G, S, autS, radius, upper_bound, solver, warmstart)
+    return Symmetrized(name, G, S, autS, halfradius, upper_bound, solver, warmstart)
 end
 
 prefix(s::Naive) = s.name
@@ -79,12 +79,13 @@ end
 ###############################################################################
 
 function warmstart(sett::Settings)
+    warmstart_fname = filename(sett, :warmstart)
     try
         ws = load(filename(sett, :warmstart), "warmstart")
-        @info "Loaded $(filename(sett, :warmstart))."
+        @info "Loaded $warmstart_fname."
         return ws
-    catch
+    catch ex
-        @info "Couldn't load $(filename(sett, :warmstart))."
+        @warn "$(ex.msg). Not providing a warmstart to the solver."
         return nothing
     end
 end
@@ -128,10 +129,12 @@ function approximate_by_SOS(sett::Symmetrized,
 
     orbit_data = try
         orbit_data = load(filename(sett, :OrbitData), "OrbitData")
-        @info "Loaded Orbit Data"
+        @info "Loaded orbit data."
         orbit_data
-    catch
+    catch ex
+        @warn ex.msg
         isdefined(parent(orderunit), :basis) || throw("You need to define basis of Group Ring to compute orbit decomposition!")
+        @info "Computing orbit and Wedderburn decomposition..."
         orbit_data = OrbitData(parent(orderunit), sett.autS)
         save(filename(sett, :OrbitData), "OrbitData", orbit_data)
         orbit_data
@@ -245,7 +248,7 @@ function print_summary(sett::Settings)
     separator = "="^76
     info_strs = [separator,
     "Running tests for $(sett.name):",
-    "Upper bound for λ: $(sett.upper_bound), on radius $(sett.radius).",
+    "Upper bound for λ: $(sett.upper_bound), on halfradius $(sett.halfradius).",
     "Warmstart: $(sett.warmstart)",
     "Results will be stored in ./$(PropertyT.prepath(sett))",
     "Solver: $(typeof(sett.solver()))",
@@ -275,26 +278,39 @@ function spectral_gap(sett::Settings)
     isdir(fp) || mkpath(fp)
 
     Δ = try
-        @info "Loading precomputed Δ..."
+        Δ = loadGRElem(filename(sett,:Δ), sett.G)
-        loadGRElem(filename(sett,:Δ), sett.G)
+        @info "Loaded precomputed Δ."
-    catch
+        Δ
-        # compute
+    catch ex
-        Δ = Laplacian(sett.S, sett.radius)
+        @warn ex.msg
+        @info "Computing Δ..."
+        Δ = Laplacian(sett.S, sett.halfradius)
         saveGRElem(filename(sett, :Δ), Δ)
         Δ
     end
 
-    λ, P = try
+    function compute(sett, Δ)
-        sett.warmstart && error()
+        @info "Computing λ and P..."
-        load(filename(sett, :solution), "λ", "P")
-    catch
         λ, P = approximate_by_SOS(sett, Δ^2, Δ;
             solverlog=filename(sett, :solverlog))
 
         save(filename(sett, :solution), "λ", λ, "P", P)
 
         λ < 0 && @warn "Solver did not produce a valid solution!"
-        λ, P
+        return λ, P
+    end
+
+    if sett.warmstart
+        λ, P = compute(sett, Δ)
+    else
+        λ, P =try
+            λ, P = load(filename(sett, :solution), "λ", "P")
+            @info "Loaded existing λ and P."
+            λ, P
+        catch ex
+            @warn ex.msg
+            compute(sett, Δ)
+        end
     end
 
     info_strs = ["Numerical metrics of matrix solution:",
@@ -307,7 +323,7 @@ function spectral_gap(sett::Settings)
         @warn "The solution matrix doesn't seem to be positive definite!"
 
     @time Q = real(sqrt(Symmetric( (P.+ P')./2 )))
-    certified_sgap = spectral_gap(Δ, λ, Q, R=sett.radius)
+    certified_sgap = spectral_gap(Δ, λ, Q, R=sett.halfradius)
 
     return certified_sgap
 end

src/1812.03456.jl

@@ -0,0 +1,26 @@
+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 EltaryMat(M::MatAlgebra, 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::MatAlgebra, n=nrows(M))
+    elts = [EltaryMat(M, i,j) for (i,j) in indexing(n)]
+    return elem_type(M)[elts; inv.(elts)]
+end
+
+include("sqadjop.jl")

src/PropertyT.jl

@@ -13,15 +13,19 @@ using GroupRings
 using JLD
 using JuMP
 
-import AbstractAlgebra: Group, Ring, perm
+import AbstractAlgebra: Group, NCRing, perm
 
 import MathProgBase.SolverInterface.AbstractMathProgSolver
 
+AbstractAlgebra.one(G::Group) = G()
+
 include("laplacians.jl")
 include("RGprojections.jl")
 include("orbitdata.jl")
 include("sos_sdps.jl")
 include("checksolution.jl")
+
 include("1712.07167.jl")
+include("1812.03456.jl")
 
 end # module Property(T)

src/RGprojections.jl

@@ -70,7 +70,7 @@ end
 
 function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int})
     result = RG(zeros(T, length(RG.basis)))
-    dim = chi(RG.group())
+    dim = chi(one(RG.group))
     ord = Int(order(RG.group))
 
     for g in RG.basis
@@ -187,8 +187,8 @@ function rankOne_projections(RBn::GroupRing{G}, T::Type=Rational{Int}) where {G<
 
     r = collect(1:N)
     for i in 1:N-1
-        first_emb = g->Bn(Generic.emb!(Bn.P(), g, view(r, 1:i)))
+        first_emb = g->Bn(Generic.emb!(one(Bn.P), g, view(r, 1:i)))
-        last_emb = g->Bn(Generic.emb!(Bn.P(), g, view(r, (i+1):N)))
+        last_emb = g->Bn(Generic.emb!(one(Bn.P), g, view(r, (i+1):N)))
 
         Sk_first = (RBn(first_emb, p) for p in Sn_rankOnePr[i])
         Sk_last = (RBn(last_emb, p) for p in Sn_rankOnePr[N-i])
@@ -238,7 +238,7 @@ end
 > The identity element `Id` and binary operation function `op` can be supplied
 > to e.g. take advantage of additive group structure.
 """
-function generateGroup(gens::Vector{T}, r=2, Id::T=parent(first(gens))(), op=*) where {T<:GroupElem}
+function generateGroup(gens::Vector{T}, r=2, Id::T=one(parent(first(gens))), op=*) where {T<:GroupElem}
     n = 0
     R = 1
     elts = gens

src/laplacians.jl

@@ -4,16 +4,7 @@
 #
 ###############################################################################
 
-function spLaplacian(RG::GroupRing, S, T::Type=Float64)
+function spLaplacian(RG::GroupRing, S::AbstractVector, T::Type=Float64)
-    result = RG(T)
-    result[RG.group()] = T(length(S))
-    for s in S
-        result[s] -= one(T)
-    end
-    return result
-end
-
-function spLaplacian(RG::GroupRing, S::Vector{REl}, T::Type=Float64) where {REl<:AbstractAlgebra.ModuleElem}
     result = RG(T)
     result[one(RG.group)] = T(length(S))
     for s in S
@@ -22,26 +13,17 @@ function spLaplacian(RG::GroupRing, S::Vector{REl}, T::Type=Float64) where {REl<
     return result
 end
 
-function Laplacian(S::Vector{E}, radius) where E<:AbstractAlgebra.ModuleElem
+function Laplacian(S::AbstractVector{REl}, halfradius) where REl<:Union{NCRingElem, GroupElem}
-    R = parent(first(S))
-    return Laplacian(S, one(R), radius)
-end
-
-function Laplacian(S::Vector{E}, radius) where E<:AbstractAlgebra.GroupElem
     G = parent(first(S))
-    return Laplacian(S, G(), radius)
+    @info "Generating metric ball of radius" radius=2halfradius
-end
+    @time E_R, sizes = Groups.generate_balls(S, one(G), radius=2halfradius)
-
-function Laplacian(S, Id, radius)
-    @info "Generating metric ball of radius" radius=2radius
-    @time E_R, sizes = Groups.generate_balls(S, Id, radius=2radius)
     @info "Generated balls:" sizes
 
     @info "Creating product matrix..."
     rdict = GroupRings.reverse_dict(E_R)
-    @time pm = GroupRings.create_pm(E_R, rdict, sizes[radius]; twisted=true)
+    @time pm = GroupRings.create_pm(E_R, rdict, sizes[halfradius]; twisted=true)
 
-    RG = GroupRing(parent(Id), E_R, rdict, pm)
+    RG = GroupRing(G, E_R, rdict, pm)
     Δ = spLaplacian(RG, S)
     return Δ
 end
@@ -56,7 +38,7 @@ function loadGRElem(fname::String, RG::GroupRing)
     return GroupRingElem(coeffs, RG)
 end
 
-function loadGRElem(fname::String, G::Group)
+function loadGRElem(fname::String, G::Union{Group, NCRing})
     pm = load(fname, "pm")
     RG = GroupRing(G, pm)
     return loadGRElem(fname, RG)

src/orbitdata.jl

@@ -28,7 +28,7 @@ function OrbitData(RG::GroupRing, autS::Group, verbose=true)
     @time Uπs = [orthSVD(matrix_repr(p, mreps)) for p in autS_mps]
 
     multiplicities = size.(Uπs,2)
-    dimensions = [Int(p[autS()]*Int(order(autS))) for p in autS_mps]
+    dimensions = [Int(p[one(autS)]*Int(order(autS))) for p in autS_mps]
     if verbose
         info_strs = ["",
         lpad("multiplicities", 14) * "  =" * join(lpad.(multiplicities, 4), ""),
@@ -213,12 +213,12 @@ function (g::perm)(y::GroupRingElem{T, <:SparseVector}) where T
     return result
 end
 
-function (p::perm)(A::MatElem)
+function (p::perm)(A::MatAlgElem)
     length(p.d) == size(A, 1) == size(A,2) || throw("Can't act via $p on matrix of size $(size(A))")
     result = similar(A)
     @inbounds for i in 1:size(A, 1)
         for j in 1:size(A, 2)
-            result[i, j] = A[p[i], p[j]] # action by permuting rows and colums/conjugation
+            result[p[i], p[j]] = A[i,j] # action by permuting rows and colums/conjugation
         end
     end
     return result
@@ -226,7 +226,7 @@ end
 
 ###############################################################################
 #
-#  WreathProductElems action on Nemo.MatElem
+#  WreathProductElems action on MatAlgElem
 #
 ###############################################################################
 
@@ -242,7 +242,7 @@ function (g::WreathProductElem)(y::GroupRingElem)
     return result
 end
 
-function (g::WreathProductElem{N})(A::MatElem) where N
+function (g::WreathProductElem{N})(A::MatAlgElem) where N
     # @assert N == size(A,1) == size(A,2)
     flips = ntuple(i->(g.n[i].d[1]==1 && g.n[i].d[2]==2 ? 1 : -1), N)
     result = similar(A)
@@ -251,11 +251,11 @@ function (g::WreathProductElem{N})(A::MatElem) where N
 
     @inbounds for i = 1:size(A,1)
         for j = 1:size(A,2)
-            x = A[g.p[i], g.p[j]]
+            x = A[i, j]
             if flips[i]*flips[j] == 1
-                result[i, j] = x
+                result[g.p[i], g.p[j]] = x
             else
-                result[i, j] = -x
+                result[g.p[i], g.p[j]] = -x
             end
             # result[i, j] = AbstractAlgebra.mul!(x, x, flips[i]*flips[j])
             # this mul! needs to be separately defined, but is 2x faster
@@ -273,8 +273,8 @@ end
 function AutFG_emb(A::AutGroup, g::WreathProductElem)
     isa(A.objectGroup, FreeGroup) || throw("Not an Aut(Fₙ)")
     parent(g).P.n == length(A.objectGroup.gens) || throw("No natural embedding of $(parent(g)) into $A")
-    elt = A()
+    elt = one(A)
-    Id = parent(g.n.elts[1])()
+    Id = one(parent(g.n.elts[1]))
     flips = Groups.AutSymbol[Groups.flip_autsymbol(i) for i in 1:length(g.p.d) if g.n.elts[i] != Id]
     Groups.r_multiply!(elt, flips, reduced=false)
     Groups.r_multiply!(elt, [Groups.perm_autsymbol(g.p)])

src/sqadjop.jl

@@ -0,0 +1,98 @@
+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{<:MatAlgebra}) = $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/1703.09680.jl

@@ -2,38 +2,50 @@
 
     @testset "SL(2,Z)" begin
         N = 2
-        G = MatrixSpace(Nemo.ZZ, N,N)
+        G = MatrixAlgebra(zz, N)
-        S = Groups.gens(G)
+        S = PropertyT.generating_set(G)
-        S = [S; inv.(S)]
 
         rm("SL($N,Z)", recursive=true, force=true)
-        sett = PropertyT.Settings("SL($N,Z)", G, S, solver(20000, accel=20); upper_bound=0.1)
+        sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(20000, accel=20); upper_bound=0.1)
-        
+
-        @test PropertyT.check_property_T(sett) == false
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ < 0.0
+        @test PropertyT.interpret_results(sett, λ) == false
     end
 
     @testset "SL(3,Z)" begin
         N = 3
-        G = MatrixSpace(Nemo.ZZ, N,N)
+        G = MatrixAlgebra(zz, N)
-        S = Groups.gens(G)
+        S = PropertyT.generating_set(G)
-        S = [S; inv.(S)]
+
-        
         rm("SL($N,Z)", recursive=true, force=true)
-        sett = PropertyT.Settings("SL($N,Z)", G, S, solver(1000, accel=20); upper_bound=0.1)
+        sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(1000, accel=20); upper_bound=0.1)
-        
+
-        @test PropertyT.check_property_T(sett) == true
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ > 0.0999
+        @test PropertyT.interpret_results(sett, λ) == true
+
+        @test PropertyT.check_property_T(sett) == true #second run should be fast
     end
-    
+
     @testset "SAut(F₂)" begin
         N = 2
         G = SAut(FreeGroup(N))
-        S = Groups.gens(G)
+        S = PropertyT.generating_set(G)
-        S = [S; inv.(S)]
+
-        
         rm("SAut(F$N)", recursive=true, force=true)
-        sett = PropertyT.Settings("SAut(F$N)", G, S, solver(20000);
+        sett = PropertyT.Settings("SAut(F$N)", G, S, with_SCS(10000);
         upper_bound=0.15, warmstart=false)
-        
+
-        @test PropertyT.check_property_T(sett) == false
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ < 0.0
+        @test PropertyT.interpret_results(sett, λ) == false
+
     end
 end

test/1712.07167.jl

@@ -2,58 +2,97 @@
 
     @testset "oSL(3,Z)" begin
         N = 3
-        G = MatrixSpace(Nemo.ZZ, N,N)
+        G = MatrixAlgebra(zz, N)
-        S = Groups.gens(G)
+        S = PropertyT.generating_set(G)
-        S = [S; inv.(S)]
         autS = WreathProduct(PermGroup(2), PermGroup(N))
-        
-        rm("oSL($N,Z)", recursive=true, force=true)
-        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, solver(2000, accel=20);
-        upper_bound=0.27, warmstart=false)
-       
-        @test PropertyT.check_property_T(sett) == false
-        #second run just checks the solution 
-        @test PropertyT.check_property_T(sett) == false 
 
-        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, solver(2000, accel=20);
+        rm("oSL($N,Z)", recursive=true, force=true)
+        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, with_SCS(2000, accel=20);
+        upper_bound=0.27, warmstart=false)
+
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ < 0.0
+        @test PropertyT.interpret_results(sett, λ) == false
+
+        # second run just checks the solution due to warmstart=false above
+        @test λ == PropertyT.spectral_gap(sett)
+        @test PropertyT.check_property_T(sett) == false
+
+        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, with_SCS(2000, accel=20);
         upper_bound=0.27, warmstart=true)
-        
+
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ > 0.269999
+        @test PropertyT.interpret_results(sett, λ) == true
+
+        # this should be very fast due to warmstarting:
+        @test λ ≈ PropertyT.spectral_gap(sett) atol=1e-5
         @test PropertyT.check_property_T(sett) == true
+
+        ##########
+        # Symmetrizing by PermGroup(3):
+
+        sett = PropertyT.Settings("SL($N,Z)", G, S, PermGroup(N), with_SCS(4000, accel=20);
+        upper_bound=0.27, warmstart=true)
+
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ > 0.269999
+        @test PropertyT.interpret_results(sett, λ) == true
     end
-    
+
     @testset "oSL(4,Z)" begin
         N = 4
-        G = MatrixSpace(Nemo.ZZ, N,N)
+        G = MatrixAlgebra(zz, N)
-        S = Groups.gens(G)
+        S = PropertyT.generating_set(G)
-        S = [S; inv.(S)]
         autS = WreathProduct(PermGroup(2), PermGroup(N))
-        
-        rm("oSL($N,Z)", recursive=true, force=true)
-        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, solver(2000, accel=20);
-        upper_bound=1.3, warmstart=false)
-        
-        @test PropertyT.check_property_T(sett) == false
-        #second run just checks the obtained solution
-        @test PropertyT.check_property_T(sett) == false 
 
-        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, solver(5000, accel=20);
+        rm("oSL($N,Z)", recursive=true, force=true)
+        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, with_SCS(2000, accel=20);
+        upper_bound=1.3, warmstart=false)
+
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ < 0.0
+        @test PropertyT.interpret_results(sett, λ) == false
+
+        # second run just checks the solution due to warmstart=false above
+        @test λ == PropertyT.spectral_gap(sett)
+        @test PropertyT.check_property_T(sett) == false
+
+        sett = PropertyT.Settings("SL($N,Z)", G, S, autS, with_SCS(5000, accel=20);
         upper_bound=1.3, warmstart=true)
-        
+
+        PropertyT.print_summary(sett)
+
+        λ = PropertyT.spectral_gap(sett)
+        @test λ > 1.2999
+        @test PropertyT.interpret_results(sett, λ) == true
+
+        # this should be very fast due to warmstarting:
+        @test λ ≈ PropertyT.spectral_gap(sett) atol=1e-5
         @test PropertyT.check_property_T(sett) == true
     end
 
     @testset "SAut(F₃)" begin
         N = 3
         G = SAut(FreeGroup(N))
-        S = Groups.gens(G)
+        S = PropertyT.generating_set(G)
-        S = [S; inv.(S)]
         autS = WreathProduct(PermGroup(2), PermGroup(N))
-        
+
         rm("oSAut(F$N)", recursive=true, force=true)
-        
+
-        sett = PropertyT.Settings("SAut(F$N)", G, S, autS, solver(5000);
+        sett = PropertyT.Settings("SAut(F$N)", G, S, autS, with_SCS(1000);
         upper_bound=0.15, warmstart=false)
-        
+
+        PropertyT.print_summary(sett)
+
         @test PropertyT.check_property_T(sett) == false
     end
 end

test/1812.03456.jl

@@ -0,0 +1,245 @@
+@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 = MatrixAlgebra(zz, N)
+
+            @test PropertyT.EltaryMat(M, 1, 2) isa MatAlgElem
+            e12 = PropertyT.EltaryMat(M, 1, 2)
+            @test e12[1,2] == 1
+            @test inv(e12)[1,2] == -1
+
+            S = PropertyT.generating_set(M)
+            @test e12 ∈ S
+
+            @test length(PropertyT.generating_set(M)) == 2N*(N-1)
+            @test all(inv(s) ∈ S for s in S)
+
+            A = SAut(FreeGroup(N))
+            @test length(PropertyT.generating_set(A)) == 4N*(N-1)
+            S = PropertyT.generating_set(A)
+            @test all(inv(s) ∈ S for s in S)
+        end
+
+        @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 = MatrixAlgebra(zz, 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.EltaryMat(M, 1,2)
+        h = PropertyT.EltaryMat(M, 1,3)
+        k = PropertyT.EltaryMat(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
+
+
+@testset "1812.03456 examples" begin
+
+    function SOS_residual(x::GroupRingElem, Q::Matrix)
+        RG = parent(x)
+        @time sos = PropertyT.compute_SOS(RG, Q);
+        return x - sos
+    end
+
+    function check_positivity(elt, Δ, orbit_data, upper_bound, warm=nothing; with_solver=with_SCS(20_000, accel=10))
+        SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=upper_bound)
+
+        status, warm = PropertyT.solve(SDP_problem, with_solver, warm);
+        Base.Libc.flush_cstdio()
+        @info "Optimization status:" status
+
+        λ = value(SDP_problem[:λ])
+        Ps = [value.(P) for P in varP]
+
+        Qs = real.(sqrt.(Ps));
+        Q = PropertyT.reconstruct(Qs, orbit_data);
+
+        b = SOS_residual(elt - λ*Δ, Q)
+        return b, λ, warm
+    end
+
+    @testset "SL(3,Z)" begin
+        N = 3
+        halfradius = 2
+        M = MatrixAlgebra(zz, N)
+        S = PropertyT.generating_set(M)
+        Δ = PropertyT.Laplacian(S, halfradius)
+        RG = parent(Δ)
+        orbit_data = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
+        orbit_data = PropertyT.decimate(orbit_data);
+
+        @testset "Sq₃ is SOS" begin
+            elt = PropertyT.Sq(RG)
+            UB = 0.05 # 0.105?
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
+            Base.Libc.flush_cstdio()
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
+
+            @test 2^2*norm(residual, 1) < λ/100
+        end
+
+        @testset "Adj₃ is SOS" begin
+            elt = PropertyT.Adj(RG)
+            UB = 0.1 # 0.157?
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
+            Base.Libc.flush_cstdio()
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) < λ
+            @test 2^2*norm(residual, 1) < λ/100
+        end
+
+        @testset "Op₃ is empty, so can not be certified" begin
+            elt = PropertyT.Op(RG)
+            UB = Inf
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
+            Base.Libc.flush_cstdio()
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) > λ
+        end
+    end
+
+    @testset "SL(4,Z)" begin
+        N = 4
+        halfradius = 2
+        M = MatrixAlgebra(zz, N)
+        S = PropertyT.generating_set(M)
+        Δ = PropertyT.Laplacian(S, halfradius)
+        RG = parent(Δ)
+        orbit_data = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
+        orbit_data = PropertyT.decimate(orbit_data);
+
+        @testset "Sq₄ is SOS" begin
+            elt = PropertyT.Sq(RG)
+            UB = 0.2 # 0.3172
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
+            Base.Libc.flush_cstdio()
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
+            @test 2^2*norm(residual, 1) < λ/100
+        end
+
+        @testset "Adj₄ is SOS" begin
+            elt = PropertyT.Adj(RG)
+            UB = 0.3 # 0.5459?
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
+            @test 2^2*norm(residual, 1) < λ/100
+        end
+
+        @testset "we can't cerify that Op₄ SOS" begin
+            elt = PropertyT.Op(RG)
+            UB = 2.0
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB,
+            with_solver=with_SCS(20_000, accel=10, eps=2e-10))
+            Base.Libc.flush_cstdio()
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) > λ # i.e. we can certify positivity
+        end
+
+        @testset "Adj₄ + Op₄ is SOS" begin
+            elt = PropertyT.Adj(RG) + PropertyT.Op(RG)
+            UB = 0.6 # 0.82005
+
+            residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
+            Base.Libc.flush_cstdio()
+            @info "obtained λ and residual" λ norm(residual, 1)
+
+            @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
+            @test 2^2*norm(residual, 1) < λ/100
+        end
+    end
+
+    # @testset "Adj₄ + 100 Op₄ ∈ ISAut(F₄) is SOS" begin
+    #     N = 4
+    #     halfradius = 2
+    #     M = SAut(FreeGroup(N))
+    #     S = PropertyT.generating_set(M)
+    #     Δ = PropertyT.Laplacian(S, halfradius)
+    #     RG = parent(Δ)
+    #     orbit_data = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
+    #     orbit_data = PropertyT.decimate(orbit_data);
+    #
+    #     @time elt = PropertyT.Adj(RG) + 100PropertyT.Op(RG)
+    #     UB = 0.05
+    #
+    #     warm = nothing
+    #
+    #     residual, λ, warm = check_positivity(elt, Δ, orbit_data, UB, warm, with_solver=with_SCS(20_000, accel=10))
+    #     @info "obtained λ and residual" λ norm(residual, 1)
+    #
+    #     @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
+    #     @test 2^2*norm(residual, 1) < λ/100
+    # end
+end

test/SOS_correctness.jl

@@ -1,5 +1,3 @@
-using PropertyT.GroupRings
-
 @testset "Correctness of HPC SOS computation" begin
 
     function prepare(G_name, λ, S_size)
@@ -35,7 +33,7 @@ using PropertyT.GroupRings
     @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
     @time sos_hpc = PropertyT.compute_SOS(pm, Q)
 
-    @test norm(sos_sqr - sos_hpc, 1) < 3e-12
+    @test norm(sos_sqr - sos_hpc, 1) < 4e-12
     @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
 
     #########################################################

test/actions.jl

@@ -0,0 +1,104 @@
+@testset "actions on Group[Rings]" begin
+    Eij = PropertyT.EltaryMat
+    ssgs(M::MatAlgebra, i, j) = (S = [Eij(M, i, j), Eij(M, j, i)];
+        S = unique([S; inv.(S)]); S)
+
+    rmul = Groups.rmul_autsymbol
+    lmul = Groups.lmul_autsymbol
+
+    function ssgs(A::AutGroup, i, j)
+        rmuls = [rmul(i,j), rmul(j,i)]
+        lmuls = [lmul(i,j), lmul(j,i)]
+        gen_set = A.([rmuls; lmuls])
+        return unique([gen_set; inv.(gen_set)])
+    end
+
+@testset "actions on SL(3,Z) and its group ring" begin
+    N = 3
+    halfradius = 2
+    M = MatrixAlgebra(zz, N)
+    S = PropertyT.generating_set(M)
+    E_R, sizes = Groups.generate_balls(S, one(M), radius=2halfradius);
+
+    rdict = GroupRings.reverse_dict(E_R)
+    pm = GroupRings.create_pm(E_R, rdict, sizes[halfradius]; twisted=false);
+    RG = GroupRing(M, E_R, rdict, pm)
+
+    @testset "correctness of actions" begin
+        Δ = length(S)*RG(1) - sum(RG(s) for s in S)
+        @test Δ == PropertyT.spLaplacian(RG, S)
+
+        elt = S[5]
+        x = RG(1) - RG(elt)
+        elt2 = E_R[rand(sizes[1]:sizes[2])]
+        y = 2RG(elt2) - RG(elt)
+
+        for G in [PermGroup(N), WreathProduct(PermGroup(2), PermGroup(N))]
+            @test all(g(one(M)) == one(M) for g in G)
+            @test all(rdict[g(m)] <= sizes[1] for g in G for m in S)
+            @test all(g(m)*g(n) == g(m*n) for g in G for m in S for n in S)
+
+            @test all(g(Δ) == Δ for g in G)
+            @test all(g(x) == RG(1) - RG(g(elt)) for g in G)
+
+            @test all(2RG(g(elt2)) - RG(g(elt)) == g(y) for g in G)
+        end
+    end
+
+    @testset "small Laplacians" begin
+        for (i,j) in PropertyT.indexing(N)
+            Sij = ssgs(M, i,j)
+            Δij= PropertyT.spLaplacian(RG, Sij)
+
+            @test all(p(Δij) == PropertyT.spLaplacian(RG, ssgs(M, p[i], p[j])) for p in PermGroup(N))
+
+            @test all(g(Δij) == PropertyT.spLaplacian(RG, ssgs(M, g.p[i], g.p[j])) for g in WreathProduct(PermGroup(2), PermGroup(N)))
+        end
+    end
+end
+
+@testset "actions on SAut(F_3) and its group ring" begin
+    N = 3
+    halfradius = 2
+    M = SAut(FreeGroup(N))
+    S = PropertyT.generating_set(M)
+    E_R, sizes = Groups.generate_balls(S, one(M), radius=2halfradius);
+
+    rdict = GroupRings.reverse_dict(E_R)
+    pm = GroupRings.create_pm(E_R, rdict, sizes[halfradius]; twisted=false);
+    RG = GroupRing(M, E_R, rdict, pm)
+
+
+    @testset "correctness of actions" begin
+
+        Δ = length(S)*RG(1) - sum(RG(s) for s in S)
+        @test Δ == PropertyT.spLaplacian(RG, S)
+
+        elt = S[5]
+        x = RG(1) - RG(elt)
+        elt2 = E_R[rand(sizes[1]:sizes[2])]
+        y = 2RG(elt2) - RG(elt)
+
+        for G in [PermGroup(N), WreathProduct(PermGroup(2), PermGroup(N))]
+            @test all(g(one(M)) == one(M) for g in G)
+            @test all(rdict[g(m)] <= sizes[1] for g in G for m in S)
+            @test all(g(m)*g(n) == g(m*n) for g in G for m in S for n in S)
+
+            @test all(g(Δ) == Δ for g in G)
+            @test all(g(x) == RG(1) - RG(g(elt)) for g in G)
+
+            @test all(2RG(g(elt2)) - RG(g(elt)) == g(y) for g in G)
+        end
+    end
+
+    for (i,j) in PropertyT.indexing(N)
+        Sij = ssgs(M, i,j)
+        Δij= PropertyT.spLaplacian(RG, Sij)
+
+        @test all(p(Δij) == PropertyT.spLaplacian(RG, ssgs(M, p[i], p[j])) for p in PermGroup(N))
+
+        @test all(g(Δij) == PropertyT.spLaplacian(RG, ssgs(M, g.p[i], g.p[j])) for g in WreathProduct(PermGroup(2), PermGroup(N)))
+    end
+end
+
+end

test/runtests.jl

@@ -1,25 +1,18 @@
-using AbstractAlgebra, Nemo, Groups, SCS
-using SparseArrays
-using JLD
-using PropertyT
 using Test
-using JuMP
+using LinearAlgebra, SparseArrays
+using AbstractAlgebra, Groups, GroupRings
+using PropertyT
+using JLD
 
-indexing(n) = [(i,j) for i in 1:n for j in (i+1):n]
+using JuMP, SCS
-function Groups.gens(M::MatSpace)
-    @assert ncols(M) == nrows(M)
-    N = ncols(M)
-    E(i,j) = begin g = M(1); g[i,j] = 1; g end
-    S = [E(i,j) for (i,j) in indexing(N)]
-    S = [S; transpose.(S)]
-    return S
-end
 
-solver(iters; accel=1) =
+with_SCS(iters; accel=1, eps=1e-10) =
     with_optimizer(SCS.Optimizer,
     linear_solver=SCS.Direct, max_iters=iters,
-    acceleration_lookback=accel, eps=1e-10, warm_start=true)
+    acceleration_lookback=accel, eps=eps, warm_start=true)
 
 include("1703.09680.jl")
+include("actions.jl")
 include("1712.07167.jl")
 include("SOS_correctness.jl")
+include("1812.03456.jl")