kalmar/PropertyT.jl
1
0
Fork 0
mirror of https://github.com/kalmarek/PropertyT.jl.git synced 2024-12-03 18:46:27 +01:00
Code Issues Releases Wiki Activity

Compare commits

base: kalmar:902be093db39c4d9d52cea5f84ad0c61227ce9be
Branches Tags
kalmar:master
kalmar:mk/pra_group
kalmar:mk/update_PermutationGroups_etal
kalmar:mk/orbit_basis
kalmar:mk/tests_chevalley
kalmar:enh/G₂_here_we_come
kalmar:enh/update_deps
kalmar:1712.07167
kalmar:1703.09680v1
..
compare: kalmar:0127d05594b83dd03b2908ade48368cf626e8620
Branches Tags
kalmar:mk/pra_group
kalmar:master
kalmar:mk/update_PermutationGroups_etal
kalmar:mk/orbit_basis
kalmar:mk/tests_chevalley
kalmar:enh/G₂_here_we_come
kalmar:enh/update_deps
kalmar:1712.07167
kalmar:1703.09680v1

No commits in common. "902be093db39c4d9d52cea5f84ad0c61227ce9be" and "0127d05594b83dd03b2908ade48368cf626e8620" have entirely different histories.
902be093db ... 0127d05594

28 changed files with 578 additions and 1438 deletions
Show Stats Expand all files Collapse all files

9
.JuliaFormatter.toml
View File

@ -1,9 +0,0 @@
# Configuration file for JuliaFormatter.jl
# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/
always_for_in = true
always_use_return = true
margin = 80
remove_extra_newlines = true
separate_kwargs_with_semicolon = true
short_to_long_function_def = true

8
Project.toml
View File

@ -1,12 +1,11 @@
name = "PropertyT" name = "PropertyT"
uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d" uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d"
authors = ["Marek Kaluba <kalmar@amu.edu.pl>"] authors = ["Marek Kaluba <kalmar@amu.edu.pl>"]
version = "0.5.0" version = "0.4.0"
[deps] [deps]
Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557" Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
@ -18,12 +17,11 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
COSMO = "0.8" COSMO = "0.8"
Groups = "0.7" Groups = "0.7"
IntervalArithmetic = "0.20" IntervalArithmetic = "0.20"
IntervalMatrices = "0.8"
JuMP = "1.3" JuMP = "1.3"
ProgressMeter = "1.7" ProgressMeter = "1.7"
SCS = "1.1" SCS = "1.1.0"
StaticArrays = "1" StaticArrays = "1"
SymbolicWedderburn = "0.3.4" SymbolicWedderburn = "0.3.2"
julia = "1.6" julia = "1.6"
[extras] [extras]

98
scripts/G₂_Adj.jl
View File

@ -1,98 +0,0 @@
using LinearAlgebra
BLAS.set_num_threads(8)
using MKL_jll
ENV["OMP_NUM_THREADS"] = 4
using Groups
import Groups.MatrixGroups
include(joinpath(@__DIR__, "../test/optimizers.jl"))
using PropertyT
using PropertyT.SymbolicWedderburn
using PropertyT.PermutationGroups
using PropertyT.StarAlgebras
include(joinpath(@__DIR__, "argparse.jl"))
include(joinpath(@__DIR__, "utils.jl"))
# const N = parsed_args["N"]
const HALFRADIUS = parsed_args["halfradius"]
const UPPER_BOUND = parsed_args["upper_bound"]
include(joinpath(@__DIR__, "./G₂_gens.jl"))
G, roots, Weyl = G₂_roots_weyl()
@info "Running Adj² - λ·Δ sum of squares decomposition for G₂"
@info "computing group algebra structure"
RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
@info "computing WedderburnDecomposition"
wd = let Σ = Weyl, RG = RG
act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
Dict(g => PermutationGroups.perm(g) for g in Σ),
)
@time SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
semisimple = false,
)
end
@info wd
function desubscriptify(symbol::Symbol)
digits = [
Int(l) - 0x2080 for
l in reverse(string(symbol)) if 0 ≤ Int(l) - 0x2080 ≤ 9
]
res = 0
for (i, d) in enumerate(digits)
res += 10^(i - 1) * d
end
return res
end
function PropertyT.grading(g::MatrixGroups.MatrixElt, roots = roots)
id = desubscriptify(g.id)
return roots[id]
end
Δ = RG(length(S)) - sum(RG(s) for s in S)
Δs = PropertyT.laplacians(
RG,
S,
x -> (gx = PropertyT.grading(x); Set([gx, -gx])),
)
elt = PropertyT.Adj(Δs)
@assert elt == Δ^2 - PropertyT.Sq(Δs)
unit = Δ
@time model, varP = PropertyT.sos_problem_primal(
elt,
unit,
wd;
upper_bound = UPPER_BOUND,
augmented = true,
show_progress = true,
)
solve_in_loop(
model,
wd,
varP;
logdir = "./log/G2/r=$HALFRADIUS/Adj-$(UPPER_BOUND)Δ",
optimizer = scs_optimizer(;
linear_solver = SCS.MKLDirectSolver,
eps = 1e-9,
max_iters = 100_000,
accel = 50,
alpha = 1.95,
),
data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
)

308
scripts/G₂_gens.jl
View File

@ -1,308 +0,0 @@
#= GAP code to generate matrices
alg := SimpleLieAlgebra("G", 2, Rationals);
root_sys := RootSystem(alg);
pos_gens := PositiveRootVectors(root_sys);
pos_rts := PositiveRoots(root_sys);
neg_gens := NegativeRootVectors(root_sys);
neg_rts := NegativeRoots(root_sys);
alg_gens := ShallowCopy(pos_gens);;
Append(alg_gens, neg_gens);
grading := ShallowCopy(pos_rts);
Append(grading, neg_rts);
mats := List(alg_gens, x->AdjointMatrix(Basis(alg), x));
W := WeylGroup(root_sys);
PW := Action(W, grading, OnRight);
=#
using LinearAlgebra
function matrix_exp(M::AbstractMatrix{<:Integer})
res = zeros(Rational{eltype(M)}, size(M))
res += I
k = 0
expM = one(M)
while !iszero(expM)
k += 1
expM *= M
@. res += 1 // factorial(k) * expM
if k == 20
@warn "matrix exponential did not converge" norm(expM - exp(M))
break
end
end
@debug "matrix_exp converged after $k iterations"
return res
end
const gap_adj_mats = [
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -2],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1],
[2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0],
[3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 1],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2],
[0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1],
[0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 0],
[0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0],
[0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -1],
[0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
],
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0],
],
]
function G₂_matrices_roots()
adj_mats = map(gap_adj_mats) do m
return hcat(m...)
end
adj_mats = filter!(!isdiag, adj_mats) # remove the ones from center
gens_mats = [convert(Matrix{Int}, matrix_exp(m')) for m in adj_mats]
#=
The roots from
G₂roots_gap = [
[2, -1], # α = e₁ - e₂
[-3, 2], # A = -α + β = -e₁ + 2e₂ - e₃
[-1, 1], # β = e₂ - e₃
[1, 0], # α + β = e₁ - e₃
[3, -1], # B = 2α + β = 2e₁ - e₂ - e₃
[0, 1], # A + B = α + 2β = e₁ + e₂ - 2e₃
[-2, 1], # -α
[3, -2], # -A
[1, -1], # -β
[-1, 0], # -α - β
[-3, 1], # -B
[0, -1], # -A - B
]
G₂roots_gap are the ones from cartan matrix. To obtain the standard
(hexagonal) picture map them by `T` defined as follows:
```julia
cartan = hcat(G₂roots_gap[1:2]...)
rot(α) = [cos(α) -sin(α); sin(α) cos(α)]
c₁ = [√2, 0]
c₂ = rot(5π / 6) * [√2, 0] * √3 # (= 1/2[√6, 1])
T = hcat(c₁, c₂) * inv(cartan)
```
By plotting one against the others (or by blind calculation) one can see
the following assignment. Here `⟨α, β⟩_ℤ = A₂` and `⟨A, B⟩_ℤ ≅ √3/√2 A₂`.
=#
e₁ = PropertyT.Roots.𝕖(3, 1)
e₂ = PropertyT.Roots.𝕖(3, 2)
e₃ = PropertyT.Roots.𝕖(3, 3)
α = e₁ - e₂
β = e₂ - e₃
A = -α + β
B = α + (α + β)
roots = [α, A, β, α + β, B, A + B, -α, -A, -β, -α - β, -B, -A - B]
return gens_mats, roots
end
function G₂_roots_weyl()
(mats, roots) = G₂_matrices_roots()
d = size(first(mats), 1)
G₂ = Groups.MatrixGroup{d}(mats)
m = Groups.gens(G₂)
σ = let w = m[1] * inv(m[7]) * m[1], m = union(m, inv.(m))
PermutationGroups.Perm([findfirst(==(inv(w) * x * w), m) for x in m])
end
τ = let w = m[2] * inv(m[8]) * m[2], m = union(m, inv.(m))
PermutationGroups.Perm([findfirst(==(inv(w) * x * w), m) for x in m])
end
W = PermGroup(σ, τ)
return G₂, roots, W
end

98
scripts/G₂_has_T.jl
View File

@ -1,98 +0,0 @@
using LinearAlgebra
BLAS.set_num_threads(8)
ENV["OMP_NUM_THREADS"] = 4
using Groups
import Groups.MatrixGroups
include(joinpath(@__DIR__, "../test/optimizers.jl"))
using PropertyT
using PropertyT.SymbolicWedderburn
using PropertyT.PermutationGroups
using PropertyT.StarAlgebras
include(joinpath(@__DIR__, "argparse.jl"))
include(joinpath(@__DIR__, "utils.jl"))
# const N = parsed_args["N"]
const HALFRADIUS = parsed_args["halfradius"]
const UPPER_BOUND = parsed_args["upper_bound"]
include(joinpath(@__DIR__, "./G₂_gens.jl"))
G, roots, Weyl = G₂_roots_weyl()
@info "Running Δ² - λ·Δ sum of squares decomposition for G₂"
@info "computing group algebra structure"
RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
@info "computing WedderburnDecomposition"
wd = let Σ = Weyl, RG = RG
act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
Dict(g => PermutationGroups.perm(g) for g in Σ),
)
@time SymbolicWedderburn.WedderburnDecomposition(
Float64,
Σ,
act,
basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
semisimple = false,
)
end
@info wd
Δ = RG(length(S)) - sum(RG(s) for s in S)
elt = Δ^2
unit = Δ
@time model, varP = PropertyT.sos_problem_primal(
elt,
unit,
wd;
upper_bound = UPPER_BOUND,
augmented = true,
show_progress = false,
)
warm = nothing
status = JuMP.OPTIMIZE_NOT_CALLED
certified, λ = false, nothing
while status ≠ JuMP.OPTIMAL
@time status, warm = PropertyT.solve(
model,
scs_optimizer(;
eps = 1e-10,
max_iters = 20_000,
accel = 50,
alpha = 1.95,
),
warm,
)
@info "reconstructing the solution"
Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP]
Qs = real.(sqrt.(Ps))
PropertyT.reconstruct(Qs, wd)
end
@info "certifying the solution"
@time certified, λ = PropertyT.certify_solution(
elt,
unit,
JuMP.objective_value(model),
Q;
halfradius = HALFRADIUS,
augmented = true,
)
end
if certified && λ > 0
Κ(λ, S) = round(sqrt(2λ / length(S)), Base.RoundDown; digits = 5)
@info "Certified result: G₂ has property (T):" N λ Κ(λ, S)
else
@info "Could NOT certify the result:" certified λ
end

30
scripts/SpN_Adj.jl
View File

@ -24,7 +24,8 @@ const GENUS = 2N
G = MatrixGroups.SymplecticGroup{GENUS}(Int8) G = MatrixGroups.SymplecticGroup{GENUS}(Int8)
RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS) RG, S, sizes =
@time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true)
wd = let RG = RG, N = N wd = let RG = RG, N = N
G = StarAlgebras.object(RG) G = StarAlgebras.object(RG)
@ -41,8 +42,6 @@ wd = let RG = RG, N = N
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
) )
@info wdfl
wdfl
end end
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -59,22 +58,23 @@ unit = Δ
@time model, varP = PropertyT.sos_problem_primal( @time model, varP = PropertyT.sos_problem_primal(
elt, elt,
unit, unit,
wd; wd,
upper_bound = UPPER_BOUND, upper_bound=UPPER_BOUND,
augmented = true, augmented=true,
show_progress = true, show_progress=true
) )
solve_in_loop( solve_in_loop(
model, model,
wd, wd,
varP; varP,
logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-$(UPPER_BOUND)Δ", logdir="./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ",
optimizer = cosmo_optimizer(; optimizer=cosmo_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 20_000, max_iters=20_000,
accel = 50, accel=50,
alpha = 1.95, alpha=1.95,
), ),
data = (elt = elt, unit = unit, halfradius = HALFRADIUS), data=(elt=elt, unit=unit, halfradius=HALFRADIUS)
) )

