From 4c4ef195e11ccc93be01df08f7735c862628b22e Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 30 Jun 2019 11:56:57 +0200 Subject: [PATCH 01/14] replace radius by halfradius where appropriate --- src/1712.07167.jl | 18 +++++++++--------- src/laplacians.jl | 16 ++++++++-------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/1712.07167.jl b/src/1712.07167.jl index 5b4afa0..308b489 100644 --- a/src/1712.07167.jl +++ b/src/1712.07167.jl @@ -12,7 +12,7 @@ struct Naive{El} <: Settings name::String G::Group S::Vector{El} - radius::Int + halfradius::Int upper_bound::Float64 solver::JuMP.OptimizerFactory @@ -24,7 +24,7 @@ struct Symmetrized{El} <: Settings G::Group S::Vector{El} autS::Group - radius::Int + halfradius::Int upper_bound::Float64 solver::JuMP.OptimizerFactory @@ -33,14 +33,14 @@ 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) + halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) + 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) + halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) + return Symmetrized(name, G, S, autS, halfradius, upper_bound, solver, warmstart) end prefix(s::Naive) = s.name @@ -245,7 +245,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()))", @@ -279,7 +279,7 @@ function spectral_gap(sett::Settings) loadGRElem(filename(sett,:Δ), sett.G) catch # compute - Δ = Laplacian(sett.S, sett.radius) + Δ = Laplacian(sett.S, sett.halfradius) saveGRElem(filename(sett, :Δ), Δ) Δ end @@ -307,7 +307,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/laplacians.jl b/src/laplacians.jl index d30bc13..c9538af 100644 --- a/src/laplacians.jl +++ b/src/laplacians.jl @@ -22,24 +22,24 @@ 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::Vector{E}, halfradius) where E<:AbstractAlgebra.ModuleElem R = parent(first(S)) - return Laplacian(S, one(R), radius) + return Laplacian(S, one(R), halfradius) end -function Laplacian(S::Vector{E}, radius) where E<:AbstractAlgebra.GroupElem +function Laplacian(S::Vector{E}, halfradius) where E<:AbstractAlgebra.GroupElem G = parent(first(S)) - return Laplacian(S, G(), radius) + return Laplacian(S, G(), halfradius) 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) +function Laplacian(S, Id, halfradius) + @info "Generating metric ball of radius" radius=2halfradius + @time E_R, sizes = Groups.generate_balls(S, Id, 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) Δ = spLaplacian(RG, S) From b9dc701f1700165a773c37a1a5d39fe7db9143cc Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 30 Jun 2019 11:58:32 +0200 Subject: [PATCH 02/14] when acting on MatElem always act on the lhs columns/rows --- src/orbitdata.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/orbitdata.jl b/src/orbitdata.jl index 815b049..116e23d 100644 --- a/src/orbitdata.jl +++ b/src/orbitdata.jl @@ -218,7 +218,7 @@ function (p::perm)(A::MatElem) 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 @@ -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 From 6c906b05cbc95cabad20ec62caac2324e874137c Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 30 Jun 2019 13:19:24 +0200 Subject: [PATCH 03/14] add sqadjop.jl and unit tests --- src/PropertyT.jl | 3 ++ src/sqadjop.jl | 125 +++++++++++++++++++++++++++++++++++++++++++++ test/1812.03456.jl | 93 +++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 4 files changed, 222 insertions(+) create mode 100644 src/sqadjop.jl create mode 100644 test/1812.03456.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 247a928..038bcc9 100644 --- a/src/PropertyT.jl +++ b/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) diff --git a/src/sqadjop.jl b/src/sqadjop.jl new file mode 100644 index 0000000..d060b57 --- /dev/null +++ b/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 diff --git a/test/1812.03456.jl b/test/1812.03456.jl new file mode 100644 index 0000000..da55d05 --- /dev/null +++ b/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 + diff --git a/test/runtests.jl b/test/runtests.jl index 773fe0a..90c6063 100644 --- a/test/runtests.jl +++ b/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") From 80e338f19130488b87e560458c49fe6e84a0d13e Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 30 Jun 2019 13:20:39 +0200 Subject: [PATCH 04/14] add 1812.03456 positivity tests in SL(n,Z) --- test/1812.03456.jl | 132 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) diff --git a/test/1812.03456.jl b/test/1812.03456.jl index da55d05..524643e 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -91,3 +91,135 @@ end end + +@testset "1812.03456 examples" begin + + with_SCS = with_optimizer(SCS.Optimizer, + linear_solver=SCS.Direct, + eps=2e-10, + max_iters=20000, + alpha=1.5, + acceleration_lookback=10, + warm_start=true) + + 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) + SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=upper_bound) + + status, warm = PropertyT.solve(SDP_problem, with_solver, warm); + @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 = MatrixSpace(Nemo.ZZ, N,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) + @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) + @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) + @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 = MatrixSpace(Nemo.ZZ, N,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) + @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) + @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) + @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 + +end From 88b460a95972a008a19a9477e2fdc8920904d66a Mon Sep 17 00:00:00 2001 From: kalmarek Date: Sun, 30 Jun 2019 19:05:46 +0200 Subject: [PATCH 05/14] use NCRing[Elem]s instead of MatSpace and use NCRing in Settings --- src/1712.07167.jl | 12 ++++++------ src/PropertyT.jl | 2 +- src/laplacians.jl | 20 ++++++-------------- 3 files changed, 13 insertions(+), 21 deletions(-) diff --git a/src/1712.07167.jl b/src/1712.07167.jl index 308b489..df0f51f 100644 --- a/src/1712.07167.jl +++ b/src/1712.07167.jl @@ -10,7 +10,7 @@ abstract type Settings end struct Naive{El} <: Settings name::String - G::Group + G::Union{Group, NCRing} S::Vector{El} halfradius::Int upper_bound::Float64 @@ -21,7 +21,7 @@ end struct Symmetrized{El} <: Settings name::String - G::Group + G::Union{Group, NCRing} S::Vector{El} autS::Group halfradius::Int @@ -32,14 +32,14 @@ struct Symmetrized{El} <: Settings end function Settings(name::String, - G::Group, S::Vector{<:GroupElem}, solver::JuMP.OptimizerFactory; - halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) + 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; - halfradius::Integer=2, upper_bound::Float64=1.0, warmstart=true) + 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 diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 038bcc9..f082f5c 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -13,7 +13,7 @@ using GroupRings using JLD using JuMP -import AbstractAlgebra: Group, Ring, perm +import AbstractAlgebra: Group, NCRing, perm import MathProgBase.SolverInterface.AbstractMathProgSolver diff --git a/src/laplacians.jl b/src/laplacians.jl index c9538af..5df4cb7 100644 --- a/src/laplacians.jl +++ b/src/laplacians.jl @@ -4,30 +4,22 @@ # ############################################################################### -function spLaplacian(RG::GroupRing, S, T::Type=Float64) +function spLaplacian(RG::GroupRing, S::AbstractVector{El}, T::Type=Float64) where El result = RG(T) - result[RG.group()] = T(length(S)) + id = (El <: AbstractAlgebra.NCRingElem ? one(RG.group) : RG.group()) + result[id] = 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 - result[s] -= one(T) - end - return result -end - -function Laplacian(S::Vector{E}, halfradius) where E<:AbstractAlgebra.ModuleElem +function Laplacian(S::AbstractVector{REl}, halfradius) where REl<:AbstractAlgebra.NCRingElem R = parent(first(S)) return Laplacian(S, one(R), halfradius) end -function Laplacian(S::Vector{E}, halfradius) where E<:AbstractAlgebra.GroupElem +function Laplacian(S::AbstractVector{E}, halfradius) where E<:AbstractAlgebra.GroupElem G = parent(first(S)) return Laplacian(S, G(), halfradius) end @@ -56,7 +48,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) From cc9d3db8464a186e125255a3e6257832b04d66ca Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 1 Jul 2019 01:37:33 +0200 Subject: [PATCH 06/14] we act on MatAlgElem now; use MatAlgebra --- src/orbitdata.jl | 6 +++--- src/sqadjop.jl | 10 ++++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/orbitdata.jl b/src/orbitdata.jl index 116e23d..6818c06 100644 --- a/src/orbitdata.jl +++ b/src/orbitdata.jl @@ -213,7 +213,7 @@ 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) @@ -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) diff --git a/src/sqadjop.jl b/src/sqadjop.jl index d060b57..2fa94e1 100644 --- a/src/sqadjop.jl +++ b/src/sqadjop.jl @@ -9,7 +9,7 @@ function generating_set(G::AutGroup{N}, n=N) where N return [gen_set; inv.(gen_set)] end -function E(M::MatSpace, i::Integer, j::Integer, val=1) +function E(M::MatAlgebra, i::Integer, j::Integer, val=1) @assert i ≠ j @assert 1 ≤ i ≤ nrows(M) @assert 1 ≤ j ≤ ncols(M) @@ -18,8 +18,7 @@ function E(M::MatSpace, i::Integer, j::Integer, val=1) return m end -function generating_set(M::MatSpace, n=nrows(M)) - @assert nrows(M) == ncols(M) +function generating_set(M::MatAlgebra, n=M.n) elts = [E(M, i,j) for (i,j) in indexing(n)] return elem_type(M)[elts; inv.(elts)] @@ -92,10 +91,13 @@ function Op(RG::GroupRing, N::Integer) return op÷factorial(N-2)^2 end +AbstractAlgebra.nrows(M::MatAlgebra) = M.n +AbstractAlgebra.ncols(M::MatAlgebra) = M.n + 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)) + $Elt(RG::GroupRing{<:MatAlgebra}) = $Elt(RG, nrows(RG.group)) end end From 75e997b2e8d9ecd7a8bf9f96d679e11db2b852cb Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 1 Jul 2019 01:38:30 +0200 Subject: [PATCH 07/14] update tests and test spectral_gap directly --- test/1703.09680.jl | 52 ++++++++++++++--------- test/1712.07167.jl | 93 ++++++++++++++++++++++++++--------------- test/1812.03456.jl | 53 ++++++++++++----------- test/SOS_correctness.jl | 2 - test/runtests.jl | 23 ++++------ 5 files changed, 125 insertions(+), 98 deletions(-) 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..b1bd59e 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -2,58 +2,85 @@ @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 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 index 524643e..f674227 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -7,26 +7,25 @@ @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) + M = MatrixAlgebra(zz, N) + + @test PropertyT.E(M, 1, 2) isa MatAlgElem + e12 = PropertyT.E(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 - 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)") @@ -40,7 +39,7 @@ @testset "Sq, Adj, Op" begin N = 4 - M = MatrixSpace(Nemo.ZZ, N,N) + M = MatrixAlgebra(zz, N) S = PropertyT.generating_set(M) Δ = PropertyT.Laplacian(S, 2) RG = parent(Δ) @@ -94,24 +93,17 @@ end @testset "1812.03456 examples" begin - with_SCS = with_optimizer(SCS.Optimizer, - linear_solver=SCS.Direct, - eps=2e-10, - max_iters=20000, - alpha=1.5, - acceleration_lookback=10, - warm_start=true) - - function SOS_residual(x::GroupRingElem, Q::Matrix) + 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) + 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[:λ]) @@ -127,7 +119,7 @@ end @testset "SL(3,Z)" begin N = 3 halfradius = 2 - M = MatrixSpace(Nemo.ZZ, N,N) + M = MatrixAlgebra(zz, N) S = PropertyT.generating_set(M) Δ = PropertyT.Laplacian(S, halfradius) RG = parent(Δ) @@ -139,6 +131,7 @@ end 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 @@ -151,6 +144,7 @@ end 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) < λ @@ -162,6 +156,7 @@ end 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) > λ @@ -171,7 +166,7 @@ end @testset "SL(4,Z)" begin N = 4 halfradius = 2 - M = MatrixSpace(Nemo.ZZ, N,N) + M = MatrixAlgebra(zz, N) S = PropertyT.generating_set(M) Δ = PropertyT.Laplacian(S, halfradius) RG = parent(Δ) @@ -183,6 +178,7 @@ end 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 @@ -204,7 +200,9 @@ end elt = PropertyT.Op(RG) UB = 2.0 - residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB) + 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 @@ -215,6 +213,7 @@ end 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 diff --git a/test/SOS_correctness.jl b/test/SOS_correctness.jl index f991ee9..e31fb77 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) diff --git a/test/runtests.jl b/test/runtests.jl index 90c6063..96bb5e0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,24 +1,15 @@ -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("1712.07167.jl") From 2ef8de7d42fdc59a30d27b0e1a9fcc5d9c22653f Mon Sep 17 00:00:00 2001 From: kalmarek Date: Mon, 1 Jul 2019 01:46:32 +0200 Subject: [PATCH 08/14] =?UTF-8?q?add=20commented=20test=20on=20Adj?= =?UTF-8?q?=E2=82=84=20+=20100=20Op=E2=82=84=20in=20SAut(F=E2=82=84)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test/1812.03456.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/1812.03456.jl b/test/1812.03456.jl index f674227..77c9191 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -221,4 +221,25 @@ end 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 From e9bb6f13ddacacf846bb65d32ca08bd4ee381c34 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 5 Jul 2019 18:57:39 +0200 Subject: [PATCH 09/14] add tests for actions --- test/1712.07167.jl | 12 ++++++ test/actions.jl | 104 +++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 3 files changed, 117 insertions(+) create mode 100644 test/actions.jl diff --git a/test/1712.07167.jl b/test/1712.07167.jl index b1bd59e..259836a 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -32,6 +32,18 @@ # 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 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 96bb5e0..1a0e85f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,6 +12,7 @@ with_SCS(iters; accel=1, eps=1e-10) = 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") From e23b6a85dc886ccf9237e21140fe1758f49ad3a9 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 4 Jul 2019 15:31:55 +0200 Subject: [PATCH 10/14] include sqadjop.jl from 1812.03456.jl file --- src/1812.03456.jl | 26 ++++++++++++++++++++++++++ src/PropertyT.jl | 3 +-- src/sqadjop.jl | 29 ----------------------------- test/1812.03456.jl | 12 ++++++------ 4 files changed, 33 insertions(+), 37 deletions(-) create mode 100644 src/1812.03456.jl 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 f082f5c..415d8a3 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -22,9 +22,8 @@ include("RGprojections.jl") include("orbitdata.jl") include("sos_sdps.jl") include("checksolution.jl") -include("sqadjop.jl") include("1712.07167.jl") - +include("1812.03456.jl") end # module Property(T) diff --git a/src/sqadjop.jl b/src/sqadjop.jl index 2fa94e1..b516d10 100644 --- a/src/sqadjop.jl +++ b/src/sqadjop.jl @@ -1,29 +1,3 @@ -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::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=M.n) - - 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] @@ -91,9 +65,6 @@ function Op(RG::GroupRing, N::Integer) return op÷factorial(N-2)^2 end -AbstractAlgebra.nrows(M::MatAlgebra) = M.n -AbstractAlgebra.ncols(M::MatAlgebra) = M.n - for Elt in [:Sq, :Adj, :Op] @eval begin $Elt(RG::GroupRing{AutGroup{N}}) where N = $Elt(RG, N) diff --git a/test/1812.03456.jl b/test/1812.03456.jl index 77c9191..7b6b121 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -9,8 +9,8 @@ for N in [3,4] M = MatrixAlgebra(zz, N) - @test PropertyT.E(M, 1, 2) isa MatAlgElem - e12 = PropertyT.E(M, 1, 2) + @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 @@ -65,9 +65,9 @@ @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) + 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 @@ -93,7 +93,7 @@ end @testset "1812.03456 examples" begin - function SOS_residual(x::GroupRingElem, Q::Matrix) + function SOS_residual(x::GroupRingElem, Q::Matrix) RG = parent(x) @time sos = PropertyT.compute_SOS(RG, Q); return x - sos From 4bce106c2d62cd87f19d38962191e1f2da2767b4 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 4 Jul 2019 22:44:56 +0200 Subject: [PATCH 11/14] better logging --- src/1712.07167.jl | 44 ++++++++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/src/1712.07167.jl b/src/1712.07167.jl index df0f51f..1d246f5 100644 --- a/src/1712.07167.jl +++ b/src/1712.07167.jl @@ -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 @@ -275,26 +278,39 @@ function spectral_gap(sett::Settings) isdir(fp) || mkpath(fp) Δ = try - @info "Loading precomputed Δ..." - loadGRElem(filename(sett,:Δ), sett.G) - catch - # compute + Δ = 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:", From 5e821286e5310cfa7d6cfdde8cc01d845171f513 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 5 Jul 2019 18:55:25 +0200 Subject: [PATCH 12/14] define and use one(::Group) --- src/PropertyT.jl | 2 ++ src/RGprojections.jl | 8 ++++---- src/laplacians.jl | 20 +++++--------------- src/orbitdata.jl | 6 +++--- 4 files changed, 14 insertions(+), 22 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 415d8a3..53199b3 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -17,6 +17,8 @@ import AbstractAlgebra: Group, NCRing, perm import MathProgBase.SolverInterface.AbstractMathProgSolver +AbstractAlgebra.one(G::Group) = G() + include("laplacians.jl") include("RGprojections.jl") include("orbitdata.jl") 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 5df4cb7..fbf6d8a 100644 --- a/src/laplacians.jl +++ b/src/laplacians.jl @@ -4,36 +4,26 @@ # ############################################################################### -function spLaplacian(RG::GroupRing, S::AbstractVector{El}, T::Type=Float64) where El +function spLaplacian(RG::GroupRing, S::AbstractVector, T::Type=Float64) result = RG(T) - id = (El <: AbstractAlgebra.NCRingElem ? one(RG.group) : RG.group()) - result[id] = T(length(S)) + result[one(RG.group)] = T(length(S)) for s in S result[s] -= one(T) end return result end -function Laplacian(S::AbstractVector{REl}, halfradius) where REl<:AbstractAlgebra.NCRingElem - R = parent(first(S)) - return Laplacian(S, one(R), halfradius) -end - -function Laplacian(S::AbstractVector{E}, halfradius) where E<:AbstractAlgebra.GroupElem +function Laplacian(S::AbstractVector{REl}, halfradius) where REl<:Union{NCRingElem, GroupElem} G = parent(first(S)) - return Laplacian(S, G(), halfradius) -end - -function Laplacian(S, Id, halfradius) @info "Generating metric ball of radius" radius=2halfradius - @time E_R, sizes = Groups.generate_balls(S, Id, 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[halfradius]; twisted=true) - RG = GroupRing(parent(Id), E_R, rdict, pm) + RG = GroupRing(G, E_R, rdict, pm) Δ = spLaplacian(RG, S) return Δ end diff --git a/src/orbitdata.jl b/src/orbitdata.jl index 6818c06..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), ""), @@ -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)]) From f39680042d5d713cf0fe84f42f9e85294146a1d2 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Fri, 5 Jul 2019 19:22:49 +0200 Subject: [PATCH 13/14] fix: avoid random numeric fails --- test/SOS_correctness.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/SOS_correctness.jl b/test/SOS_correctness.jl index e31fb77..f76f360 100644 --- a/test/SOS_correctness.jl +++ b/test/SOS_correctness.jl @@ -33,7 +33,7 @@ @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) ######################################################### From 874881b1abab25d5fea9971e7e46b000ad15e6c9 Mon Sep 17 00:00:00 2001 From: kalmarek Date: Thu, 4 Jul 2019 17:46:35 +0200 Subject: [PATCH 14/14] remove Nemo, Manifest update, version 0.3 --- .travis.yml | 1 + Manifest.toml | 28 +++++++++++----------------- Project.toml | 5 ++--- 3 files changed, 14 insertions(+), 20 deletions(-) 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"