diff --git a/.travis.yml b/.travis.yml index 4c7dc0b..8084d03 100644 --- a/.travis.yml +++ b/.travis.yml @@ -6,6 +6,7 @@ os: julia: - 1.0 - 1.1 + - 1.2 - nightly notifications: email: true diff --git a/Manifest.toml b/Manifest.toml index 199472e..c7cfdbe 100644 --- a/Manifest.toml +++ b/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"] -git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648" +deps = ["Libdl", "Logging", "SHA"] +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" diff --git a/Project.toml b/Project.toml index 21135bd..ed05ad8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PropertyT" uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d" authors = ["Marek Kaluba "] -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" diff --git a/src/1712.07167.jl b/src/1712.07167.jl index 5b4afa0..1d246f5 100644 --- a/src/1712.07167.jl +++ b/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; - radius::Integer=2, upper_bound::Float64=1.0, warmstart=true) - return Naive(name, G, S, radius, upper_bound, solver, warmstart) + G::Union{Group, NCRing}, S::AbstractVector{El}, solver::JuMP.OptimizerFactory; + halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) where El <: Union{GroupElem, NCRingElem} + 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; - radius::Integer=2, upper_bound::Float64=1.0, warmstart=true) - return Symmetrized(name, G, S, autS, radius, upper_bound, solver, warmstart) + G::Union{Group, NCRing}, S::AbstractVector{El}, autS::Group, solver::JuMP.OptimizerFactory; + halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) where El <: Union{GroupElem, NCRingElem} + 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 - @info "Couldn't load $(filename(sett, :warmstart))." + catch ex + @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) - catch - # compute - Δ = Laplacian(sett.S, sett.radius) + Δ = loadGRElem(filename(sett,:Δ), sett.G) + @info "Loaded precomputed Δ." + Δ + catch ex + @warn ex.msg + @info "Computing Δ..." + Δ = Laplacian(sett.S, sett.halfradius) saveGRElem(filename(sett, :Δ), Δ) Δ end - λ, P = try - sett.warmstart && error() - load(filename(sett, :solution), "λ", "P") - catch + function compute(sett, Δ) + @info "Computing λ and P..." λ, 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 diff --git a/src/1812.03456.jl b/src/1812.03456.jl new file mode 100644 index 0000000..bf36e78 --- /dev/null +++ b/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") diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 247a928..53199b3 100644 --- a/src/PropertyT.jl +++ b/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) diff --git a/src/RGprojections.jl b/src/RGprojections.jl index 4a470ab..ffa3f11 100644 --- a/src/RGprojections.jl +++ b/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))) - last_emb = g->Bn(Generic.emb!(Bn.P(), g, view(r, (i+1):N))) + first_emb = g->Bn(Generic.emb!(one(Bn.P), g, view(r, 1:i))) + 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 diff --git a/src/laplacians.jl b/src/laplacians.jl index d30bc13..fbf6d8a 100644 --- a/src/laplacians.jl +++ b/src/laplacians.jl @@ -4,16 +4,7 @@ # ############################################################################### -function spLaplacian(RG::GroupRing, S, 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} +function spLaplacian(RG::GroupRing, S::AbstractVector, T::Type=Float64) 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 - R = parent(first(S)) - return Laplacian(S, one(R), radius) -end - -function Laplacian(S::Vector{E}, radius) where E<:AbstractAlgebra.GroupElem +function Laplacian(S::AbstractVector{REl}, halfradius) where REl<:Union{NCRingElem, GroupElem} G = parent(first(S)) - return Laplacian(S, G(), radius) -end - -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 "Generating metric ball of radius" radius=2halfradius + @time E_R, sizes = Groups.generate_balls(S, one(G), radius=2halfradius) @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) diff --git a/src/orbitdata.jl b/src/orbitdata.jl index 815b049..1e8c189 100644 --- a/src/orbitdata.jl +++ b/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() - Id = parent(g.n.elts[1])() + elt = one(A) + 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)]) diff --git a/src/sqadjop.jl b/src/sqadjop.jl new file mode 100644 index 0000000..b516d10 --- /dev/null +++ b/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 diff --git a/test/1703.09680.jl b/test/1703.09680.jl index ddc8fc3..f8c66af 100644 --- a/test/1703.09680.jl +++ b/test/1703.09680.jl @@ -2,38 +2,50 @@ @testset "SL(2,Z)" begin N = 2 - G = MatrixSpace(Nemo.ZZ, N,N) - S = Groups.gens(G) - S = [S; inv.(S)] + G = MatrixAlgebra(zz, N) + S = PropertyT.generating_set(G) rm("SL($N,Z)", recursive=true, force=true) - sett = PropertyT.Settings("SL($N,Z)", G, S, solver(20000, accel=20); upper_bound=0.1) - - @test PropertyT.check_property_T(sett) == false + sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(20000, accel=20); upper_bound=0.1) + + 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) - S = Groups.gens(G) - S = [S; inv.(S)] - + G = MatrixAlgebra(zz, N) + S = PropertyT.generating_set(G) + rm("SL($N,Z)", recursive=true, force=true) - sett = PropertyT.Settings("SL($N,Z)", G, S, solver(1000, accel=20); upper_bound=0.1) - - @test PropertyT.check_property_T(sett) == true + sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(1000, accel=20); upper_bound=0.1) + + 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 = [S; inv.(S)] - + S = PropertyT.generating_set(G) + 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 diff --git a/test/1712.07167.jl b/test/1712.07167.jl index 2c08b05..259836a 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -2,58 +2,97 @@ @testset "oSL(3,Z)" begin N = 3 - G = MatrixSpace(Nemo.ZZ, N,N) - S = Groups.gens(G) - S = [S; inv.(S)] + G = MatrixAlgebra(zz, N) + S = PropertyT.generating_set(G) 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) - S = Groups.gens(G) - S = [S; inv.(S)] + G = MatrixAlgebra(zz, N) + S = PropertyT.generating_set(G) 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 = [S; inv.(S)] + S = PropertyT.generating_set(G) 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 diff --git a/test/1812.03456.jl b/test/1812.03456.jl new file mode 100644 index 0000000..7b6b121 --- /dev/null +++ b/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 diff --git a/test/SOS_correctness.jl b/test/SOS_correctness.jl index f991ee9..f76f360 100644 --- a/test/SOS_correctness.jl +++ b/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) ######################################################### diff --git a/test/actions.jl b/test/actions.jl new file mode 100644 index 0000000..fd9ce04 --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 773fe0a..1a0e85f 100644 --- a/test/runtests.jl +++ b/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] -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 +using JuMP, SCS -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")