47
scripts/utils.jl
View File

@ -11,19 +11,12 @@ function get_solution(model)
return solution return solution
end end
function get_solution(model, wd, varP, eps = 1e-10) function get_solution(model, wd, varP; logdir)
λ = JuMP.value(model[:λ]) λ = JuMP.value(model[:λ])
@info "reconstructing the solution" Qs = [real.(sqrt(JuMP.value.(P))) for P in varP]
Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP], eps = eps Q = PropertyT.reconstruct(Qs, wd)
PropertyT.__droptol!.(Ps, 100eps)
Qs = real.(sqrt.(Ps))
PropertyT.__droptol!.(Qs, eps)
PropertyT.reconstruct(Qs, wd)
end
solution = Dict(:λ => λ, :Q => Q) solution = Dict(:λ => λ, :Q => Q)
return solution return solution
end end
@ -49,23 +42,22 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
# logstream = current_logger().logger.stream # logstream = current_logger().logger.stream
# v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint # v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint
# @warn v # @warn v
status, warm = status, warm = @time PropertyT.solve(log_file, model, optimizer, warm)
@time PropertyT.solve(log_file, model, optimizer, warm)
solution = get_solution(model, args...; logdir=logdir)
solution[:warm] = warm
solution = get_solution(model, args...)
serialize(joinpath(logdir, "solution_$date.sjl"), solution) serialize(joinpath(logdir, "solution_$date.sjl"), solution)
serialize(joinpath(logdir, "solution.sjl"), solution) serialize(joinpath(logdir, "solution.sjl"), solution)
solution[:warm] = warm flag, λ_cert = open(log_file, append=true) do io
flag, λ_cert = open(log_file; append = true) do io
with_logger(SimpleLogger(io)) do with_logger(SimpleLogger(io)) do
return PropertyT.certify_solution( PropertyT.certify_solution(
data.elt, data.elt,
data.unit, data.unit,
solution[:λ], solution[:λ],
solution[:Q]; solution[:Q],
halfradius = data.halfradius, halfradius=data.halfradius,
) )
end end
end end
@ -74,21 +66,20 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
end end
if flag == true && certified_λ ≥ 0 if flag == true && certified_λ ≥ 0
@info "Certification done with λ = $certified_λ" certified_λ status @info "Certification done with λ = $certified_λ"
return certified_λ return certified_λ
else else
rel_change = rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda))
abs(certified_λ - old_lambda) / @info "Certification failed with λ = $λ" certified_λ rel_change
(abs(certified_λ) + abs(old_lambda)) end
@info "Certification failed with λ = $λ" certified_λ rel_change status
old_lambda = certified_λ
if rel_change < 1e-9 if rel_change < 1e-9
@info "No progress detected, breaking" certified_λ rel_change status @info "No progress detected, breaking"
break break
end end
end end
old_lambda = certified_λ
end
return status == JuMP.OPTIMAL ? old_lambda : NaN return status == JuMP.OPTIMAL ? old_lambda : NaN
end end

11
src/PropertyT.jl
View File

@ -3,6 +3,7 @@ module PropertyT
using LinearAlgebra using LinearAlgebra
using SparseArrays using SparseArrays
using IntervalArithmetic
using JuMP using JuMP
using Groups using Groups
@ -25,17 +26,17 @@ include("gradings.jl")
include("actions/actions.jl") include("actions/actions.jl")
function group_algebra(G::Groups.Group, S = gens(G); halfradius::Integer) function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool)
S = union!(S, inv.(S)) S = union!(S, inv.(S))
@info "generating wl-metric ball of radius $(2halfradius)" @info "generating wl-metric ball of radius $(2halfradius)"
@time E, sizes = Groups.wlmetric_ball(S; radius = 2halfradius) @time E, sizes = Groups.wlmetric_ball(S, radius=2halfradius)
@info "sizes = $(sizes)" @info "sizes = $(sizes)"
@info "computing the *-algebra structure for G" @info "computing the *-algebra structure for G"
@time RG = StarAlgebras.StarAlgebra( @time RG = StarAlgebras.StarAlgebra{twisted}(
G, G,
StarAlgebras.Basis{UInt32}(E), StarAlgebras.Basis{UInt32}(E),
(sizes[halfradius], sizes[halfradius]); (sizes[halfradius], sizes[halfradius]),
precompute = false, precompute=false,
) )
return RG, S, sizes return RG, S, sizes
end end

3
src/actions/actions.jl
View File

@ -1,4 +1,5 @@
import SymbolicWedderburn.action import SymbolicWedderburn.action
StarAlgebras.star(g::Groups.GroupElement) = inv(g)
include("alphabet_permutation.jl") include("alphabet_permutation.jl")
@ -9,7 +10,7 @@ include("autfn_conjugation.jl")
function SymbolicWedderburn.action( function SymbolicWedderburn.action(
act::SymbolicWedderburn.ByPermutations, act::SymbolicWedderburn.ByPermutations,
g::Groups.GroupElement, g::Groups.GroupElement,
α::StarAlgebras.AlgebraElement, α::StarAlgebras.AlgebraElement
) )
res = StarAlgebras.zero!(similar(α)) res = StarAlgebras.zero!(similar(α))
B = basis(parent(α)) B = basis(parent(α))

80
src/certify.jl
View File

