From 9511e34de4db6388765884a98edd6c93e1f59873 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 15:34:30 +0100 Subject: [PATCH 01/19] replace BlockDecomposition by SymbolicWedderburn --- Project.toml | 9 +- src/PropertyT.jl | 8 +- src/RGprojections.jl | 251 --------------------------------- src/blockdecomposition.jl | 287 -------------------------------------- 4 files changed, 2 insertions(+), 553 deletions(-) delete mode 100644 src/RGprojections.jl delete mode 100644 src/blockdecomposition.jl diff --git a/Project.toml b/Project.toml index 40d907f..155508e 100644 --- a/Project.toml +++ b/Project.toml @@ -4,23 +4,16 @@ authors = ["Marek Kaluba "] version = "0.3.2" [deps] -AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" -GroupRings = "0befed6a-bd73-11e8-1e41-a1190947c2f5" -Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" -JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2" [compat] -AbstractAlgebra = "^0.10.0" -GroupRings = "^0.3.2" -Groups = "^0.5.0" IntervalArithmetic = "^0.16.0" -JLD = "^0.9.0" JuMP = "^0.20.0" SCS = "^0.7.0" julia = "^1.3.0" diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 2ea9c30..8cba26d 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -1,22 +1,16 @@ __precompile__() module PropertyT -using AbstractAlgebra using LinearAlgebra using SparseArrays using Dates using Groups -using GroupRings +using SymbolicWedderburn -using JLD using JuMP -import AbstractAlgebra: Group, NCRing - include("laplacians.jl") -include("RGprojections.jl") -include("blockdecomposition.jl") include("sos_sdps.jl") include("checksolution.jl") diff --git a/src/RGprojections.jl b/src/RGprojections.jl deleted file mode 100644 index 50ceb86..0000000 --- a/src/RGprojections.jl +++ /dev/null @@ -1,251 +0,0 @@ -module Projections - -using AbstractAlgebra -using Groups -using GroupRings - -export PermCharacter, DirectProdCharacter, rankOne_projections - -############################################################################### -# -# Characters of Symmetric Group and DirectProduct -# -############################################################################### - -abstract type AbstractCharacter end - -struct PermCharacter <: AbstractCharacter - p::Generic.Partition -end - -struct DirectProdCharacter{N, T<:AbstractCharacter} <: AbstractCharacter - chars::NTuple{N, T} -end - -function (chi::DirectProdCharacter)(g::DirectPowerGroupElem) - res = 1 - for (χ, elt) in zip(chi.chars, g.elts) - res *= χ(elt) - end - return res -end - -function (chi::PermCharacter)(g::Generic.Perm) - R = AbstractAlgebra.partitionseq(chi.p) - p = Partition(Generic.permtype(g)) - return Int(Generic.MN1inner(R, p, 1, Generic._charvalsTable)) -end - -AbstractAlgebra.dim(χ::PermCharacter) = dim(YoungTableau(χ.p)) - -for T in [PermCharacter, DirectProdCharacter] - @eval begin - function (chi::$T)(X::GroupRingElem) - RG = parent(X) - z = zero(eltype(X)) - result = z - for i in 1:length(X.coeffs) - if X.coeffs[i] != z - result += chi(RG.basis[i])*X.coeffs[i] - end - end - return result - end - end -end - -characters(G::Generic.SymmetricGroup) = (PermCharacter(p) for p in AllParts(G.n)) - -function characters(G::DirectPowerGroup{N}) where N - nfold_chars = Iterators.repeated(characters(G.group), N) - return (DirectProdCharacter(idx) for idx in Iterators.product(nfold_chars...)) -end - -############################################################################### -# -# Projections -# -############################################################################### - -function central_projection(RG::GroupRing, chi::AbstractCharacter, T::Type=Rational{Int}) - result = RG(zeros(T, length(RG.basis))) - dim = chi(one(RG.group)) - ord = Int(order(RG.group)) - - for g in RG.basis - result[g] = convert(T, (dim//ord)*chi(g)) - end - - return result -end - -function alternating_emb(RG::GroupRing{Gr,T}, V::Vector{T}, S::Type=Rational{Int}) where {Gr<:AbstractAlgebra.AbstractPermutationGroup, T<:GroupElem} - res = RG(S) - for g in V - res[g] += sign(g) - end - return res -end - -function idempotents(RG::GroupRing{Generic.SymmetricGroup{S}}, T::Type=Rational{Int}) where S<:Integer - if RG.group.n == 1 - return GroupRingElem{T}[one(RG,T)] - elseif RG.group.n == 2 - Id = one(RG,T) - transp = RG(perm"(1,2)", T) - return GroupRingElem{T}[1//2*(Id + transp), 1//2*(Id - transp)] - end - - projs = Vector{Vector{Generic.Perm{S}}}() - for l in 2:RG.group.n - u = RG.group([circshift([i for i in 1:l], -1); [i for i in l+1:RG.group.n]]) - i = 0 - while (l-1)*i <= RG.group.n - v = RG.group(circshift(collect(1:RG.group.n), i)) - k = inv(v)*u*v - push!(projs, generateGroup([k], RG.group.n)) - i += 1 - end - end - - idems = Vector{GroupRingElem{T}}() - - for p in projs - append!(idems, [RG(p, T)//length(p), alternating_emb(RG, p, T)//length(p)]) - end - - return unique(idems) -end - -function rankOne_projection(chi::PermCharacter, idems::Vector{T}) where {T<:GroupRingElem} - - RG = parent(first(idems)) - S = eltype(first(idems)) - - ids = [one(RG, S); idems] - zzz = zero(S) - - for (i,j,k) in Base.product(ids, ids, ids) - if chi(i) == zzz || chi(j) == zzz || chi(k) == zzz - continue - else - elt = i*j*k - if elt^2 != elt - continue - elseif chi(elt) == one(S) - return elt - # return (i,j,k) - end - end - end - throw("Couldn't find rank-one projection for $chi") -end - -function rankOne_projections(RG::GroupRing{<:Generic.SymmetricGroup}, T::Type=Rational{Int}) - if RG.group.n == 1 - return [GroupRingElem([one(T)], RG)] - end - - RGidems = idempotents(RG, T) - - min_projs = [central_projection(RG,chi)*rankOne_projection(chi,RGidems) for chi in characters(RG.group)] - - return min_projs -end - -function ifelsetuple(a,b, k, n) - x = [repeat([a], k); repeat([b], n-k)] - return tuple(x...) -end - -function orbit_selector(n::Integer, k::Integer, - chi::AbstractCharacter, psi::AbstractCharacter) - return Projections.DirectProdCharacter(ifelsetuple(chi, psi, k, n)) -end - -function rankOne_projections(RBn::GroupRing{G}, T::Type=Rational{Int}) where {G<:WreathProduct} - - Bn = RBn.group - N = Bn.P.n - # projections as elements of the group rings RSₙ - Sn_rankOnePr = [rankOne_projections( - GroupRing(SymmetricGroup(i), collect(SymmetricGroup(i)))) - for i in typeof(N)(1):N] - - # embedding into group ring of BN - RN = GroupRing(Bn.N, collect(Bn.N)) - - sign, id = collect(characters(Bn.N.group)) - # Bn.N = (Z/2Z)ⁿ characters corresponding to the first k coordinates: - BnN_orbits = Dict(i => orbit_selector(N, i, sign, id) for i in 0:N) - - Q = Dict(i => RBn(g -> Bn(g), central_projection(RN, BnN_orbits[i], T)) for i in 0:N) - Q = Dict(key => GroupRings.dense(val) for (key, val) in Q) - - all_projs = [Q[0]*RBn(g->Bn(g), p) for p in Sn_rankOnePr[N]] - - r = collect(1:N) - for i in 1:N-1 - 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]) - - append!(all_projs, - [Q[i]*p1*p2 for (p1,p2) in Base.product(Sk_first,Sk_last)]) - end - - append!(all_projs, [Q[N]*RBn(g->Bn(g), p) for p in Sn_rankOnePr[N]]) - - return all_projs -end - -############################################################################## -# -# General Groups Misc -# -############################################################################## - -""" - products(X::Vector{GroupElem}, Y::Vector{GroupElem}[, op=*]) -Return a vector of all possible products (or `op(x,y)`), where `x ∈ X` and -`y ∈ Y`. You may change which operation is used by specifying `op` argument. -""" -function products(X::AbstractVector{T}, Y::AbstractVector{T}, op=*) where {T<:GroupElem} - result = Vector{T}() - seen = Set{T}() - sizehint!(result, length(X)*length(Y)) - sizehint!(seen, length(X)*length(Y)) - for x in X - for y in Y - z = op(x,y) - if !in(z, seen) - push!(seen, z) - push!(result, z) - end - end - end - return result -end - -""" - generateGroup(gens::Vector{GroupElem}, r, [op=*]) -Produce all elements of a group generated by elements in `gens` in ball of -radius `r` (word-length metric induced by `gens`). The binary operation can -be optionally specified. -""" -function generateGroup(gens::Vector{T}, r, op=*) where {T<:GroupElem} - n = 0 - R = 1 - elts = [one(first(gens)); gens] - while n ≠ length(elts) && R < r - # @show elts - R += 1 - n = length(elts) - elts = products(gens, elts, op) - end - return elts -end - -end # of module Projections diff --git a/src/blockdecomposition.jl b/src/blockdecomposition.jl deleted file mode 100644 index d4af161..0000000 --- a/src/blockdecomposition.jl +++ /dev/null @@ -1,287 +0,0 @@ -############################################################################### -# -# BlockDecomposition -# -############################################################################### - -struct BlockDecomposition{T<:AbstractArray{Float64, 2}, GEl<:GroupElem, P<:Generic.Perm} - orbits::Vector{Vector{Int}} - preps::Dict{GEl, P} - Uπs::Vector{T} - dims::Vector{Int} -end - -function BlockDecomposition(RG::GroupRing, autS::Group; verbose=true) - verbose && @info "Decomposing basis of RG into orbits of" autS - @time orbs = orbit_decomposition(autS, RG.basis, RG.basis_dict) - @assert sum(length(o) for o in orbs) == length(RG.basis) - verbose && @info "The action has $(length(orbs)) orbits" - - verbose && @info "Finding projections in the Group Ring of" autS - @time autS_mps = Projections.rankOne_projections(GroupRing(autS, collect(autS))) - - verbose && @info "Finding AutS-action matrix representation" - @time preps = perm_reps(autS, RG.basis[1:size(RG.pm,1)], RG.basis_dict) - @time mreps = matrix_reps(preps) - - verbose && @info "Computing the projection matrices Uπs" - @time Uπs = [orthSVD(matrix_repr(p, mreps)) for p in autS_mps] - - multiplicities = size.(Uπs,2) - dimensions = [Int(p[one(autS)]*Int(order(autS))) for p in autS_mps] - if verbose - info_strs = ["", - lpad("multiplicities", 14) * " =" * join(lpad.(multiplicities, 4), ""), - lpad("dimensions", 14) * " =" * join(lpad.(dimensions, 4), "") - ] - @info join(info_strs, "\n") - end - @assert dot(multiplicities, dimensions) == size(RG.pm,1) - - return BlockDecomposition(orbs, preps, Uπs, dimensions) -end - -function decimate(od::BlockDecomposition; verbose=true) - nzros = [i for i in 1:length(od.Uπs) if !isempty(od.Uπs[i])] - - Us = sparsify!.(od.Uπs, eps(Float64) * 1e4, verbose = verbose)[nzros] - #dimensions of the corresponding Uπs: - dims = od.dims[nzros] - - return BlockDecomposition(od.orbits, od.preps, Array{Float64}.(Us), dims) -end - -function orthSVD(M::AbstractMatrix{T}) where {T<:AbstractFloat} - fact = svd(convert(Matrix{T}, M)) - M_rank = sum(fact.S .> maximum(size(M)) * eps(T)) - return fact.U[:, 1:M_rank] -end - -orbit_decomposition( - G::Group, - E::AbstractVector, - rdict = GroupRings.reverse_dict(E); - op = ^, -) = orbit_decomposition(collect(G), E, rdict; op=op) - -function orbit_decomposition(elts::AbstractVector{<:GroupElem}, E::AbstractVector, rdict=GroupRings.reverse_dict(E); op=^) - - tovisit = trues(size(E)); - orbits = Vector{Vector{Int}}() - - orbit = zeros(Int, length(elts)) - - for i in eachindex(E) - if tovisit[i] - g = E[i] - Threads.@threads for j in eachindex(elts) - orbit[j] = rdict[op(g, elts[j])] - end - tovisit[orbit] .= false - push!(orbits, unique(orbit)) - end - end - return orbits -end - -############################################################################### -# -# Sparsification -# -############################################################################### - -dens(M::SparseMatrixCSC) = nnz(M)/length(M) -dens(M::AbstractArray) = count(!iszero, M)/length(M) - -function sparsify!(M::SparseMatrixCSC{Tv,Ti}, tol=eps(Tv); verbose=false) where {Tv,Ti} - - densM = dens(M) - droptol!(M, tol) - - verbose && @info( - "Sparsified density:", - rpad(densM, 20), - " → ", - rpad(dens(M), 20), - " ($(nnz(M)) non-zeros)" - ) - - return M -end - -function sparsify!(M::AbstractArray{T}, tol=eps(T); verbose=false) where T - densM = dens(M) - clamp_small!(M, tol) - - if verbose - @info("Sparsifying $(size(M))-matrix... \n $(rpad(densM, 20)) → $(rpad(dens(M),20))), ($(count(!iszero, M)) non-zeros)") - end - - return sparse(M) -end - -function clamp_small!(M::AbstractArray{T}, tol=eps(T)) where T - for n in eachindex(M) - if abs(M[n]) < tol - M[n] = zero(T) - end - end - return M -end - -function sparsify(U::AbstractArray{T}, tol=eps(T); verbose=false) where T - return sparsify!(deepcopy(U), tol, verbose=verbose) -end - -############################################################################### -# -# perm-, matrix-, representations -# -############################################################################### - -function perm_repr(g::GroupElem, E::Vector, E_dict) - p = Vector{Int}(undef, length(E)) - for (i,elt) in enumerate(E) - p[i] = E_dict[elt^g] - end - return p -end - -function perm_reps(G::Group, E::Vector, E_rdict=GroupRings.reverse_dict(E)) - elts = collect(G) - l = length(elts) - preps = Vector{Generic.Perm}(undef, l) - - permG = SymmetricGroup(length(E)) - - Threads.@threads for i in 1:l - preps[i] = permG(PropertyT.perm_repr(elts[i], E, E_rdict), false) - end - - return Dict(elts[i]=>preps[i] for i in 1:l) -end - -function matrix_repr(x::GroupRingElem, mreps::Dict) - nzeros = findall(!iszero, x.coeffs) - return sum(x[i].*mreps[parent(x).basis[i]] for i in nzeros) -end - -function matrix_reps(preps::Dict{T,Generic.Perm{I}}) where {T<:GroupElem, I<:Integer} - kk = collect(keys(preps)) - mreps = Vector{SparseMatrixCSC{Float64, Int}}(undef, length(kk)) - Threads.@threads for i in 1:length(kk) - mreps[i] = AbstractAlgebra.matrix_repr(preps[kk[i]]) - end - return Dict(kk[i] => mreps[i] for i in 1:length(kk)) -end - -############################################################################### -# -# actions -# -############################################################################### - -function Base.:^(y::GroupRingElem, g::GroupRingElem, op = ^) - res = parent(y)() - for elt in GroupRings.supp(g) - res += g[elt] * ^(y, elt, op) - end - return res -end - -function Base.:^(y::GroupRingElem, g::GroupElem, op = ^) - RG = parent(y) - result = zero(RG, eltype(y.coeffs)) - - for (idx, c) in enumerate(y.coeffs) - if !iszero(c) - result[op(RG.basis[idx], g)] = c - end - end - return result -end - -function Base.:^( - y::GroupRingElem{T,<:SparseVector}, - g::GroupElem, - op = ^, -) where {T} - RG = parent(y) - index = [RG.basis_dict[op(RG.basis[idx], g)] for idx in y.coeffs.nzind] - - result = GroupRingElem(sparsevec(index, y.coeffs.nzval, y.coeffs.n), RG) - - return result -end - -############################################################################### -# -# perm && WreathProductElems actions: MatAlgElem -# -############################################################################### - -function Base.:^(A::MatAlgElem, p::Generic.Perm) - 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 = 1:size(A, 1) - for j = 1:size(A, 2) - result[p[i], p[j]] = A[i, j] # action by permuting rows and colums/conjugation - end - end - return result -end - -function Base.:^(A::MatAlgElem, g::WreathProductElem{N}) 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) - R = base_ring(parent(A)) - tmp = R(1) - - @inbounds for i = 1:size(A, 1) - for j = 1:size(A, 2) - x = A[i, j] - if flips[i] * flips[j] == 1 - result[g.p[i], g.p[j]] = x - else - result[g.p[i], g.p[j]] = -x - end - end - end - return result -end - -############################################################################### -# -# perm && WreathProductElems actions: Automorphism -# -############################################################################### - -function Base.:^(a::Automorphism, g::GroupElem) - Ag = parent(a)(g) - Ag_inv = inv(Ag) - res = append!(Ag, a, Ag_inv) - return Groups.freereduce!(res) -end - -(A::AutGroup)(p::Generic.Perm) = A(Groups.AutSymbol(p)) - -function (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 = one(A) - Id = one(parent(g.n.elts[1])) - for i = 1:length(g.p.d) - if g.n.elts[i] != Id - push!(elt, Groups.flip(i)) - end - end - push!(elt, Groups.AutSymbol(g.p)) - return elt -end - - -# fallback: -Base.one(p::Generic.Perm) = Perm(length(p.d)) From bb0354d3a05308c0062de44d2ad41695d51ea9ec Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 15:42:10 +0100 Subject: [PATCH 02/19] rewrite sos_spds.jl This includes changes related to: * SymbolicWedderburn * new formulation of constraints * general warmstart for JuMP-^1.3 --- src/PropertyT.jl | 7 +- src/sos_sdps.jl | 525 ++++++++++++++++++++++++++++------------------- 2 files changed, 316 insertions(+), 216 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 8cba26d..da53cfa 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -5,11 +5,12 @@ using LinearAlgebra using SparseArrays using Dates -using Groups -using SymbolicWedderburn - using JuMP +using Groups +using StarAlgebras +using SymbolicWedderburn + include("laplacians.jl") include("sos_sdps.jl") include("checksolution.jl") diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index ea2e20b..5c77de4 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -1,260 +1,359 @@ -############################################################################### -# -# Constraints -# -############################################################################### +""" + sos_problem_dual(X, [u = zero(X); upper_bound=Inf]) +Formulate the dual to the sum of squares decomposition problem for `X - λ·u`. -function constraints(pm::Matrix{I}, total_length=maximum(pm)) where {I<:Integer} - cnstrs = [Vector{I}() for _ in 1:total_length] - for i in eachindex(pm) - push!(cnstrs[pm[i]], i) +See also [sos_problem_primal](@ref). +""" +function sos_problem_dual( + elt::AlgebraElement, + order_unit::AlgebraElement=zero(elt); + lower_bound=-Inf +) + @assert parent(elt) == parent(order_unit) + algebra = parent(elt) + mstructure = if StarAlgebras._istwisted(algebra.mstructure) + algebra.mstructure + else + StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure)) + end + + # 1 variable for every primal constraint + # 1 dual variable for every element of basis + # Symmetrized: + # 1 dual variable for every orbit of G acting on basis + model = Model() + @variable model y[1:length(basis(algebra))] + @constraint model λ_dual dot(order_unit, y) == 1 + @constraint(model, psd, y[mstructure] in PSDCone()) + + if !isinf(lower_bound) + throw("Not Implemented yet") + @variable model λ_ub_dual + @objective model Min dot(elt, y) + lower_bound * λ_ub_dual + else + @objective model Min dot(elt, y) + end + + return model +end + +function constraints( + mstr::AbstractMatrix{<:Integer}, + total_length=maximum(mstr), +) + cnstrs = [Vector{eltype(mstr)}() for _ = 1:total_length] + li = LinearIndices(mstr) + + for (idx, val) in pairs(mstr) + push!(cnstrs[val], li[idx]) end return cnstrs end -function orbit_constraint!(result::SparseMatrixCSC, cnstrs, orbit; val=1.0/length(orbit)) +function constraints( + basis::StarAlgebras.AbstractBasis, + mstr::AbstractMatrix{<:Integer}; + augmented::Bool=false, + table_size=size(mstr) +) + cnstrs = [signed(eltype(mstr))[] for _ in basis] + LI = LinearIndices(table_size) + + for ci in CartesianIndices(table_size) + k = LI[ci] + a_star_b = basis[mstr[k]] + push!(cnstrs[basis[a_star_b]], k) + if augmented + # (1-a_star)(1-b) = 1 - a_star - b + a_star_b + + i, j = Tuple(ci) + a, b = basis[i], basis[j] + + push!(cnstrs[basis[one(a)]], k) + push!(cnstrs[basis[star(a)]], -k) + push!(cnstrs[basis[b]], -k) + end + end + + return Dict( + basis[i] => ConstraintMatrix(c, table_size..., 1) for (i, c) in pairs(cnstrs) + ) +end + +function constraints(A::StarAlgebra; augmented::Bool, twisted::Bool) + mstructure = if StarAlgebras._istwisted(A.mstructure) == twisted + A.mstructure + else + StarAlgebras.MTable{twisted}(basis(A), table_size=size(A.mstructure)) + end + return constraints(basis(A), mstructure, augmented=augmented) +end + +""" + sos_problem_primal(X, [u = zero(X); upper_bound=Inf]) +Formulate sum of squares decomposition problem for `X - λ·u`. + +Given element `X` of a star-algebra `A` and an order unit `u` of `Σ²A` the sum +of squares cone in `A`, fomulate sum of squares decomposition problem: + +``` +max: λ +subject to: X - λ·u ∈ Σ²A +``` + +If `upper_bound` is given a finite value, additionally `λ ≤ upper_bound` will +be added to the model. This may improve the accuracy of the solution if +`upper_bound` is less than the optimal `λ`. + +The default `u = zero(X)` formulates a simple feasibility problem. +""" +function sos_problem_primal( + elt::AlgebraElement, + order_unit::AlgebraElement=zero(elt); + upper_bound=Inf, + augmented::Bool=iszero(aug(elt)) && iszero(aug(order_unit)) +) + @assert parent(elt) === parent(order_unit) + + N = LinearAlgebra.checksquare(parent(elt).mstructure) + model = JuMP.Model() + P = JuMP.@variable(model, P[1:N, 1:N], Symmetric) + JuMP.@constraint(model, psd, P in PSDCone()) + + if iszero(order_unit) && isfinite(upper_bound) + @warn "Setting `upper_bound` together with zero `order_unit` has no effect" + end + + A = constraints(parent(elt), augmented=augmented, twisted=true) + + if !iszero(order_unit) + λ = JuMP.@variable(model, λ) + if isfinite(upper_bound) + JuMP.@constraint model λ <= upper_bound + end + JuMP.@objective(model, Max, λ) + + for b in basis(parent(elt)) + JuMP.@constraint(model, elt(b) - λ * order_unit(b) == dot(A[b], P)) + end + else + for b in basis(parent(elt)) + JuMP.@constraint(model, elt(b) == dot(A[b], P)) + end + end + + return model +end + +function invariant_constraint!( + result::AbstractMatrix, + basis::StarAlgebras.AbstractBasis, + cnstrs::AbstractDict{K,CM}, + invariant_vec::SparseVector, +) where {K,CM<:ConstraintMatrix} result .= zero(eltype(result)) - dropzeros!(result) - for constraint in cnstrs[orbit] - for idx in constraint - result[idx] = val + for i in SparseArrays.nonzeroinds(invariant_vec) + g = basis[i] + A = cnstrs[g] + for (idx, v) in nzpairs(A) + result[idx] += invariant_vec[i] * v end end return result end -function orbit_spvector(vect::AbstractVector, orbits) - orb_vector = spzeros(length(orbits)) - - for (i,o) in enumerate(orbits) - k = vect[collect(o)] - val = k[1] - @assert all(k .== val) - orb_vector[i] = val - end - - return orb_vector +function isorth_projection(ds::SymbolicWedderburn.DirectSummand) + U = SymbolicWedderburn.image_basis(ds) + return isapprox(U * U', I) end -############################################################################### -# -# Naive SDP -# -############################################################################### +sos_problem_primal( + elt::AlgebraElement, + wedderburn::WedderburnDecomposition; + kwargs... +) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...) -function SOS_problem_dual(elt::GroupRingElem, order_unit::GroupRingElem; - lower_bound::Float64=Inf) - @assert parent(elt) == parent(order_unit) +function sos_problem_primal( + elt::AlgebraElement, + orderunit::AlgebraElement, + wedderburn::WedderburnDecomposition; + upper_bound=Inf, + augmented=iszero(aug(elt)) && iszero(aug(orderunit)) +) - RG = parent(elt) - m = Model() - - y = @variable(m, y[1:length(elt.coeffs)]) - - @constraint(m, λ_dual, dot(order_unit.coeffs, y) == 1) - @constraint(m, psd, [y[i] for i in RG.pm] in PSDCone()) - - if !isinf(lower_bound) - @variable(m, λ_ub_dual) - expr = dot(elt.coeffs, y) + lower_bound*λ_ub_dual - # @constraint m expr >= lower_bound - @objective m Min expr - else - @objective m Min dot(elt.coeffs, y) + @assert parent(elt) === parent(orderunit) + if any(!isorth_projection, direct_summands(wedderburn)) + error("Wedderburn decomposition contains a non-orthogonal projection") end - return m -end + feasibility_problem = iszero(orderunit) -function SOS_problem_primal(X::GroupRingElem, orderunit::GroupRingElem; - upper_bound::Float64=Inf) + model = JuMP.Model() + if !feasibility_problem # add λ or not? + λ = JuMP.@variable(model, λ) + JuMP.@objective(model, Max, λ) - N = size(parent(X).pm, 1) - m = JuMP.Model(); - - JuMP.@variable(m, P[1:N, 1:N]) - # SP = Symmetric(P) - JuMP.@SDconstraint(m, sdp, P >= 0) - - if iszero(aug(X)) && iszero(aug(orderunit)) - JuMP.@constraint(m, augmentation, sum(P) == 0) - end - - if upper_bound < Inf - λ = JuMP.@variable(m, λ <= upper_bound) - else - λ = JuMP.@variable(m, λ) - end - - cnstrs = constraints(parent(X).pm) - @assert length(cnstrs) == length(X.coeffs) == length(orderunit.coeffs) - - x, u = X.coeffs, orderunit.coeffs - JuMP.@constraint(m, lincnstr[i=1:length(cnstrs)], - x[i] - λ*u[i] == sum(P[cnstrs[i]])) - - JuMP.@objective(m, Max, λ) - - return m -end - -############################################################################### -# -# Symmetrized SDP -# -############################################################################### - -function SOS_problem_primal(X::GroupRingElem, orderunit::GroupRingElem, data::BlockDecomposition; upper_bound::Float64=Inf) - Ns = size.(data.Uπs, 2) - m = JuMP.Model(); - - Ps = Vector{Matrix{JuMP.VariableRef}}(undef, length(Ns)) - - for (k,s) in enumerate(Ns) - Ps[k] = JuMP.@variable(m, [1:s, 1:s]) - JuMP.@SDconstraint(m, Ps[k] >= 0) - end - - if upper_bound < Inf - λ = JuMP.@variable(m, λ <= upper_bound) - else - λ = JuMP.@variable(m, λ) - end - - @info "Adding $(length(data.orbits)) constraints..." - @time addconstraints!(m, Ps, X, orderunit, data) - - JuMP.@objective(m, Max, λ) - - return m, Ps -end - -function constraintLHS!(M, cnstr, Us, Ust, dims, eps=1000*eps(eltype(first(M)))) - for π in 1:lastindex(Us) - M[π] = dims[π].*PropertyT.clamp_small!(Ust[π]*(cnstr*Us[π]), eps) - end -end - -function addconstraints!(m::JuMP.Model, - P::Vector{Matrix{JuMP.VariableRef}}, - X::GroupRingElem, orderunit::GroupRingElem, data::BlockDecomposition) - - orderunit_orb = orbit_spvector(orderunit.coeffs, data.orbits) - X_orb = orbit_spvector(X.coeffs, data.orbits) - UπsT = [U' for U in data.Uπs] - - cnstrs = constraints(parent(X).pm) - orb_cnstr = spzeros(Float64, size(parent(X).pm)...) - - M = [Array{Float64}(undef, n,n) for n in size.(UπsT,1)] - - λ = m[:λ] - - for (t, orbit) in enumerate(data.orbits) - orbit_constraint!(orb_cnstr, cnstrs, orbit) - constraintLHS!(M, orb_cnstr, data.Uπs, UπsT, data.dims) - - x, u = X_orb[t], orderunit_orb[t] - - JuMP.@constraints m begin - x - λ*u == sum(dot(M[π], P[π]) for π in eachindex(data.Uπs)) - end - end - return m -end - -function reconstruct(Ps::Vector{Matrix{F}}, data::BlockDecomposition) where F - return reconstruct(Ps, data.preps, data.Uπs, data.dims) -end - -function reconstruct(Ps::Vector{M}, - preps::Dict{GEl, P}, Uπs::Vector{U}, dims::Vector{Int}) where - {M<:AbstractMatrix, GEl<:GroupElem, P<:Generic.Perm, U<:AbstractMatrix} - - lU = length(Uπs) - transfP = [dims[π].*Uπs[π]*Ps[π]*Uπs[π]' for π in 1:lU] - tmp = [zeros(Float64, size(first(transfP))) for _ in 1:lU] - - Threads.@threads for π in 1:lU - tmp[π] = perm_avg!(tmp[π], transfP[π], values(preps)) - end - - recP = sum(tmp)./length(preps) - - return recP -end - -function perm_avg!(result, P, perms) - lp = length(first(perms).d) - for p in perms - # result .+= view(P, p.d, p.d) - @inbounds for j in 1:lp - k = p[j] - for i in 1:lp - result[i,j] += P[p[i], k] + if isfinite(upper_bound) + JuMP.@constraint(model, λ <= upper_bound) + if feasibility_problem + @warn "setting `upper_bound` with zero `orderunit` has no effect" end end end - return result -end -############################################################################### -# -# Low-level solve -# -############################################################################### - -function setwarmstart!(m::JuMP.Model, warmstart) - if solver_name(m) == "SCS" - primal, dual, slack = warmstart - m.moi_backend.optimizer.model.optimizer.data.primal = primal - m.moi_backend.optimizer.model.optimizer.data.dual = dual - m.moi_backend.optimizer.model.optimizer.data.slack = slack - else - @warn "Setting warmstart for $(solver_name(m)) is not implemented! Ignoring..." + 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()) + P end - return m -end -function getwarmstart(m::JuMP.Model) - if solver_name(m) == "SCS" - warmstart = ( - primal = m.moi_backend.optimizer.model.optimizer.data.primal, - dual = m.moi_backend.optimizer.model.optimizer.data.dual, - slack = m.moi_backend.optimizer.model.optimizer.data.slack - ) - else - @warn "Saving warmstart for $(solver_name(m)) is not implemented!" - return (primal=Float64[], dual=Float64[], slack=Float64[]) + begin # preallocating + T = eltype(wedderburn) + M = zeros.(T, size.(P)) + M_orb = zeros(T, size(parent(elt).mstructure)...) + tmps = SymbolicWedderburn._tmps(wedderburn) end - return warmstart + + X = convert(Vector{T}, coeffs(elt)) + U = convert(Vector{T}, coeffs(orderunit)) + + # defining constraints based on the multiplicative structure + cnstrs = constraints(parent(elt), augmented=augmented, twisted=true) + + @info "Adding $(length(invariant_vectors(wedderburn))) constraints" + + for iv in invariant_vectors(wedderburn) + 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, wedderburn, tmps) + SymbolicWedderburn.zerotol!.(M, atol=1e-12) + + M_dot_P = sum(dot(M[π], P[π]) for π in eachindex(M) if !iszero(M[π])) + + if feasibility_problem + JuMP.@constraint(model, x == M_dot_P) + else + JuMP.@constraint(model, x - λ * u == M_dot_P) + end + end + return model, P end -function solve(m::JuMP.Model, with_optimizer::JuMP.OptimizerFactory, warmstart=nothing) +function reconstruct(Ps, wd::WedderburnDecomposition) + N = size(first(direct_summands(wd)), 2) + P = zeros(eltype(wd), N, N) + return reconstruct!(P, Ps, wd) +end - set_optimizer(m, with_optimizer) +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 .*= 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 ./= order(Int, G) + + 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) - if warmstart != nothing - setwarmstart!(m, warmstart) - end + m = setwarmstart!(m, warmstart) - optimize!(m) + JuMP.optimize!(m) Base.Libc.flush_cstdio() - status = termination_status(m) + status = JuMP.termination_status(m) return status, getwarmstart(m) end -function solve(solverlog::String, m::JuMP.Model, with_optimizer::JuMP.OptimizerFactory, warmstart=nothing) +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 = PropertyT.solve(m, with_optimizer, warmstart) + status, warmstart = solve(m, optimizer, warmstart) status, warmstart end end From ecea3dfbcba5289322b944666956e9df92cbcc1f Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 15:45:18 +0100 Subject: [PATCH 03/19] add ConstraintMatrix --- src/PropertyT.jl | 1 + src/constraint_matrix.jl | 129 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 130 insertions(+) create mode 100644 src/constraint_matrix.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index da53cfa..72e9247 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -12,6 +12,7 @@ using StarAlgebras using SymbolicWedderburn include("laplacians.jl") +include("constraint_matrix.jl") include("sos_sdps.jl") include("checksolution.jl") diff --git a/src/constraint_matrix.jl b/src/constraint_matrix.jl new file mode 100644 index 0000000..8fbd36a --- /dev/null +++ b/src/constraint_matrix.jl @@ -0,0 +1,129 @@ +""" + ConstraintMatrix{T,I} <: AbstractMatrix{T} + +Special type of sparse matrix used to store constraints in SOS problems. + +This matrix has in general very few non-zero values which also are multiples of each other. + +The constructor accepts +* `nzeros` - a vector of non-zero indices; negative values are used to signify + negative values; repetitions are allowed +* `n`, `m` - the size of matrix +* `val` - the greatest common factor of the values + +To iterate efficiently over `A::ConstraintMatrix` use [`nzpairs(A)`](@ref). + +# Examples +```julia-repl +julia> ConstraintMatrix{Float64}([-1,2,-1,1,4,2,6], 3,2, π) +3×2 ConstraintMatrix{Float64, Int64}: + -3.14159 3.14159 + 6.28319 0.0 + 0.0 3.14159 +``` +""" +struct ConstraintMatrix{T,I} <: AbstractMatrix{T} + pos::Vector{I} # list of positive indices + neg::Vector{I} # list of negative indices + size::Tuple{Int,Int} + val::T + + function ConstraintMatrix{T}(nzeros::AbstractArray{<:Integer}, n, m, val) where {T} + @assert n ≥ 1 + @assert m ≥ 1 + + if !isempty(nzeros) + sort!(nzeros) + a, b = first(nzeros), last(nzeros) + @assert 1 ≤ abs(a) ≤ n * m + @assert 1 ≤ abs(b) ≤ n * m + end + k = searchsortedlast(nzeros, 0) + neg = @view nzeros[begin:k] + pos = @view nzeros[k+1:end] + return new{T,eltype(nzeros)}(pos, -neg, (n, m), val) + end +end + +ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} = + ConstraintMatrix{T}(nzeros, n, m, val) + +Base.size(cm::ConstraintMatrix) = cm.size + +__get_positive(cm::ConstraintMatrix, idx::Integer) = + convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx))) +__get_negative(cm::ConstraintMatrix, idx::Integer) = + convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx))) + +Base.@propagate_inbounds function Base.getindex( + cm::ConstraintMatrix, + i::Integer, + j::Integer, +) + li = LinearIndices(cm) + idx = li[i, j] + pos = __get_positive(cm, idx) + neg = __get_negative(cm, idx) + return pos - neg +end + +struct NZPairsIter{T} + m::ConstraintMatrix{T} +end + +Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T} +Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown() + +# TODO: iterate over (idx=>val) pairs combining vals +function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 1)) + k = iterate(itr.m.pos, state[1]) + isnothing(k) && return iterate(itr, state[2]) + idx, st = k + return idx => itr.m.val, (st, 1) +end + +function Base.iterate(itr::NZPairsIter, state::Int) + k = iterate(itr.m.neg, state[1]) + isnothing(k) && return nothing + idx, st = k + return idx => -itr.m.val, st +end + +""" + nzpairs(cm::ConstraintMatrix) +Efficiently iterate over non-zero `(idx=>value)` pairs. + +If the `cm` was created with repetitions (or contains negative values) there will +be repetitions in the returned sequence of pairs. + +# Examples +```julia +julia> ConstraintMatrix{Float64}([-1,2,-1,1,4,2,6], 3,2, π) +3×2 ConstraintMatrix{Float64, Int64}: + -3.14159 3.14159 + 6.28319 0.0 + 0.0 3.14159 + +julia> collect(nzpairs(M)) +7-element Vector{Pair{Int64, Float64}}: + 1 => 3.141592653589793 + 2 => 3.141592653589793 + 2 => 3.141592653589793 + 4 => 3.141592653589793 + 6 => 3.141592653589793 + 1 => -3.141592653589793 + 1 => -3.141592653589793 +``` +""" +nzpairs(cm::ConstraintMatrix) = NZPairsIter(cm) + +function LinearAlgebra.dot(cm::ConstraintMatrix, m::AbstractMatrix{T}) where {T} + if isempty(cm.pos) && isempty(cm.neg) + isempty(m) && return zero(T) + return zero(first(m) + first(m)) + end + + pos = isempty(cm.pos) ? zero(first(m)) : sum(@view m[cm.pos]) + neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg]) + return convert(eltype(cm), cm.val) * (pos - neg) +end From 085f6bce3c5561049934c57f4a6fc577ea3a1cdb Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 16:01:35 +0100 Subject: [PATCH 04/19] replace checksoltuion by certify --- src/PropertyT.jl | 3 +- src/certify.jl | 168 +++++++++++++++++++++++++++++++++++++++++++ src/checksolution.jl | 77 -------------------- 3 files changed, 170 insertions(+), 78 deletions(-) create mode 100644 src/certify.jl delete mode 100644 src/checksolution.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 72e9247..2e3d9bf 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -5,6 +5,7 @@ using LinearAlgebra using SparseArrays using Dates +using IntervalArithmetic using JuMP using Groups @@ -14,7 +15,7 @@ using SymbolicWedderburn include("laplacians.jl") include("constraint_matrix.jl") include("sos_sdps.jl") -include("checksolution.jl") +include("certify.jl") include("1712.07167.jl") include("1812.03456.jl") diff --git a/src/certify.jl b/src/certify.jl new file mode 100644 index 0000000..bca18c6 --- /dev/null +++ b/src/certify.jl @@ -0,0 +1,168 @@ +function augment_columns!(Q::AbstractMatrix) + for c in eachcol(Q) + c .-= sum(c) ./ length(c) + end + return Q +end + +function _fma_SOS_thr!( + result::AbstractVector{T}, + mstructure::AbstractMatrix{<:Integer}, + Q::AbstractMatrix{T}, + acc_matrix=zeros(T, size(mstructure)...), +) where {T} + + 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 + end + end + end + + @inbounds for j = 1:s2 + for i = 1:s1 + result[mstructure[i, j]] += acc_matrix[i, j] + end + end + + return result +end + +function _cnstr_sos!(res::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::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::StarAlgebra, Q::AbstractMatrix; augmented::Bool) + if augmented + z = zeros(eltype(Q), length(basis(A))) + res = 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 AlgebraElement(z, A) + end +end + +function sufficient_λ(residual::AlgebraElement, λ; halfradius) + L1_norm = norm(residual, 1) + suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm + + eq_sign = let T = eltype(residual) + if T <: Interval + "∈" + elseif T <: Union{Rational,Integer} + "=" + else # if T <: AbstractFloat + "≈" + end + end + + info_strs = [ + "Numerical metrics of the obtained SOS:", + "ɛ(elt - λu - ∑ξᵢ*ξᵢ) $eq_sign $(aug(residual))", + "‖elt - λu - ∑ξᵢ*ξᵢ‖₁ $eq_sign $(L1_norm)", + " λ $eq_sign $suff_λ", + ] + @info join(info_strs, "\n") + + return suff_λ +end + +function sufficient_λ( + elt::AlgebraElement, + order_unit::AlgebraElement, + λ, + sos::AlgebraElement; + halfradius +) + + @assert parent(elt) === parent(order_unit) == parent(sos) + residual = (elt - λ * order_unit) - sos + + return sufficient_λ(residual, λ; halfradius=halfradius) +end + +function certify_solution( + elt::AlgebraElement, + orderunit::AlgebraElement, + λ, + Q::AbstractMatrix{<:AbstractFloat}; + halfradius, + augmented=iszero(aug(elt)) && iszero(aug(orderunit)) +) + + should_we_augment = !augmented && aug(elt) == aug(orderunit) == 0 + + Q = should_we_augment ? augment_columns!(Q) : Q + @time sos = compute_sos(parent(elt), Q, augmented=augmented) + + @info "Checking in $(eltype(sos)) arithmetic with" λ + + λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius) + + if λ_flpoint ≤ 0 + return false, λ_flpoint + end + + λ_int = @interval(λ) + Q_int = [@interval(q) for q in Q] + + check, sos_int = @time if should_we_augment + @info("Projecting columns of Q to the augmentation ideal...") + Q_int = augment_columns!(Q_int) + @info "Checking that sum of every column contains 0.0..." + check_augmented = all(0 ∈ sum(c) for c in eachcol(Q_int)) + check_augmented || @error( + "Augmentation failed! The following numbers are not certified!" + ) + sos_int = compute_sos(parent(elt), Q_int; augmented=augmented) + check_augmented, sos_int + else + true, compute_sos(parent(elt), Q_int, augmented=augmented) + end + + @info "Checking in $(eltype(sos_int)) arithmetic with" λ + + λ_certified = + sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius) + + return check && inf(λ_certified) > 0.0, inf(λ_certified) +end diff --git a/src/checksolution.jl b/src/checksolution.jl deleted file mode 100644 index 8a4a166..0000000 --- a/src/checksolution.jl +++ /dev/null @@ -1,77 +0,0 @@ -using IntervalArithmetic - -IntervalArithmetic.setrounding(Interval, :tight) -IntervalArithmetic.setformat(sigfigs=12) - -function fma_SOS_thr!(result::AbstractVector{T}, pm::AbstractMatrix{<:Integer}, - Q::AbstractMatrix{T}, acc_matrix=zeros(T, size(pm)...)) where T - - s1, s2 = size(pm) - - @inbounds for k in 1:s2 - let k=k, s1=s1, s2=s2, Q=Q, acc_matrix=acc_matrix - Threads.@threads for j in 1:s2 - for i in 1:s1 - @inbounds acc_matrix[i,j] = muladd(Q[i, k], Q[j, k], acc_matrix[i,j]) - end - end - end - end - - @inbounds for j in 1:s2 - for i in 1:s1 - result[pm[i,j]] += acc_matrix[i,j] - end - end - - return result -end - -function compute_SOS(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix) - result = zeros(eltype(Q), maximum(pm)); - return fma_SOS_thr!(result, pm, Q) -end - -function compute_SOS(RG::GroupRing, Q::AbstractMatrix{<:Real}) - result = compute_SOS(RG.pm, Q) - return GroupRingElem(result, RG) -end - -function compute_SOS_square(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix) - result = zeros(eltype(Q), maximum(pm)); - - for i in 1:size(Q,2) - GroupRings.fmac!(result, view(Q,:,i), view(Q,:,i), pm) - end - - return result -end - -function compute_SOS_square(RG::GroupRing, Q::AbstractMatrix) - return GroupRingElem(compute_SOS_square(RG.pm, Q), RG) -end - -function augIdproj(Q::AbstractMatrix{T}) where T - result = zeros(T, size(Q)) - l = size(Q, 2) - Threads.@threads for j in 1:l - col = sum(view(Q, :,j))/l - for i in 1:size(Q, 1) - result[i,j] = Q[i,j] - col - end - end - return result -end - -function augIdproj(::Type{Interval}, Q::AbstractMatrix{T}) where {T<:Real} - result = zeros(Interval{T}, size(Q)) - l = size(Q, 2) - Threads.@threads for j in 1:l - col = sum(view(Q, :,j))/l - for i in 1:size(Q, 1) - result[i,j] = @interval(Q[i,j] - col) - end - end - check = all([zero(T) in sum(view(result, :, i)) for i in 1:size(result, 2)]) - return result, check -end From 4b8efd2a40849e5db087979aa62d943945e7b94d Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 16:09:02 +0100 Subject: [PATCH 05/19] adjust sq, adj, op to the new world --- src/PropertyT.jl | 2 + src/sqadjop.jl | 134 ++++++++++++++++++++--------------------------- 2 files changed, 60 insertions(+), 76 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 2e3d9bf..10d74ba 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -17,6 +17,8 @@ include("constraint_matrix.jl") include("sos_sdps.jl") include("certify.jl") +include("sqadjop.jl") + include("1712.07167.jl") include("1812.03456.jl") diff --git a/src/sqadjop.jl b/src/sqadjop.jl index 33aa056..2915660 100644 --- a/src/sqadjop.jl +++ b/src/sqadjop.jl @@ -1,98 +1,80 @@ -isopposite(σ::Generic.Perm, τ::Generic.Perm, i=1, j=2) = - σ[i] ≠ τ[i] && σ[i] ≠ τ[j] && - σ[j] ≠ τ[i] && σ[j] ≠ τ[j] +import PermutationGroups.AbstractPerm -isadjacent(σ::Generic.Perm, τ::Generic.Perm, i=1, j=2) = - (σ[i] == τ[i] && σ[j] ≠ τ[j]) || # first equal, second differ - (σ[j] == τ[j] && σ[i] ≠ τ[i]) || # second equal, first differ - (σ[i] == τ[j] && σ[j] ≠ τ[i]) || # first σ equal to second τ - (σ[j] == τ[i] && σ[i] ≠ τ[j]) # second σ equal to first τ +# move to Groups +Base.keys(a::Alphabet) = keys(1:length(a)) -Base.div(X::GroupRingElem, x::Number) = parent(X)(X.coeffs.÷x) +## the old 1812.03456 definitions -function Sq(RG::GroupRing, N::Integer) - S₂ = generating_set(RG.group, 2) - ℤ = Int64 - Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ); +isopposite(σ::AbstractPerm, τ::AbstractPerm, i=1, j=2) = + i^σ ≠ i^τ && i^σ ≠ j^τ && j^σ ≠ i^τ && j^σ ≠ j^τ - Alt_N = [g for g in SymmetricGroup(N) if parity(g) == 0] +isadjacent(σ::AbstractPerm, τ::AbstractPerm, i=1, j=2) = + (i^σ == i^τ && j^σ ≠ j^τ) || # first equal, second differ + (j^σ == j^τ && i^σ ≠ i^τ) || # second equal, first differ + (i^σ == j^τ && j^σ ≠ i^τ) || # first σ equal to second τ + (j^σ == i^τ && i^σ ≠ j^τ) # second σ equal to first τ - sq = RG() - for σ in Alt_N - GroupRings.addeq!(sq, *(Δ₂^σ, Δ₂^σ, false)) +function _ncycle(start, length, n=start + length - 1) + p = Perm(Int8(n)) + @assert n ≥ start + length - 1 + for k in start:start+length-2 + p[k] = k + 1 end - return sq÷factorial(N-2) + p[start+length-1] = start + return p end -function Adj(RG::GroupRing, N::Integer) - S₂ = generating_set(RG.group, 2) - ℤ = Int64 - Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ); +alternating_group(n::Integer) = PermGroup([_ncycle(i, 3) for i in 1:n-2]) - Alt_N = [g for g in SymmetricGroup(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 +function small_gens(G::MatrixGroups.SpecialLinearGroup) + A = alphabet(G) + S = map([(1, 2), (2, 1)]) do (i, j) + k = findfirst(l -> (l.i == i && l.j == j), A) + @assert !isnothing(k) + G(Groups.word_type(G)([k])) end - return adj÷factorial(N-2)^2 + return union!(S, inv.(S)) end -function Op(RG::GroupRing, N::Integer) - if N < 4 - return RG() +function small_gens(G::Groups.AutomorphismGroup{<:FreeGroup}) + A = alphabet(G) + S = map([(:ϱ, 1, 2), (:ϱ, 2, 1), (:λ, 1, 2), (:λ, 2, 1)]) do (id, i, j) + k = findfirst(l -> (l.id == id && l.i == i && l.j == j), A) + @assert !isnothing(k) + G(Groups.word_type(G)([k])) end - S₂ = generating_set(RG.group, 2) - ℤ = Int64 - Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ); + return union!(S, inv.(S)) +end - Alt_N = [g for g in SymmetricGroup(N) if parity(g) == 0] - Δ₂s = Dict(σ=>Δ₂^σ for σ in Alt_N) - op = RG() +function small_laplacian(RG::StarAlgebra) + G = StarAlgebras.object(RG) + S₂ = small_gens(G) + return length(S₂) * one(RG) - sum(RG(s) for s in S₂) +end - for σ in Alt_N - for τ in Alt_N +function SqAdjOp(A::StarAlgebra, n::Integer, Δ₂=small_laplacian(A)) + @assert parent(Δ₂) === A + + alt_n = alternating_group(n) + G = StarAlgebras.object(A) + act = action_by_conjugation(G, alt_n) + + Δ₂s = Dict(σ => SymbolicWedderburn.action(act, σ, Δ₂) for σ in alt_n) + + sq, adj, op = zero(A), zero(A), zero(A) + tmp = zero(A) + + for σ in alt_n + StarAlgebras.add!(sq, sq, StarAlgebras.mul!(tmp, Δ₂s[σ], Δ₂s[σ])) + for τ in alt_n if isopposite(σ, τ) - GroupRings.addeq!(op, *(Δ₂s[σ], Δ₂s[τ], false)) - end - end - end - return op÷factorial(N-2)^2 -end - -for Elt in [:Sq, :Adj, :Op] - @eval begin - $Elt(RG::GroupRing{AutGroup{N}}) where N = $Elt(RG, N) - $Elt(RG::GroupRing{<:MatAlgebra}) = $Elt(RG, nrows(RG.group)) - end -end - -function SqAdjOp(RG::GroupRing, N::Integer) - S₂ = generating_set(RG.group, 2) - ℤ = Int64 - Δ₂ = length(S₂)*one(RG, ℤ) - RG(S₂, ℤ); - - Alt_N = [σ for σ in SymmetricGroup(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)) + StarAlgebras.add!(op, op, StarAlgebras.mul!(tmp, Δ₂s[σ], Δ₂s[τ])) elseif isadjacent(σ, τ) - GroupRings.addeq!(adj, *(Δ₂s[σ], Δ₂s[τ], false)) + StarAlgebras.add!(adj, adj, StarAlgebras.mul!(tmp, Δ₂s[σ], Δ₂s[τ])) end end end - k = factorial(N-2) - return sq÷k, adj÷k^2, op÷k^2 + k = factorial(n - 2) + return sq ÷ k, adj ÷ k^2, op ÷ k^2 end From 147211ea7af4a79d4ab5729e657b68f3407a4e8f Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 16:13:39 +0100 Subject: [PATCH 06/19] add graded-by-root-system Adj --- src/PropertyT.jl | 4 ++ src/gradings.jl | 80 +++++++++++++++++++++ src/roots.jl | 183 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 267 insertions(+) create mode 100644 src/gradings.jl create mode 100644 src/roots.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 10d74ba..c53e69f 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -19,6 +19,10 @@ include("certify.jl") include("sqadjop.jl") +include("roots.jl") +import .Roots +include("gradings.jl") + include("1712.07167.jl") include("1812.03456.jl") diff --git a/src/gradings.jl b/src/gradings.jl new file mode 100644 index 0000000..d1b2e9b --- /dev/null +++ b/src/gradings.jl @@ -0,0 +1,80 @@ +## something about roots + +Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} = + Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j) + +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 + 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)) + end +end + +function Roots.positive( + generating_set::AbstractVector{<:MatrixGroups.ElementarySymplectic}, +) + r = Roots._positive_direction(Roots.Root(first(generating_set))) + pos_gens = [ + s for s in generating_set if s.val > 0.0 && dot(Roots.Root(s), r) ≥ 0.0 + ] + return pos_gens +end + +grading(s::MatrixGroups.ElementarySymplectic) = Roots.Root(s) +grading(e::MatrixGroups.ElementaryMatrix) = Roots.Root(e) +grading(t::Groups.Transvection) = grading(Groups._abelianize(t)) + +function grading(g::FPGroupElement) + if length(word(g)) == 1 + A = alphabet(parent(g)) + return grading(A[first(word(g))]) + else + throw("Grading is implemented only for generators") + end +end + +_groupby(f, iter::AbstractVector) = _groupby(f.(iter), iter) +function _groupby(keys::AbstractVector{K}, vals::AbstractVector{V}) where {K,V} + @assert length(keys) == length(vals) + d = Dict(k => V[] for k in keys) + for (k, v) in zip(keys, vals) + push!(d[k], v) + end + return d +end + +function laplacians(RG::StarAlgebra, S, grading) + d = _groupby(grading, S) + Δs = Dict(α => RG(length(Sα)) - sum(RG(s) for s in Sα) for (α, Sα) in d) + return Δs +end + +function Adj(rootsystem::AbstractDict, subtype::Symbol) + roots = let W = mapreduce(collect, union, keys(rootsystem)) + W = union!(W, -1 .* W) + end + + return reduce( + +, + ( + Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if + PropertyT_new.Roots.classify_sub_root_system( + roots, + first(α), + first(β), + ) == subtype + ), + init=zero(first(values(rootsystem))), + ) +end + +function level(rootsystem, level::Integer) + 1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}") + level == 1 && return Adj(rootsystem, :C₁) # always positive + level == 2 && return Adj(rootsystem, :A₁) + Adj(rootsystem, Symbol("C₁×C₁")) + Adj(rootsystem, :C₂) # C₂ is not positive + level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁")) + level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive +end diff --git a/src/roots.jl b/src/roots.jl new file mode 100644 index 0000000..080190a --- /dev/null +++ b/src/roots.jl @@ -0,0 +1,183 @@ +module Roots + +using StaticArrays +using LinearAlgebra + +export Root, isproportional, isorthogonal, ~, ⟂ + +abstract type AbstractRoot{N,T} end + +struct Root{N,T} <: AbstractRoot{N,T} + coord::SVector{N,T} +end + +Root(a) = Root(SVector(a...)) + +function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N} + M == N || return false + r.coord == s.coord || return false + return true +end + +Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h)) + +Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord) +Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord) +Base.:-(r::Root{N}) where {N} = Root(-r.coord) + +Base.:*(a::Number, r::Root) = Root(a * r.coord) +Base.:*(r::Root, a::Number) = a * r + +Base.length(r::AbstractRoot) = norm(r, 2) + +LinearAlgebra.norm(r::Root, p::Real=2) = norm(r.coord, p) +LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord) + +cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b)) + +function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M} + N == M || return false + val = abs(cos_angle(α, β)) + return isapprox(val, one(val), atol=eps(one(val))) +end + +function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M} + N == M || return false + val = cos_angle(α, β) + return isapprox(val, zero(val), atol=eps(one(val))) +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)), + ) +end + +function positive(roots::AbstractVector{<:Root{N}}) where {N} + # return those roots for which dot(α, Root([½, ¼, …])) > 0.0 + pd = _positive_direction(first(roots)) + return filter(α -> dot(α, pd) > 0.0, roots) +end + +Base.:~(α::AbstractRoot, β::AbstractRoot) = isproportional(α, β) +⟂(α::AbstractRoot, β::AbstractRoot) = isorthogonal(α, β) + +function Base.show(io::IO, r::Root{N}) where {N} + print(io, "Root$(r.coord)") +end + +function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N} + lngth² = sum(x -> x^2, r.coord) + l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)" + print(io, "Root in ℝ^$N of length $l\n", r.coord) +end + +E(N, i::Integer) = Root(ntuple(k -> k == i ? 1 : 0, N)) +𝕖(N, i) = E(N, i) +𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N)) + +""" + classify_root_system(α, β) +Return the symbol of smallest system generated by roots `α` and `β`. + +The classification is based only on roots length and +proportionality/orthogonality. +""" +function classify_root_system(α::AbstractRoot, β::AbstractRoot) + lα, lβ = length(α), length(β) + if isproportional(α, β) + if lα ≈ lβ ≈ √2 + return :A₁ + elseif lα ≈ lβ ≈ 2.0 + return :C₁ + else + error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") + end + elseif isorthogonal(α, β) + if lα ≈ lβ ≈ √2 + return Symbol("A₁×A₁") + elseif lα ≈ lβ ≈ 2.0 + return Symbol("C₁×C₁") + elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2) + return Symbol("A₁×C₁") + else + error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") + end + else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal + if lα ≈ lβ ≈ √2 + return :A₂ + elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2) + return :C₂ + else + error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") + end + end +end + +function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root) + k = findfirst(v -> isproportional(α, v), Ω) + if isnothing(k) + error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω") + end + return Ω[k] +end + +struct Plane{R<:Root} + v1::R + v2::R + vectors::Vector{R} +end + +Plane(α::R, β::R) where {R<:Root} = + Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3]) + +function Base.in(r::R, plane::Plane{R}) where {R} + return any(isproportional(r, v) for v in plane.vectors) +end + +function classify_sub_root_system( + Ω::AbstractVector{<:Root{N}}, + α::Root{N}, + β::Root{N}, +) where {N} + + v = proportional_root_from_system(Ω, α) + w = proportional_root_from_system(Ω, β) + + subsystem = filter(ω -> ω in Plane(v, w), Ω) + @assert length(subsystem) > 0 + subsystem = positive(union(subsystem, -1 .* subsystem)) + + l = length(subsystem) + if l == 1 + x = first(subsystem) + return classify_root_system(x, x) + elseif l == 2 + return classify_root_system(subsystem...) + elseif l == 3 + a = classify_root_system(subsystem[1], subsystem[2]) + b = classify_root_system(subsystem[2], subsystem[3]) + c = classify_root_system(subsystem[1], subsystem[3]) + + if a == b == c # it's only A₂ + return a + end + + C = (:C₂, Symbol("C₁×C₁")) + if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c)) + return :C₂ + end + elseif l == 4 + for i = 1:l + for j = (i+1):l + T = classify_root_system(subsystem[i], subsystem[j]) + T == :C₂ && return :C₂ + end + end + end + @error "Unknown root subsystem generated by" α β + throw("Unknown root system: $subsystem") +end + +end # of module Roots From 0e5799862b5d586ab2e9e5d19cdb98345f0181dd Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 16:21:58 +0100 Subject: [PATCH 07/19] add action by alphabet permutation --- src/PropertyT.jl | 2 ++ src/actions/actions.jl | 29 ++++++++++++++++++++ src/actions/alphabet_permutation.jl | 42 +++++++++++++++++++++++++++++ src/actions/autfn_conjugation.jl | 26 ++++++++++++++++++ src/actions/sln_conjugation.jl | 20 ++++++++++++++ src/actions/spn_conjugation.jl | 27 +++++++++++++++++++ 6 files changed, 146 insertions(+) create mode 100644 src/actions/actions.jl create mode 100644 src/actions/alphabet_permutation.jl create mode 100644 src/actions/autfn_conjugation.jl create mode 100644 src/actions/sln_conjugation.jl create mode 100644 src/actions/spn_conjugation.jl diff --git a/src/PropertyT.jl b/src/PropertyT.jl index c53e69f..02d2536 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -23,6 +23,8 @@ include("roots.jl") import .Roots include("gradings.jl") +include("actions/actions.jl") + include("1712.07167.jl") include("1812.03456.jl") diff --git a/src/actions/actions.jl b/src/actions/actions.jl new file mode 100644 index 0000000..29a2fb3 --- /dev/null +++ b/src/actions/actions.jl @@ -0,0 +1,29 @@ +import SymbolicWedderburn.action +StarAlgebras.star(g::GroupElement) = inv(g) + +include("alphabet_permutation.jl") + +include("sln_conjugation.jl") +include("autfn_conjugation.jl") + +function SymbolicWedderburn.action( + act::SymbolicWedderburn.ByPermutations, + g::Groups.GroupElement, + α::StarAlgebras.AlgebraElement +) + res = StarAlgebras.zero!(similar(α)) + B = basis(parent(α)) + for (idx, val) in StarAlgebras._nzpairs(StarAlgebras.coeffs(α)) + a = B[idx] + a_g = SymbolicWedderburn.action(act, g, a) + res[a_g] += val + end + return res +end + +function Base.:^( + w::W, + p::PermutationGroups.AbstractPerm, +) where {W<:Groups.AbstractWord} + return W([l^p for l in w]) +end diff --git a/src/actions/alphabet_permutation.jl b/src/actions/alphabet_permutation.jl new file mode 100644 index 0000000..2f2655c --- /dev/null +++ b/src/actions/alphabet_permutation.jl @@ -0,0 +1,42 @@ +## action induced from permuting letters of an alphabet + +struct AlphabetPermutation{GEl,I} <: SymbolicWedderburn.ByPermutations + perms::Dict{GEl,Perm{I}} +end + +function AlphabetPermutation( + A::Alphabet, + Γ::PermutationGroups.AbstractPermutationGroup, + op, +) + return AlphabetPermutation( + Dict(γ => inv(Perm([A[op(l, γ)] for l in A])) for γ in Γ), + ) +end + +function AlphabetPermutation(A::Alphabet, W::Constructions.WreathProduct, op) + return AlphabetPermutation( + Dict( + w => inv(Perm([A[op(op(l, w.p), w.n)] for l in A])) for + w in W + ), + ) +end + +function SymbolicWedderburn.action( + act::AlphabetPermutation, + γ::GroupElement, + w::Groups.AbstractWord, +) + return w^(act.perms[γ]) +end + +function SymbolicWedderburn.action( + act::AlphabetPermutation, + γ::GroupElement, + g::Groups.AbstractFPGroupElement, +) + G = parent(g) + w = word(g)^(act.perms[γ]) + return G(w) +end diff --git a/src/actions/autfn_conjugation.jl b/src/actions/autfn_conjugation.jl new file mode 100644 index 0000000..99559d5 --- /dev/null +++ b/src/actions/autfn_conjugation.jl @@ -0,0 +1,26 @@ +## Particular definitions for actions on Aut(F_n) + +function _conj( + t::Groups.Transvection, + σ::PermutationGroups.AbstractPerm, +) + return Groups.Transvection(t.id, t.i^inv(σ), t.j^inv(σ), t.inv) +end + +function _flip(t::Groups.Transvection, g::Groups.GroupElement) + isone(g) && return t + return Groups.Transvection(t.id === :ϱ ? :λ : :ϱ, t.i, t.j, t.inv) +end + +function _conj( + t::Groups.Transvection, + x::Groups.Constructions.DirectPowerElement, +) + @assert Groups.order(Int, parent(x).group) == 2 + i, j = t.i, t.j + t = ifelse(isone(x.elts[i] * x.elts[j]), t, inv(t)) + return _flip(t, x.elts[i]) +end + +action_by_conjugation(sautfn::Groups.AutomorphismGroup{<:Groups.FreeGroup}, Σ::Groups.Group) = + AlphabetPermutation(alphabet(sautfn), Σ, _conj) diff --git a/src/actions/sln_conjugation.jl b/src/actions/sln_conjugation.jl new file mode 100644 index 0000000..bdb016f --- /dev/null +++ b/src/actions/sln_conjugation.jl @@ -0,0 +1,20 @@ +## Particular definitions for actions on SL(n,ℤ) + +function _conj( + t::MatrixGroups.ElementaryMatrix{N}, + σ::PermutationGroups.AbstractPerm, +) where {N} + return MatrixGroups.ElementaryMatrix{N}(t.i^inv(σ), t.j^inv(σ), t.val) +end + +function _conj( + t::MatrixGroups.ElementaryMatrix{N}, + x::Groups.Constructions.DirectPowerElement, +) where {N} + @assert Groups.order(Int, parent(x).group) == 2 + just_one_flips = xor(isone(x.elts[t.i]), isone(x.elts[t.j])) + return ifelse(just_one_flips, inv(t), t) +end + +action_by_conjugation(sln::Groups.MatrixGroups.SpecialLinearGroup, Σ::Groups.Group) = + AlphabetPermutation(alphabet(sln), Σ, _conj) diff --git a/src/actions/spn_conjugation.jl b/src/actions/spn_conjugation.jl new file mode 100644 index 0000000..25b8bbc --- /dev/null +++ b/src/actions/spn_conjugation.jl @@ -0,0 +1,27 @@ +## Particular definitions for actions on Sp(n,ℤ) + +function _conj( + t::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) +end + +function _conj( + t::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) +end + +action_by_conjugation(sln::Groups.MatrixGroups.SymplecticGroup, Σ::Groups.Group) = + AlphabetPermutation(alphabet(sln), Σ, _conj) From f00bfb7ca909adffa472862153b6d219afece131 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 16:29:26 +0100 Subject: [PATCH 08/19] remove old laplacians.jl, 1712.* 1812.* --- src/1712.07167.jl | 329 ---------------------------------------------- src/1812.03456.jl | 26 ---- src/PropertyT.jl | 1 - src/laplacians.jl | 51 ------- 4 files changed, 407 deletions(-) delete mode 100644 src/1712.07167.jl delete mode 100644 src/1812.03456.jl delete mode 100644 src/laplacians.jl diff --git a/src/1712.07167.jl b/src/1712.07167.jl deleted file mode 100644 index 7e2cdba..0000000 --- a/src/1712.07167.jl +++ /dev/null @@ -1,329 +0,0 @@ -############################################################################### -# Settings and filenames - -abstract type Settings end - -struct Naive{El} <: Settings - name::String - G::Union{Group, NCRing} - S::Vector{El} - halfradius::Int - upper_bound::Float64 - - solver::JuMP.OptimizerFactory - force_compute::Bool -end - -struct Symmetrized{El} <: Settings - name::String - G::Union{Group, NCRing} - S::Vector{El} - autS::Group - halfradius::Int - upper_bound::Float64 - - solver::JuMP.OptimizerFactory - force_compute::Bool -end - -function Settings(name::String, - G::Union{Group, NCRing}, S::AbstractVector{El}, solver::JuMP.OptimizerFactory; - halfradius=2, upper_bound=1.0, force_compute=false) where El <: Union{GroupElem, NCRingElem} - return Naive(name, G, S, halfradius, upper_bound, solver, force_compute) -end - -function Settings(name::String, - G::Union{Group, NCRing}, S::AbstractVector{El}, autS::Group, solver::JuMP.OptimizerFactory; - halfradius=2, upper_bound=1.0, force_compute=false) where El <: Union{GroupElem, NCRingElem} - return Symmetrized(name, G, S, autS, halfradius, upper_bound, solver, force_compute) -end - -suffix(s::Settings) = "$(s.upper_bound)" -prepath(s::Settings) = s.name -fullpath(s::Settings) = joinpath(prepath(s), suffix(s)) - -filename(sett::Settings, s::Symbol; kwargs...) = filename(sett, Val{s}; kwargs...) - -filename(sett::Settings, ::Type{Val{:fulllog}}; kwargs...) = - filename(fullpath(sett), "full", "log", suffix=Dates.now(); kwargs...) -filename(sett::Settings, ::Type{Val{:solverlog}}; kwargs...) = - filename(fullpath(sett), "solver", "log", suffix=Dates.now(); kwargs...) - -filename(sett::Settings, ::Type{Val{:Δ}}; kwargs...) = - filename(prepath(sett), "delta", "jld"; kwargs...) -filename(sett::Settings, ::Type{Val{:BlockDecomposition}}; kwargs...) = - filename(prepath(sett), "BlockDecomposition", "jld"; kwargs...) - -filename(sett::Settings, ::Type{Val{:solution}}; kwargs...) = - filename(fullpath(sett), "solution", "jld"; kwargs...) - -function filename(sett::Settings, ::Type{Val{:warmstart}}; kwargs...) - filename(fullpath(sett), "warmstart", "jld"; kwargs...) -end - -function filename(path::String, name, extension; prefix=nothing, suffix=nothing) - pre = isnothing(prefix) ? "" : "$(prefix)_" - suf = isnothing(suffix) ? "" : "_$(suffix)" - return joinpath(path, "$pre$name$suf.$extension") -end - -############################################################################### -# Approximation by SOS (logged & warmstarted) - -function warmstart(sett::Settings) - warmstart_fname = filename(sett, :warmstart) - try - ws = load(warmstart_fname, "warmstart") - @info "Loaded $warmstart_fname." - return ws - catch ex - @warn "$(ex.msg). Could not provide a warmstart to the solver." - return nothing - end -end - -function approximate_by_SOS(sett::Naive, - elt::GroupRingElem, orderunit::GroupRingElem; - solverlog=tempname()*".log") - - isdir(fullpath(sett)) || mkpath(fullpath(sett)) - - @info "Creating SDP problem..." - SDP_problem = SOS_problem_primal(elt, orderunit, upper_bound=sett.upper_bound) - @info Base.repr(SDP_problem) - - @info "Logging solver's progress into $solverlog" - - ws = warmstart(sett) - @time status, ws = PropertyT.solve(solverlog, SDP_problem, sett.solver, ws) - @info "Optimization finished:" status - - P = value.(SDP_problem[:P]) - λ = value(SDP_problem[:λ]) - - if any(isnan, P) - @warn "The solution seems to contain NaNs. Not overriding warmstart.jld" - else - save(filename(sett, :warmstart), - "warmstart", (ws.primal, ws.dual, ws.slack), - "P", P, - "λ", λ) - end - - save(filename(sett, :warmstart, suffix=Dates.now()), - "warmstart", (ws.primal, ws.dual, ws.slack), - "P", P, - "λ", λ) - - return λ, P -end - -function approximate_by_SOS(sett::Symmetrized, - elt::GroupRingElem, orderunit::GroupRingElem; - solverlog=tempname()*".log") - - isdir(fullpath(sett)) || mkpath(fullpath(sett)) - - orbit_data = try - orbit_data = load(filename(sett, :BlockDecomposition), "BlockDecomposition") - @info "Loaded orbit data." - orbit_data - catch ex - @warn ex.msg - GroupRings.hasbasis(parent(orderunit)) || - throw("You need to define basis of Group Ring to compute orbit decomposition!") - @info "Computing orbit and Wedderburn decomposition..." - orbit_data = BlockDecomposition(parent(orderunit), sett.autS) - save(filename(sett, :BlockDecomposition), "BlockDecomposition", orbit_data) - orbit_data - end - - orbit_data = decimate(orbit_data) - - @info "Creating SDP problem..." - SDP_problem, varP = SOS_problem_primal(elt, orderunit, orbit_data, upper_bound=sett.upper_bound) - @info Base.repr(SDP_problem) - - @info "Logging solver's progress into $solverlog" - - ws = warmstart(sett) - @time status, ws = PropertyT.solve(solverlog, SDP_problem, sett.solver, ws) - @info "Optimization finished:" status - - λ = value(SDP_problem[:λ]) - Ps = [value.(P) for P in varP] - - if any(any(isnan, P) for P in Ps) - @warn "The solution seems to contain NaNs. Not overriding warmstart.jld" - else - save(filename(sett, :warmstart), - "warmstart", (ws.primal, ws.dual, ws.slack), - "Ps", Ps, - "λ", λ) - end - - save(filename(sett, :warmstart, suffix=Dates.now()), - "warmstart", (ws.primal, ws.dual, ws.slack), - "Ps", Ps, - "λ", λ) - - @info "Reconstructing P..." - @time P = reconstruct(Ps, orbit_data) - - return λ, P -end - -############################################################################### -# Checking solution - -function certify_SOS_decomposition(elt::GroupRingElem, orderunit::GroupRingElem, - λ::Number, Q::AbstractMatrix; R::Int=2) - separator = "-"^76 - @info "$separator\nChecking in floating-point arithmetic..." λ - eoi = elt - λ*orderunit - - @info("Computing sum of squares decomposition...") - @time residual = eoi - compute_SOS(parent(eoi), augIdproj(Q)) - - L1_norm = norm(residual,1) - floatingpoint_λ = λ - 2.0^(2ceil(log2(R)))*L1_norm - - info_strs = ["Numerical metrics of the obtained SOS:", - "ɛ(elt - λu - ∑ξᵢ*ξᵢ) ≈ $(aug(residual))", - "‖elt - λu - ∑ξᵢ*ξᵢ‖₁ ≈ $(L1_norm)", - "Floating point (NOT certified) λ ≈"] - @info join(info_strs, "\n") floatingpoint_λ - - if floatingpoint_λ ≤ 0 - return floatingpoint_λ - end - - λ = @interval(λ) - info_strs = [separator, - "Checking in interval arithmetic...", - "λ ∈ $λ"] - @info(join(info_strs, "\n")) - eoi = elt - λ*orderunit - - @info("Projecting columns of Q to the augmentation ideal...") - @time Q, check = augIdproj(Interval, Q) - @info "Checking that sum of every column contains 0.0..." check_augmented=check - check || @error("The following numbers are meaningless!") - - @info("Computing sum of squares decomposition...") - @time residual = eoi - compute_SOS(parent(eoi), Q) - - L1_norm = norm(residual,1) - certified_λ = λ - 2.0^(2ceil(log2(R)))*L1_norm - - info_strs = ["Numerical metrics of the obtained SOS:", - "ɛ(elt - λu - ∑ξᵢ*ξᵢ) ∈ $(aug(residual))", - "‖elt - λu - ∑ξᵢ*ξᵢ‖₁ ∈ $(L1_norm)", - "Interval aritmetic (certified) λ ∈"] - @info join(info_strs, "\n") certified_λ - - return certified_λ -end - -function spectral_gap(Δ::GroupRingElem, λ::Number, Q::AbstractMatrix; R::Int=2) - @info "elt = Δ², u = Δ" - return certify_SOS_decomposition(Δ^2, Δ, λ, Q, R=R) -end - -############################################################################### -# Interpreting the numerical results - -Kazhdan_constant(λ::Number, N::Integer) = sqrt(2*λ/N) -Kazhdan_constant(λ::Interval, N::Integer) = IntervalArithmetic.inf(sqrt(2*λ/N)) - -function check_property_T(sett::Settings) - @info sett - certified_sgap = spectral_gap(sett) - return interpret_results(sett, certified_sgap) -end - -function Base.show(io::IO, sett::Settings) - info_strs = ["PropertyT Settings:", - "Group: $(sett.name)", - "Upper bound for λ: $(sett.upper_bound), on halfradius $(sett.halfradius).", - "Force computations: $(sett.force_compute);", - "Results will be stored in ./$(PropertyT.prepath(sett));", - "Solver: $(typeof(sett.solver()))", - "Solvers options: "] - append!(info_strs, [rpad(" $k", 30)* "→ \t$v" for (k,v) in sett.solver().options]) - push!(info_strs, "="^76) - print(io, join(info_strs, "\n")) -end - -function interpret_results(name::String, sgap::Number, N::Integer) - if sgap > 0 - κ = Kazhdan_constant(sgap, N) - @info "κ($name, S) ≥ $κ: Group HAS property (T)!" - return true - end - info_strs = [ - "The certified lower bound on the spectral gap is negative:", - "λ($name, S) ≥ 0.0 > $sgap", - "This tells us nothing about property (T)", - ] - @info join(info_strs, "\n") - return false -end - -interpret_results(sett::Settings, sgap::Number) = - interpret_results(sett.name, sgap, length(sett.S)) - -function spectral_gap(sett::Settings) - fp = PropertyT.fullpath(sett) - isdir(fp) || mkpath(fp) - - Δ = try - Δ = loadGRElem(filename(sett,:Δ), sett.G) - @info "Loaded precomputed Δ." - Δ - catch ex - @warn ex.msg - @info "Computing Δ..." - Δ = Laplacian(sett.S, sett.halfradius) - saveGRElem(filename(sett, :Δ), Δ) - Δ - end - - 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!" - return λ, P - end - - if sett.force_compute - λ, 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:", - "sum(P) = $(sum(P))", - "maximum(P) = $(maximum(P))", - "minimum(P) = $(minimum(P))"] - @info join(info_strs, "\n") - - isapprox(eigvals(P), abs.(eigvals(P))) || - @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.halfradius) - - return certified_sgap -end diff --git a/src/1812.03456.jl b/src/1812.03456.jl deleted file mode 100644 index ff3a71b..0000000 --- a/src/1812.03456.jl +++ /dev/null @@ -1,26 +0,0 @@ -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.transvection_R(i,j) for (i,j) in indexing(n)] - lmuls = [Groups.transvection_L(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 02d2536..1a7ae44 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -12,7 +12,6 @@ using Groups using StarAlgebras using SymbolicWedderburn -include("laplacians.jl") include("constraint_matrix.jl") include("sos_sdps.jl") include("certify.jl") diff --git a/src/laplacians.jl b/src/laplacians.jl deleted file mode 100644 index 7138f8f..0000000 --- a/src/laplacians.jl +++ /dev/null @@ -1,51 +0,0 @@ -############################################################################### -# -# Laplacians -# -############################################################################### - -function spLaplacian(RG::GroupRing, S::AbstractVector, T::Type=Float64) - 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::AbstractVector{REl}, halfradius) where REl<:Union{NCRingElem, GroupElem} - G = parent(first(S)) - @info "Generating metric ball of radius" radius=2halfradius - @time E_R, sizes = Groups.wlmetric_ball(S, 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(G, E_R, rdict, pm) - Δ = spLaplacian(RG, S) - return Δ -end - -function saveGRElem(fname::String, g::GroupRingElem) - RG = parent(g) - JLD.save(fname, "coeffs", g.coeffs, "pm", RG.pm, "G", RG.group) -end - -function loadGRElem(fname::String, RG::GroupRing) - coeffs = load(fname, "coeffs") - return GroupRingElem(coeffs, RG) -end - -function loadGRElem(fname::String, G::Union{Group, NCRing}) - pm = load(fname, "pm") - RG = GroupRing(G, pm) - return loadGRElem(fname, RG) -end - -function loadGRElem(fname::String) - pm, G = load(fname, "pm", "G") - RG = GroupRing(G, pm) - return loadGRElem(fname, RG) -end From 633f06548873c3fe88feecf274997eb5c6c16772 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 17:01:06 +0100 Subject: [PATCH 09/19] cleanup Project.toml and fix imports --- Project.toml | 23 ++++++++++++------ src/PropertyT.jl | 5 ++-- src/actions/actions.jl | 2 +- src/actions/alphabet_permutation.jl | 12 ++++++---- src/certify.jl | 28 +++++++++++----------- src/gradings.jl | 4 ++-- src/sos_sdps.jl | 37 ++++++++++++++++------------- src/sqadjop.jl | 10 ++++---- 8 files changed, 68 insertions(+), 53 deletions(-) diff --git a/Project.toml b/Project.toml index 155508e..a14a559 100644 --- a/Project.toml +++ b/Project.toml @@ -4,23 +4,32 @@ authors = ["Marek Kaluba "] version = "0.3.2" [deps] -Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2" [compat] -IntervalArithmetic = "^0.16.0" -JuMP = "^0.20.0" -SCS = "^0.7.0" -julia = "^1.3.0" +COSMO = "0.8" +Groups = "0.7" +IntervalArithmetic = "0.20" +JuMP = "1.3" +SCS = "1.1.0" +StaticArrays = "1" +SymbolicWedderburn = "0.3.1" +julia = "1.6" [extras] +COSMO = "1e616198-aa4e-51ec-90a2-23f7fbd31d8d" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" +Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "SCS"] +scripts = ["Dates", "Logging", "Serialization", "SCS", "COSMO"] +test = ["Test", "SCS", "COSMO"] diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 1a7ae44..838f40f 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -3,14 +3,15 @@ module PropertyT using LinearAlgebra using SparseArrays -using Dates using IntervalArithmetic using JuMP using Groups -using StarAlgebras +import Groups.GroupsCore using SymbolicWedderburn +import SymbolicWedderburn.StarAlgebras +import SymbolicWedderburn.PermutationGroups include("constraint_matrix.jl") include("sos_sdps.jl") diff --git a/src/actions/actions.jl b/src/actions/actions.jl index 29a2fb3..aa02b15 100644 --- a/src/actions/actions.jl +++ b/src/actions/actions.jl @@ -1,5 +1,5 @@ import SymbolicWedderburn.action -StarAlgebras.star(g::GroupElement) = inv(g) +StarAlgebras.star(g::Groups.GroupElement) = inv(g) include("alphabet_permutation.jl") diff --git a/src/actions/alphabet_permutation.jl b/src/actions/alphabet_permutation.jl index 2f2655c..a3e15e0 100644 --- a/src/actions/alphabet_permutation.jl +++ b/src/actions/alphabet_permutation.jl @@ -1,7 +1,9 @@ ## action induced from permuting letters of an alphabet +import Groups: Constructions + struct AlphabetPermutation{GEl,I} <: SymbolicWedderburn.ByPermutations - perms::Dict{GEl,Perm{I}} + perms::Dict{GEl,PermutationGroups.Perm{I}} end function AlphabetPermutation( @@ -10,14 +12,14 @@ function AlphabetPermutation( op, ) return AlphabetPermutation( - Dict(γ => inv(Perm([A[op(l, γ)] for l in A])) for γ in Γ), + Dict(γ => inv(PermutationGroups.Perm([A[op(l, γ)] for l in A])) for γ in Γ), ) end function AlphabetPermutation(A::Alphabet, W::Constructions.WreathProduct, op) return AlphabetPermutation( Dict( - w => inv(Perm([A[op(op(l, w.p), w.n)] for l in A])) for + w => inv(PermutationGroups.Perm([A[op(op(l, w.p), w.n)] for l in A])) for w in W ), ) @@ -25,7 +27,7 @@ end function SymbolicWedderburn.action( act::AlphabetPermutation, - γ::GroupElement, + γ::Groups.GroupElement, w::Groups.AbstractWord, ) return w^(act.perms[γ]) @@ -33,7 +35,7 @@ end function SymbolicWedderburn.action( act::AlphabetPermutation, - γ::GroupElement, + γ::Groups.GroupElement, g::Groups.AbstractFPGroupElement, ) G = parent(g) diff --git a/src/certify.jl b/src/certify.jl index bca18c6..089621d 100644 --- a/src/certify.jl +++ b/src/certify.jl @@ -34,7 +34,7 @@ function _fma_SOS_thr!( return result end -function _cnstr_sos!(res::AlgebraElement, Q::AbstractMatrix, cnstrs) +function _cnstr_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix, cnstrs) StarAlgebras.zero!(res) Q² = Q' * Q for (g, A_g) in cnstrs @@ -43,7 +43,7 @@ function _cnstr_sos!(res::AlgebraElement, Q::AbstractMatrix, cnstrs) return res end -function _augmented_sos!(res::AlgebraElement, Q::AbstractMatrix) +function _augmented_sos!(res::StarAlgebras.AlgebraElement, Q::AbstractMatrix) A = parent(res) StarAlgebras.zero!(res) Q² = Q' * Q @@ -64,10 +64,10 @@ function _augmented_sos!(res::AlgebraElement, Q::AbstractMatrix) return res end -function compute_sos(A::StarAlgebra, Q::AbstractMatrix; augmented::Bool) +function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool) if augmented z = zeros(eltype(Q), length(basis(A))) - res = AlgebraElement(z, 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) @@ -77,11 +77,11 @@ function compute_sos(A::StarAlgebra, Q::AbstractMatrix; augmented::Bool) _fma_SOS_thr!(z, A.mstructure, Q) - return AlgebraElement(z, A) + return StarAlgebras.AlgebraElement(z, A) end end -function sufficient_λ(residual::AlgebraElement, λ; halfradius) +function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius) L1_norm = norm(residual, 1) suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm @@ -97,7 +97,7 @@ function sufficient_λ(residual::AlgebraElement, λ; halfradius) info_strs = [ "Numerical metrics of the obtained SOS:", - "ɛ(elt - λu - ∑ξᵢ*ξᵢ) $eq_sign $(aug(residual))", + "ɛ(elt - λu - ∑ξᵢ*ξᵢ) $eq_sign $(StarAlgebras.aug(residual))", "‖elt - λu - ∑ξᵢ*ξᵢ‖₁ $eq_sign $(L1_norm)", " λ $eq_sign $suff_λ", ] @@ -107,10 +107,10 @@ function sufficient_λ(residual::AlgebraElement, λ; halfradius) end function sufficient_λ( - elt::AlgebraElement, - order_unit::AlgebraElement, + elt::StarAlgebras.AlgebraElement, + order_unit::StarAlgebras.AlgebraElement, λ, - sos::AlgebraElement; + sos::StarAlgebras.AlgebraElement; halfradius ) @@ -121,15 +121,15 @@ function sufficient_λ( end function certify_solution( - elt::AlgebraElement, - orderunit::AlgebraElement, + elt::StarAlgebras.AlgebraElement, + orderunit::StarAlgebras.AlgebraElement, λ, Q::AbstractMatrix{<:AbstractFloat}; halfradius, - augmented=iszero(aug(elt)) && iszero(aug(orderunit)) + augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)) ) - should_we_augment = !augmented && aug(elt) == aug(orderunit) == 0 + should_we_augment = !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0 Q = should_we_augment ? augment_columns!(Q) : Q @time sos = compute_sos(parent(elt), Q, augmented=augmented) diff --git a/src/gradings.jl b/src/gradings.jl index d1b2e9b..20edd29 100644 --- a/src/gradings.jl +++ b/src/gradings.jl @@ -46,7 +46,7 @@ function _groupby(keys::AbstractVector{K}, vals::AbstractVector{V}) where {K,V} return d end -function laplacians(RG::StarAlgebra, S, grading) +function laplacians(RG::StarAlgebras.StarAlgebra, S, grading) d = _groupby(grading, S) Δs = Dict(α => RG(length(Sα)) - sum(RG(s) for s in Sα) for (α, Sα) in d) return Δs @@ -61,7 +61,7 @@ function Adj(rootsystem::AbstractDict, subtype::Symbol) +, ( Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if - PropertyT_new.Roots.classify_sub_root_system( + Roots.classify_sub_root_system( roots, first(α), first(β), diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index 5c77de4..250388c 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -5,8 +5,8 @@ Formulate the dual to the sum of squares decomposition problem for `X - λ·u`. See also [sos_problem_primal](@ref). """ function sos_problem_dual( - elt::AlgebraElement, - order_unit::AlgebraElement=zero(elt); + elt::StarAlgebras.AlgebraElement, + order_unit::StarAlgebras.AlgebraElement=zero(elt); lower_bound=-Inf ) @assert parent(elt) == parent(order_unit) @@ -70,7 +70,7 @@ function constraints( a, b = basis[i], basis[j] push!(cnstrs[basis[one(a)]], k) - push!(cnstrs[basis[star(a)]], -k) + push!(cnstrs[basis[StarAlgebras.star(a)]], -k) push!(cnstrs[basis[b]], -k) end end @@ -80,7 +80,7 @@ function constraints( ) end -function constraints(A::StarAlgebra; augmented::Bool, twisted::Bool) +function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool, twisted::Bool) mstructure = if StarAlgebras._istwisted(A.mstructure) == twisted A.mstructure else @@ -108,10 +108,10 @@ be added to the model. This may improve the accuracy of the solution if The default `u = zero(X)` formulates a simple feasibility problem. """ function sos_problem_primal( - elt::AlgebraElement, - order_unit::AlgebraElement=zero(elt); + elt::StarAlgebras.AlgebraElement, + order_unit::StarAlgebras.AlgebraElement=zero(elt); upper_bound=Inf, - augmented::Bool=iszero(aug(elt)) && iszero(aug(order_unit)) + augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit)) ) @assert parent(elt) === parent(order_unit) @@ -168,22 +168,25 @@ function isorth_projection(ds::SymbolicWedderburn.DirectSummand) end sos_problem_primal( - elt::AlgebraElement, + elt::StarAlgebras.AlgebraElement, wedderburn::WedderburnDecomposition; kwargs... ) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...) function sos_problem_primal( - elt::AlgebraElement, - orderunit::AlgebraElement, + elt::StarAlgebras.AlgebraElement, + orderunit::StarAlgebras.AlgebraElement, wedderburn::WedderburnDecomposition; upper_bound=Inf, - augmented=iszero(aug(elt)) && iszero(aug(orderunit)) + augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)), + check_orthogonality=true ) @assert parent(elt) === parent(orderunit) - if any(!isorth_projection, direct_summands(wedderburn)) - error("Wedderburn decomposition contains a non-orthogonal projection") + if check_orthogonality + if any(!isorth_projection, direct_summands(wedderburn)) + error("Wedderburn decomposition contains a non-orthogonal projection") + end end feasibility_problem = iszero(orderunit) @@ -215,8 +218,8 @@ function sos_problem_primal( tmps = SymbolicWedderburn._tmps(wedderburn) end - X = convert(Vector{T}, coeffs(elt)) - U = convert(Vector{T}, coeffs(orderunit)) + X = convert(Vector{T}, StarAlgebras.coeffs(elt)) + U = convert(Vector{T}, StarAlgebras.coeffs(orderunit)) # defining constraints based on the multiplicative structure cnstrs = constraints(parent(elt), augmented=augmented, twisted=true) @@ -278,7 +281,7 @@ function reconstruct!( if eltype(res) <: AbstractFloat SymbolicWedderburn.zerotol!(tmp2, atol=1e-12) end - tmp2 .*= degree(ds) + tmp2 .*= SymbolicWedderburn.degree(ds) @assert size(tmp2) == size(res) @@ -291,7 +294,7 @@ function reconstruct!( end end end - res ./= order(Int, G) + res ./= Groups.order(Int, G) return res end diff --git a/src/sqadjop.jl b/src/sqadjop.jl index 2915660..9df0a06 100644 --- a/src/sqadjop.jl +++ b/src/sqadjop.jl @@ -1,4 +1,4 @@ -import PermutationGroups.AbstractPerm +import SymbolicWedderburn.PermutationGroups.AbstractPerm # move to Groups Base.keys(a::Alphabet) = keys(1:length(a)) @@ -15,7 +15,7 @@ isadjacent(σ::AbstractPerm, τ::AbstractPerm, i=1, j=2) = (j^σ == i^τ && i^σ ≠ j^τ) # second σ equal to first τ function _ncycle(start, length, n=start + length - 1) - p = Perm(Int8(n)) + p = PermutationGroups.Perm(Int8(n)) @assert n ≥ start + length - 1 for k in start:start+length-2 p[k] = k + 1 @@ -24,7 +24,7 @@ function _ncycle(start, length, n=start + length - 1) return p end -alternating_group(n::Integer) = PermGroup([_ncycle(i, 3) for i in 1:n-2]) +alternating_group(n::Integer) = PermutationGroups.PermGroup([_ncycle(i, 3) for i in 1:n-2]) function small_gens(G::MatrixGroups.SpecialLinearGroup) A = alphabet(G) @@ -46,13 +46,13 @@ function small_gens(G::Groups.AutomorphismGroup{<:FreeGroup}) return union!(S, inv.(S)) end -function small_laplacian(RG::StarAlgebra) +function small_laplacian(RG::StarAlgebras.StarAlgebra) G = StarAlgebras.object(RG) S₂ = small_gens(G) return length(S₂) * one(RG) - sum(RG(s) for s in S₂) end -function SqAdjOp(A::StarAlgebra, n::Integer, Δ₂=small_laplacian(A)) +function SqAdjOp(A::StarAlgebras.StarAlgebra, n::Integer, Δ₂=small_laplacian(A)) @assert parent(Δ₂) === A alt_n = alternating_group(n) From 7907137fb5b6ca3c3ca5ff572f16331102baf6f2 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 18:44:42 +0100 Subject: [PATCH 10/19] define group_algebra --- src/PropertyT.jl | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/PropertyT.jl b/src/PropertyT.jl index 838f40f..c9adba4 100644 --- a/src/PropertyT.jl +++ b/src/PropertyT.jl @@ -25,7 +25,19 @@ include("gradings.jl") include("actions/actions.jl") -include("1712.07167.jl") -include("1812.03456.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) + @info "sizes = $(sizes)" + @info "computing the *-algebra structure for G" + @time RG = StarAlgebras.StarAlgebra{twisted}( + G, + StarAlgebras.Basis{UInt32}(E), + (sizes[halfradius], sizes[halfradius]), + precompute=false, + ) + return RG, S, sizes +end end # module Property(T) From a080ae74c311209521ab604dc7a4a375c26fd4d3 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 18:45:12 +0100 Subject: [PATCH 11/19] reorganize/rewrite tests! --- test/1703.09680.jl | 206 ++++++++++++++++++++++---- test/1712.07167.jl | 307 +++++++++++++++++++++++++++------------ test/1812.03456.jl | 355 ++++++++++++++++++++++----------------------- test/actions.jl | 261 +++++++++++++++++++++++---------- test/graded_adj.jl | 40 +++++ test/optimizers.jl | 54 +++++++ test/runtests.jl | 35 +++-- 7 files changed, 871 insertions(+), 387 deletions(-) create mode 100644 test/graded_adj.jl create mode 100644 test/optimizers.jl diff --git a/test/1703.09680.jl b/test/1703.09680.jl index cbed108..75e76c1 100644 --- a/test/1703.09680.jl +++ b/test/1703.09680.jl @@ -1,51 +1,203 @@ +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 N = 2 - G = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(G) + G = MatrixGroups.SpecialLinearGroup{N}(Int8) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) - rm("SL($N,Z)", recursive=true, force=true) - sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(20000, accel=20); upper_bound=0.1) + Δ = let RG = RG, S = S + RG(length(S)) - sum(RG(s) for s in S) + end - @info sett + elt = Δ^2 + unit = Δ + ub = 0.1 - λ = PropertyT.spectral_gap(sett) - @test λ < 0.0 - @test PropertyT.interpret_results(sett, λ) == false + status, certified, λ = check_positivity( + elt, + unit, + upper_bound=ub, + halfradius=2, + optimizer=scs_optimizer( + eps=1e-10, + max_iters=5_000, + accel=50, + alpha=1.9, + ) + ) + + @test status == JuMP.ALMOST_OPTIMAL + @test !certified + @test λ < 0 end - @testset "SL(3,Z)" begin + @testset "SL(3,F₅)" begin N = 3 - G = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(G) + G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{5}) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) - rm("SL($N,Z)", recursive=true, force=true) - sett = PropertyT.Settings("SL($N,Z)", G, S, with_SCS(1000, accel=20); upper_bound=0.1) + Δ = let RG = RG, S = S + RG(length(S)) - sum(RG(s) for s in S) + end - @info sett + elt = Δ^2 + unit = Δ + ub = 1.01 # 1.5 - λ = PropertyT.spectral_gap(sett) - @test λ > 0.099 - @test PropertyT.interpret_results(sett, λ) == true + status, certified, λ = check_positivity( + elt, + unit, + upper_bound=ub, + halfradius=2, + optimizer=scs_optimizer( + eps=1e-10, + max_iters=5_000, + accel=50, + alpha=1.9, + ) + ) - @test PropertyT.check_property_T(sett) == true #second run should be fast + @test status == JuMP.OPTIMAL + @test certified + @test λ > 1 end @testset "SAut(F₂)" begin N = 2 - G = SAut(FreeGroup(N)) - S = PropertyT.generating_set(G) + G = SpecialAutomorphismGroup(FreeGroup(N)) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) - rm("SAut(F$N)", recursive=true, force=true) - sett = PropertyT.Settings("SAut(F$N)", G, S, with_SCS(10000); - upper_bound=0.15) + Δ = let RG = RG, S = S + RG(length(S)) - sum(RG(s) for s in S) + end - @info sett + elt = Δ^2 + unit = Δ + ub = 0.1 - λ = PropertyT.spectral_gap(sett) - @test λ < 0.0 - @test PropertyT.interpret_results(sett, λ) == false + status, certified, λ = check_positivity( + elt, + unit, + upper_bound=0.1, + halfradius=2, + optimizer=scs_optimizer( + eps=1e-10, + max_iters=5_000, + accel=50, + alpha=1.9, + ) + ) + @test status == JuMP.ALMOST_OPTIMAL + @test λ < 0 + @test !certified + end + + @testset "SL(3,Z) has (T)" begin + n = 3 + + SL = MatrixGroups.SpecialLinearGroup{n}(Int8) + RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true) + + Δ = RSL(length(S)) - sum(RSL(s) for s in S) + + @testset "basic formulation" begin + elt = Δ^2 + unit = Δ + ub = 0.1 + + opt_problem = PropertyT.sos_problem_primal( + elt, + unit, + upper_bound=ub, + augmented=false, + ) + + status, _ = PropertyT.solve( + opt_problem, + cosmo_optimizer( + eps=1e-10, + max_iters=10_000, + accel=0, + alpha=1.5, + ), + ) + + @test status == JuMP.OPTIMAL + + λ = JuMP.value(opt_problem[:λ]) + @test λ > 0.09 + Q = real.(sqrt(JuMP.value.(opt_problem[:P]))) + + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + λ, + Q, + halfradius=2, + augmented=false, + ) + + @test certified + @test isapprox(λ_cert, λ, rtol=1e-5) + end + + @testset "augmented formulation" begin + elt = Δ^2 + unit = Δ + ub = 0.20001 # Inf + + opt_problem = PropertyT.sos_problem_primal( + elt, + unit, + upper_bound=ub, + augmented=true, + ) + + status, _ = PropertyT.solve( + opt_problem, + scs_optimizer( + eps=1e-10, + max_iters=10_000, + accel=-10, + alpha=1.5, + ), + ) + + @test status == JuMP.OPTIMAL + + λ = JuMP.value(opt_problem[:λ]) + Q = real.(sqrt(JuMP.value.(opt_problem[:P]))) + + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + λ, + Q, + halfradius=2, + augmented=true, + ) + + @test certified + @test isapprox(λ_cert, λ, rtol=1e-5) + @test λ_cert > 2 // 10 + end end end diff --git a/test/1712.07167.jl b/test/1712.07167.jl index 8c9dbd0..96ea904 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -1,104 +1,225 @@ +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 + flPs = [real.(sqrt(JuMP.value.(P))) for P in Ps] + PropertyT.reconstruct(flPs, 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) +end + @testset "1712.07167 Examples" begin - @testset "SL(3,Z)" begin - N = 3 - G = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(G) - autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N)) - - NAME = "SL($N,Z)_orbit" - - rm(NAME, recursive=true, force=true) - sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(1000, accel=20); - upper_bound=0.27, force_compute=false) - - @info sett - - λ = PropertyT.spectral_gap(sett) - @test λ < 0.0 - @test PropertyT.interpret_results(sett, λ) == false - - # second run just checks the solution due to force_compute=false above - @test λ == PropertyT.spectral_gap(sett) - @test PropertyT.check_property_T(sett) == false - - sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(4000, accel=20); - upper_bound=0.27, force_compute=true) - - @info sett - - λ = PropertyT.spectral_gap(sett) - @test λ > 0.269999 - @test PropertyT.interpret_results(sett, λ) == true - - # this should be very fast due to warmstarting: - @test λ ≈ PropertyT.spectral_gap(sett) atol=1e-5 - @test PropertyT.check_property_T(sett) == true - - ########## - # Symmetrizing by SymmetricGroup(3): - - sett = PropertyT.Settings(NAME, G, S, SymmetricGroup(N), with_SCS(4000, accel=20, warm_start=false); - upper_bound=0.27, force_compute=true) - - @info sett - - λ = PropertyT.spectral_gap(sett) - @test λ > 0.269999 - @test PropertyT.interpret_results(sett, λ) == true - end - - @testset "SL(4,Z)" begin - N = 4 - G = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(G) - autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N)) - - NAME = "SL($N,Z)_orbit" - - rm(NAME, recursive=true, force=true) - sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(1000, accel=20); - upper_bound=1.3, force_compute=false) - - @info sett - - λ = PropertyT.spectral_gap(sett) - @test λ < 0.0 - @test PropertyT.interpret_results(sett, λ) == false - - # second run just checks the solution due to force_compute=false above - @test λ == PropertyT.spectral_gap(sett) - @test PropertyT.check_property_T(sett) == false - - sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(7000, accel=20, warm_start=true); - upper_bound=1.3, force_compute=true) - - @info sett - - @time λ = PropertyT.spectral_gap(sett) - @test λ > 1.2999 - @test PropertyT.interpret_results(sett, λ) == true - - # this should be very fast due to warmstarting: - @time @test λ ≈ PropertyT.spectral_gap(sett) atol=1e-5 - @time @test PropertyT.check_property_T(sett) == true - end - @testset "SAut(F₃)" begin N = 3 - G = SAut(FreeGroup(N)) - S = PropertyT.generating_set(G) - autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N)) + G = SpecialAutomorphismGroup(FreeGroup(N)) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) - NAME = "SAut(F$N)_orbit" + 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[2]]), + ) - rm(NAME, recursive=true, force=true) + Δ = let RG = RG, S = S + RG(length(S)) - sum(RG(s) for s in S) + end - sett = PropertyT.Settings(NAME, G, S, autS, with_SCS(1000); - upper_bound=0.15) + elt = Δ^2 + unit = Δ + ub = Inf - @info sett + 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 PropertyT.check_property_T(sett) == false + @test status == JuMP.OPTIMAL + @test !certified + @test λ_cert < 0 + end + + @testset "SL(3,Z) has (T)" begin + n = 3 + + SL = MatrixGroups.SpecialLinearGroup{n}(Int8) + RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true) + + Δ = RSL(length(S)) - sum(RSL(s) for s in S) + + @testset "Wedderburn formulation" begin + P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) + Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + act = PropertyT.action_by_conjugation(SL, Σ) + wd = WedderburnDecomposition( + Rational{Int}, + Σ, + act, + basis(RSL), + StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), + ) + + elt = Δ^2 + unit = Δ + ub = 0.2801 + + @test_throws ErrorException PropertyT.sos_problem_primal( + elt, + unit, + wd, + upper_bound=ub, + augmented=false, + ) + + wdfl = SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RSL), + StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), + ) + + model, varP = PropertyT.sos_problem_primal( + elt, + unit, + wdfl, + upper_bound=ub, + augmented=false, + ) + + status, warm = PropertyT.solve( + model, + scs_optimizer( + eps=1e-10, + max_iters=20_000, + accel=-20, + alpha=1.2, + ), + ) + + Q = @time let varP = varP + Qs = map(varP) do P + real.(sqrt(JuMP.value.(P))) + end + PropertyT.reconstruct(Qs, wdfl) + end + λ = JuMP.value(model[:λ]) + + sos = PropertyT.compute_sos(parent(elt), Q; augmented=false) + + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + λ, + Q, + halfradius=2, + augmented=false, + ) + + @test certified + @test λ_cert >= 28 // 100 + end + + @testset "augmented Wedderburn formulation" begin + elt = Δ^2 + unit = Δ + ub = Inf + + P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) + Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + act = PropertyT.action_by_conjugation(SL, Σ) + + wdfl = SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RSL), + StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), + ) + + opt_problem, varP = PropertyT.sos_problem_primal( + elt, + unit, + wdfl, + upper_bound=ub, + # augmented = true # since both elt and unit are augmented + ) + + status, _ = PropertyT.solve( + opt_problem, + scs_optimizer( + eps=1e-8, + max_iters=20_000, + accel=0, + alpha=1.9, + ), + ) + + @test status == JuMP.OPTIMAL + + Q = @time let varP = varP + Qs = map(varP) do P + real.(sqrt(JuMP.value.(P))) + end + PropertyT.reconstruct(Qs, wdfl) + end + + certified, λ_cert = PropertyT.certify_solution( + elt, + unit, + JuMP.objective_value(opt_problem), + Q, + halfradius=2, + # augmented = true # since both elt and unit are augmented + ) + + @test certified + @test λ_cert > 28 // 100 + end end end diff --git a/test/1812.03456.jl b/test/1812.03456.jl index bc31c74..c3dbe5c 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -1,3 +1,5 @@ +using SparseArrays + @testset "Sq, Adj, Op" begin function isconstant_on_orbit(v, orb) isempty(orb) && return true @@ -6,240 +8,237 @@ end @testset "unit tests" begin - for N in [3,4] - M = MatrixAlgebra(zz, N) - - @test PropertyT.EltaryMat(M, 1, 2) isa MatAlgElem - e12 = PropertyT.EltaryMat(M, 1, 2) - @test e12[1,2] == 1 - @test inv(e12)[1,2] == -1 - - S = PropertyT.generating_set(M) - @test e12 ∈ S - - @test length(PropertyT.generating_set(M)) == 2N*(N-1) - @test all(inv(s) ∈ S for s in S) - - A = SAut(FreeGroup(N)) - @test length(PropertyT.generating_set(A)) == 4N*(N-1) - S = PropertyT.generating_set(A) - @test all(inv(s) ∈ S for s in S) - end - @test PropertyT.isopposite(perm"(1,2,3)(4)", perm"(1,4,2)") @test PropertyT.isadjacent(perm"(1,2,3)", perm"(1,2)(3)") @test !PropertyT.isopposite(perm"(1,2,3)", perm"(1,2)(3)") @test !PropertyT.isadjacent(perm"(1,4)", perm"(2,3)(4)") - @test isconstant_on_orbit([1,1,1,2,2], [2,3]) - @test !isconstant_on_orbit([1,1,1,2,2], [2,3,4]) + @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 - + @testset "Sq, Adj, Op in SL(4,Z)" begin N = 4 - M = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(M) - Δ = PropertyT.Laplacian(S, 2) - RG = parent(Δ) + G = MatrixGroups.SpecialLinearGroup{N}(Int8) - autS = WreathProduct(SymmetricGroup(2), SymmetricGroup(N)) - orbits = PropertyT.orbit_decomposition(autS, RG.basis) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) - @test PropertyT.Sq(RG) isa GroupRingElem - sq = PropertyT.Sq(RG) - @test all(isconstant_on_orbit(sq, orb) for orb in orbits) + Δ = let RG = RG, S = S + RG(length(S)) - sum(RG(s) for s in S) + end - @test PropertyT.Adj(RG) isa GroupRingElem - adj = PropertyT.Adj(RG) - @test all(isconstant_on_orbit(adj, orb) for orb in orbits) + P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) + Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + act = PropertyT.action_by_conjugation(G, Σ) - @test PropertyT.Op(RG) isa GroupRingElem - op = PropertyT.Op(RG) - @test all(isconstant_on_orbit(op, orb) for orb in orbits) + wd = WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), + ) + ivs = SymbolicWedderburn.invariant_vectors(wd) sq, adj, op = PropertyT.SqAdjOp(RG, N) - @test sq == PropertyT.Sq(RG) - @test adj == PropertyT.Adj(RG) - @test op == PropertyT.Op(RG) - e = one(M) - g = PropertyT.EltaryMat(M, 1,2) - h = PropertyT.EltaryMat(M, 1,3) - k = PropertyT.EltaryMat(M, 3,4) + @test all( + isconstant_on_orbit(sq, SparseArrays.nonzeroinds(iv)) for iv in ivs + ) - edges = N*(N-1)÷2 - @test sq[e] == 20*edges + @test all( + isconstant_on_orbit(adj, SparseArrays.nonzeroinds(iv)) for iv in ivs + ) + + @test all( + isconstant_on_orbit(op, SparseArrays.nonzeroinds(iv)) for iv in ivs + ) + + e = one(G) + g = G([alphabet(G)[MatrixGroups.ElementaryMatrix{N}(1, 2, Int8(1))]]) + h = G([alphabet(G)[MatrixGroups.ElementaryMatrix{N}(1, 3, Int8(1))]]) + k = G([alphabet(G)[MatrixGroups.ElementaryMatrix{N}(3, 4, Int8(1))]]) + + @test sq[e] == 120 @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[e] == 384 + @test adj[g] == adj[h] == -32 @test adj[g^2] == adj[h^2] == 0 - @test adj[g*h] == adj[h*g] # == ... + @test adj[g*h] == adj[h*g] == 2 + @test adj[k*h] == adj[h*k] == 1 - - # @test op[e] == ... - @test op[g] == op[h] # == ... + @test op[e] == 96 + @test op[g] == op[h] == -8 @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[g*k] == op[k*g] == 2 @test op[h*k] == op[k*h] == 0 end + + @testset "SAut(F₃)" begin + n = 3 + G = SpecialAutomorphismGroup(FreeGroup(n)) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) + sq, adj, op = PropertyT.SqAdjOp(RG, n) + + @test sq(one(G)) == 216 + @test all(sq(g) == -16 for g in gens(G)) + @test adj(one(G)) == 384 + @test all(adj(g) == -32 for g in gens(G)) + @test iszero(op) + end end - @testset "1812.03456 examples" begin - - function SOS_residual(x::GroupRingElem, Q::Matrix) - RG = parent(x) - @time sos = PropertyT.compute_SOS(RG, Q); - return x - sos - end - - function check_positivity(elt, Δ, orbit_data, upper_bound, warm=nothing; with_solver=with_SCS(20_000, accel=10)) - SDP_problem, varP = PropertyT.SOS_problem_primal(elt, Δ, orbit_data; upper_bound=upper_bound) - - status, warm = PropertyT.solve(SDP_problem, with_solver, warm); - Base.Libc.flush_cstdio() - @info "Optimization status:" status - - λ = value(SDP_problem[:λ]) - Ps = [value.(P) for P in varP] - - Qs = real.(sqrt.(Ps)); - Q = PropertyT.reconstruct(Qs, orbit_data); - - b = SOS_residual(elt - λ*Δ, Q) - return b, λ, warm - end - @testset "SL(3,Z)" begin - N = 3 - halfradius = 2 - M = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(M) - Δ = PropertyT.Laplacian(S, halfradius) - RG = parent(Δ) - orbit_data = PropertyT.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N))) - orbit_data = PropertyT.decimate(orbit_data); + n = 3 + + G = MatrixGroups.SpecialLinearGroup{n}(Int8) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) + + Δ = RG(length(S)) - sum(RG(s) for s in S) + + 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 = SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), + ) + + sq, adj, op = PropertyT.SqAdjOp(RG, n) @testset "Sq₃ is SOS" begin - elt = PropertyT.Sq(RG) - UB = 0.05 # 0.105? + elt = sq + UB = Inf # λ ≈ 0.1040844 - residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB) - Base.Libc.flush_cstdio() - @info "obtained λ and residual" λ norm(residual, 1) - - @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity - - @test 2^2*norm(residual, 1) < 2λ/100 + status, certified, λ_cert = check_positivity( + elt, + Δ, + wd, + upper_bound=UB, + halfradius=2, + optimizer=cosmo_optimizer(accel=50, alpha=1.9) + ) + @test status == JuMP.OPTIMAL + @test certified + @test λ_cert > 104 // 1000 end @testset "Adj₃ is SOS" begin - elt = PropertyT.Adj(RG) - UB = 0.1 # 0.157? + elt = adj + UB = Inf # λ ≈ 0.15858018 - residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB) - Base.Libc.flush_cstdio() - @info "obtained λ and residual" λ norm(residual, 1) - - @test 2^2*norm(residual, 1) < λ - @test 2^2*norm(residual, 1) < λ/100 + status, certified, λ_cert = check_positivity( + elt, + Δ, + wd, + upper_bound=UB, + halfradius=2, + optimizer=cosmo_optimizer(accel=50, alpha=1.9) + ) + @test status == JuMP.OPTIMAL + @test certified + @test λ_cert > 1585 // 10000 end @testset "Op₃ is empty, so can not be certified" begin - elt = PropertyT.Op(RG) + elt = op 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) > λ + status, certified, λ_cert = check_positivity( + elt, + Δ, + wd, + upper_bound=UB, + halfradius=2, + optimizer=cosmo_optimizer(accel=50, alpha=1.9) + ) + @test status == JuMP.OPTIMAL + @test !certified + @test λ_cert < 0 end end @testset "SL(4,Z)" begin - N = 4 - halfradius = 2 - M = MatrixAlgebra(zz, N) - S = PropertyT.generating_set(M) - Δ = PropertyT.Laplacian(S, halfradius) - RG = parent(Δ) - orbit_data = PropertyT.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N))) - orbit_data = PropertyT.decimate(orbit_data); + n = 4 - @testset "Sq₄ is SOS" begin - elt = PropertyT.Sq(RG) - UB = 0.2 # 0.3172 + G = MatrixGroups.SpecialLinearGroup{n}(Int8) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true) - residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB) - Base.Libc.flush_cstdio() - @info "obtained λ and residual" λ norm(residual, 1) + Δ = RG(length(S)) - sum(RG(s) for s in S) - @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity - @test 2^2*norm(residual, 1) < λ/100 + 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 = SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), + ) + + sq, adj, op = PropertyT.SqAdjOp(RG, n) + + @testset "Sq is SOS" begin + elt = sq + UB = Inf # λ ≈ 0.31670 + + status, certified, λ_cert = check_positivity( + elt, + Δ, + wd, + upper_bound=UB, + halfradius=2, + optimizer=cosmo_optimizer(accel=50, alpha=1.9) + ) + @test status == JuMP.OPTIMAL + @test certified + @test λ_cert > 316 // 1000 end - @testset "Adj₄ is SOS" begin - elt = PropertyT.Adj(RG) - UB = 0.3 # 0.5459? + @testset "Adj is SOS" begin + elt = adj + UB = 0.541 # λ ≈ 0.545710 - 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 + status, certified, λ_cert = check_positivity( + elt, + Δ, + wd, + upper_bound=UB, + halfradius=2, + optimizer=cosmo_optimizer(accel=50, alpha=1.9) + ) + @test status == JuMP.OPTIMAL + @test certified + @test λ_cert > 54 // 100 end - @testset "we can't cerify that Op₄ SOS" begin - elt = PropertyT.Op(RG) - UB = 2.0 + @testset "Op is a sum of squares, but not an order unit" begin + elt = op + UB = Inf - 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't certify positivity - end - - @testset "Adj₄ + Op₄ is SOS" begin - elt = PropertyT.Adj(RG) + PropertyT.Op(RG) - UB = 0.6 # 0.82005 - - residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB) - Base.Libc.flush_cstdio() - @info "obtained λ and residual" λ norm(residual, 1) - - @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity - @test 2^2*norm(residual, 1) < λ/100 + status, certified, λ_cert = check_positivity( + elt, + Δ, + wd, + upper_bound=UB, + halfradius=2, + optimizer=cosmo_optimizer(accel=50, alpha=1.9) + ) + @test status == JuMP.OPTIMAL + @test !certified + @test -1e-2 < λ_cert < 0 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.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N))) - # orbit_data = PropertyT.decimate(orbit_data); - # - # @time elt = PropertyT.Adj(RG) + 100PropertyT.Op(RG) - # UB = 0.05 - # - # warm = nothing - # - # residual, λ, warm = check_positivity(elt, Δ, orbit_data, UB, warm, with_solver=with_SCS(20_000, accel=10)) - # @info "obtained λ and residual" λ norm(residual, 1) - # - # @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity - # @test 2^2*norm(residual, 1) < λ/100 - # end end diff --git a/test/actions.jl b/test/actions.jl index 2cee0a7..af7d123 100644 --- a/test/actions.jl +++ b/test/actions.jl @@ -1,101 +1,210 @@ -@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) - - function ssgs(A::AutGroup, i, j) - rmuls = [Groups.transvection_R(i,j), Groups.transvection_R(j,i)] - lmuls = [Groups.transvection_L(i,j), Groups.transvection_L(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.wlmetric_ball(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 [SymmetricGroup(N), WreathProduct(SymmetricGroup(2), SymmetricGroup(N))] - @test all(one(M)^g == one(M) for g in G) - @test all(rdict[m^g] <= sizes[1] for g in G for m in S) - @test all(m^g*n^g == (m*n)^g for g in G for m in S for n in S) - - @test all(Δ^g == Δ for g in G) - @test all(x^g == RG(1) - RG(elt^g) for g in G) - - @test all(2RG(elt2^g) - RG(elt^g) == y^g for g in G) +function test_action(basis, group, act) + action = SymbolicWedderburn.action + return @testset "action definition" begin + @test all(basis) do b + e = one(group) + action(act, e, b) == b end - end - @testset "small Laplacians" begin - for (i,j) in PropertyT.indexing(N) - Sij = ssgs(M, i,j) - Δij= PropertyT.spLaplacian(RG, Sij) + a = let a = rand(basis) + while isone(a) + a = rand(basis) + end + @assert !isone(a) + a + end - @test all(Δij^p == PropertyT.spLaplacian(RG, ssgs(M, p[i], p[j])) for p in SymmetricGroup(N)) + g, h = let g_h = rand(group, 2) + while any(isone, g_h) + g_h = rand(group, 2) + end + @assert all(!isone, g_h) + g_h + end - @test all(Δij^g == PropertyT.spLaplacian(RG, ssgs(M, g.p[i], g.p[j])) for g in WreathProduct(SymmetricGroup(2), SymmetricGroup(N))) + action = SymbolicWedderburn.action + @test action(act, g, a) in basis + @test action(act, h, a) in basis + @test action(act, h, action(act, g, a)) == action(act, g * h, a) + + @test all([(g, h) for g in group for h in group]) do (g, h) + x = action(act, h, action(act, g, a)) + y = action(act, g * h, a) + x == y + end + + if act isa SymbolicWedderburn.ByPermutations + @test all(basis) do b + action(act, g, b) ∈ basis && action(act, h, b) ∈ basis + end 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.wlmetric_ball(S, one(M), radius=2halfradius); +## Testing - 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 "Actions on SL(3,ℤ)" begin + n = 3 - @testset "correctness of actions" begin + SL = MatrixGroups.SpecialLinearGroup{n}(Int8) + RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true) - Δ = length(S)*RG(1) - sum(RG(s) for s in S) - @test Δ == PropertyT.spLaplacian(RG, S) + @testset "Permutation action" begin - elt = S[5] - x = RG(1) - RG(elt) - elt2 = E_R[rand(sizes[1]:sizes[2])] - y = 2RG(elt2) - RG(elt) + Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) + ΓpA = PropertyT.action_by_conjugation(SL, Γ) - for G in [SymmetricGroup(N), WreathProduct(SymmetricGroup(2), SymmetricGroup(N))] - @test all(one(M)^g == one(M) for g in G) - @test all(rdict[m^g] <= sizes[1] for g in G for m in S) - @test all(m^g*n^g == (m*n)^g for g in G for m in S for n in S) + test_action(basis(RSL), Γ, ΓpA) - @test all(Δ^g == Δ for g in G) - @test all(x^g == RG(1) - RG(elt^g) for g in G) + @testset "mps is successful" begin - @test all(2RG(elt2^g) - RG(elt^g) == y^g for g in G) + charsΓ = + SymbolicWedderburn.Character{ + Rational{Int}, + }.(SymbolicWedderburn.irreducible_characters(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @test all(simple) + end + + @testset "Wedderburn decomposition" begin + wd = SymbolicWedderburn.WedderburnDecomposition( + Rational{Int}, + Γ, + ΓpA, + basis(RSL), + StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]) + ) + + @test length(invariant_vectors(wd)) == 918 + @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [40, 23, 18] + @test all(issimple, direct_summands(wd)) end end - for (i,j) in PropertyT.indexing(N) - Sij = ssgs(M, i,j) - Δij= PropertyT.spLaplacian(RG, Sij) + @testset "Wreath action" begin - @test all(Δij^p == PropertyT.spLaplacian(RG, ssgs(M, p[i], p[j])) for p in SymmetricGroup(N)) + Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) + PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + end - @test all(Δij^g == PropertyT.spLaplacian(RG, ssgs(M, g.p[i], g.p[j])) for g in WreathProduct(SymmetricGroup(2), SymmetricGroup(N))) + ΓpA = PropertyT.action_by_conjugation(SL, Γ) + + test_action(basis(RSL), Γ, ΓpA) + + @testset "mps is successful" begin + + charsΓ = + SymbolicWedderburn.Character{ + Rational{Int}, + }.(SymbolicWedderburn.irreducible_characters(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @test all(simple) + end + + @testset "Wedderburn decomposition" begin + wd = SymbolicWedderburn.WedderburnDecomposition( + Rational{Int}, + Γ, + ΓpA, + basis(RSL), + StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]) + ) + + @test length(invariant_vectors(wd)) == 247 + @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [14, 9, 6, 14, 12] + @test all(issimple, direct_summands(wd)) + end end end +@testset "Actions on SAut(F4)" begin + + n = 4 + + SAutFn = SpecialAutomorphismGroup(FreeGroup(n)) + RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn, halfradius=1, twisted=true) + + @testset "Permutation action" begin + + Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) + ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ) + + test_action(basis(RSAutFn), Γ, ΓpA) + + @testset "mps is successful" begin + + charsΓ = + SymbolicWedderburn.Character{ + Rational{Int}, + }.(SymbolicWedderburn.irreducible_characters(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @test all(simple) + end + + @testset "Wedderburn decomposition" begin + wd = SymbolicWedderburn.WedderburnDecomposition( + Rational{Int}, + Γ, + ΓpA, + basis(RSAutFn), + StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]) + ) + + @test length(invariant_vectors(wd)) == 93 + @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [4, 8, 5, 4] + @test all(issimple, direct_summands(wd)) + end + end + + @testset "Wreath action" begin + + Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) + PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + end + + ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ) + + test_action(basis(RSAutFn), Γ, ΓpA) + + @testset "mps is successful" begin + + charsΓ = + SymbolicWedderburn.Character{ + Rational{Int}, + }.(SymbolicWedderburn.irreducible_characters(Γ)) + + RΓ = SymbolicWedderburn._group_algebra(Γ) + + @time mps, simple = + SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) + @test all(simple) + end + + @testset "Wedderburn decomposition" begin + wd = SymbolicWedderburn.WedderburnDecomposition( + Rational{Int}, + Γ, + ΓpA, + basis(RSAutFn), + StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]) + ) + + @test length(invariant_vectors(wd)) == 18 + @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [1, 1, 2, 2, 1, 2, 2, 1] + @test all(issimple, direct_summands(wd)) + end + end end diff --git a/test/graded_adj.jl b/test/graded_adj.jl new file mode 100644 index 0000000..117b1dd --- /dev/null +++ b/test/graded_adj.jl @@ -0,0 +1,40 @@ +@testset "Adj for SpN via grading" begin + + genus = 3 + halfradius = 2 + + 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 + + @testset "Adj correctness: genus=$genus" begin + + all_subtypes = ( + :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂ + ) + + @test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384 + @test iszero(PropertyT.Adj(Δs, Symbol("A₁×A₁"))) + @test iszero(PropertyT.Adj(Δs, Symbol("C₁×C₁"))) + + @testset "divisibility by 16" begin + for subtype in all_subtypes + subtype in (:A₁, :C₁) && continue + @test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16) + end + end + @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2 + end + +end + diff --git a/test/optimizers.jl b/test/optimizers.jl new file mode 100644 index 0000000..f415de7 --- /dev/null +++ b/test/optimizers.jl @@ -0,0 +1,54 @@ +## Optimizers + +import JuMP +import SCS + +function scs_optimizer(; + accel=10, + alpha=1.5, + eps=1e-9, + max_iters=100_000, + verbose=true, + linear_solver=SCS.DirectSolver +) + return JuMP.optimizer_with_attributes( + SCS.Optimizer, + "acceleration_lookback" => accel, + "acceleration_interval" => 10, + "alpha" => alpha, + "eps_abs" => eps, + "eps_rel" => eps, + "linear_solver" => linear_solver, + "max_iters" => max_iters, + "warm_start" => true, + "verbose" => verbose, + ) +end + +import COSMO + +function cosmo_optimizer(; + accel=15, + alpha=1.6, + eps=1e-9, + max_iters=100_000, + verbose=true, + verbose_timing=verbose, + decompose=false +) + return JuMP.optimizer_with_attributes( + COSMO.Optimizer, + "accelerator" => COSMO.with_options( + COSMO.AndersonAccelerator, + mem=max(accel, 2) + ), + "alpha" => alpha, + "decompose" => decompose, + "eps_abs" => eps, + "eps_rel" => eps, + "max_iter" => max_iters, + "verbose" => verbose, + "verbose_timing" => verbose_timing, + "check_termination" => 200, + ) +end diff --git a/test/runtests.jl b/test/runtests.jl index dfb0413..daf0aa3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,18 +1,27 @@ using Test -using LinearAlgebra, SparseArrays -using AbstractAlgebra, Groups, GroupRings +using LinearAlgebra +using SparseArrays +BLAS.set_num_threads(1) +ENV["OMP_NUM_THREADS"] = 4 + +using Groups +using Groups.GroupsCore +import Groups.MatrixGroups + using PropertyT -using JLD +using SymbolicWedderburn +using SymbolicWedderburn.StarAlgebras +using SymbolicWedderburn.PermutationGroups -using JuMP, SCS +include("optimizers.jl") -with_SCS(iters; accel=0, eps=1e-10, warm_start=true) = - with_optimizer(SCS.Optimizer, - linear_solver=SCS.DirectSolver, max_iters=iters, - acceleration_lookback=accel, eps=eps, warm_start=warm_start) +@testset "PropertyT" begin -include("1703.09680.jl") -include("actions.jl") -include("1712.07167.jl") -include("SOS_correctness.jl") -include("1812.03456.jl") + include("actions.jl") + + include("1703.09680.jl") + include("1712.07167.jl") + include("1812.03456.jl") + + include("graded_adj.jl") +end From abf2d014cc0aec4d6dff086defb56dcc3d328393 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Mon, 7 Nov 2022 18:48:26 +0100 Subject: [PATCH 12/19] add/update workflows --- .github/workflows/CompatHelper.yml | 16 +++++++++++++ .github/workflows/TagBot.yml | 10 +++++++-- .github/workflows/{runtests.yml => ci.yml} | 26 +++++++++++++--------- 3 files changed, 39 insertions(+), 13 deletions(-) create mode 100644 .github/workflows/CompatHelper.yml rename .github/workflows/{runtests.yml => ci.yml} (62%) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..bcdb51a --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,16 @@ +name: CompatHelper +on: + schedule: + - cron: '00 00 * * *' + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.COMPATHELPER_PRIV }} # optional + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index d77d3a0..33fd52d 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -1,11 +1,17 @@ name: TagBot on: - schedule: - - cron: 0 * * * * + issue_comment: # THIS BIT IS NEW + types: + - created + workflow_dispatch: jobs: TagBot: + # THIS 'if' LINE IS NEW + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + # NOTHING BELOW HAS CHANGED runs-on: ubuntu-latest steps: - uses: JuliaRegistries/TagBot@v1 with: token: ${{ secrets.GITHUB_TOKEN }} + # ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.github/workflows/runtests.yml b/.github/workflows/ci.yml similarity index 62% rename from .github/workflows/runtests.yml rename to .github/workflows/ci.yml index bd20140..db65c50 100644 --- a/.github/workflows/runtests.yml +++ b/.github/workflows/ci.yml @@ -1,11 +1,6 @@ name: CI on: - push: - branches: - - master - pull_request: - branches: - - master + - pull_request jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} @@ -13,22 +8,32 @@ jobs: strategy: matrix: version: - - '1.3' - - '1.5' + - '1.6' + - '1' - 'nightly' os: - ubuntu-latest - - macOS-latest arch: - x64 + allow_failures: + - julia: nightly fail-fast: false - steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest - uses: julia-actions/julia-processcoverage@v1 @@ -38,4 +43,3 @@ jobs: flags: unittests name: codecov-umbrella fail_ci_if_error: false - token: ${{ secrets.CODECOV_TOKEN }} From 6e9cc8e0f68400969ac61188c7513ea9ee04def0 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 09:59:07 +0100 Subject: [PATCH 13/19] action on FPGroupElement via action on the underlying word --- src/actions/alphabet_permutation.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/actions/alphabet_permutation.jl b/src/actions/alphabet_permutation.jl index a3e15e0..8dfe89d 100644 --- a/src/actions/alphabet_permutation.jl +++ b/src/actions/alphabet_permutation.jl @@ -28,17 +28,17 @@ end function SymbolicWedderburn.action( act::AlphabetPermutation, γ::Groups.GroupElement, - w::Groups.AbstractWord, + g::Groups.AbstractFPGroupElement, ) - return w^(act.perms[γ]) + G = parent(g) + w = SymbolicWedderburn.action(act, γ, word(g)) + return G(w) end function SymbolicWedderburn.action( act::AlphabetPermutation, γ::Groups.GroupElement, - g::Groups.AbstractFPGroupElement, + w::Groups.AbstractWord, ) - G = parent(g) - w = word(g)^(act.perms[γ]) - return G(w) + return w^(act.perms[γ]) end From 903c1683ff0d592ee22e835dea77e153dfe00e66 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 09:59:26 +0100 Subject: [PATCH 14/19] remove unused (old) constraints function --- src/sos_sdps.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl index 250388c..d2387ed 100644 --- a/src/sos_sdps.jl +++ b/src/sos_sdps.jl @@ -37,19 +37,6 @@ function sos_problem_dual( return model end -function constraints( - mstr::AbstractMatrix{<:Integer}, - total_length=maximum(mstr), -) - cnstrs = [Vector{eltype(mstr)}() for _ = 1:total_length] - li = LinearIndices(mstr) - - for (idx, val) in pairs(mstr) - push!(cnstrs[val], li[idx]) - end - return cnstrs -end - function constraints( basis::StarAlgebras.AbstractBasis, mstr::AbstractMatrix{<:Integer}; From 518c661dc8323b85ee55f81eea7cd3dec69bd90c Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 10:01:31 +0100 Subject: [PATCH 15/19] add feasibility/warmstarting tests --- test/1703.09680.jl | 39 ++++++++++++++++++++++++++++++++++++++- test/1712.07167.jl | 17 ++++++++++++++++- test/1812.03456.jl | 11 +++++++++++ 3 files changed, 65 insertions(+), 2 deletions(-) diff --git a/test/1703.09680.jl b/test/1703.09680.jl index 75e76c1..074d263 100644 --- a/test/1703.09680.jl +++ b/test/1703.09680.jl @@ -77,6 +77,17 @@ end @test status == JuMP.OPTIMAL @test certified @test λ > 1 + + m = PropertyT.sos_problem_dual(elt, unit) + PropertyT.solve(m, scs_optimizer( + eps=1e-10, + 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 end @testset "SAut(F₂)" begin @@ -95,7 +106,7 @@ end status, certified, λ = check_positivity( elt, unit, - upper_bound=0.1, + upper_bound=ub, halfradius=2, optimizer=scs_optimizer( eps=1e-10, @@ -108,6 +119,32 @@ end @test status == JuMP.ALMOST_OPTIMAL @test λ < 0 @test !certified + + @time sos_problem = + PropertyT.sos_problem_primal(elt, upper_bound=ub) + + status, _ = PropertyT.solve( + sos_problem, + cosmo_optimizer( + eps=1e-7, + max_iters=10_000, + accel=0, + alpha=1.9, + ) + ) + @test status == JuMP.OPTIMAL + P = JuMP.value.(sos_problem[:P]) + Q = real.(sqrt(P)) + certified, λ_cert = PropertyT.certify_solution( + elt, + zero(elt), + 0.0, + Q, + halfradius=2, + ) + @test !certified + @test λ_cert < 0 + end @testset "SL(3,Z) has (T)" begin diff --git a/test/1712.07167.jl b/test/1712.07167.jl index 96ea904..982806e 100644 --- a/test/1712.07167.jl +++ b/test/1712.07167.jl @@ -134,14 +134,29 @@ end status, warm = PropertyT.solve( model, - scs_optimizer( + cosmo_optimizer( eps=1e-10, max_iters=20_000, + accel=50, + alpha=1.9, + ), + ) + + @test status == JuMP.OPTIMAL + + status, _ = PropertyT.solve( + model, + scs_optimizer( + eps=1e-10, + max_iters=100, accel=-20, alpha=1.2, ), + warm ) + @test status == JuMP.OPTIMAL + Q = @time let varP = varP Qs = map(varP) do P real.(sqrt(JuMP.value.(P))) diff --git a/test/1812.03456.jl b/test/1812.03456.jl index c3dbe5c..b0d9482 100644 --- a/test/1812.03456.jl +++ b/test/1812.03456.jl @@ -148,10 +148,21 @@ end @test status == JuMP.OPTIMAL @test certified @test λ_cert > 1585 // 10000 + + m, _ = PropertyT.sos_problem_primal(elt, wd) + PropertyT.solve( + m, + scs_optimizer(max_iters=5000, accel=50, alpha=1.9) + ) + + @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT) + @test abs(JuMP.objective_value(m)) < 1e-3 end @testset "Op₃ is empty, so can not be certified" begin elt = op + @test iszero(op) + UB = Inf status, certified, λ_cert = check_positivity( From b2b08e68c17a1444409286102e2ca3703ddbc3ac Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 10:02:21 +0100 Subject: [PATCH 16/19] include tests for constraint Matrices --- test/constratint_matrices.jl | 21 +++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 22 insertions(+) create mode 100644 test/constratint_matrices.jl diff --git a/test/constratint_matrices.jl b/test/constratint_matrices.jl new file mode 100644 index 0000000..afdee61 --- /dev/null +++ b/test/constratint_matrices.jl @@ -0,0 +1,21 @@ +@testset "ConstraintMatrix" begin + @test PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) isa AbstractMatrix + + cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) + + @test cm == Float64[ + -π π + 2π 0.0 + 0.0 π + ] + + @test collect(PropertyT.nzpairs(cm)) == [ + 1 => 3.141592653589793 + 2 => 3.141592653589793 + 2 => 3.141592653589793 + 4 => 3.141592653589793 + 6 => 3.141592653589793 + 1 => -3.141592653589793 + 1 => -3.141592653589793 + ] +end diff --git a/test/runtests.jl b/test/runtests.jl index daf0aa3..cdeac69 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,7 @@ include("optimizers.jl") @testset "PropertyT" begin + include("constratint_matrices.jl") include("actions.jl") include("1703.09680.jl") From 4e93e78a9b131930324fcebac9f4f2c8a63a78fb Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 11:49:04 +0100 Subject: [PATCH 17/19] test: old sq,adj,op = the graded Adj for SL/SAut --- src/gradings.jl | 1 - src/roots.jl | 8 +--- test/graded_adj.jl | 110 ++++++++++++++++++++++++++++++++------------- 3 files changed, 80 insertions(+), 39 deletions(-) diff --git a/src/gradings.jl b/src/gradings.jl index 20edd29..f6de948 100644 --- a/src/gradings.jl +++ b/src/gradings.jl @@ -25,7 +25,6 @@ end grading(s::MatrixGroups.ElementarySymplectic) = Roots.Root(s) grading(e::MatrixGroups.ElementaryMatrix) = Roots.Root(e) -grading(t::Groups.Transvection) = grading(Groups._abelianize(t)) function grading(g::FPGroupElement) if length(word(g)) == 1 diff --git a/src/roots.jl b/src/roots.jl index 080190a..94beb1d 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -60,10 +60,7 @@ function positive(roots::AbstractVector{<:Root{N}}) where {N} return filter(α -> dot(α, pd) > 0.0, roots) end -Base.:~(α::AbstractRoot, β::AbstractRoot) = isproportional(α, β) -⟂(α::AbstractRoot, β::AbstractRoot) = isorthogonal(α, β) - -function Base.show(io::IO, r::Root{N}) where {N} +function Base.show(io::IO, r::Root) print(io, "Root$(r.coord)") end @@ -73,8 +70,7 @@ function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N} print(io, "Root in ℝ^$N of length $l\n", r.coord) end -E(N, i::Integer) = Root(ntuple(k -> k == i ? 1 : 0, N)) -𝕖(N, i) = E(N, i) +𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N)) 𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N)) """ diff --git a/test/graded_adj.jl b/test/graded_adj.jl index 117b1dd..7c65e30 100644 --- a/test/graded_adj.jl +++ b/test/graded_adj.jl @@ -1,40 +1,86 @@ -@testset "Adj for SpN via grading" begin +@testset "Adj via grading" begin - genus = 3 - halfradius = 2 + @testset "SL(n,Z) & Aut(F₄)" begin + n = 4 + halfradius = 1 + SL = MatrixGroups.SpecialLinearGroup{n}(Int8) + RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=halfradius, twisted=true) - SpN = MatrixGroups.SymplecticGroup{2genus}(Int8) + Δ = RSL(length(S)) - sum(RSL(s) for s in S) - 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 - - @testset "Adj correctness: genus=$genus" begin - - all_subtypes = ( - :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂ - ) - - @test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384 - @test iszero(PropertyT.Adj(Δs, Symbol("A₁×A₁"))) - @test iszero(PropertyT.Adj(Δs, Symbol("C₁×C₁"))) - - @testset "divisibility by 16" begin - for subtype in all_subtypes - subtype in (:A₁, :C₁) && continue - @test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16) - end + Δs = let ψ = identity + PropertyT.laplacians( + RSL, + S, + x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])), + ) end - @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2 + + sq, adj, op = PropertyT.SqAdjOp(RSL, n) + + @test PropertyT.Adj(Δs, :A₁) == sq + @test PropertyT.Adj(Δs, :A₂) == adj + @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op + + + halfradius = 1 + G = SpecialAutomorphismGroup(FreeGroup(n)) + RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true) + + Δ = RG(length(S)) - sum(RG(s) for s in S) + + Δs = let ψ = Groups.Homomorphism(Groups._abelianize, G, SL) + PropertyT.laplacians( + RG, + S, + x -> (gx = PropertyT.grading(ψ(x)); Set([gx, -gx])), + ) + end + + sq, adj, op = PropertyT.SqAdjOp(RG, n) + + @test PropertyT.Adj(Δs, :A₁) == sq + @test PropertyT.Adj(Δs, :A₂) == adj + @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op end + + @testset "Symplectic group" begin + genus = 3 + 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 + + @testset "Adj correctness: genus=$genus" begin + + all_subtypes = ( + :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂ + ) + + @test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384 + @test iszero(PropertyT.Adj(Δs, Symbol("A₁×A₁"))) + @test iszero(PropertyT.Adj(Δs, Symbol("C₁×C₁"))) + + @testset "divisibility by 16" begin + for subtype in all_subtypes + subtype in (:A₁, :C₁) && continue + @test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16) + end + end + @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2 + end + end end From fee9c38537a56662ae419aeb98373d5edffb1066 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 12:00:20 +0100 Subject: [PATCH 18/19] remove old unsed definition for Positive --- src/gradings.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/gradings.jl b/src/gradings.jl index f6de948..1f20d6b 100644 --- a/src/gradings.jl +++ b/src/gradings.jl @@ -13,16 +13,6 @@ function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N} end end -function Roots.positive( - generating_set::AbstractVector{<:MatrixGroups.ElementarySymplectic}, -) - r = Roots._positive_direction(Roots.Root(first(generating_set))) - pos_gens = [ - s for s in generating_set if s.val > 0.0 && dot(Roots.Root(s), r) ≥ 0.0 - ] - return pos_gens -end - grading(s::MatrixGroups.ElementarySymplectic) = Roots.Root(s) grading(e::MatrixGroups.ElementaryMatrix) = Roots.Root(e) From ccffa9e7cf4680c410123a55ccc7e574278eb52b Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 8 Nov 2022 12:00:55 +0100 Subject: [PATCH 19/19] bump to 0.4 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a14a559..ade1715 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PropertyT" uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d" authors = ["Marek Kaluba "] -version = "0.3.2" +version = "0.4.0" [deps] Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"