From 560476b4b375e00fd531846631d388dea07b314d Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Thu, 10 Nov 2022 19:38:20 +0100 Subject: [PATCH 01/17] remove old ci files --- .codecov.yml | 1 - .gitignore | 4 ---- .travis.yml | 28 ---------------------------- appveyor.yml | 34 ---------------------------------- 4 files changed, 67 deletions(-) delete mode 100644 .codecov.yml delete mode 100644 .travis.yml delete mode 100644 appveyor.yml diff --git a/.codecov.yml b/.codecov.yml deleted file mode 100644 index 69cb760..0000000 --- a/.codecov.yml +++ /dev/null @@ -1 +0,0 @@ -comment: false diff --git a/.gitignore b/.gitignore index d87b441..3f02ca7 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,3 @@ *.jl.*.cov *.jl.mem Manifest.toml -test/SL* -test/oSL* -test/SAut* -test/oSAut* diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index a686713..0000000 --- a/.travis.yml +++ /dev/null @@ -1,28 +0,0 @@ -# Documentation: http://docs.travis-ci.com/user/languages/julia/ -language: julia -os: - - linux - - osx -julia: - - 1.1 - - 1.2 - - 1.3 - - nightly -notifications: - email: true -matrix: - fast_finish: true - allow_failures: - - julia: nightly - - os: osx - -addons: - apt: - packages: - - hdf5-tools - -## uncomment the following lines to override the default test -# script: - # - julia -e 'using Pkg; Pkg.build(); Pkg.test(coverage=true);' - -codecov: true diff --git a/appveyor.yml b/appveyor.yml deleted file mode 100644 index 475c68c..0000000 --- a/appveyor.yml +++ /dev/null @@ -1,34 +0,0 @@ -environment: - matrix: - - JULIAVERSION: "julialang/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe" - - JULIAVERSION: "julialang/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe" - - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe" - - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe" - -branches: - only: - - master - - /release-.*/ - -notifications: - - provider: Email - on_build_success: false - on_build_failure: false - on_build_status_changed: false - -install: -# Download most recent Julia Windows binary - - ps: (new-object net.webclient).DownloadFile( - $("http://s3.amazonaws.com/"+$env:JULIAVERSION), - "C:\projects\julia-binary.exe") -# Run installer silently, output to C:\projects\julia - - C:\projects\julia-binary.exe /S /D=C:\projects\julia - -build_script: -# Need to convert from shallow to complete for Pkg.clone to work - - IF EXIST .git\shallow (git fetch --unshallow) - - C:\projects\julia\bin\julia -e "versioninfo(); - Pkg.clone(pwd(), \"Property(T)\"); Pkg.build(\"Property(T)\")" - -test_script: - - C:\projects\julia\bin\julia -e "Pkg.test(\"Property(T)\")" From aa9ccb882a220c7359860363998d86d818110b95 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Thu, 10 Nov 2022 19:39:06 +0100 Subject: [PATCH 02/17] tiny changes to the README --- README.md | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 58b8a55..4e7ea5b 100644 --- a/README.md +++ b/README.md @@ -83,7 +83,10 @@ julia> include("test/optimizers.jl"); Now we have everything what we need to solve the problem! ```julia -julia> status, warmstart = PropertyT.solve(opt_problem, scs_optimizer(max_iters=5_000, accel=50, alpha=1.9)); +julia> status, warmstart = PropertyT.solve( + opt_problem, + scs_optimizer(max_iters=5_000, accel=50, alpha=1.9), +); ------------------------------------------------------------------ SCS v3.2.1 - Splitting Conic Solver (c) Brendan O'Donoghue, Stanford University, 2012 @@ -119,8 +122,12 @@ julia> status ALMOST_OPTIMAL::TerminationStatusCode = 7 ``` The solver didn't manage to solve the problem but it got quite close! (duality gap is ~`1.63e-6`). Let's try once again this time warmstarting the solver: -``` -julia> status, warmstart = PropertyT.solve(opt_problem, scs_optimizer(max_iters=10_000, accel=50, alpha=1.9), warmstart); +```julia +julia> status, warmstart = PropertyT.solve( + opt_problem, + scs_optimizer(max_iters=10_000, accel=50, alpha=1.9), + warmstart, +); ------------------------------------------------------------------ SCS v3.2.1 - Splitting Conic Solver (c) Brendan O'Donoghue, Stanford University, 2012 From 65286e09d26f3db440f8201ab2f585366bb5b663 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Sun, 13 Nov 2022 14:27:58 +0100 Subject: [PATCH 03/17] fix _positive_direction so that Adj works for N=2 --- src/roots.jl | 6 ++---- test/graded_adj.jl | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+), 4 deletions(-) diff --git a/src/roots.jl b/src/roots.jl index 94beb1d..3c5d6f3 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -48,10 +48,8 @@ function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M} end function _positive_direction(α::Root{N}) where {N} - last = -1 / √2^(N - 1) - return Root{N,Float64}( - SVector(ntuple(i -> ifelse(i == N, last, (√2)^-i), N)), - ) + v = α.coord + 1 / (N * 100) * rand(N) + return Root{N,Float64}(v / norm(v, 2)) end function positive(roots::AbstractVector{<:Root{N}}) where {N} diff --git a/test/graded_adj.jl b/test/graded_adj.jl index 7c65e30..2c20506 100644 --- a/test/graded_adj.jl +++ b/test/graded_adj.jl @@ -46,6 +46,28 @@ @testset "Symplectic group" begin + @testset "Sp2(ℤ)" begin + genus = 2 + halfradius = 1 + + SpN = MatrixGroups.SymplecticGroup{2genus}(Int8) + + RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true) + + Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity + Δ = RG(length(S)) - sum(RG(s) for s in S) + Δs = PropertyT.laplacians( + RG, + S, + x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])), + ) + Δ, Δs + end + + sq = sum(Δᵢ^2 for Δᵢ in values(Δs)) + @test PropertyT.Adj(Δs, :C₂) + sq == Δ^2 + end + genus = 3 halfradius = 1 From 227e82d5516c58a70f09f8442d7ee87bb1815bb1 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:45:41 +0100 Subject: [PATCH 04/17] make augmented compute_sos fast --- src/certify.jl | 90 ++++++++++++++++------------------------------ test/1712.07167.jl | 37 ++++++------------- 2 files changed, 40 insertions(+), 87 deletions(-) diff --git a/src/certify.jl b/src/certify.jl index 089621d..cdebcc9 100644 --- a/src/certify.jl +++ b/src/certify.jl @@ -5,80 +5,50 @@ function augment_columns!(Q::AbstractMatrix) return Q end -function _fma_SOS_thr!( - result::AbstractVector{T}, - mstructure::AbstractMatrix{<:Integer}, - Q::AbstractMatrix{T}, - acc_matrix=zeros(T, size(mstructure)...), -) where {T} +function __sos_via_sqr!( + res::StarAlgebras.AlgebraElement, + P::AbstractMatrix; + augmented::Bool +) + StarAlgebras.zero!(res) + A = parent(res) + b = basis(A) + @assert size(A.mstructure) == size(P) + e = b[one(b[1])] - s1, s2 = size(mstructure) - - @inbounds for k = 1:s2 - let k = k, s1 = s1, s2 = s2, Q = Q, acc_matrix = acc_matrix - Threads.@threads for j = 1:s2 - for i = 1:s1 - @inbounds acc_matrix[i, j] = - muladd(Q[i, k], Q[j, k], acc_matrix[i, j]) - end + for i in axes(A.mstructure, 1) + x = StarAlgebras._istwisted(A.mstructure) ? StarAlgebras.star(b[i]) : b[i] + for j in axes(A.mstructure, 2) + p = P[i, j] + xy = b[A.mstructure[i, j]] + # either result += P[x,y]*(x*y) + res[xy] += p + if augmented + # or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy) + y = b[j] + res[e] += p + res[x] -= p + res[y] -= p end end end - @inbounds for j = 1:s2 - for i = 1:s1 - result[mstructure[i, j]] += acc_matrix[i, j] - end - end - - return result + return res end -function _cnstr_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix, cnstrs) +function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix, cnstrs) StarAlgebras.zero!(res) - Q² = Q' * Q for (g, A_g) in cnstrs res[g] = dot(A_g, Q²) end return res end -function _augmented_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix) - A = parent(res) - StarAlgebras.zero!(res) - Q² = Q' * Q - - N = LinearAlgebra.checksquare(A.mstructure) - augmented_basis = [A(1) - A(b) for b in @view basis(A)[1:N]] - tmp = zero(res) - - for (j, y) in enumerate(augmented_basis) - for (i, x) in enumerate(augmented_basis) - # res += Q²[i, j] * x * y - - StarAlgebras.mul!(tmp, x, y) - StarAlgebras.mul!(tmp, tmp, Q²[i, j]) - StarAlgebras.add!(res, res, tmp) - end - end - return res -end - function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool) - if augmented - z = zeros(eltype(Q), length(basis(A))) - res = StarAlgebras.AlgebraElement(z, A) - return _augmented_sos!(res, Q) - cnstrs = constraints(basis(A), A.mstructure; augmented=true) - return _cnstr_sos!(res, Q, cnstrs) - else - @assert size(A.mstructure) == size(Q) - z = zeros(eltype(Q), length(basis(A))) - - _fma_SOS_thr!(z, A.mstructure, Q) - - return StarAlgebras.AlgebraElement(z, A) - end + Q² = Q' * Q + res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A) + res = __sos_via_sqr!(res, Q², augmented=augmented) + return res end function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius) @@ -159,7 +129,7 @@ function certify_solution( true, compute_sos(parent(elt), Q_int, augmented=augmented) end - @info "Checking in $(eltype(sos_int)) arithmetic with" λ + @info "Checking in $(eltype(sos_int)) arithmetic with" λ_int λ_certified = sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius) diff --git a/test/1712.07167.jl b/test/1712.07167.jl index 982806e..9a1adcc 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -6,37 +6,20 @@ function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimize @time status, _ = PropertyT.solve(sos_problem, optimizer) Q = let Ps = Ps - flPs = [real.(sqrt(JuMP.value.(P))) for P in Ps] - PropertyT.reconstruct(flPs, wd) + Qs = [real.(sqrt(JuMP.value.(P))) for P in Ps] + PropertyT.reconstruct(Qs, wd) end λ = JuMP.value(sos_problem[:λ]) - sos = let RG = parent(elt), Q = Q - z = zeros(eltype(Q), length(basis(RG))) - res = AlgebraElement(z, RG) - cnstrs = PropertyT.constraints(basis(RG), RG.mstructure, augmented=true) - PropertyT._cnstr_sos!(res, Q, cnstrs) - end - - residual = elt - λ * unit - sos - λ_fl = PropertyT.sufficient_λ(residual, λ, halfradius=2) - - λ_fl < 0 && return status, false, λ_fl - - sos = let RG = parent(elt), Q = [PropertyT.IntervalArithmetic.@interval(q) for q in Q] - z = zeros(eltype(Q), length(basis(RG))) - res = AlgebraElement(z, RG) - cnstrs = PropertyT.constraints(basis(RG), RG.mstructure, augmented=true) - PropertyT._cnstr_sos!(res, Q, cnstrs) - end - - λ_int = PropertyT.IntervalArithmetic.@interval(λ) - - residual_int = elt - λ_int * unit - sos - λ_int = PropertyT.sufficient_λ(residual_int, λ_int, halfradius=2) - - return status, λ_int > 0, PropertyT.IntervalArithmetic.inf(λ_int) + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + λ, + Q, + halfradius=2 + ) + return status, certified, λ_cert end @testset "1712.07167 Examples" begin From 2f89538eb092e8dc8ed40dfd9c2fd355df54d83e Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:47:38 +0100 Subject: [PATCH 05/17] add __fast_recursive_dot to speed up constraints --- src/sos_sdps.jl | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index d2387ed..8f086c7 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -160,6 +160,28 @@ sos_problem_primal( kwargs... ) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...) +function __fast_recursive_dot!( + res::JuMP.AffExpr, + Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}}, + Ms::AbstractArray{<:AbstractSparseMatrix}; +) + @assert length(Ps) == length(Ms) + + for (A, P) in zip(Ms, Ps) + iszero(Ms) && continue + rows = rowvals(A) + vals = nonzeros(A) + for cidx in axes(A, 2) + for i in nzrange(A, cidx) + ridx = rows[i] + val = vals[i] + JuMP.add_to_expression!(res, P[ridx, cidx], val) + end + end + end + return res +end + function sos_problem_primal( elt::StarAlgebras.AlgebraElement, orderunit::StarAlgebras.AlgebraElement, @@ -225,9 +247,9 @@ function sos_problem_primal( M_dot_P = sum(dot(M[π], P[π]) for π in eachindex(M) if !iszero(M[π])) if feasibility_problem - JuMP.@constraint(model, x == M_dot_P) + JuMP.@constraint(model, x == __fast_recursive_dot!(JuMP.AffExpr(), P, M)) else - JuMP.@constraint(model, x - λ * u == M_dot_P) + JuMP.@constraint(model, x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, M)) end end return model, P From 971e07b81964902f75c73cff5d5d26e382557118 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:50:09 +0100 Subject: [PATCH 06/17] move to using sparse matrices in symmetrized sdp dense are faster for small sizes only --- src/sos_sdps.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index 8f086c7..bbcba67 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -216,15 +216,14 @@ function sos_problem_primal( P = map(direct_summands(wedderburn)) do ds dim = size(ds, 1) P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric) - @constraint(model, P in PSDCone()) + JuMP.@constraint(model, P in PSDCone()) P end begin # preallocating T = eltype(wedderburn) - M = zeros.(T, size.(P)) + M = spzeros.(T, size.(P)) M_orb = zeros(T, size(parent(elt).mstructure)...) - tmps = SymbolicWedderburn._tmps(wedderburn) end X = convert(Vector{T}, StarAlgebras.coeffs(elt)) @@ -235,16 +234,27 @@ function sos_problem_primal( @info "Adding $(length(invariant_vectors(wedderburn))) constraints" - for iv in invariant_vectors(wedderburn) + ds = SymbolicWedderburn.direct_summands(wedderburn) + Uπs = SymbolicWedderburn.image_basis.(ds) + T = eltype(first(Uπs)) + degrees = SymbolicWedderburn.degree.(ds) + + for (i, iv) in enumerate(invariant_vectors(wedderburn)) + # @debug i + i % 100 == 0 && print('.') + i % 10000 === 0 && print('\n') + x = dot(X, iv) u = dot(U, iv) M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv) + M = SymbolicWedderburn.diagonalize!(M, M_orb, Uπs, degrees) + SparseArrays.droptol!.(M, 10 * eps(T) * max(size(M_orb)...)) - M = SymbolicWedderburn.diagonalize!(M, M_orb, wedderburn, tmps) - SymbolicWedderburn.zerotol!.(M, atol=1e-12) - - M_dot_P = sum(dot(M[π], P[π]) for π in eachindex(M) if !iszero(M[π])) + # @debug [nnz(m) / length(m) for m in M] + # spM = sparse.(M) + # @time M_dot_P = sum(dot(spM[π], P[π]) for π in eachindex(spM) if !iszero(spM[π])) + # @info density = [count(!iszero, m) / sum(length, m) for m in M] if feasibility_problem JuMP.@constraint(model, x == __fast_recursive_dot!(JuMP.AffExpr(), P, M)) From 703e69fc6274059511bc6e8b2a828e5fbcc3f266 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:50:27 +0100 Subject: [PATCH 07/17] fix actions on SpNs --- src/actions/actions.jl | 1 + src/actions/spn_conjugation.jl | 37 +++++++++++++++++++++++++--------- 2 files changed, 28 insertions(+), 10 deletions(-) diff --git a/src/actions/actions.jl b/src/actions/actions.jl index aa02b15..9cc6872 100644 --- a/src/actions/actions.jl +++ b/src/actions/actions.jl @@ -4,6 +4,7 @@ StarAlgebras.star(g::Groups.GroupElement) = inv(g) include("alphabet_permutation.jl") include("sln_conjugation.jl") +include("spn_conjugation.jl") include("autfn_conjugation.jl") function SymbolicWedderburn.action( diff --git a/src/actions/spn_conjugation.jl b/src/actions/spn_conjugation.jl index 25b8bbc..ec760e9 100644 --- a/src/actions/spn_conjugation.jl +++ b/src/actions/spn_conjugation.jl @@ -1,26 +1,43 @@ ## Particular definitions for actions on Sp(n,ℤ) function _conj( - t::MatrixGroups.ElementarySymplectic{N,T}, + s::MatrixGroups.ElementarySymplectic{N,T}, σ::PermutationGroups.AbstractPerm, ) where {N,T} @assert iseven(N) - @assert degree(σ) == N ÷ 2 "Got degree = $(degree(σ)); N = $N" - i = mod1(t.i, N ÷ 2) - ib = i == t.i ? 0 : N ÷ 2 - j = mod1(t.j, N ÷ 2) - jb = j == t.j ? 0 : N ÷ 2 - return MatrixGroups.ElementarySymplectic{N}(t.symbol, i^inv(σ) + ib, j^inv(σ) + jb, t.val) + @assert PermutationGroups.degree(σ) == N ÷ 2 "Got degree = $(PermutationGroups.degree(σ)); N = $N" + n = N ÷ 2 + @assert 1 ≤ s.i ≤ N + @assert 1 ≤ s.j ≤ N + if s.symbol == :A + @assert 1 ≤ s.i ≤ n + @assert 1 ≤ s.j ≤ n + i = s.i^inv(σ) + j = s.j^inv(σ) + else + @assert s.symbol == :B + @assert xor(s.i > n, s.j > n) + if s.i > n + i = (s.i - n)^inv(σ) + n + j = s.j^inv(σ) + elseif s.j > n + i = s.i^inv(σ) + j = (s.j - n)^inv(σ) + n + end + end + return MatrixGroups.ElementarySymplectic{N}(s.symbol, i, j, s.val) end function _conj( - t::MatrixGroups.ElementarySymplectic{N,T}, + s::MatrixGroups.ElementarySymplectic{N,T}, x::Groups.Constructions.DirectPowerElement, ) where {N,T} @assert Groups.order(Int, parent(x).group) == 2 @assert iseven(N) - just_one_flips = xor(isone(x.elts[mod1(t.i, N ÷ 2)]), isone(x.elts[mod1(t.j, N ÷ 2)])) - return ifelse(just_one_flips, inv(t), t) + n = N ÷ 2 + i, j = ifelse(s.i <= n, s.i, s.i - n), ifelse(s.j <= n, s.j, s.j - n) + just_one_flips = xor(isone(x.elts[i]), isone(x.elts[j])) + return ifelse(just_one_flips, inv(s), s) end action_by_conjugation(sln::Groups.MatrixGroups.SymplecticGroup, Σ::Groups.Group) = From e6bd862e7a46deaef872f61ceb23da483ba65685 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:50:50 +0100 Subject: [PATCH 08/17] small tweaks --- src/gradings.jl | 3 ++- test/1703.09680.jl | 6 +++--- test/graded_adj.jl | 2 +- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/gradings.jl b/src/gradings.jl index 1f20d6b..8f34efa 100644 --- a/src/gradings.jl +++ b/src/gradings.jl @@ -6,7 +6,8 @@ Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} = function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N} if s.symbol === :A return Roots.𝕖(N ÷ 2, s.i) - Roots.𝕖(N ÷ 2, s.j) - else#if s.symbol === :B + else + @assert s.symbol === :B n = N ÷ 2 i, j = ifelse(s.i <= n, s.i, s.i - n), ifelse(s.j <= n, s.j, s.j - n) return (-1)^(s.i > s.j) * (Roots.𝕖(n, i) + Roots.𝕖(n, j)) diff --git a/test/1703.09680.jl b/test/1703.09680.jl index 074d263..edf0306 100644 --- a/test/1703.09680.jl +++ b/test/1703.09680.jl @@ -79,15 +79,15 @@ end @test λ > 1 m = PropertyT.sos_problem_dual(elt, unit) - PropertyT.solve(m, scs_optimizer( - eps=1e-10, + PropertyT.solve(m, cosmo_optimizer( + eps=1e-6, max_iters=5_000, accel=50, alpha=1.9, )) @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL) - @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-3 + @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2 end @testset "SAut(F₂)" begin diff --git a/test/graded_adj.jl b/test/graded_adj.jl index 2c20506..bf5a804 100644 --- a/test/graded_adj.jl +++ b/test/graded_adj.jl @@ -85,7 +85,7 @@ Δ, Δs end - @testset "Adj correctness: genus=$genus" begin + @testset "Adj numerics for genus=$genus" begin all_subtypes = ( :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂ From fca7733280f356493926ed30fda30d5d88313b3b Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 14 Nov 2022 19:51:38 +0100 Subject: [PATCH 09/17] remove old (unported) SOS_correctness tests --- test/SOS_correctness.jl | 59 ----------------------------------------- 1 file changed, 59 deletions(-) delete mode 100644 test/SOS_correctness.jl diff --git a/test/SOS_correctness.jl b/test/SOS_correctness.jl deleted file mode 100644 index c13abe2..0000000 --- a/test/SOS_correctness.jl +++ /dev/null @@ -1,59 +0,0 @@ -@testset "Correctness of HPC SOS computation" begin - - function prepare(G_name, λ, S_size) - pm = load("$G_name/delta.jld", "pm") - P = load("$G_name/$λ/solution.jld", "P") - @time Q = real(sqrt(P)) - - Δ_coeff = SparseVector(maximum(pm), collect(1:1+S_size), [S_size; ((-1.0) for i in 1:S_size)...]) - - Δ²_coeff = GroupRings.GRmul!(spzeros(length(Δ_coeff)), Δ_coeff, Δ_coeff, pm) - - eoi = Δ²_coeff - λ*Δ_coeff - - Q = PropertyT.augIdproj(Q) - - return eoi, pm, Q - end - - ######################################################### - NAME = "SL(3,Z)" - eoi, pm, Q = prepare(NAME, 0.1, 3*2*2) - - @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 - @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) - - ######################################################### - NAME = "SL(3,Z)_orbit" - eoi, pm, Q = prepare(NAME, 0.27, 3*2*2) - - @time sos_sqr = PropertyT.compute_SOS_square(pm, Q) - @time sos_hpc = PropertyT.compute_SOS(pm, Q) - - @test norm(sos_sqr - sos_hpc, 1) < 5e-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) - - ######################################################### - NAME = "SL(4,Z)_orbit" - eoi, pm, Q = prepare(NAME, 1.3, 4*3*2) - - @time sos_sqr = PropertyT.compute_SOS_square(pm, Q) - @time sos_hpc = PropertyT.compute_SOS(pm, Q) - - @test norm(sos_sqr - sos_hpc, 1) < 2e-10 - @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) - - ######################################################### - NAME = "SAut(F3)_orbit" - eoi, pm, Q = prepare(NAME, 0.15, 4*3*2*2) - - @time sos_sqr = PropertyT.compute_SOS_square(pm, Q) - @time sos_hpc = PropertyT.compute_SOS(pm, Q) - - @test norm(sos_sqr - sos_hpc, 1) < 6e-11 - @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) - -end From ca6a17accad860b7dfed554b86af88349e2a4fbd Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 15 Nov 2022 00:27:09 +0100 Subject: [PATCH 10/17] use ProgressMeter instead of poormans progressbar --- Project.toml | 2 ++ src/sos_sdps.jl | 46 +++++++++++++++++++++++++++------------------- 2 files changed, 29 insertions(+), 19 deletions(-) diff --git a/Project.toml b/Project.toml index ade1715..cad256d 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2" @@ -17,6 +18,7 @@ COSMO = "0.8" Groups = "0.7" IntervalArithmetic = "0.20" JuMP = "1.3" +ProgressMeter = "1.7" SCS = "1.1.0" StaticArrays = "1" SymbolicWedderburn = "0.3.1" diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index bbcba67..96c71c3 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -182,13 +182,17 @@ function __fast_recursive_dot!( return res end +import ProgressMeter +__show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")] + function sos_problem_primal( elt::StarAlgebras.AlgebraElement, orderunit::StarAlgebras.AlgebraElement, wedderburn::WedderburnDecomposition; upper_bound=Inf, augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)), - check_orthogonality=true + check_orthogonality=true, + show_progress=false ) @assert parent(elt) === parent(orderunit) @@ -222,7 +226,7 @@ function sos_problem_primal( begin # preallocating T = eltype(wedderburn) - M = spzeros.(T, size.(P)) + Ms = spzeros.(T, size.(P)) M_orb = zeros(T, size(parent(elt).mstructure)...) end @@ -232,36 +236,40 @@ function sos_problem_primal( # defining constraints based on the multiplicative structure cnstrs = constraints(parent(elt), augmented=augmented, twisted=true) - @info "Adding $(length(invariant_vectors(wedderburn))) constraints" - - ds = SymbolicWedderburn.direct_summands(wedderburn) - Uπs = SymbolicWedderburn.image_basis.(ds) - T = eltype(first(Uπs)) - degrees = SymbolicWedderburn.degree.(ds) + prog = ProgressMeter.Progress( + length(invariant_vectors(wedderburn)), + dt=1, + desc="Adding constraints... ", + enabled=show_progress, + barlen=60, + showspeed=true + ) for (i, iv) in enumerate(invariant_vectors(wedderburn)) - # @debug i - i % 100 == 0 && print('.') - i % 10000 === 0 && print('\n') + ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n)) x = dot(X, iv) u = dot(U, iv) M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv) - M = SymbolicWedderburn.diagonalize!(M, M_orb, Uπs, degrees) - SparseArrays.droptol!.(M, 10 * eps(T) * max(size(M_orb)...)) + Ms = SymbolicWedderburn.diagonalize!(Ms, M_orb, wedderburn) + SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...)) - # @debug [nnz(m) / length(m) for m in M] - # spM = sparse.(M) - # @time M_dot_P = sum(dot(spM[π], P[π]) for π in eachindex(spM) if !iszero(spM[π])) - # @info density = [count(!iszero, m) / sum(length, m) for m in M] + # @info [nnz(m) / length(m) for m in Ms] if feasibility_problem - JuMP.@constraint(model, x == __fast_recursive_dot!(JuMP.AffExpr(), P, M)) + JuMP.@constraint( + model, + x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms) + ) else - JuMP.@constraint(model, x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, M)) + JuMP.@constraint( + model, + x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms) + ) end end + ProgressMeter.finish!(prog) return model, P end From 35c5110a3749289957ca3d86f760f2dc0afce3cd Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 15 Nov 2022 18:51:43 +0100 Subject: [PATCH 11/17] move solve to a separate file --- src/PropertyT.jl | 1 + src/solve.jl | 63 +++++++++++++++++++++++++++++++++++++++++++++++ src/sos_sdps.jl | 64 ------------------------------------------------ 3 files changed, 64 insertions(+), 64 deletions(-) create mode 100644 src/solve.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index c9adba4..50495a5 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -15,6 +15,7 @@ import SymbolicWedderburn.PermutationGroups include("constraint_matrix.jl") include("sos_sdps.jl") +include("solve.jl") include("certify.jl") include("sqadjop.jl") diff --git a/src/solve.jl b/src/solve.jl new file mode 100644 index 0000000..9d76194 --- /dev/null +++ b/src/solve.jl @@ -0,0 +1,63 @@ +## Low-level solve + +setwarmstart!(model::JuMP.Model, ::Nothing) = model + +function setwarmstart!(model::JuMP.Model, warmstart) + constraint_map = Dict( + ct => JuMP.all_constraints(model, ct...) for + ct in JuMP.list_of_constraint_types(model) + ) + + JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal) + + for (ct, idx) in pairs(constraint_map) + JuMP.set_start_value.(idx, warmstart.slack[ct]) + JuMP.set_dual_start_value.(idx, warmstart.dual[ct]) + end + return model +end + +function getwarmstart(model::JuMP.Model) + constraint_map = Dict( + ct => JuMP.all_constraints(model, ct...) for + ct in JuMP.list_of_constraint_types(model) + ) + + primal = value.(JuMP.all_variables(model)) + + slack = Dict(k => value.(v) for (k, v) in constraint_map) + duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map) + + return (primal=primal, dual=duals, slack=slack) +end + +function solve(m::JuMP.Model, optimizer, warmstart=nothing) + + JuMP.set_optimizer(m, optimizer) + MOIU.attach_optimizer(m) + + m = setwarmstart!(m, warmstart) + + JuMP.optimize!(m) + Base.Libc.flush_cstdio() + + status = JuMP.termination_status(m) + + return status, getwarmstart(m) +end + +function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing) + + isdir(dirname(solverlog)) || mkpath(dirname(solverlog)) + + Base.flush(Base.stdout) + Base.Libc.flush_cstdio() + status, warmstart = open(solverlog, "a+") do logfile + redirect_stdout(logfile) do + status, warmstart = solve(m, optimizer, warmstart) + status, warmstart + end + end + + return status, warmstart +end diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index 96c71c3..795ffac 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -326,67 +326,3 @@ function reconstruct!( return res end -## -# Low-level solve - -setwarmstart!(model::JuMP.Model, ::Nothing) = model - -function setwarmstart!(model::JuMP.Model, warmstart) - constraint_map = Dict( - ct => JuMP.all_constraints(model, ct...) for - ct in JuMP.list_of_constraint_types(model) - ) - - JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal) - - for (ct, idx) in pairs(constraint_map) - JuMP.set_start_value.(idx, warmstart.slack[ct]) - JuMP.set_dual_start_value.(idx, warmstart.dual[ct]) - end - return model -end - -function getwarmstart(model::JuMP.Model) - constraint_map = Dict( - ct => JuMP.all_constraints(model, ct...) for - ct in JuMP.list_of_constraint_types(model) - ) - - primal = value.(JuMP.all_variables(model)) - - slack = Dict(k => value.(v) for (k, v) in constraint_map) - duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map) - - return (primal=primal, dual=duals, slack=slack) -end - -function solve(m::JuMP.Model, optimizer, warmstart=nothing) - - JuMP.set_optimizer(m, optimizer) - MOIU.attach_optimizer(m) - - m = setwarmstart!(m, warmstart) - - JuMP.optimize!(m) - Base.Libc.flush_cstdio() - - status = JuMP.termination_status(m) - - return status, getwarmstart(m) -end - -function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing) - - isdir(dirname(solverlog)) || mkpath(dirname(solverlog)) - - Base.flush(Base.stdout) - Base.Libc.flush_cstdio() - status, warmstart = open(solverlog, "a+") do logfile - redirect_stdout(logfile) do - status, warmstart = solve(m, optimizer, warmstart) - status, warmstart - end - end - - return status, warmstart -end From f053bffefe1f3ef857358cccae6e891f13d1fda3 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 15 Nov 2022 18:53:03 +0100 Subject: [PATCH 12/17] rewrite reconstruct! with a better architecture --- src/PropertyT.jl | 2 +- src/reconstruct.jl | 69 ++++++++++++++++++++++++++++++++++++++++++++++ src/sos_sdps.jl | 54 ------------------------------------ 3 files changed, 70 insertions(+), 55 deletions(-) create mode 100644 src/reconstruct.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 50495a5..ac32da9 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -1,4 +1,3 @@ -__precompile__() module PropertyT using LinearAlgebra @@ -16,6 +15,7 @@ import SymbolicWedderburn.PermutationGroups include("constraint_matrix.jl") include("sos_sdps.jl") include("solve.jl") +include("reconstruct.jl") include("certify.jl") include("sqadjop.jl") diff --git a/src/reconstruct.jl b/src/reconstruct.jl new file mode 100644 index 0000000..75451ee --- /dev/null +++ b/src/reconstruct.jl @@ -0,0 +1,69 @@ +__outer_dim(wd::WedderburnDecomposition) = size(first(direct_summands(wd)), 2) + +function __group_of(wd::WedderburnDecomposition) + # this is veeeery hacky... ;) + return parent(first(keys(wd.hom.cache))) +end + +function reconstruct( + Ms::AbstractVector{<:AbstractMatrix}, + wbdec::WedderburnDecomposition; + atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec) +) + n = __outer_dim(wbdec) + res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds) + res = similar(M, n, n) + reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol) + end + return res +end + +function reconstruct!( + res::AbstractMatrix, + M::AbstractMatrix, + ds::SymbolicWedderburn.DirectSummand, + G, + hom; + atol=eps(real(eltype(ds))) * 10max(size(res)...) +) + res .= zero(eltype(res)) + U = SymbolicWedderburn.image_basis(ds) + d = SymbolicWedderburn.degree(ds) + tmp = (U' * M * U) .* d + + res = average!(res, tmp, G, hom) + if eltype(res) <: AbstractFloat + __droptol!(res, atol) # TODO: is this really necessary?! + end + return res +end + +function __droptol!(M::AbstractMatrix, tol) + for i in eachindex(M) + if abs(M[i]) < tol + M[i] = zero(M[i]) + end + end + return M +end + +# implement average! for other actions when needed +function average!( + res::AbstractMatrix, + M::AbstractMatrix, + G::Groups.Group, + hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations} +) + @assert size(M) == size(res) + for g in G + gext = SymbolicWedderburn.induce(hom, g) + Threads.@threads for c in axes(res, 2) + for r in axes(res, 1) + res[r, c] += M[r^gext, c^gext] + end + end + end + o = Groups.order(Int, G) + res ./= o + return res +end diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index 795ffac..f6135db 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -272,57 +272,3 @@ function sos_problem_primal( ProgressMeter.finish!(prog) return model, P end - -function reconstruct(Ps, wd::WedderburnDecomposition) - N = size(first(direct_summands(wd)), 2) - P = zeros(eltype(wd), N, N) - return reconstruct!(P, Ps, wd) -end - -function group_of(wd::WedderburnDecomposition) - # this is veeeery hacky... ;) - return parent(first(keys(wd.hom.cache))) -end - -# TODO: move to SymbolicWedderburn -SymbolicWedderburn.action(wd::WedderburnDecomposition) = - SymbolicWedderburn.action(wd.hom) - -function reconstruct!( - res::AbstractMatrix, - Ps, - wedderburn::WedderburnDecomposition, -) - G = group_of(wedderburn) - - act = SymbolicWedderburn.action(wedderburn) - - @assert act isa SymbolicWedderburn.ByPermutations - - for (π, ds) in pairs(direct_summands(wedderburn)) - Uπ = SymbolicWedderburn.image_basis(ds) - - # LinearAlgebra.mul!(tmp, Uπ', P[π]) - # LinearAlgebra.mul!(tmp2, tmp, Uπ) - tmp2 = Uπ' * Ps[π] * Uπ - if eltype(res) <: AbstractFloat - SymbolicWedderburn.zerotol!(tmp2, atol=1e-12) - end - tmp2 .*= SymbolicWedderburn.degree(ds) - - @assert size(tmp2) == size(res) - - for g in G - p = SymbolicWedderburn.induce(wedderburn.hom, g) - for c in axes(res, 2) - for r in axes(res, 1) - res[r, c] += tmp2[r^p, c^p] - end - end - end - end - res ./= Groups.order(Int, G) - - return res -end - From b92f12107f9f66fca6104b9079defbe7db22382d Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 15 Nov 2022 18:53:19 +0100 Subject: [PATCH 13/17] update SymbolicWedderburn to 0.3.2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index cad256d..8255f05 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ JuMP = "1.3" ProgressMeter = "1.7" SCS = "1.1.0" StaticArrays = "1" -SymbolicWedderburn = "0.3.1" +SymbolicWedderburn = "0.3.2" julia = "1.6" [extras] From d1d46d13ef2fa9445d4a2a698072ab1d839f442d Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Wed, 16 Nov 2022 00:14:06 +0100 Subject: [PATCH 14/17] use spzeros constructor compatible with julia-1.6 --- src/sos_sdps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index f6135db..2f6faea 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -226,7 +226,7 @@ function sos_problem_primal( begin # preallocating T = eltype(wedderburn) - Ms = spzeros.(T, size.(P)) + Ms = [spzeros.(T, size(p)...) for p in P] M_orb = zeros(T, size(parent(elt).mstructure)...) end From 18f681b12b030cff9ad0e776d423240ae197e15d Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Wed, 16 Nov 2022 00:16:01 +0100 Subject: [PATCH 15/17] add really quick tests to be run first tests run really fast, but compilation time is killing it... --- test/1703.09680.jl | 17 --------- test/1712.07167.jl | 24 ------------ test/check_positivity.jl | 41 ++++++++++++++++++++ test/quick_tests.jl | 80 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 21 ++++++----- 5 files changed, 132 insertions(+), 51 deletions(-) create mode 100644 test/check_positivity.jl create mode 100644 test/quick_tests.jl diff --git a/test/1703.09680.jl b/test/1703.09680.jl index edf0306..6ed0f79 100644 --- a/test/1703.09680.jl +++ b/test/1703.09680.jl @@ -1,20 +1,3 @@ -function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer) - @time sos_problem = - PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound) - - status, _ = PropertyT.solve(sos_problem, optimizer) - P = JuMP.value.(sos_problem[:P]) - Q = real.(sqrt(P)) - certified, λ_cert = PropertyT.certify_solution( - elt, - unit, - JuMP.objective_value(sos_problem), - Q, - halfradius=halfradius, - ) - return status, certified, λ_cert -end - @testset "1703.09680 Examples" begin @testset "SL(2,Z)" begin diff --git a/test/1712.07167.jl b/test/1712.07167.jl index 9a1adcc..be479d9 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -1,27 +1,3 @@ -function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer) - @assert aug(elt) == aug(unit) == 0 - @time sos_problem, Ps = - PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound) - - @time status, _ = PropertyT.solve(sos_problem, optimizer) - - Q = let Ps = Ps - Qs = [real.(sqrt(JuMP.value.(P))) for P in Ps] - PropertyT.reconstruct(Qs, wd) - end - - λ = JuMP.value(sos_problem[:λ]) - - certified, λ_cert = PropertyT.certify_solution( - elt, - unit, - λ, - Q, - halfradius=2 - ) - return status, certified, λ_cert -end - @testset "1712.07167 Examples" begin @testset "SAut(F₃)" begin diff --git a/test/check_positivity.jl b/test/check_positivity.jl new file mode 100644 index 0000000..974077b --- /dev/null +++ b/test/check_positivity.jl @@ -0,0 +1,41 @@ +function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer) + @time sos_problem = + PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound) + + status, _ = PropertyT.solve(sos_problem, optimizer) + P = JuMP.value.(sos_problem[:P]) + Q = real.(sqrt(P)) + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + JuMP.objective_value(sos_problem), + Q, + halfradius=halfradius, + ) + return status, certified, λ_cert +end + +function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer) + @assert aug(elt) == aug(unit) == 0 + @time sos_problem, Ps = + PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound) + + @time status, _ = PropertyT.solve(sos_problem, optimizer) + + Q = let Ps = Ps + Qs = [real.(sqrt(JuMP.value.(P))) for P in Ps] + PropertyT.reconstruct(Qs, wd) + end + + λ = JuMP.value(sos_problem[:λ]) + + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + λ, + Q, + halfradius=halfradius + ) + return status, certified, λ_cert +end + diff --git a/test/quick_tests.jl b/test/quick_tests.jl new file mode 100644 index 0000000..9a38b5c --- /dev/null +++ b/test/quick_tests.jl @@ -0,0 +1,80 @@ +@testset "Quick tests" begin + + @testset "SL(2,F₇)" begin + N = 2 + p = 7 + halfradius = 3 + G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p}) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true) + + Δ = let RG = RG, S = S + RG(length(S)) - sum(RG(s) for s in S) + end + + elt = Δ^2 + unit = Δ + ub = 0.58578# Inf# 1.5 + + @testset "standard formulation" begin + status, certified, λ_cert = check_positivity( + elt, + unit, + upper_bound=ub, + halfradius=2, + optimizer=cosmo_optimizer( + eps=1e-7, + max_iters=5_000, + accel=50, + alpha=1.95, + ) + ) + + @test status == JuMP.OPTIMAL + @test certified + @test λ_cert > 5857 // 10000 + + m = PropertyT.sos_problem_dual(elt, unit) + PropertyT.solve(m, cosmo_optimizer( + eps=1e-7, + max_iters=10_000, + accel=50, + alpha=1.95, + )) + + @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL) + @test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2 + end + + @testset "Wedderburn decomposition" begin + P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) + Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + act = PropertyT.action_by_conjugation(G, Σ) + + wd = WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[halfradius]]), + ) + + status, certified, λ_cert = check_positivity( + elt, + unit, + wd, + upper_bound=ub, + halfradius=2, + optimizer=cosmo_optimizer( + eps=1e-7, + max_iters=10_000, + accel=50, + alpha=1.9, + ), + ) + + @test status == JuMP.OPTIMAL + @test certified + @test λ_cert > 5857 // 10000 + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index cdeac69..7003c1e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,6 @@ using Test using LinearAlgebra using SparseArrays -BLAS.set_num_threads(1) -ENV["OMP_NUM_THREADS"] = 4 using Groups using Groups.GroupsCore @@ -14,15 +12,18 @@ using SymbolicWedderburn.StarAlgebras using SymbolicWedderburn.PermutationGroups include("optimizers.jl") +include("check_positivity.jl") +include("quick_tests.jl") -@testset "PropertyT" begin +if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI") + @testset "PropertyT" begin + include("constratint_matrices.jl") + include("actions.jl") - include("constratint_matrices.jl") - include("actions.jl") + include("1703.09680.jl") + include("1712.07167.jl") + include("1812.03456.jl") - include("1703.09680.jl") - include("1712.07167.jl") - include("1812.03456.jl") - - include("graded_adj.jl") + include("graded_adj.jl") + end end From 5349c21034a401da6548f38302cbac8e042949be Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Thu, 17 Nov 2022 02:50:48 +0100 Subject: [PATCH 16/17] add script for SpN_Adj.jl --- scripts/SpN_Adj.jl | 80 ++++++++++++++++++++++++++++++++++++++++++ scripts/argparse.jl | 19 ++++++++++ scripts/utils.jl | 85 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 184 insertions(+) create mode 100644 scripts/SpN_Adj.jl create mode 100644 scripts/argparse.jl create mode 100644 scripts/utils.jl diff --git a/scripts/SpN_Adj.jl b/scripts/SpN_Adj.jl new file mode 100644 index 0000000..80b9b65 --- /dev/null +++ b/scripts/SpN_Adj.jl @@ -0,0 +1,80 @@ +using LinearAlgebra +BLAS.set_num_threads(8) + +ENV["OMP_NUM_THREADS"] = 1 + +using Groups +import Groups.MatrixGroups + +include(joinpath(@__DIR__, "../test/optimizers.jl")) +using PropertyT + +using PropertyT.SymbolicWedderburn +using PropertyT.PermutationGroups +using PropertyT.StarAlgebras + +include(joinpath(@__DIR__, "argparse.jl")) +include(joinpath(@__DIR__, "utils.jl")) + +const N = parsed_args["N"] +const HALFRADIUS = parsed_args["halfradius"] +const UPPER_BOUND = parsed_args["upper_bound"] + +const GENUS = 2N + +G = MatrixGroups.SymplecticGroup{GENUS}(Int8) + +RG, S, sizes = + @time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true) + +wd = let RG = RG, N = N + G = StarAlgebras.object(RG) + P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) + Σ = Groups.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + # Σ = P + act = PropertyT.action_by_conjugation(G, Σ) + @info "Computing WedderburnDecomposition" + + wdfl = @time SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]), + ) +end + +Δ = RG(length(S)) - sum(RG(s) for s in S) +Δs = PropertyT.laplacians( + RG, + S, + x -> (gx = PropertyT.grading(x); Set([gx, -gx])), +) + +# elt = Δ^2 +elt = PropertyT.Adj(Δs, :C₂) +unit = Δ + +@time model, varP = PropertyT.sos_problem_primal( + elt, + unit, + wd, + upper_bound=UPPER_BOUND, + augmented=true, + show_progress=true +) + +solve_in_loop( + model, + wd, + varP, + logdir="./log/r=$HALFRADIUS/Sp($N,Z)/Adj_C₂-InfΔ", + optimizer=cosmo_optimizer( + eps=1e-10, + max_iters=20_000, + accel=50, + alpha=1.95, + ), + data=(elt=elt, unit=unit, halfradius=HALFRADIUS) +) + diff --git a/scripts/argparse.jl b/scripts/argparse.jl new file mode 100644 index 0000000..6bc414b --- /dev/null +++ b/scripts/argparse.jl @@ -0,0 +1,19 @@ +using ArgParse + +args_settings = ArgParseSettings() +@add_arg_table! args_settings begin + "-N" + help = "the degree/genus/etc. parameter for a group" + arg_type = Int + default = 3 + "--halfradius", "-R" + help = "the halfradius on which perform the sum of squares decomposition" + arg_type = Int + default = 2 + "--upper_bound", "-u" + help = "set upper bound for the optimization problem to speed-up the convergence" + arg_type = Float64 + default = Inf +end + +parsed_args = parse_args(ARGS, args_settings) diff --git a/scripts/utils.jl b/scripts/utils.jl new file mode 100644 index 0000000..df891b3 --- /dev/null +++ b/scripts/utils.jl @@ -0,0 +1,85 @@ +using Dates +using Serialization +using Logging + +import JuMP + +function get_solution(model) + λ = JuMP.value(model[:λ]) + Q = real.(sqrt(JuMP.value.(model[:P]))) + solution = Dict(:λ => λ, :Q => Q) + return solution +end + +function get_solution(model, wd, varP; logdir) + λ = JuMP.value(model[:λ]) + + Qs = [real.(sqrt(JuMP.value.(P))) for P in varP] + Q = PropertyT.reconstruct(Qs, wd) + solution = Dict(:λ => λ, :Q => Q) + return solution +end + +function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data) + @info "logging to $logdir" + status = JuMP.UNKNOWN_RESULT_STATUS + warm = try + solution = deserialize(joinpath(logdir, "solution.sjl")) + warm = solution[:warm] + @info "trying to warm-start model with λ=$(solution[:λ])..." + warm + catch + nothing + end + old_lambda = 0.0 + while status != JuMP.OPTIMAL + date = now() + log_file = joinpath(logdir, "solver_$date.log") + @info "Current logfile is $log_file." + isdir(dirname(log_file)) || mkpath(dirname(log_file)) + + λ, flag, certified_λ = let + # logstream = current_logger().logger.stream + # v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint + # @warn v + status, warm = @time PropertyT.solve(log_file, model, optimizer, warm) + + solution = get_solution(model, args...; logdir=logdir) + solution[:warm] = warm + + serialize(joinpath(logdir, "solution_$date.sjl"), solution) + serialize(joinpath(logdir, "solution.sjl"), solution) + + flag, λ_cert = open(log_file, append=true) do io + with_logger(SimpleLogger(io)) do + PropertyT.certify_solution( + data.elt, + data.unit, + solution[:λ], + solution[:Q], + halfradius=data.halfradius, + ) + end + end + + solution[:λ], flag, λ_cert + end + + if flag == true && certified_λ ≥ 0 + @info "Certification done with λ = $certified_λ" + return certified_λ + else + rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda)) + @info "Certification failed with λ = $λ" certified_λ rel_change + end + + old_lambda = certified_λ + + if rel_change < 1e-9 + @info "No progress detected, breaking" + break + end + end + + return status == JuMP.OPTIMAL ? old_lambda : NaN +end From 2cc944466780357baa7b66fd4782f704852cde87 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Thu, 17 Nov 2022 11:22:17 +0100 Subject: [PATCH 17/17] small fixes --- scripts/SpN_Adj.jl | 2 +- src/PropertyT.jl | 2 +- src/solve.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/SpN_Adj.jl b/scripts/SpN_Adj.jl index 80b9b65..c5922de 100644 --- a/scripts/SpN_Adj.jl +++ b/scripts/SpN_Adj.jl @@ -68,7 +68,7 @@ solve_in_loop( model, wd, varP, - logdir="./log/r=$HALFRADIUS/Sp($N,Z)/Adj_C₂-InfΔ", + logdir="./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ", optimizer=cosmo_optimizer( eps=1e-10, max_iters=20_000, diff --git a/src/PropertyT.jl b/src/PropertyT.jl index ac32da9..6b97d5d 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -29,7 +29,7 @@ include("actions/actions.jl") function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool) S = union!(S, inv.(S)) @info "generating wl-metric ball of radius $(2halfradius)" - @time E, sizes = Groups.wlmetric_ball_serial(S, radius=2halfradius) + @time E, sizes = Groups.wlmetric_ball(S, radius=2halfradius) @info "sizes = $(sizes)" @info "computing the *-algebra structure for G" @time RG = StarAlgebras.StarAlgebra{twisted}( diff --git a/src/solve.jl b/src/solve.jl index 9d76194..e7df555 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -39,7 +39,6 @@ function solve(m::JuMP.Model, optimizer, warmstart=nothing) m = setwarmstart!(m, warmstart) JuMP.optimize!(m) - Base.Libc.flush_cstdio() status = JuMP.termination_status(m) @@ -55,6 +54,7 @@ function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing) status, warmstart = open(solverlog, "a+") do logfile redirect_stdout(logfile) do status, warmstart = solve(m, optimizer, warmstart) + Base.Libc.flush_cstdio() status, warmstart end end