@ -1,6 +1,3 @@
import IntervalArithmetic
import IntervalMatrices
function augment_columns!(Q::AbstractMatrix) function augment_columns!(Q::AbstractMatrix)
for c in eachcol(Q) for c in eachcol(Q)
c .-= sum(c) ./ length(c) c .-= sum(c) ./ length(c)
@ -11,25 +8,26 @@ end
function __sos_via_sqr!( function __sos_via_sqr!(
res::StarAlgebras.AlgebraElement, res::StarAlgebras.AlgebraElement,
P::AbstractMatrix; P::AbstractMatrix;
augmented::Bool, augmented::Bool
id = (b = basis(parent(res)); b[one(first(b))]),
) )
A = parent(res)
mstr = A.mstructure
@assert size(mstr) == size(P)
StarAlgebras.zero!(res) StarAlgebras.zero!(res)
for j in axes(mstr, 2) A = parent(res)
for i in axes(mstr, 1) b = basis(A)
@assert size(A.mstructure) == size(P)
e = b[one(b[1])]
for i in axes(A.mstructure, 1)
x = StarAlgebras._istwisted(A.mstructure) ? StarAlgebras.star(b[i]) : b[i]
for j in axes(A.mstructure, 2)
p = P[i, j] p = P[i, j]
x_star_y = mstr[-i, j] xy = b[A.mstructure[i, j]]
res[x_star_y] += p # either result += P[x,y]*(x*y)
# either result += P[x,y]*(x'*y) res[xy] += p
if augmented if augmented
# or result += P[x,y]*(1-x)'*(1-y) == P[x,y]*(1-x'-y+x'y) # or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy)
res[id] += p y = b[j]
x_star, y = mstr[-i, id], j res[e] += p
res[x_star] -= p res[x] -= p
res[y] -= p res[y] -= p
end end
end end
@ -38,11 +36,7 @@ function __sos_via_sqr!(
return res return res
end end
function __sos_via_cnstr!( function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix, cnstrs)
res::StarAlgebras.AlgebraElement,
Q²::AbstractMatrix,
cnstrs,
)
StarAlgebras.zero!(res) StarAlgebras.zero!(res)
for (g, A_g) in cnstrs for (g, A_g) in cnstrs
res[g] = dot(A_g, Q²) res[g] = dot(A_g, Q²)
@ -50,14 +44,10 @@ function __sos_via_cnstr!(
return res return res
end end
function compute_sos( function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool)
A::StarAlgebras.StarAlgebra,
Q::AbstractMatrix;
augmented::Bool,
)
Q² = Q' * Q Q² = Q' * Q
res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A) res = StarAlgebras.AlgebraElement(zeros(eltype(Q²), length(basis(A))), A)
res = __sos_via_sqr!(res, Q²; augmented = augmented) res = __sos_via_sqr!(res, Q², augmented=augmented)
return res return res
end end
@ -66,7 +56,7 @@ function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
eq_sign = let T = eltype(residual) eq_sign = let T = eltype(residual)
if T <: IntervalArithmetic.Interval if T <: Interval
"∈" "∈"
elseif T <: Union{Rational,Integer} elseif T <: Union{Rational,Integer}
"=" "="
@ -91,12 +81,13 @@ function sufficient_λ(
order_unit::StarAlgebras.AlgebraElement, order_unit::StarAlgebras.AlgebraElement,
λ, λ,
sos::StarAlgebras.AlgebraElement; sos::StarAlgebras.AlgebraElement;
halfradius, halfradius
) )
@assert parent(elt) === parent(order_unit) == parent(sos) @assert parent(elt) === parent(order_unit) == parent(sos)
residual = (elt - λ * order_unit) - sos residual = (elt - λ * order_unit) - sos
return sufficient_λ(residual, λ; halfradius = halfradius) return sufficient_λ(residual, λ; halfradius=halfradius)
end end
function certify_solution( function certify_solution(
@ -105,27 +96,24 @@ function certify_solution(
λ, λ,
Q::AbstractMatrix{<:AbstractFloat}; Q::AbstractMatrix{<:AbstractFloat};
halfradius, halfradius,
augmented = iszero(StarAlgebras.aug(elt)) && augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit))
iszero(StarAlgebras.aug(orderunit)),
) )
should_we_augment =
!augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0 should_we_augment = !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0
Q = should_we_augment ? augment_columns!(Q) : Q Q = should_we_augment ? augment_columns!(Q) : Q
@time sos = compute_sos(parent(elt), Q; augmented = augmented) @time sos = compute_sos(parent(elt), Q, augmented=augmented)
@info "Checking in $(eltype(sos)) arithmetic with" λ @info "Checking in $(eltype(sos)) arithmetic with" λ
λ_flpoint = sufficient_λ(elt, orderunit, λ, sos; halfradius = halfradius) λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius)
if λ_flpoint ≤ 0 if λ_flpoint ≤ 0
return false, λ_flpoint return false, λ_flpoint
end end
λ_int = IntervalArithmetic.@interval(λ) λ_int = @interval(λ)
Q_int = IntervalMatrices.IntervalMatrix([ Q_int = [@interval(q) for q in Q]
IntervalArithmetic.@interval(q) for q in Q
])
check, sos_int = @time if should_we_augment check, sos_int = @time if should_we_augment
@info("Projecting columns of Q to the augmentation ideal...") @info("Projecting columns of Q to the augmentation ideal...")
@ -135,16 +123,16 @@ function certify_solution(
check_augmented || @error( check_augmented || @error(
"Augmentation failed! The following numbers are not certified!" "Augmentation failed! The following numbers are not certified!"
) )
sos_int = compute_sos(parent(elt), Q_int; augmented = augmented) sos_int = compute_sos(parent(elt), Q_int; augmented=augmented)
check_augmented, sos_int check_augmented, sos_int
else else
true, compute_sos(parent(elt), Q_int; augmented = augmented) true, compute_sos(parent(elt), Q_int, augmented=augmented)
end end
@info "Checking in $(eltype(sos_int)) arithmetic with" λ_int @info "Checking in $(eltype(sos_int)) arithmetic with" λ_int
λ_certified = λ_certified =
sufficient_λ(elt, orderunit, λ_int, sos_int; halfradius = halfradius) sufficient_λ(elt, orderunit, λ_int, sos_int, halfradius=halfradius)
return check && IntervalArithmetic.inf(λ_certified) > 0.0, λ_certified return check && inf(λ_certified) > 0.0, inf(λ_certified)
end end

96
src/constraint_matrix.jl
View File

@ -28,12 +28,7 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
size::Tuple{Int,Int} size::Tuple{Int,Int}
val::T val::T
function ConstraintMatrix{T}( function ConstraintMatrix{T}(nzeros::AbstractArray{<:Integer}, n, m, val) where {T}
nzeros::AbstractArray{<:Integer},
n,
m,
val,
) where {T}
@assert n ≥ 1 @assert n ≥ 1
@assert m ≥ 1 @assert m ≥ 1
@ -50,26 +45,15 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
end end
end end
function ConstraintMatrix( ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} =
nzeros::AbstractArray{<:Integer}, ConstraintMatrix{T}(nzeros, n, m, val)
n,
m,
val::T,
) where {T}
return ConstraintMatrix{T}(nzeros, n, m, val)
end
Base.size(cm::ConstraintMatrix) = cm.size Base.size(cm::ConstraintMatrix) = cm.size
function __get_positive(cm::ConstraintMatrix, idx::Integer) __get_positive(cm::ConstraintMatrix, idx::Integer) =
return convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx))) convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
end __get_negative(cm::ConstraintMatrix, idx::Integer) =
function __get_negative(cm::ConstraintMatrix, idx::Integer) convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx)))
return convert(
eltype(cm),
cm.val * length(searchsorted(cm.neg, idx; rev = true)),
)
end
Base.@propagate_inbounds function Base.getindex( Base.@propagate_inbounds function Base.getindex(
cm::ConstraintMatrix, cm::ConstraintMatrix,
@ -83,29 +67,26 @@ Base.@propagate_inbounds function Base.getindex(
return pos - neg return pos - neg
end end
struct NZPairsIter{T,I} struct NZPairsIter{T}
m::ConstraintMatrix{T,I} m::ConstraintMatrix{T}
end end
Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T} Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown() Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
# TODO: iterate over (idx=>val) pairs combining vals # TODO: iterate over (idx=>val) pairs combining vals
function Base.iterate( function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 1))
itr::NZPairsIter,
state::Tuple{Int,Nothing} = (1, nothing),
)
k = iterate(itr.m.pos, state[1]) k = iterate(itr.m.pos, state[1])
isnothing(k) && return iterate(itr, (nothing, 1)) isnothing(k) && return iterate(itr, state[2])
idx, st = k idx, st = k
return idx => itr.m.val, (st, nothing) return idx => itr.m.val, (st, 1)
end end
function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int}) function Base.iterate(itr::NZPairsIter, state::Int)
k = iterate(itr.m.neg, state[2]) k = iterate(itr.m.neg, state[1])
isnothing(k) && return nothing isnothing(k) && return nothing
idx, st = k idx, st = k
return idx => -itr.m.val, (nothing, st) return idx => -itr.m.val, st
end end
""" """
@ -146,50 +127,3 @@ function LinearAlgebra.dot(cm::ConstraintMatrix, m::AbstractMatrix{T}) where {T}
neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg]) neg = isempty(cm.neg) ? zero(first(m)) : sum(@view m[cm.neg])
return convert(eltype(cm), cm.val) * (pos - neg) return convert(eltype(cm), cm.val) * (pos - neg)
end end
function constraints(A::StarAlgebras.StarAlgebra; augmented::Bool)
return constraints(basis(A), A.mstructure; augmented = augmented)
end
function constraints(
basis::StarAlgebras.AbstractBasis,
mstr::StarAlgebras.MultiplicativeStructure;
augmented = false,
)
cnstrs = _constraints(
mstr;
augmented = augmented,
num_constraints = length(basis),
id = basis[one(first(basis))],
)
return Dict(
basis[i] => ConstraintMatrix(c, size(mstr)..., 1) for
(i, c) in pairs(cnstrs)
)
end
function _constraints(
mstr::StarAlgebras.MultiplicativeStructure;
augmented::Bool = false,
num_constraints = maximum(mstr),
id,
)
cnstrs = [signed(eltype(mstr))[] for _ in 1:num_constraints]
LI = LinearIndices(size(mstr))
for ci in CartesianIndices(size(mstr))
k = LI[ci]
i, j = Tuple(ci)
a_star_b = mstr[-i, j]
push!(cnstrs[a_star_b], k)
if augmented
# (1-a)'(1-b) = 1 - a' - b + a'b
push!(cnstrs[id], k)
a_star, b = mstr[-i, id], j
push!(cnstrs[a_star], -k)
push!(cnstrs[b], -k)
end
end
return cnstrs
end

29
src/gradings.jl
View File

@ -1,8 +1,7 @@
## something about roots ## something about roots
function Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
return Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j) Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
end
function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N} function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
if s.symbol === :A if s.symbol === :A
@ -52,28 +51,20 @@ function Adj(rootsystem::AbstractDict, subtype::Symbol)
+, +,
( (
Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if
Roots.classify_sub_root_system(roots, first(α), first(β)) == subtype Roots.classify_sub_root_system(
); roots,
init = zero(first(values(rootsystem))), first(α),
) first(β),
end ) == subtype
),
Adj(rootsystem::AbstractDict) = sum(values(rootsystem))^2 - Sq(rootsystem) init=zero(first(values(rootsystem))),
function Sq(rootsystem::AbstractDict)
return reduce(
+,
Δα^2 for (_, Δα) in rootsystem;
init = zero(first(values(rootsystem))),
) )
end end
function level(rootsystem, level::Integer) function level(rootsystem, level::Integer)
1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}") 1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}")
level == 1 && return Adj(rootsystem, :C₁) # always positive level == 1 && return Adj(rootsystem, :C₁) # always positive
level == 2 && return Adj(rootsystem, :A₁) + level == 2 && return Adj(rootsystem, :A₁) + Adj(rootsystem, Symbol("C₁×C₁")) + Adj(rootsystem, :C₂) # C₂ is not positive
Adj(rootsystem, Symbol("C₁×C₁")) +
Adj(rootsystem, :C₂) # C₂ is not positive
level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁")) level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁"))
level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive
end end

33
src/reconstruct.jl
View File

@ -7,28 +7,33 @@ end
function reconstruct( function reconstruct(
Ms::AbstractVector{<:AbstractMatrix}, Ms::AbstractVector{<:AbstractMatrix},
wbdec::WedderburnDecomposition, wbdec::WedderburnDecomposition;
atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
) )
n = __outer_dim(wbdec) n = __outer_dim(wbdec)
res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds) res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
res = similar(M, n, n) res = similar(M, n, n)
res = _reconstruct!(res, M, ds) reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol)
return res
end end
res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
return res return res
end end
function _reconstruct!( function reconstruct!(
res::AbstractMatrix, res::AbstractMatrix,
M::AbstractMatrix, M::AbstractMatrix,
ds::SymbolicWedderburn.DirectSummand, ds::SymbolicWedderburn.DirectSummand,
G,
hom;
atol=eps(real(eltype(ds))) * 10max(size(res)...)
) )
res .= zero(eltype(res)) res .= zero(eltype(res))
if !iszero(M)
U = SymbolicWedderburn.image_basis(ds) U = SymbolicWedderburn.image_basis(ds)
d = SymbolicWedderburn.degree(ds) d = SymbolicWedderburn.degree(ds)
res = (U' * M * U) .* d tmp = (U' * M * U) .* d
res = average!(res, tmp, G, hom)
if eltype(res) <: AbstractFloat
__droptol!(res, atol) # TODO: is this really necessary?!
end end
return res return res
end end
@ -47,22 +52,18 @@ function average!(
res::AbstractMatrix, res::AbstractMatrix,
M::AbstractMatrix, M::AbstractMatrix,
G::Groups.Group, G::Groups.Group,
hom::SymbolicWedderburn.InducedActionHomomorphism{ hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations}
<:SymbolicWedderburn.ByPermutations,
},
) )
res .= zero(eltype(res))
@assert size(M) == size(res) @assert size(M) == size(res)
o = Groups.order(Int, G)
for g in G for g in G
p = SymbolicWedderburn.induce(hom, g) gext = SymbolicWedderburn.induce(hom, g)
Threads.@threads for c in axes(res, 2) Threads.@threads for c in axes(res, 2)
for r in axes(res, 1) for r in axes(res, 1)
if !iszero(M[r, c]) res[r, c] += M[r^gext, c^gext]
res[r^p, c^p] += M[r, c] / o
end
end end
end end
end end
o = Groups.order(Int, G)
res ./= o
return res return res
end end

214
src/roots.jl
View File

@ -5,104 +5,111 @@ using LinearAlgebra
export Root, isproportional, isorthogonal, ~, ⟂ export Root, isproportional, isorthogonal, ~, ⟂
abstract type AbstractRoot{N,T} end # <: AbstractVector{T} ? abstract type AbstractRoot{N,T} end
ℓ₂length(r::AbstractRoot) = norm(r, 2) struct Root{N,T} <: AbstractRoot{N,T}
ambient_dim(r::AbstractRoot) = length(r) coord::SVector{N,T}
Base.:*(r::AbstractRoot, a::Number) = a * r 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)) cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
function isproportional(α::AbstractRoot, β::AbstractRoot) function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
ambient_dim(α) == ambient_dim(β) || return false N == M || return false
val = abs(cos_angle(α, β)) val = abs(cos_angle(α, β))
return isapprox(val, one(val); atol = eps(one(val))) return isapprox(val, one(val), atol=eps(one(val)))
end end
function isorthogonal(α::AbstractRoot, β::AbstractRoot) function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
ambient_dim(α) == ambient_dim(β) || return false N == M || return false
val = cos_angle(α, β) val = cos_angle(α, β)
return isapprox(val, zero(val); atol = eps(one(val))) return isapprox(val, zero(val), atol=eps(one(val)))
end end
function positive(roots::AbstractVector{<:AbstractRoot}) function _positive_direction(α::Root{N}) where {N}
isempty(roots) && return empty(roots) v = α.coord + 1 / (N * 100) * rand(N)
return Root{N,Float64}(v / norm(v, 2))
end
function positive(roots::AbstractVector{<:Root{N}}) where {N}
# return those roots for which dot(α, Root([½, ¼, …])) > 0.0
pd = _positive_direction(first(roots)) pd = _positive_direction(first(roots))
return filter(α -> dot(α, pd) > 0.0, roots) return filter(α -> dot(α, pd) > 0.0, roots)
end end
function Base.show(io::IO, r::AbstractRoot) function Base.show(io::IO, r::Root)
return print(io, "Root $(r.coord)") print(io, "Root$(r.coord)")
end end
function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot) function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
l₂l = ℓ₂length(r) lngth² = sum(x -> x^2, r.coord)
l = round(Int, l₂l) ≈ l₂l ? "$(round(Int, l₂l))" : "√$(round(Int, l₂l^2))" l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
return print(io, "Root in ℝ^$(length(r)) of length $l\n", r.coord) print(io, "Root in ℝ^$N of length $l\n", r.coord)
end end
function reflection(α::AbstractRoot, β::AbstractRoot) 𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
return β - Int(2dot(α, β) // dot(α, α)) * α 𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
end
function cartan(α::AbstractRoot, β::AbstractRoot)
ambient_dim(α) == ambient_dim(β) || throw("incompatible ambient dimensions")
return [
ℓ₂length(reflection(a, b) - b) / ℓ₂length(a) for a in (α, β),
b in (α, β)
]
end
""" """
classify_root_system(α, β) classify_root_system(α, β)
Return the symbol of smallest system generated by roots `α` and `β`. Return the symbol of smallest system generated by roots `α` and `β`.
The classification is based only on roots length, The classification is based only on roots length and
proportionality/orthogonality and Cartan matrix. proportionality/orthogonality.
""" """
function classify_root_system( function classify_root_system(α::AbstractRoot, β::AbstractRoot)
α::AbstractRoot, lα, lβ = length(α), length(β)
β::AbstractRoot,
long::Tuple{Bool,Bool},
)
if isproportional(α, β) if isproportional(α, β)
if all(long) if lα ≈ lβ ≈ √2
return :C₁
elseif all(.!long) # both short
return :A₁ return :A₁
elseif lα ≈ lβ ≈ 2.0
return :C₁
else else
@error "Proportional roots of different length"
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end end
elseif isorthogonal(α, β) elseif isorthogonal(α, β)
if all(long) if lα ≈ lβ ≈ √2
return Symbol("C₁×C₁")
elseif all(.!long) # both short
return Symbol("A₁×A₁") return Symbol("A₁×A₁")
elseif any(long) elseif lα ≈ lβ ≈ 2.0
return Symbol("C₁×C₁")
elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
return Symbol("A₁×C₁") return Symbol("A₁×C₁")
else
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end end
else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
a, b, c, d = abs.(cartan(α, β)) if lα ≈ lβ ≈ √2
@assert a == d == 2
b, c = b < c ? (b, c) : (c, b)
if b == c == 1
return :A₂ return :A₂
elseif b == 1 && c == 2 elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
return :C₂ return :C₂
elseif b == 1 && c == 3
@warn ":G₂? really?"
return :G₂
else else
@error a, b, c, d
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β") error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
end end
end end
end end
function proportional_root_from_system( function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
Ω::AbstractVector{<:AbstractRoot},
α::AbstractRoot,
)
k = findfirst(v -> isproportional(α, v), Ω) k = findfirst(v -> isproportional(α, v), Ω)
if isnothing(k) if isnothing(k)
error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω") error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
@ -110,31 +117,25 @@ function proportional_root_from_system(
return Ω[k] return Ω[k]
end end
struct Plane{R<:AbstractRoot} struct Plane{R<:Root}
v1::R v1::R
v2::R v2::R
vectors::Vector{R} vectors::Vector{R}
end end
function Plane(α::AbstractRoot, β::AbstractRoot) Plane(α::R, β::R) where {R<:Root} =
return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3]) Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
end
function Base.in(r::AbstractRoot, plane::Plane) function Base.in(r::R, plane::Plane{R}) where {R}
return any(isproportional(r, v) for v in plane.vectors) return any(isproportional(r, v) for v in plane.vectors)
end end
function _islong(α::AbstractRoot, Ω)
lα = ℓ₂length(α)
return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
end
function classify_sub_root_system( function classify_sub_root_system(
Ω::AbstractVector{<:AbstractRoot{N}}, Ω::AbstractVector{<:Root{N}},
α::AbstractRoot{N}, α::Root{N},
β::AbstractRoot{N}, β::Root{N},
) where {N} ) where {N}
@assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
v = proportional_root_from_system(Ω, α) v = proportional_root_from_system(Ω, α)
w = proportional_root_from_system(Ω, β) w = proportional_root_from_system(Ω, β)
@ -145,75 +146,32 @@ function classify_sub_root_system(
l = length(subsystem) l = length(subsystem)
if l == 1 if l == 1
x = first(subsystem) x = first(subsystem)
long = _islong(x, Ω) return classify_root_system(x, x)
return classify_root_system(x, -x, (long, long))
elseif l == 2 elseif l == 2
x, y = subsystem return classify_root_system(subsystem...)
return classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω)))
elseif l == 3 elseif l == 3
x, y, z = subsystem a = classify_root_system(subsystem[1], subsystem[2])
l1, l2, l3 = _islong(x, Ω), _islong(y, Ω), _islong(z, Ω) b = classify_root_system(subsystem[2], subsystem[3])
a = classify_root_system(x, y, (l1, l2)) c = classify_root_system(subsystem[1], subsystem[3])
b = classify_root_system(y, z, (l2, l3))
c = classify_root_system(x, z, (l1, l3))
if :A₂ == a == b == c # it's only A₂ if a == b == c # it's only A₂
return a return a
end end
throw("Unknown subroot system! $((x,y,z))") C = (:C₂, Symbol("C₁×C₁"))
elseif l == 4 if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
subtypes = [
classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω))) for
x in subsystem for y in subsystem if x ≠ y
]
if :C₂ in subtypes
return :C₂ return :C₂
end 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 end
@error "Unknown root subsystem generated by" α β @error "Unknown root subsystem generated by" α β
throw("Unknown root system: $subsystem") throw("Unknown root system: $subsystem")
end end
## concrete implementation:
struct Root{N,T} <: AbstractRoot{N,T}
coord::SVector{N,T}
end
Root(a) = Root(SVector(a...))
# convienience constructors
𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
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))
function Base.:+(r::Root, s::Root)
ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
return Root(r.coord + s.coord)
end
function Base.:-(r::Root, s::Root)
ambient_dim(r) == ambient_dim(s) || throw("incompatible ambient dimensions")
return Root(r.coord - s.coord)
end
Base.:-(r::Root) = Root(-r.coord)
Base.:*(a::Number, r::Root) = Root(a * r.coord)
Base.length(r::Root) = length(r.coord)
LinearAlgebra.norm(r::Root, p::Real = 2) = norm(r.coord, p)
LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
function _positive_direction(α::Root{N}) where {N}
v = α.coord + 1 / (N * 100) * rand(N)
return Root{N,Float64}(v / norm(v, 2))
end
end # of module Roots end # of module Roots

15
src/solve.jl
View File

@ -7,16 +7,13 @@ function setwarmstart!(model::JuMP.Model, warmstart)
ct => JuMP.all_constraints(model, ct...) for ct => JuMP.all_constraints(model, ct...) for
ct in JuMP.list_of_constraint_types(model) ct in JuMP.list_of_constraint_types(model)
) )
try
JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal) JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
for (ct, idx) in pairs(constraint_map) for (ct, idx) in pairs(constraint_map)
JuMP.set_start_value.(idx, warmstart.slack[ct]) JuMP.set_start_value.(idx, warmstart.slack[ct])
JuMP.set_dual_start_value.(idx, warmstart.dual[ct]) JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
end end
catch e
@warn "Warmstarting was not succesfull!" e
end
return model return model
end end
@ -31,10 +28,11 @@ function getwarmstart(model::JuMP.Model)
slack = Dict(k => value.(v) for (k, v) in constraint_map) slack = Dict(k => value.(v) for (k, v) in constraint_map)
duals = Dict(k => JuMP.dual.(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) return (primal=primal, dual=duals, slack=slack)
end end
function solve(m::JuMP.Model, optimizer, warmstart = nothing) function solve(m::JuMP.Model, optimizer, warmstart=nothing)
JuMP.set_optimizer(m, optimizer) JuMP.set_optimizer(m, optimizer)
MOIU.attach_optimizer(m) MOIU.attach_optimizer(m)
@ -47,7 +45,8 @@ function solve(m::JuMP.Model, optimizer, warmstart = nothing)
return status, getwarmstart(m) return status, getwarmstart(m)
end end
function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart = nothing) function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart=nothing)
isdir(dirname(solverlog)) || mkpath(dirname(solverlog)) isdir(dirname(solverlog)) || mkpath(dirname(solverlog))
Base.flush(Base.stdout) Base.flush(Base.stdout)
@ -56,7 +55,7 @@ function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart = nothing)
redirect_stdout(logfile) do redirect_stdout(logfile) do
status, warmstart = solve(m, optimizer, warmstart) status, warmstart = solve(m, optimizer, warmstart)
Base.Libc.flush_cstdio() Base.Libc.flush_cstdio()
return status, warmstart status, warmstart
end end
end end

205
src/sos_sdps.jl
View File

@ -6,13 +6,15 @@ See also [sos_problem_primal](@ref).
""" """
function sos_problem_dual( function sos_problem_dual(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement = zero(elt); order_unit::StarAlgebras.AlgebraElement=zero(elt);
lower_bound = -Inf, lower_bound=-Inf
) )
@assert parent(elt) == parent(order_unit) @assert parent(elt) == parent(order_unit)
algebra = parent(elt) algebra = parent(elt)
moment_matrix = let m = algebra.mstructure mstructure = if StarAlgebras._istwisted(algebra.mstructure)
[m[-i, j] for i in axes(m, 1), j in axes(m, 2)] algebra.mstructure
else
StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
end end
# 1 variable for every primal constraint # 1 variable for every primal constraint
@ -20,21 +22,60 @@ function sos_problem_dual(
# Symmetrized: # Symmetrized:
# 1 dual variable for every orbit of G acting on basis # 1 dual variable for every orbit of G acting on basis
model = Model() model = Model()
JuMP.@variable(model, y[1:length(basis(algebra))]) @variable model y[1:length(basis(algebra))]
JuMP.@constraint(model, λ_dual, dot(order_unit, y) == 1) @constraint model λ_dual dot(order_unit, y) == 1
JuMP.@constraint(model, psd, y[moment_matrix] in PSDCone()) @constraint(model, psd, y[mstructure] in PSDCone())
if !isinf(lower_bound) if !isinf(lower_bound)
throw("Not Implemented yet") throw("Not Implemented yet")
JuMP.@variable(model, λ_ub_dual) @variable model λ_ub_dual
JuMP.@objective(model, Min, dot(elt, y) + lower_bound * λ_ub_dual) @objective model Min dot(elt, y) + lower_bound * λ_ub_dual
else else
JuMP.@objective(model, Min, dot(elt, y)) @objective model Min dot(elt, y)
end end
return model return model
end end
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[StarAlgebras.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::StarAlgebras.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]) sos_problem_primal(X, [u = zero(X); upper_bound=Inf])
Formulate sum of squares decomposition problem for `X - λ·u`. Formulate sum of squares decomposition problem for `X - λ·u`.
@ -55,10 +96,9 @@ The default `u = zero(X)` formulates a simple feasibility problem.
""" """
function sos_problem_primal( function sos_problem_primal(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
order_unit::StarAlgebras.AlgebraElement = zero(elt); order_unit::StarAlgebras.AlgebraElement=zero(elt);
upper_bound = Inf, upper_bound=Inf,
augmented::Bool = iszero(StarAlgebras.aug(elt)) && augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit))
iszero(StarAlgebras.aug(order_unit)),
) )
@assert parent(elt) === parent(order_unit) @assert parent(elt) === parent(order_unit)
@ -71,7 +111,7 @@ function sos_problem_primal(
@warn "Setting `upper_bound` together with zero `order_unit` has no effect" @warn "Setting `upper_bound` together with zero `order_unit` has no effect"
end end
A = constraints(parent(elt); augmented = augmented) A = constraints(parent(elt), augmented=augmented, twisted=true)
if !iszero(order_unit) if !iszero(order_unit)
λ = JuMP.@variable(model, λ) λ = JuMP.@variable(model, λ)
@ -95,11 +135,11 @@ end
function invariant_constraint!( function invariant_constraint!(
result::AbstractMatrix, result::AbstractMatrix,
basis::StarAlgebras.AbstractBasis, basis::StarAlgebras.AbstractBasis,
cnstrs::AbstractDict{K,<:ConstraintMatrix}, cnstrs::AbstractDict{K,CM},
invariant_vec::SparseVector, invariant_vec::SparseVector,
) where {K} ) where {K,CM<:ConstraintMatrix}
result .= zero(eltype(result)) result .= zero(eltype(result))
@inbounds for i in SparseArrays.nonzeroinds(invariant_vec) for i in SparseArrays.nonzeroinds(invariant_vec)
g = basis[i] g = basis[i]
A = cnstrs[g] A = cnstrs[g]
for (idx, v) in nzpairs(A) for (idx, v) in nzpairs(A)
@ -109,47 +149,25 @@ function invariant_constraint!(
return result return result
end end
function invariant_constraint(basis, cnstrs, invariant_vec)
I = UInt32[]
J = UInt32[]
V = Float64[]
_M = first(values(cnstrs))
CI = CartesianIndices(_M)
@inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
g = basis[i]
A = cnstrs[g]
for (idx, v) in nzpairs(A)
ci = CI[idx]
push!(I, ci[1])
push!(J, ci[2])
push!(V, invariant_vec[i] * v)
end
end
return sparse(I, J, V, size(_M)...)
end
function isorth_projection(ds::SymbolicWedderburn.DirectSummand) function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
U = SymbolicWedderburn.image_basis(ds) U = SymbolicWedderburn.image_basis(ds)
return isapprox(U * U', I) return isapprox(U * U', I)
end end
function sos_problem_primal( sos_problem_primal(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
wedderburn::WedderburnDecomposition; wedderburn::WedderburnDecomposition;
kwargs..., kwargs...
) ) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
return sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
end
function __fast_recursive_dot!( function __fast_recursive_dot!(
res::JuMP.AffExpr, res::JuMP.AffExpr,
Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}}, Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
Ms::AbstractArray{<:AbstractSparseMatrix}, Ms::AbstractArray{<:AbstractSparseMatrix};
weights,
) )
@assert length(Ps) == length(Ms) @assert length(Ps) == length(Ms)
for (w, A, P) in zip(weights, Ms, Ps) for (A, P) in zip(Ms, Ps)
iszero(Ms) && continue iszero(Ms) && continue
rows = rowvals(A) rows = rowvals(A)
vals = nonzeros(A) vals = nonzeros(A)
@ -157,21 +175,13 @@ function __fast_recursive_dot!(
for i in nzrange(A, cidx) for i in nzrange(A, cidx)
ridx = rows[i] ridx = rows[i]
val = vals[i] val = vals[i]
JuMP.add_to_expression!(res, P[ridx, cidx], w * val) JuMP.add_to_expression!(res, P[ridx, cidx], val)
end end
end end
end end
return res return res
end end
function _dot(
Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
Ms::AbstractArray{<:AbstractMatrix{T}},
weights = Iterators.repeated(one(T), length(Ms)),
) where {T}
return __fast_recursive_dot!(JuMP.AffExpr(), Ps, Ms, weights)
end
import ProgressMeter import ProgressMeter
__show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")] __show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")]
@ -179,38 +189,21 @@ function sos_problem_primal(
elt::StarAlgebras.AlgebraElement, elt::StarAlgebras.AlgebraElement,
orderunit::StarAlgebras.AlgebraElement, orderunit::StarAlgebras.AlgebraElement,
wedderburn::WedderburnDecomposition; wedderburn::WedderburnDecomposition;
upper_bound = Inf, upper_bound=Inf,
augmented = iszero(StarAlgebras.aug(elt)) && augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
iszero(StarAlgebras.aug(orderunit)), check_orthogonality=true,
check_orthogonality = true, show_progress=false
show_progress = false,
) )
@assert parent(elt) === parent(orderunit) @assert parent(elt) === parent(orderunit)
if check_orthogonality if check_orthogonality
if any(!isorth_projection, direct_summands(wedderburn)) if any(!isorth_projection, direct_summands(wedderburn))
error( error("Wedderburn decomposition contains a non-orthogonal projection")
"Wedderburn decomposition contains a non-orthogonal projection",
)
end end
end end
id_one = findfirst(invariant_vectors(wedderburn)) do v
b = basis(parent(elt))
return sparsevec([b[one(first(b))]], [1 // 1], length(v)) == v
end
prog = ProgressMeter.Progress(
length(invariant_vectors(wedderburn));
dt = 1,
desc = "Adding constraints: ",
enabled = show_progress,
barlen = 60,
showspeed = true,
)
feasibility_problem = iszero(orderunit) feasibility_problem = iszero(orderunit)
# problem creation starts here
model = JuMP.Model() model = JuMP.Model()
if !feasibility_problem # add λ or not? if !feasibility_problem # add λ or not?
λ = JuMP.@variable(model, λ) λ = JuMP.@variable(model, λ)
@ -224,52 +217,58 @@ function sos_problem_primal(
end end
end end
# semidefinite constraints as described by wedderburn P = map(direct_summands(wedderburn)) do ds
Ps = map(direct_summands(wedderburn)) do ds
dim = size(ds, 1) dim = size(ds, 1)
P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric) P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
JuMP.@constraint(model, P in PSDCone()) JuMP.@constraint(model, P in PSDCone())
return P P
end end
begin # Ms are preallocated for the constraints loop begin # preallocating
T = eltype(wedderburn) T = eltype(wedderburn)
Ms = [spzeros.(T, size(p)...) for p in Ps] Ms = [spzeros.(T, size(p)...) for p in P]
_eps = 10 * eps(T) * max(size(parent(elt).mstructure)...) M_orb = zeros(T, size(parent(elt).mstructure)...)
end end
X = StarAlgebras.coeffs(elt) X = convert(Vector{T}, StarAlgebras.coeffs(elt))
U = StarAlgebras.coeffs(orderunit) U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
# defining constraints based on the multiplicative structure # defining constraints based on the multiplicative structure
cnstrs = constraints(parent(elt); augmented = augmented) cnstrs = constraints(parent(elt), augmented=augmented, twisted=true)
prog = ProgressMeter.Progress(
length(invariant_vectors(wedderburn)),
dt=1,
desc="Adding constraints... ",
enabled=show_progress,
barlen=60,
showspeed=true
)
# adding linear constraints: one per orbit
for (i, iv) in enumerate(invariant_vectors(wedderburn)) for (i, iv) in enumerate(invariant_vectors(wedderburn))
ProgressMeter.next!(prog; showvalues = __show_itrs(i, prog.n)) ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n))
augmented && i == id_one && continue
# i == 500 && break
x = dot(X, iv) x = dot(X, iv)
u = dot(U, iv) u = dot(U, iv)
spM_orb = invariant_constraint(basis(parent(elt)), cnstrs, iv) M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
Ms = SymbolicWedderburn.diagonalize!(Ms, M_orb, wedderburn)
SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
# @info [nnz(m) / length(m) for m in Ms]
Ms = SymbolicWedderburn.diagonalize!(
Ms,
spM_orb,
wedderburn;
trace_preserving = true,
)
for M in Ms
SparseArrays.droptol!(M, _eps)
end
if feasibility_problem if feasibility_problem
JuMP.@constraint(model, x == _dot(Ps, Ms)) JuMP.@constraint(
model,
x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
)
else else
JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms)) JuMP.@constraint(
model,
x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
)
end end
end end
ProgressMeter.finish!(prog) ProgressMeter.finish!(prog)
return model, Ps return model, P
end end

155
test/1703.09680.jl
View File

@ -1,8 +1,9 @@
@testset "1703.09680 Examples" begin @testset "1703.09680 Examples" begin
@testset "SL(2,Z)" begin @testset "SL(2,Z)" begin
N = 2 N = 2
G = MatrixGroups.SpecialLinearGroup{N}(Int8) G = MatrixGroups.SpecialLinearGroup{N}(Int8)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2) RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -14,15 +15,15 @@
status, certified, λ = check_positivity( status, certified, λ = check_positivity(
elt, elt,
unit; unit,
upper_bound = ub, upper_bound=ub,
halfradius = 2, halfradius=2,
optimizer = scs_optimizer(; optimizer=scs_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 5_000, max_iters=5_000,
accel = 50, accel=50,
alpha = 1.9, alpha=1.9,
), )
) )
@test status == JuMP.ALMOST_OPTIMAL @test status == JuMP.ALMOST_OPTIMAL
@ -32,10 +33,8 @@
@testset "SL(3,F₅)" begin @testset "SL(3,F₅)" begin
N = 3 N = 3
G = MatrixGroups.SpecialLinearGroup{N}( G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{5})
SymbolicWedderburn.Characters.FiniteFields.GF{5}, RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -47,15 +46,15 @@
status, certified, λ = check_positivity( status, certified, λ = check_positivity(
elt, elt,
unit; unit,
upper_bound = ub, upper_bound=ub,
halfradius = 2, halfradius=2,
optimizer = scs_optimizer(; optimizer=scs_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 5_000, max_iters=5_000,
accel = 50, accel=50,
alpha = 1.9, alpha=1.9,
), )
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@ -63,15 +62,12 @@
@test λ > 1 @test λ > 1
m = PropertyT.sos_problem_dual(elt, unit) m = PropertyT.sos_problem_dual(elt, unit)
PropertyT.solve( PropertyT.solve(m, cosmo_optimizer(
m, eps=1e-6,
cosmo_optimizer(; max_iters=5_000,
eps = 1e-6, accel=50,
max_iters = 5_000, alpha=1.9,
accel = 50, ))
alpha = 1.9,
),
)
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL) @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
@test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2 @test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2
@ -80,7 +76,7 @@
@testset "SAut(F₂)" begin @testset "SAut(F₂)" begin
N = 2 N = 2
G = SpecialAutomorphismGroup(FreeGroup(N)) G = SpecialAutomorphismGroup(FreeGroup(N))
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2) RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -92,46 +88,53 @@
status, certified, λ = check_positivity( status, certified, λ = check_positivity(
elt, elt,
unit; unit,
upper_bound = ub, upper_bound=ub,
halfradius = 2, halfradius=2,
optimizer = scs_optimizer(; optimizer=scs_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 5_000, max_iters=5_000,
accel = 50, accel=50,
alpha = 1.9, alpha=1.9,
), )
) )
@test status == JuMP.ALMOST_OPTIMAL @test status == JuMP.ALMOST_OPTIMAL
@test λ < 0 @test λ < 0
@test !certified @test !certified
@time sos_problem = PropertyT.sos_problem_primal(elt; upper_bound = ub) @time sos_problem =
PropertyT.sos_problem_primal(elt, upper_bound=ub)
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
sos_problem, sos_problem,
cosmo_optimizer(; cosmo_optimizer(
eps = 1e-7, eps=1e-7,
max_iters = 10_000, max_iters=10_000,
accel = 0, accel=0,
alpha = 1.9, alpha=1.9,
), )
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
P = JuMP.value.(sos_problem[:P]) P = JuMP.value.(sos_problem[:P])
Q = real.(sqrt(P)) Q = real.(sqrt(P))
certified, λ_cert = certified, λ_cert = PropertyT.certify_solution(
PropertyT.certify_solution(elt, zero(elt), 0.0, Q; halfradius = 2) elt,
zero(elt),
0.0,
Q,
halfradius=2,
)
@test !certified @test !certified
@test λ_cert < 0 @test λ_cert < 0
end end
@testset "SL(3,Z) has (T)" begin @testset "SL(3,Z) has (T)" begin
n = 3 n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2) RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
Δ = RSL(length(S)) - sum(RSL(s) for s in S) Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@ -142,18 +145,18 @@
opt_problem = PropertyT.sos_problem_primal( opt_problem = PropertyT.sos_problem_primal(
elt, elt,
unit; unit,
upper_bound = ub, upper_bound=ub,
augmented = false, augmented=false,
) )
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
opt_problem, opt_problem,
cosmo_optimizer(; cosmo_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 10_000, max_iters=10_000,
accel = 0, accel=0,
alpha = 1.5, alpha=1.5,
), ),
) )
@ -167,13 +170,13 @@
elt, elt,
unit, unit,
λ, λ,
Q; Q,
halfradius = 2, halfradius=2,
augmented = false, augmented=false,
) )
@test certified @test certified
@test isapprox(λ_cert, λ, rtol = 1e-5) @test isapprox(λ_cert, λ, rtol=1e-5)
end end
@testset "augmented formulation" begin @testset "augmented formulation" begin
@ -183,18 +186,18 @@
opt_problem = PropertyT.sos_problem_primal( opt_problem = PropertyT.sos_problem_primal(
elt, elt,
unit; unit,
upper_bound = ub, upper_bound=ub,
augmented = true, augmented=true,
) )
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
opt_problem, opt_problem,
scs_optimizer(; scs_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 10_000, max_iters=10_000,
accel = -10, accel=-10,
alpha = 1.5, alpha=1.5,
), ),
) )
@ -207,13 +210,13 @@
elt, elt,
unit, unit,
λ, λ,
Q; Q,
halfradius = 2, halfradius=2,
augmented = true, augmented=true,
) )
@test certified @test certified
@test isapprox(λ_cert, λ, rtol = 1e-5) @test isapprox(λ_cert, λ, rtol=1e-5)
@test λ_cert > 2 // 10 @test λ_cert > 2 // 10
end end
end end

88
test/1712.07167.jl
View File

@ -1,9 +1,9 @@
@testset "1712.07167 Examples" begin @testset "1712.07167 Examples" begin
@testset "SAut(F₃)" begin @testset "SAut(F₃)" begin
N = 3 N = 3
G = SpecialAutomorphismGroup(FreeGroup(N)) G = SpecialAutomorphismGroup(FreeGroup(N))
@info "running tests for" G RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
@ -15,7 +15,6 @@
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -28,14 +27,14 @@
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
unit, unit,
wd; wd,
upper_bound = ub, upper_bound=ub,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; optimizer=cosmo_optimizer(
eps = 1e-7, eps=1e-7,
max_iters = 10_000, max_iters=10_000,
accel = 50, accel=50,
alpha = 1.9, alpha=1.9,
), ),
) )
@ -48,8 +47,7 @@
n = 3 n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
@info "running tests for" SL RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
Δ = RSL(length(S)) - sum(RSL(s) for s in S) Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@ -64,7 +62,6 @@
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
) )
@info wd
elt = Δ^2 elt = Δ^2
unit = Δ unit = Δ
@ -74,8 +71,8 @@
elt, elt,
unit, unit,
wd, wd,
upper_bound = ub, upper_bound=ub,
augmented = false, augmented=false,
) )
wdfl = SymbolicWedderburn.WedderburnDecomposition( wdfl = SymbolicWedderburn.WedderburnDecomposition(
@ -89,18 +86,18 @@
model, varP = PropertyT.sos_problem_primal( model, varP = PropertyT.sos_problem_primal(
elt, elt,
unit, unit,
wdfl; wdfl,
upper_bound = ub, upper_bound=ub,
augmented = false, augmented=false,
) )
status, warm = PropertyT.solve( status, warm = PropertyT.solve(
model, model,
cosmo_optimizer(; cosmo_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 20_000, max_iters=20_000,
accel = 50, accel=50,
alpha = 1.9, alpha=1.9,
), ),
) )
@ -108,34 +105,34 @@
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
model, model,
scs_optimizer(; scs_optimizer(
eps = 1e-10, eps=1e-10,
max_iters = 100, max_iters=100,
accel = -20, accel=-20,
alpha = 1.2, alpha=1.2,
), ),
warm, warm
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
Q = @time let varP = varP Q = @time let varP = varP
Qs = map(varP) do P Qs = map(varP) do P
return real.(sqrt(JuMP.value.(P))) real.(sqrt(JuMP.value.(P)))
end end
PropertyT.reconstruct(Qs, wdfl) PropertyT.reconstruct(Qs, wdfl)
end end
λ = JuMP.value(model[:λ]) λ = JuMP.value(model[:λ])
sos = PropertyT.compute_sos(parent(elt), Q; augmented = false) sos = PropertyT.compute_sos(parent(elt), Q; augmented=false)
certified, λ_cert = PropertyT.certify_solution( certified, λ_cert = PropertyT.certify_solution(
elt, elt,
unit, unit,
λ, λ,
Q; Q,
halfradius = 2, halfradius=2,
augmented = false, augmented=false,
) )
@test certified @test certified
@ -158,23 +155,22 @@
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
) )
@info wdfl
opt_problem, varP = PropertyT.sos_problem_primal( opt_problem, varP = PropertyT.sos_problem_primal(
elt, elt,
unit, unit,
wdfl; wdfl,
upper_bound = ub, upper_bound=ub,
# augmented = true # since both elt and unit are augmented # augmented = true # since both elt and unit are augmented
) )
status, _ = PropertyT.solve( status, _ = PropertyT.solve(
opt_problem, opt_problem,
scs_optimizer(; scs_optimizer(
eps = 1e-8, eps=1e-8,
max_iters = 20_000, max_iters=20_000,
accel = 0, accel=0,
alpha = 1.9, alpha=1.9,
), ),
) )
@ -182,7 +178,7 @@
Q = @time let varP = varP Q = @time let varP = varP
Qs = map(varP) do P Qs = map(varP) do P
return real.(sqrt(JuMP.value.(P))) real.(sqrt(JuMP.value.(P)))
end end
PropertyT.reconstruct(Qs, wdfl) PropertyT.reconstruct(Qs, wdfl)
end end
@ -191,8 +187,8 @@
elt, elt,
unit, unit,
JuMP.objective_value(opt_problem), JuMP.objective_value(opt_problem),
Q; Q,
halfradius = 2, halfradius=2,
# augmented = true # since both elt and unit are augmented # augmented = true # since both elt and unit are augmented
) )

67
test/1812.03456.jl
View File

@ -21,9 +21,8 @@ using SparseArrays
@testset "Sq, Adj, Op in SL(4,Z)" begin @testset "Sq, Adj, Op in SL(4,Z)" begin
N = 4 N = 4
G = MatrixGroups.SpecialLinearGroup{N}(Int8) G = MatrixGroups.SpecialLinearGroup{N}(Int8)
@info "running tests for" G
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2) RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -40,7 +39,6 @@ using SparseArrays
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
ivs = SymbolicWedderburn.invariant_vectors(wd) ivs = SymbolicWedderburn.invariant_vectors(wd)
sq, adj, op = PropertyT.SqAdjOp(RG, N) sq, adj, op = PropertyT.SqAdjOp(RG, N)
@ -84,7 +82,7 @@ using SparseArrays
@testset "SAut(F₃)" begin @testset "SAut(F₃)" begin
n = 3 n = 3
G = SpecialAutomorphismGroup(FreeGroup(n)) G = SpecialAutomorphismGroup(FreeGroup(n))
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2) RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
sq, adj, op = PropertyT.SqAdjOp(RG, n) sq, adj, op = PropertyT.SqAdjOp(RG, n)
@test sq(one(G)) == 216 @test sq(one(G)) == 216
@ -100,8 +98,7 @@ end
n = 3 n = 3
G = MatrixGroups.SpecialLinearGroup{n}(Int8) G = MatrixGroups.SpecialLinearGroup{n}(Int8)
@info "running tests for" G RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -116,7 +113,6 @@ end
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
sq, adj, op = PropertyT.SqAdjOp(RG, n) sq, adj, op = PropertyT.SqAdjOp(RG, n)
@ -127,10 +123,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd; wd,
upper_bound = UB, upper_bound=UB,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9), optimizer=cosmo_optimizer(accel=50, alpha=1.9)
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -144,10 +140,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd; wd,
upper_bound = UB, upper_bound=UB,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9), optimizer=cosmo_optimizer(accel=50, alpha=1.9)
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -156,11 +152,10 @@ end
m, _ = PropertyT.sos_problem_primal(elt, wd) m, _ = PropertyT.sos_problem_primal(elt, wd)
PropertyT.solve( PropertyT.solve(
m, m,
scs_optimizer(; max_iters = 5000, accel = 50, alpha = 1.9), scs_optimizer(max_iters=5000, accel=50, alpha=1.9)
) )
@test JuMP.termination_status(m) in @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
(JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
@test abs(JuMP.objective_value(m)) < 1e-3 @test abs(JuMP.objective_value(m)) < 1e-3
end end
@ -173,10 +168,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd; wd,
upper_bound = UB, upper_bound=UB,
halfradius = 2, halfradius=2,
optimizer = scs_optimizer(; accel = 50, alpha = 1.9), optimizer=cosmo_optimizer(accel=50, alpha=1.9)
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test !certified @test !certified
@ -188,8 +183,7 @@ end
n = 4 n = 4
G = MatrixGroups.SpecialLinearGroup{n}(Int8) G = MatrixGroups.SpecialLinearGroup{n}(Int8)
@info "running tests for" G RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -204,7 +198,6 @@ end
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
) )
@info wd
sq, adj, op = PropertyT.SqAdjOp(RG, n) sq, adj, op = PropertyT.SqAdjOp(RG, n)
@ -215,10 +208,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd; wd,
upper_bound = UB, upper_bound=UB,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9), optimizer=cosmo_optimizer(accel=50, alpha=1.9)
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -232,10 +225,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd; wd,
upper_bound = UB, upper_bound=UB,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9), optimizer=cosmo_optimizer(accel=50, alpha=1.9)
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test certified @test certified
@ -249,10 +242,10 @@ end
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
Δ, Δ,
wd; wd,
upper_bound = UB, upper_bound=UB,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9), optimizer=cosmo_optimizer(accel=50, alpha=1.9)
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@test !certified @test !certified

168
test/Chevalley.jl
View File

@ -1,168 +0,0 @@
countmap(v) = countmap(identity, v)
function countmap(f, v)
counts = Dict{eltype(f(first(v))),Int}()
for x in v
fx = f(x)
counts[fx] = get!(counts, fx, 0) + 1
end
return counts
end
@testset "classify_root_system" begin
α = PropertyT.Roots.Root([1, -1, 0])
β = PropertyT.Roots.Root([0, 1, -1])
γ = PropertyT.Roots.Root([2, 0, 0])
@test PropertyT.Roots.classify_root_system(α, β, (false, false)) == :A₂
@test PropertyT.Roots.classify_root_system(α, γ, (false, true)) == :C₂
@test PropertyT.Roots.classify_root_system(β, γ, (false, true)) ==
Symbol("A₁×C₁")
end
@testset "Exceptional root systems" begin
@testset "F4" begin
F4 = let Σ = PermutationGroups.PermGroup(perm"(1,2,3,4)", perm"(1,2)")
long = let x = (1, 1, 0, 0) .// 1
PropertyT.Roots.Root.(
union(
(x^g for g in Σ),
((x .* (-1, 1, 1, 1))^g for g in Σ),
((-1 .* x)^g for g in Σ),
),
)
end
short = let x = (1, 0, 0, 0) .// 1
PropertyT.Roots.Root.(
union((x^g for g in Σ), ((-1 .* x)^g for g in Σ))
)
end
signs = collect(Iterators.product(fill([-1, +1], 4)...))
halfs = let x = (1, 1, 1, 1) .// 2
PropertyT.Roots.Root.(union(x .* sgn for sgn in signs))
end
union(long, short, halfs)
end
@test length(F4) == 48
a = F4[1]
@test isapprox(PropertyT.Roots.ℓ₂length(a), sqrt(2))
b = F4[6]
@test isapprox(PropertyT.Roots.ℓ₂length(b), sqrt(2))
c = a + b
@test isapprox(PropertyT.Roots.ℓ₂length(c), 2.0)
@test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
long = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == sqrt(2), F4)]
short = F4[findfirst(r -> PropertyT.Roots.ℓ₂length(r) == 1.0, F4)]
subtypes = Set([:C₂, :A₂, Symbol("A₁×C₁")])
let Ω = F4, α = long
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:C₂], 6)
@test r == 0 && d == 3
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 4
end
let Ω = F4, α = short
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:C₂], 6)
@test r == 0 && d == 3
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 4
end
end
@testset "E6-7-8 exceptional root systems" begin
E8 =
let Σ = PermutationGroups.PermGroup(
perm"(1,2,3,4,5,6,7,8)",
perm"(1,2)",
)
long = let x = (1, 1, 0, 0, 0, 0, 0, 0) .// 1
PropertyT.Roots.Root.(
union(
(x^g for g in Σ),
((x .* (-1, 1, 1, 1, 1, 1, 1, 1))^g for g in Σ),
((-1 .* x)^g for g in Σ),
),
)
end
signs = collect(
p for p in Iterators.product(fill([-1, +1], 8)...) if
iseven(count(==(-1), p))
)
halfs = let x = (1, 1, 1, 1, 1, 1, 1, 1) .// 2
rts = unique(PropertyT.Roots.Root(x .* sgn) for sgn in signs)
end
union(long, halfs)
end
subtypes = Set([:A₂, Symbol("A₁×A₁")])
@testset "E8" begin
@test length(E8) == 240
@test all(r -> PropertyT.Roots.ℓ₂length(r) ≈ sqrt(2), E8)
let Ω = E8, α = first(Ω)
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 28
end
end
@testset "E7" begin
E7 = filter(r -> iszero(sum(r.coord)), E8)
@test length(E7) == 126
let Ω = E7, α = first(Ω)
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:A₂], 4)
@test r == 0 && d == 16
end
end
@testset "E6" begin
E6 = filter(
r -> r.coord[end] == r.coord[end-1] == r.coord[end-2],
E8,
)
@test length(E6) == 72
let Ω = E6, α = first(Ω)
counts = countmap([
PropertyT.Roots.classify_sub_root_system(Ω, α, γ) for
γ in Ω if !PropertyT.Roots.isproportional(α, γ)
])
@test Set(keys(counts)) == subtypes
d, r = divrem(counts[:A₂], 4)
@info d, r
@test r == 0 && d == 10
end
end
end
end

56
test/actions.jl
View File

@ -3,7 +3,7 @@ function test_action(basis, group, act)
return @testset "action definition" begin return @testset "action definition" begin
@test all(basis) do b @test all(basis) do b
e = one(group) e = one(group)
return action(act, e, b) == b action(act, e, b) == b
end end
a = let a = rand(basis) a = let a = rand(basis)
@ -30,12 +30,12 @@ function test_action(basis, group, act)
@test all([(g, h) for g in group for h in group]) do (g, h) @test all([(g, h) for g in group for h in group]) do (g, h)
x = action(act, h, action(act, g, a)) x = action(act, h, action(act, g, a))
y = action(act, g * h, a) y = action(act, g * h, a)
return x == y x == y
end end
if act isa SymbolicWedderburn.ByPermutations if act isa SymbolicWedderburn.ByPermutations
@test all(basis) do b @test all(basis) do b
return action(act, g, b) ∈ basis && action(act, h, b) ∈ basis action(act, g, b) ∈ basis && action(act, h, b) ∈ basis
end end
end end
end end
@ -44,18 +44,21 @@ end
## Testing ## Testing
@testset "Actions on SL(3,ℤ)" begin @testset "Actions on SL(3,ℤ)" begin
n = 3 n = 3
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2) RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
@testset "Permutation action" begin @testset "Permutation action" begin
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
ΓpA = PropertyT.action_by_conjugation(SL, Γ) ΓpA = PropertyT.action_by_conjugation(SL, Γ)
test_action(basis(RSL), Γ, ΓpA) test_action(basis(RSL), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -63,9 +66,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, ranks = @time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(isone, ranks) @test all(simple)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -74,17 +77,17 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]])
) )
@test length(invariant_vectors(wd)) == 918 @test length(invariant_vectors(wd)) == 918
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [40, 23, 18]
[40, 23, 18]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end
@testset "Wreath action" begin @testset "Wreath action" begin
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
end end
@ -94,6 +97,7 @@ end
test_action(basis(RSL), Γ, ΓpA) test_action(basis(RSL), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -101,9 +105,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, ranks = @time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(isone, ranks) @test all(simple)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -112,30 +116,32 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSL), basis(RSL),
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]), StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]])
) )
@test length(invariant_vectors(wd)) == 247 @test length(invariant_vectors(wd)) == 247
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [14, 9, 6, 14, 12]
[14, 9, 6, 14, 12]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end
end end
@testset "Actions on SAut(F4)" begin @testset "Actions on SAut(F4)" begin
n = 4 n = 4
SAutFn = SpecialAutomorphismGroup(FreeGroup(n)) SAutFn = SpecialAutomorphismGroup(FreeGroup(n))
RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn; halfradius = 1) RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn, halfradius=1, twisted=true)
@testset "Permutation action" begin @testset "Permutation action" begin
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ) ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ)
test_action(basis(RSAutFn), Γ, ΓpA) test_action(basis(RSAutFn), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -143,9 +149,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, ranks = @time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(isone, ranks) @test all(simple)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -154,17 +160,17 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSAutFn), basis(RSAutFn),
StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]), StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]])
) )
@test length(invariant_vectors(wd)) == 93 @test length(invariant_vectors(wd)) == 93
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [4, 8, 5, 4]
[4, 8, 5, 4]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end
@testset "Wreath action" begin @testset "Wreath action" begin
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1))) Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
end end
@ -174,6 +180,7 @@ end
test_action(basis(RSAutFn), Γ, ΓpA) test_action(basis(RSAutFn), Γ, ΓpA)
@testset "mps is successful" begin @testset "mps is successful" begin
charsΓ = charsΓ =
SymbolicWedderburn.Character{ SymbolicWedderburn.Character{
Rational{Int}, Rational{Int},
@ -181,9 +188,9 @@ end
RΓ = SymbolicWedderburn._group_algebra(Γ) RΓ = SymbolicWedderburn._group_algebra(Γ)
@time mps, ranks = @time mps, simple =
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ) SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
@test all(isone, ranks) @test all(simple)
end end
@testset "Wedderburn decomposition" begin @testset "Wedderburn decomposition" begin
@ -192,12 +199,11 @@ end
Γ, Γ,
ΓpA, ΓpA,
basis(RSAutFn), basis(RSAutFn),
StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]]), StarAlgebras.Basis{UInt16}(@view basis(RSAutFn)[1:sizes[1]])
) )
@test length(invariant_vectors(wd)) == 18 @test length(invariant_vectors(wd)) == 18
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == @test SymbolicWedderburn.size.(direct_summands(wd), 1) == [1, 1, 2, 2, 1, 2, 2, 1]
[1, 1, 2, 2, 1, 2, 2, 1]
@test all(issimple, direct_summands(wd)) @test all(issimple, direct_summands(wd))
end end
end end

35
test/check_positivity.jl
View File

@ -1,12 +1,6 @@
function check_positivity( function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
elt,
unit;
upper_bound = Inf,
halfradius = 2,
optimizer,
)
@time sos_problem = @time sos_problem =
PropertyT.sos_problem_primal(elt, unit; upper_bound = upper_bound) PropertyT.sos_problem_primal(elt, unit, upper_bound=upper_bound)
status, _ = PropertyT.solve(sos_problem, optimizer) status, _ = PropertyT.solve(sos_problem, optimizer)
P = JuMP.value.(sos_problem[:P]) P = JuMP.value.(sos_problem[:P])
@ -15,23 +9,16 @@ function check_positivity(
elt, elt,
unit, unit,
JuMP.objective_value(sos_problem), JuMP.objective_value(sos_problem),
Q; Q,
halfradius = halfradius, halfradius=halfradius,
) )
return status, certified, λ_cert return status, certified, λ_cert
end end
function check_positivity( function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer)
elt,
unit,
wd;
upper_bound = Inf,
halfradius = 2,
optimizer,
)
@assert aug(elt) == aug(unit) == 0 @assert aug(elt) == aug(unit) == 0
@time sos_problem, Ps = @time sos_problem, Ps =
PropertyT.sos_problem_primal(elt, unit, wd; upper_bound = upper_bound) PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound)
@time status, _ = PropertyT.solve(sos_problem, optimizer) @time status, _ = PropertyT.solve(sos_problem, optimizer)
@ -42,7 +29,13 @@ function check_positivity(
λ = JuMP.value(sos_problem[:λ]) λ = JuMP.value(sos_problem[:λ])
certified, λ_cert = certified, λ_cert = PropertyT.certify_solution(
PropertyT.certify_solution(elt, unit, λ, Q; halfradius = halfradius) elt,
unit,
λ,
Q,
halfradius=halfradius
)
return status, certified, λ_cert return status, certified, λ_cert
end end

14
test/constratint_matrices.jl
View File

@ -1,17 +1,12 @@
@testset "ConstraintMatrix" begin @testset "ConstraintMatrix" begin
@test PropertyT.ConstraintMatrix{Float64}( @test PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) isa AbstractMatrix
[-1, 2, -1, 1, 4, 2, 6],
3,
2,
π,
) isa AbstractMatrix
cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π)
@test cm == Float64[ @test cm == Float64[
-π π -π π
2π 0 2π 0.0
0 π 0.0 π
] ]
@test collect(PropertyT.nzpairs(cm)) == [ @test collect(PropertyT.nzpairs(cm)) == [
@ -23,7 +18,4 @@
1 => -3.141592653589793 1 => -3.141592653589793
1 => -3.141592653589793 1 => -3.141592653589793
] ]
@test PropertyT.ConstraintMatrix{Float64}([-9:-1; 1:9], 3, 3, 1.0) ==
zeros(3, 3)
end end

26
test/graded_adj.jl
View File

@ -1,9 +1,10 @@
@testset "Adj via grading" begin @testset "Adj via grading" begin
@testset "SL(n,Z) & Aut(F₄)" begin @testset "SL(n,Z) & Aut(F₄)" begin
n = 4 n = 4
halfradius = 1 halfradius = 1
SL = MatrixGroups.SpecialLinearGroup{n}(Int8) SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = halfradius) RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=halfradius, twisted=true)
Δ = RSL(length(S)) - sum(RSL(s) for s in S) Δ = RSL(length(S)) - sum(RSL(s) for s in S)
@ -21,9 +22,10 @@
@test PropertyT.Adj(Δs, :A₂) == adj @test PropertyT.Adj(Δs, :A₂) == adj
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
halfradius = 1 halfradius = 1
G = SpecialAutomorphismGroup(FreeGroup(n)) G = SpecialAutomorphismGroup(FreeGroup(n))
RG, S, sizes = PropertyT.group_algebra(G; halfradius = halfradius) RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true)
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -42,6 +44,7 @@
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
end end
@testset "Symplectic group" begin @testset "Symplectic group" begin
@testset "Sp2(ℤ)" begin @testset "Sp2(ℤ)" begin
genus = 2 genus = 2
@ -49,8 +52,7 @@
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8) SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
RSpN, S_sp, sizes_sp = RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
PropertyT.group_algebra(SpN; halfradius = halfradius)
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -71,8 +73,7 @@
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8) SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
RSpN, S_sp, sizes_sp = RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
PropertyT.group_algebra(SpN; halfradius = halfradius)
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
Δ = RG(length(S)) - sum(RG(s) for s in S) Δ = RG(length(S)) - sum(RG(s) for s in S)
@ -85,14 +86,9 @@
end end
@testset "Adj numerics for genus=$genus" begin @testset "Adj numerics for genus=$genus" begin
all_subtypes = ( all_subtypes = (
:A₁, :A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂
:C₁,
Symbol("A₁×A₁"),
Symbol("C₁×C₁"),
Symbol("A₁×C₁"),
:A₂,
:C₂,
) )
@test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384 @test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384
@ -105,8 +101,8 @@
@test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16) @test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16)
end end
end end
@test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == @test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2
Δ^2
end end
end end
end end

36
test/optimizers.jl
View File

@ -4,12 +4,12 @@ import JuMP
import SCS import SCS
function scs_optimizer(; function scs_optimizer(;
accel = 10, accel=10,
alpha = 1.5, alpha=1.5,
eps = 1e-9, eps=1e-9,
max_iters = 100_000, max_iters=100_000,
verbose = true, verbose=true,
linear_solver = SCS.DirectSolver, linear_solver=SCS.DirectSolver
) )
return JuMP.optimizer_with_attributes( return JuMP.optimizer_with_attributes(
SCS.Optimizer, SCS.Optimizer,
@ -28,27 +28,27 @@ end
import COSMO import COSMO
function cosmo_optimizer(; function cosmo_optimizer(;
accel = 15, accel=15,
alpha = 1.6, alpha=1.6,
eps = 1e-9, eps=1e-9,
max_iters = 100_000, max_iters=100_000,
verbose = true, verbose=true,
verbose_timing = verbose, verbose_timing=verbose,
decompose = false, decompose=false
) )
return JuMP.optimizer_with_attributes( return JuMP.optimizer_with_attributes(
COSMO.Optimizer, COSMO.Optimizer,
"accelerator" => "accelerator" => COSMO.with_options(
COSMO.with_options(COSMO.AndersonAccelerator; mem = max(accel, 2)), COSMO.AndersonAccelerator,
mem=max(accel, 2)
),
"alpha" => alpha, "alpha" => alpha,
"decompose" => decompose, "decompose" => decompose,
"eps_abs" => eps, "eps_abs" => eps,
"eps_rel" => eps, "eps_rel" => eps,
"eps_prim_inf" => eps,
"eps_dual_inf" => eps,
"max_iter" => max_iters, "max_iter" => max_iters,
"verbose" => verbose, "verbose" => verbose,
"verbose_timing" => verbose_timing, "verbose_timing" => verbose_timing,
"check_termination" => 250, "check_termination" => 200,
) )
end end

63
test/quick_tests.jl
View File

@ -1,12 +1,11 @@
@testset "Quick tests" begin @testset "Quick tests" begin
@testset "SL(2,F₇)" begin @testset "SL(2,F₇)" begin
N = 2 N = 2
p = 7 p = 7
halfradius = 3 halfradius = 3
G = MatrixGroups.SpecialLinearGroup{N}( G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p})
SymbolicWedderburn.Characters.FiniteFields.GF{p}, RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true)
)
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 3)
Δ = let RG = RG, S = S Δ = let RG = RG, S = S
RG(length(S)) - sum(RG(s) for s in S) RG(length(S)) - sum(RG(s) for s in S)
@ -19,15 +18,15 @@
@testset "standard formulation" begin @testset "standard formulation" begin
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
unit; unit,
upper_bound = ub, upper_bound=ub,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; optimizer=cosmo_optimizer(
eps = 1e-7, eps=1e-7,
max_iters = 5_000, max_iters=5_000,
accel = 50, accel=50,
alpha = 1.95, alpha=1.95,
), )
) )
@test status == JuMP.OPTIMAL @test status == JuMP.OPTIMAL
@ -35,18 +34,14 @@
@test λ_cert > 5857 // 10000 @test λ_cert > 5857 // 10000
m = PropertyT.sos_problem_dual(elt, unit) m = PropertyT.sos_problem_dual(elt, unit)
PropertyT.solve( PropertyT.solve(m, cosmo_optimizer(
m, eps=1e-7,
cosmo_optimizer(; max_iters=10_000,
eps = 1e-7, accel=50,
max_iters = 10_000, alpha=1.95,
accel = 50, ))
alpha = 1.95,
),
)
@test JuMP.termination_status(m) in @test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
(JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
@test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2 @test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2
end end
@ -60,22 +55,20 @@
Σ, Σ,
act, act,
basis(RG), basis(RG),
StarAlgebras.Basis{UInt16}( StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[halfradius]]),
@view basis(RG)[1:sizes[halfradius]]
),
) )
status, certified, λ_cert = check_positivity( status, certified, λ_cert = check_positivity(
elt, elt,
unit, unit,
wd; wd,
upper_bound = ub, upper_bound=ub,
halfradius = 2, halfradius=2,
optimizer = cosmo_optimizer(; optimizer=cosmo_optimizer(
eps = 1e-7, eps=1e-7,
max_iters = 10_000, max_iters=10_000,
accel = 50, accel=50,
alpha = 1.9, alpha=1.9,
), ),
) )

10
test/roots.jl
View File

@ -1,10 +0,0 @@
using PropertyT.Roots
@testset "Roots" begin
@test Roots.Root{3,Int}([1, 2, 3]) isa Roots.AbstractRoot{}
@test Roots.Root([1, 2, 3]) isa Roots.AbstractRoot{3,Int}
# io
r = Roots.Root{3,Int}([1, 2, 3])
@test contains(sprint(show, MIME"text/plain"(), r), "of length √14\n")
r = Roots.Root{3,Int}([1, 2, 2])
@test contains(sprint(show, MIME"text/plain"(), r), "of length 3\n")
end

2
test/runtests.jl
View File

@ -24,8 +24,6 @@ if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI")
include("1712.07167.jl") include("1712.07167.jl")
include("1812.03456.jl") include("1812.03456.jl")
include("roots.jl")
include("graded_adj.jl") include("graded_adj.jl")
include("Chevalley.jl")
end end
end end