mirror of
https://github.com/kalmarek/PropertyT.jl.git
synced 2024-11-30 18:05:27 +01:00
Compare commits
No commits in common. "902be093db39c4d9d52cea5f84ad0c61227ce9be" and "0127d05594b83dd03b2908ade48368cf626e8620" have entirely different histories.
902be093db
...
0127d05594
@ -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
|
||||
@ -1,12 +1,11 @@
|
||||
name = "PropertyT"
|
||||
uuid = "03b72c93-0167-51e2-8a1e-eb4ff1fb940d"
|
||||
authors = ["Marek Kaluba <kalmar@amu.edu.pl>"]
|
||||
version = "0.5.0"
|
||||
version = "0.4.0"
|
||||
|
||||
[deps]
|
||||
Groups = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
|
||||
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
|
||||
IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
|
||||
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
|
||||
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
|
||||
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
|
||||
@ -18,12 +17,11 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
|
||||
COSMO = "0.8"
|
||||
Groups = "0.7"
|
||||
IntervalArithmetic = "0.20"
|
||||
IntervalMatrices = "0.8"
|
||||
JuMP = "1.3"
|
||||
ProgressMeter = "1.7"
|
||||
SCS = "1.1"
|
||||
SCS = "1.1.0"
|
||||
StaticArrays = "1"
|
||||
SymbolicWedderburn = "0.3.4"
|
||||
SymbolicWedderburn = "0.3.2"
|
||||
julia = "1.6"
|
||||
|
||||
[extras]
|
||||
|
||||
@ -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),
|
||||
)
|
||||
@ -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
|
||||
@ -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
|
||||
@ -24,7 +24,8 @@ const GENUS = 2N
|
||||
|
||||
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
|
||||
G = StarAlgebras.object(RG)
|
||||
@ -41,8 +42,6 @@ wd = let RG = RG, N = N
|
||||
basis(RG),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
|
||||
)
|
||||
@info wdfl
|
||||
wdfl
|
||||
end
|
||||
|
||||
Δ = RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -59,22 +58,23 @@ unit = Δ
|
||||
@time model, varP = PropertyT.sos_problem_primal(
|
||||
elt,
|
||||
unit,
|
||||
wd;
|
||||
upper_bound = UPPER_BOUND,
|
||||
augmented = true,
|
||||
show_progress = true,
|
||||
wd,
|
||||
upper_bound=UPPER_BOUND,
|
||||
augmented=true,
|
||||
show_progress=true
|
||||
)
|
||||
|
||||
solve_in_loop(
|
||||
model,
|
||||
wd,
|
||||
varP;
|
||||
logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-$(UPPER_BOUND)Δ",
|
||||
optimizer = cosmo_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 20_000,
|
||||
accel = 50,
|
||||
alpha = 1.95,
|
||||
varP,
|
||||
logdir="./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ",
|
||||
optimizer=cosmo_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=20_000,
|
||||
accel=50,
|
||||
alpha=1.95,
|
||||
),
|
||||
data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
|
||||
data=(elt=elt, unit=unit, halfradius=HALFRADIUS)
|
||||
)
|
||||
|
||||
|
||||
@ -11,19 +11,12 @@ function get_solution(model)
|
||||
return solution
|
||||
end
|
||||
|
||||
function get_solution(model, wd, varP, eps = 1e-10)
|
||||
function get_solution(model, wd, varP; logdir)
|
||||
λ = JuMP.value(model[:λ])
|
||||
|
||||
@info "reconstructing the solution"
|
||||
Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP], eps = eps
|
||||
PropertyT.__droptol!.(Ps, 100eps)
|
||||
Qs = real.(sqrt.(Ps))
|
||||
PropertyT.__droptol!.(Qs, eps)
|
||||
PropertyT.reconstruct(Qs, wd)
|
||||
end
|
||||
|
||||
Qs = [real.(sqrt(JuMP.value.(P))) for P in varP]
|
||||
Q = PropertyT.reconstruct(Qs, wd)
|
||||
solution = Dict(:λ => λ, :Q => Q)
|
||||
|
||||
return solution
|
||||
end
|
||||
|
||||
@ -49,23 +42,22 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
|
||||
# logstream = current_logger().logger.stream
|
||||
# v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint
|
||||
# @warn v
|
||||
status, warm =
|
||||
@time PropertyT.solve(log_file, model, optimizer, warm)
|
||||
status, 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.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
|
||||
return PropertyT.certify_solution(
|
||||
PropertyT.certify_solution(
|
||||
data.elt,
|
||||
data.unit,
|
||||
solution[:λ],
|
||||
solution[:Q];
|
||||
halfradius = data.halfradius,
|
||||
solution[:Q],
|
||||
halfradius=data.halfradius,
|
||||
)
|
||||
end
|
||||
end
|
||||
@ -74,20 +66,19 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
|
||||
end
|
||||
|
||||
if flag == true && certified_λ ≥ 0
|
||||
@info "Certification done with λ = $certified_λ" certified_λ status
|
||||
@info "Certification done with λ = $certified_λ"
|
||||
return certified_λ
|
||||
else
|
||||
rel_change =
|
||||
abs(certified_λ - old_lambda) /
|
||||
(abs(certified_λ) + abs(old_lambda))
|
||||
@info "Certification failed with λ = $λ" certified_λ rel_change status
|
||||
if rel_change < 1e-9
|
||||
@info "No progress detected, breaking" certified_λ rel_change status
|
||||
break
|
||||
end
|
||||
rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda))
|
||||
@info "Certification failed with λ = $λ" certified_λ rel_change
|
||||
end
|
||||
|
||||
old_lambda = certified_λ
|
||||
|
||||
if rel_change < 1e-9
|
||||
@info "No progress detected, breaking"
|
||||
break
|
||||
end
|
||||
end
|
||||
|
||||
return status == JuMP.OPTIMAL ? old_lambda : NaN
|
||||
|
||||
@ -3,6 +3,7 @@ module PropertyT
|
||||
using LinearAlgebra
|
||||
using SparseArrays
|
||||
|
||||
using IntervalArithmetic
|
||||
using JuMP
|
||||
|
||||
using Groups
|
||||
@ -25,17 +26,17 @@ include("gradings.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))
|
||||
@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 "computing the *-algebra structure for G"
|
||||
@time RG = StarAlgebras.StarAlgebra(
|
||||
@time RG = StarAlgebras.StarAlgebra{twisted}(
|
||||
G,
|
||||
StarAlgebras.Basis{UInt32}(E),
|
||||
(sizes[halfradius], sizes[halfradius]);
|
||||
precompute = false,
|
||||
(sizes[halfradius], sizes[halfradius]),
|
||||
precompute=false,
|
||||
)
|
||||
return RG, S, sizes
|
||||
end
|
||||
|
||||
@ -1,4 +1,5 @@
|
||||
import SymbolicWedderburn.action
|
||||
StarAlgebras.star(g::Groups.GroupElement) = inv(g)
|
||||
|
||||
include("alphabet_permutation.jl")
|
||||
|
||||
@ -9,7 +10,7 @@ include("autfn_conjugation.jl")
|
||||
function SymbolicWedderburn.action(
|
||||
act::SymbolicWedderburn.ByPermutations,
|
||||
g::Groups.GroupElement,
|
||||
α::StarAlgebras.AlgebraElement,
|
||||
α::StarAlgebras.AlgebraElement
|
||||
)
|
||||
res = StarAlgebras.zero!(similar(α))
|
||||
B = basis(parent(α))
|
||||
|
||||
@ -1,6 +1,3 @@
|
||||
import IntervalArithmetic
|
||||
import IntervalMatrices
|
||||
|
||||
function augment_columns!(Q::AbstractMatrix)
|
||||
for c in eachcol(Q)
|
||||
c .-= sum(c) ./ length(c)
|
||||
@ -11,25 +8,26 @@ end
|
||||
function __sos_via_sqr!(
|
||||
res::StarAlgebras.AlgebraElement,
|
||||
P::AbstractMatrix;
|
||||
augmented::Bool,
|
||||
id = (b = basis(parent(res)); b[one(first(b))]),
|
||||
augmented::Bool
|
||||
)
|
||||
A = parent(res)
|
||||
mstr = A.mstructure
|
||||
@assert size(mstr) == size(P)
|
||||
|
||||
StarAlgebras.zero!(res)
|
||||
for j in axes(mstr, 2)
|
||||
for i in axes(mstr, 1)
|
||||
A = parent(res)
|
||||
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]
|
||||
x_star_y = mstr[-i, j]
|
||||
res[x_star_y] += p
|
||||
# either result += P[x,y]*(x'*y)
|
||||
xy = b[A.mstructure[i, j]]
|
||||
# either result += P[x,y]*(x*y)
|
||||
res[xy] += p
|
||||
if augmented
|
||||
# or result += P[x,y]*(1-x)'*(1-y) == P[x,y]*(1-x'-y+x'y)
|
||||
res[id] += p
|
||||
x_star, y = mstr[-i, id], j
|
||||
res[x_star] -= p
|
||||
# or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy)
|
||||
y = b[j]
|
||||
res[e] += p
|
||||
res[x] -= p
|
||||
res[y] -= p
|
||||
end
|
||||
end
|
||||
@ -38,11 +36,7 @@ function __sos_via_sqr!(
|
||||
return res
|
||||
end
|
||||
|
||||
function __sos_via_cnstr!(
|
||||
res::StarAlgebras.AlgebraElement,
|
||||
Q²::AbstractMatrix,
|
||||
cnstrs,
|
||||
)
|
||||
function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix, cnstrs)
|
||||
StarAlgebras.zero!(res)
|
||||
for (g, A_g) in cnstrs
|
||||
res[g] = dot(A_g, Q²)
|
||||
@ -50,14 +44,10 @@ function __sos_via_cnstr!(
|
||||
return res
|
||||
end
|
||||
|
||||
function compute_sos(
|
||||
A::StarAlgebras.StarAlgebra,
|
||||
Q::AbstractMatrix;
|
||||
augmented::Bool,
|
||||
)
|
||||
function compute_sos(A::StarAlgebras.StarAlgebra, Q::AbstractMatrix; augmented::Bool)
|
||||
Q² = Q' * Q
|
||||
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
|
||||
end
|
||||
|
||||
@ -66,7 +56,7 @@ function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
|
||||
suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
|
||||
|
||||
eq_sign = let T = eltype(residual)
|
||||
if T <: IntervalArithmetic.Interval
|
||||
if T <: Interval
|
||||
"∈"
|
||||
elseif T <: Union{Rational,Integer}
|
||||
"="
|
||||
@ -91,12 +81,13 @@ function sufficient_λ(
|
||||
order_unit::StarAlgebras.AlgebraElement,
|
||||
λ,
|
||||
sos::StarAlgebras.AlgebraElement;
|
||||
halfradius,
|
||||
halfradius
|
||||
)
|
||||
|
||||
@assert parent(elt) === parent(order_unit) == parent(sos)
|
||||
residual = (elt - λ * order_unit) - sos
|
||||
|
||||
return sufficient_λ(residual, λ; halfradius = halfradius)
|
||||
return sufficient_λ(residual, λ; halfradius=halfradius)
|
||||
end
|
||||
|
||||
function certify_solution(
|
||||
@ -105,27 +96,24 @@ function certify_solution(
|
||||
λ,
|
||||
Q::AbstractMatrix{<:AbstractFloat};
|
||||
halfradius,
|
||||
augmented = iszero(StarAlgebras.aug(elt)) &&
|
||||
iszero(StarAlgebras.aug(orderunit)),
|
||||
augmented=iszero(StarAlgebras.aug(elt)) && 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
|
||||
@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" λ
|
||||
|
||||
λ_flpoint = sufficient_λ(elt, orderunit, λ, sos; halfradius = halfradius)
|
||||
λ_flpoint = sufficient_λ(elt, orderunit, λ, sos, halfradius=halfradius)
|
||||
|
||||
if λ_flpoint ≤ 0
|
||||
return false, λ_flpoint
|
||||
end
|
||||
|
||||
λ_int = IntervalArithmetic.@interval(λ)
|
||||
Q_int = IntervalMatrices.IntervalMatrix([
|
||||
IntervalArithmetic.@interval(q) for q in Q
|
||||
])
|
||||
λ_int = @interval(λ)
|
||||
Q_int = [@interval(q) for q in Q]
|
||||
|
||||
check, sos_int = @time if should_we_augment
|
||||
@info("Projecting columns of Q to the augmentation ideal...")
|
||||
@ -135,16 +123,16 @@ function certify_solution(
|
||||
check_augmented || @error(
|
||||
"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
|
||||
else
|
||||
true, compute_sos(parent(elt), Q_int; augmented = augmented)
|
||||
true, compute_sos(parent(elt), Q_int, augmented=augmented)
|
||||
end
|
||||
|
||||
@info "Checking in $(eltype(sos_int)) arithmetic with" λ_int
|
||||
|
||||
λ_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
|
||||
|
||||
@ -28,12 +28,7 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
|
||||
size::Tuple{Int,Int}
|
||||
val::T
|
||||
|
||||
function ConstraintMatrix{T}(
|
||||
nzeros::AbstractArray{<:Integer},
|
||||
n,
|
||||
m,
|
||||
val,
|
||||
) where {T}
|
||||
function ConstraintMatrix{T}(nzeros::AbstractArray{<:Integer}, n, m, val) where {T}
|
||||
@assert n ≥ 1
|
||||
@assert m ≥ 1
|
||||
|
||||
@ -50,26 +45,15 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
|
||||
end
|
||||
end
|
||||
|
||||
function ConstraintMatrix(
|
||||
nzeros::AbstractArray{<:Integer},
|
||||
n,
|
||||
m,
|
||||
val::T,
|
||||
) where {T}
|
||||
return ConstraintMatrix{T}(nzeros, n, m, val)
|
||||
end
|
||||
ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} =
|
||||
ConstraintMatrix{T}(nzeros, n, m, val)
|
||||
|
||||
Base.size(cm::ConstraintMatrix) = cm.size
|
||||
|
||||
function __get_positive(cm::ConstraintMatrix, idx::Integer)
|
||||
return convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
|
||||
end
|
||||
function __get_negative(cm::ConstraintMatrix, idx::Integer)
|
||||
return convert(
|
||||
eltype(cm),
|
||||
cm.val * length(searchsorted(cm.neg, idx; rev = true)),
|
||||
)
|
||||
end
|
||||
__get_positive(cm::ConstraintMatrix, idx::Integer) =
|
||||
convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
|
||||
__get_negative(cm::ConstraintMatrix, idx::Integer) =
|
||||
convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx)))
|
||||
|
||||
Base.@propagate_inbounds function Base.getindex(
|
||||
cm::ConstraintMatrix,
|
||||
@ -83,29 +67,26 @@ Base.@propagate_inbounds function Base.getindex(
|
||||
return pos - neg
|
||||
end
|
||||
|
||||
struct NZPairsIter{T,I}
|
||||
m::ConstraintMatrix{T,I}
|
||||
struct NZPairsIter{T}
|
||||
m::ConstraintMatrix{T}
|
||||
end
|
||||
|
||||
Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
|
||||
Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
|
||||
|
||||
# TODO: iterate over (idx=>val) pairs combining vals
|
||||
function Base.iterate(
|
||||
itr::NZPairsIter,
|
||||
state::Tuple{Int,Nothing} = (1, nothing),
|
||||
)
|
||||
function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 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
|
||||
return idx => itr.m.val, (st, nothing)
|
||||
return idx => itr.m.val, (st, 1)
|
||||
end
|
||||
|
||||
function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int})
|
||||
k = iterate(itr.m.neg, state[2])
|
||||
function Base.iterate(itr::NZPairsIter, state::Int)
|
||||
k = iterate(itr.m.neg, state[1])
|
||||
isnothing(k) && return nothing
|
||||
idx, st = k
|
||||
return idx => -itr.m.val, (nothing, st)
|
||||
return idx => -itr.m.val, st
|
||||
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])
|
||||
return convert(eltype(cm), cm.val) * (pos - neg)
|
||||
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
|
||||
|
||||
@ -1,8 +1,7 @@
|
||||
## something about roots
|
||||
|
||||
function Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N}
|
||||
return Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
|
||||
end
|
||||
Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
|
||||
Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
|
||||
|
||||
function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
|
||||
if s.symbol === :A
|
||||
@ -52,28 +51,20 @@ function Adj(rootsystem::AbstractDict, subtype::Symbol)
|
||||
+,
|
||||
(
|
||||
Δα * Δβ for (α, Δα) in rootsystem for (β, Δβ) in rootsystem if
|
||||
Roots.classify_sub_root_system(roots, first(α), first(β)) == subtype
|
||||
);
|
||||
init = zero(first(values(rootsystem))),
|
||||
)
|
||||
end
|
||||
|
||||
Adj(rootsystem::AbstractDict) = sum(values(rootsystem))^2 - Sq(rootsystem)
|
||||
|
||||
function Sq(rootsystem::AbstractDict)
|
||||
return reduce(
|
||||
+,
|
||||
Δα^2 for (_, Δα) in rootsystem;
|
||||
init = zero(first(values(rootsystem))),
|
||||
Roots.classify_sub_root_system(
|
||||
roots,
|
||||
first(α),
|
||||
first(β),
|
||||
) == subtype
|
||||
),
|
||||
init=zero(first(values(rootsystem))),
|
||||
)
|
||||
end
|
||||
|
||||
function level(rootsystem, level::Integer)
|
||||
1 ≤ level ≤ 4 || throw("level is implemented only for i ∈{1,2,3,4}")
|
||||
level == 1 && return Adj(rootsystem, :C₁) # always positive
|
||||
level == 2 && return Adj(rootsystem, :A₁) +
|
||||
Adj(rootsystem, Symbol("C₁×C₁")) +
|
||||
Adj(rootsystem, :C₂) # C₂ is not positive
|
||||
level == 2 && return Adj(rootsystem, :A₁) + Adj(rootsystem, Symbol("C₁×C₁")) + Adj(rootsystem, :C₂) # C₂ is not positive
|
||||
level == 3 && return Adj(rootsystem, :A₂) + Adj(rootsystem, Symbol("A₁×C₁"))
|
||||
level == 4 && return Adj(rootsystem, Symbol("A₁×A₁")) # positive
|
||||
end
|
||||
|
||||
@ -7,28 +7,33 @@ end
|
||||
|
||||
function reconstruct(
|
||||
Ms::AbstractVector{<:AbstractMatrix},
|
||||
wbdec::WedderburnDecomposition,
|
||||
wbdec::WedderburnDecomposition;
|
||||
atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
|
||||
)
|
||||
n = __outer_dim(wbdec)
|
||||
res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
|
||||
res = similar(M, n, n)
|
||||
res = _reconstruct!(res, M, ds)
|
||||
return res
|
||||
reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol)
|
||||
end
|
||||
res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
|
||||
return res
|
||||
end
|
||||
|
||||
function _reconstruct!(
|
||||
function reconstruct!(
|
||||
res::AbstractMatrix,
|
||||
M::AbstractMatrix,
|
||||
ds::SymbolicWedderburn.DirectSummand,
|
||||
G,
|
||||
hom;
|
||||
atol=eps(real(eltype(ds))) * 10max(size(res)...)
|
||||
)
|
||||
res .= zero(eltype(res))
|
||||
if !iszero(M)
|
||||
U = SymbolicWedderburn.image_basis(ds)
|
||||
d = SymbolicWedderburn.degree(ds)
|
||||
res = (U' * M * U) .* d
|
||||
U = SymbolicWedderburn.image_basis(ds)
|
||||
d = SymbolicWedderburn.degree(ds)
|
||||
tmp = (U' * M * U) .* d
|
||||
|
||||
res = average!(res, tmp, G, hom)
|
||||
if eltype(res) <: AbstractFloat
|
||||
__droptol!(res, atol) # TODO: is this really necessary?!
|
||||
end
|
||||
return res
|
||||
end
|
||||
@ -47,22 +52,18 @@ function average!(
|
||||
res::AbstractMatrix,
|
||||
M::AbstractMatrix,
|
||||
G::Groups.Group,
|
||||
hom::SymbolicWedderburn.InducedActionHomomorphism{
|
||||
<:SymbolicWedderburn.ByPermutations,
|
||||
},
|
||||
hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations}
|
||||
)
|
||||
res .= zero(eltype(res))
|
||||
@assert size(M) == size(res)
|
||||
o = Groups.order(Int, G)
|
||||
for g in G
|
||||
p = SymbolicWedderburn.induce(hom, g)
|
||||
gext = SymbolicWedderburn.induce(hom, g)
|
||||
Threads.@threads for c in axes(res, 2)
|
||||
for r in axes(res, 1)
|
||||
if !iszero(M[r, c])
|
||||
res[r^p, c^p] += M[r, c] / o
|
||||
end
|
||||
res[r, c] += M[r^gext, c^gext]
|
||||
end
|
||||
end
|
||||
end
|
||||
o = Groups.order(Int, G)
|
||||
res ./= o
|
||||
return res
|
||||
end
|
||||
|
||||
214
src/roots.jl
214
src/roots.jl
@ -5,104 +5,111 @@ using LinearAlgebra
|
||||
|
||||
export Root, isproportional, isorthogonal, ~, ⟂
|
||||
|
||||
abstract type AbstractRoot{N,T} end # <: AbstractVector{T} ?
|
||||
abstract type AbstractRoot{N,T} end
|
||||
|
||||
ℓ₂length(r::AbstractRoot) = norm(r, 2)
|
||||
ambient_dim(r::AbstractRoot) = length(r)
|
||||
Base.:*(r::AbstractRoot, a::Number) = a * r
|
||||
struct Root{N,T} <: AbstractRoot{N,T}
|
||||
coord::SVector{N,T}
|
||||
end
|
||||
|
||||
Root(a) = Root(SVector(a...))
|
||||
|
||||
function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
|
||||
M == N || return false
|
||||
r.coord == s.coord || return false
|
||||
return true
|
||||
end
|
||||
|
||||
Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
|
||||
|
||||
Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord)
|
||||
Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord)
|
||||
Base.:-(r::Root{N}) where {N} = Root(-r.coord)
|
||||
|
||||
Base.:*(a::Number, r::Root) = Root(a * r.coord)
|
||||
Base.:*(r::Root, a::Number) = a * r
|
||||
|
||||
Base.length(r::AbstractRoot) = norm(r, 2)
|
||||
|
||||
LinearAlgebra.norm(r::Root, p::Real=2) = norm(r.coord, p)
|
||||
LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
|
||||
|
||||
cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
|
||||
|
||||
function isproportional(α::AbstractRoot, β::AbstractRoot)
|
||||
ambient_dim(α) == ambient_dim(β) || return false
|
||||
function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
|
||||
N == M || return false
|
||||
val = abs(cos_angle(α, β))
|
||||
return isapprox(val, one(val); atol = eps(one(val)))
|
||||
return isapprox(val, one(val), atol=eps(one(val)))
|
||||
end
|
||||
|
||||
function isorthogonal(α::AbstractRoot, β::AbstractRoot)
|
||||
ambient_dim(α) == ambient_dim(β) || return false
|
||||
function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
|
||||
N == M || return false
|
||||
val = cos_angle(α, β)
|
||||
return isapprox(val, zero(val); atol = eps(one(val)))
|
||||
return isapprox(val, zero(val), atol=eps(one(val)))
|
||||
end
|
||||
|
||||
function positive(roots::AbstractVector{<:AbstractRoot})
|
||||
isempty(roots) && return empty(roots)
|
||||
function _positive_direction(α::Root{N}) where {N}
|
||||
v = α.coord + 1 / (N * 100) * rand(N)
|
||||
return Root{N,Float64}(v / norm(v, 2))
|
||||
end
|
||||
|
||||
function positive(roots::AbstractVector{<:Root{N}}) where {N}
|
||||
# return those roots for which dot(α, Root([½, ¼, …])) > 0.0
|
||||
pd = _positive_direction(first(roots))
|
||||
return filter(α -> dot(α, pd) > 0.0, roots)
|
||||
end
|
||||
|
||||
function Base.show(io::IO, r::AbstractRoot)
|
||||
return print(io, "Root $(r.coord)")
|
||||
function Base.show(io::IO, r::Root)
|
||||
print(io, "Root$(r.coord)")
|
||||
end
|
||||
|
||||
function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
|
||||
l₂l = ℓ₂length(r)
|
||||
l = round(Int, l₂l) ≈ l₂l ? "$(round(Int, l₂l))" : "√$(round(Int, l₂l^2))"
|
||||
return print(io, "Root in ℝ^$(length(r)) of length $l\n", r.coord)
|
||||
function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
|
||||
lngth² = sum(x -> x^2, r.coord)
|
||||
l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
|
||||
print(io, "Root in ℝ^$N of length $l\n", r.coord)
|
||||
end
|
||||
|
||||
function reflection(α::AbstractRoot, β::AbstractRoot)
|
||||
return β - Int(2dot(α, β) // dot(α, α)) * α
|
||||
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
|
||||
𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
|
||||
𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
|
||||
|
||||
"""
|
||||
classify_root_system(α, β)
|
||||
Return the symbol of smallest system generated by roots `α` and `β`.
|
||||
|
||||
The classification is based only on roots length,
|
||||
proportionality/orthogonality and Cartan matrix.
|
||||
The classification is based only on roots length and
|
||||
proportionality/orthogonality.
|
||||
"""
|
||||
function classify_root_system(
|
||||
α::AbstractRoot,
|
||||
β::AbstractRoot,
|
||||
long::Tuple{Bool,Bool},
|
||||
)
|
||||
function classify_root_system(α::AbstractRoot, β::AbstractRoot)
|
||||
lα, lβ = length(α), length(β)
|
||||
if isproportional(α, β)
|
||||
if all(long)
|
||||
return :C₁
|
||||
elseif all(.!long) # both short
|
||||
if lα ≈ lβ ≈ √2
|
||||
return :A₁
|
||||
elseif lα ≈ lβ ≈ 2.0
|
||||
return :C₁
|
||||
else
|
||||
@error "Proportional roots of different length"
|
||||
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
|
||||
end
|
||||
elseif isorthogonal(α, β)
|
||||
if all(long)
|
||||
return Symbol("C₁×C₁")
|
||||
elseif all(.!long) # both short
|
||||
if lα ≈ lβ ≈ √2
|
||||
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₁")
|
||||
else
|
||||
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
|
||||
end
|
||||
else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
|
||||
a, b, c, d = abs.(cartan(α, β))
|
||||
@assert a == d == 2
|
||||
b, c = b < c ? (b, c) : (c, b)
|
||||
if b == c == 1
|
||||
if lα ≈ lβ ≈ √2
|
||||
return :A₂
|
||||
elseif b == 1 && c == 2
|
||||
elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
|
||||
return :C₂
|
||||
elseif b == 1 && c == 3
|
||||
@warn ":G₂? really?"
|
||||
return :G₂
|
||||
else
|
||||
@error a, b, c, d
|
||||
error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
function proportional_root_from_system(
|
||||
Ω::AbstractVector{<:AbstractRoot},
|
||||
α::AbstractRoot,
|
||||
)
|
||||
function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
|
||||
k = findfirst(v -> isproportional(α, v), Ω)
|
||||
if isnothing(k)
|
||||
error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
|
||||
@ -110,31 +117,25 @@ function proportional_root_from_system(
|
||||
return Ω[k]
|
||||
end
|
||||
|
||||
struct Plane{R<:AbstractRoot}
|
||||
struct Plane{R<:Root}
|
||||
v1::R
|
||||
v2::R
|
||||
vectors::Vector{R}
|
||||
end
|
||||
|
||||
function Plane(α::AbstractRoot, β::AbstractRoot)
|
||||
return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
|
||||
end
|
||||
Plane(α::R, β::R) where {R<:Root} =
|
||||
Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
|
||||
|
||||
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)
|
||||
end
|
||||
|
||||
function _islong(α::AbstractRoot, Ω)
|
||||
lα = ℓ₂length(α)
|
||||
return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
|
||||
end
|
||||
|
||||
function classify_sub_root_system(
|
||||
Ω::AbstractVector{<:AbstractRoot{N}},
|
||||
α::AbstractRoot{N},
|
||||
β::AbstractRoot{N},
|
||||
Ω::AbstractVector{<:Root{N}},
|
||||
α::Root{N},
|
||||
β::Root{N},
|
||||
) where {N}
|
||||
@assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
|
||||
|
||||
v = proportional_root_from_system(Ω, α)
|
||||
w = proportional_root_from_system(Ω, β)
|
||||
|
||||
@ -145,75 +146,32 @@ function classify_sub_root_system(
|
||||
l = length(subsystem)
|
||||
if l == 1
|
||||
x = first(subsystem)
|
||||
long = _islong(x, Ω)
|
||||
return classify_root_system(x, -x, (long, long))
|
||||
return classify_root_system(x, x)
|
||||
elseif l == 2
|
||||
x, y = subsystem
|
||||
return classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω)))
|
||||
return classify_root_system(subsystem...)
|
||||
elseif l == 3
|
||||
x, y, z = subsystem
|
||||
l1, l2, l3 = _islong(x, Ω), _islong(y, Ω), _islong(z, Ω)
|
||||
a = classify_root_system(x, y, (l1, l2))
|
||||
b = classify_root_system(y, z, (l2, l3))
|
||||
c = classify_root_system(x, z, (l1, l3))
|
||||
a = classify_root_system(subsystem[1], subsystem[2])
|
||||
b = classify_root_system(subsystem[2], subsystem[3])
|
||||
c = classify_root_system(subsystem[1], subsystem[3])
|
||||
|
||||
if :A₂ == a == b == c # it's only A₂
|
||||
if a == b == c # it's only A₂
|
||||
return a
|
||||
end
|
||||
|
||||
throw("Unknown subroot system! $((x,y,z))")
|
||||
elseif l == 4
|
||||
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
|
||||
C = (:C₂, Symbol("C₁×C₁"))
|
||||
if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
|
||||
return :C₂
|
||||
end
|
||||
elseif l == 4
|
||||
for i = 1:l
|
||||
for j = (i+1):l
|
||||
T = classify_root_system(subsystem[i], subsystem[j])
|
||||
T == :C₂ && return :C₂
|
||||
end
|
||||
end
|
||||
end
|
||||
@error "Unknown root subsystem generated by" α β
|
||||
throw("Unknown root system: $subsystem")
|
||||
end
|
||||
|
||||
## 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
|
||||
|
||||
23
src/solve.jl
23
src/solve.jl
@ -7,15 +7,12 @@ function setwarmstart!(model::JuMP.Model, warmstart)
|
||||
ct => JuMP.all_constraints(model, ct...) for
|
||||
ct in JuMP.list_of_constraint_types(model)
|
||||
)
|
||||
try
|
||||
JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
|
||||
|
||||
for (ct, idx) in pairs(constraint_map)
|
||||
JuMP.set_start_value.(idx, warmstart.slack[ct])
|
||||
JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
|
||||
end
|
||||
catch e
|
||||
@warn "Warmstarting was not succesfull!" e
|
||||
JuMP.set_start_value.(JuMP.all_variables(model), warmstart.primal)
|
||||
|
||||
for (ct, idx) in pairs(constraint_map)
|
||||
JuMP.set_start_value.(idx, warmstart.slack[ct])
|
||||
JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
|
||||
end
|
||||
return model
|
||||
end
|
||||
@ -31,10 +28,11 @@ function getwarmstart(model::JuMP.Model)
|
||||
slack = Dict(k => value.(v) for (k, v) in constraint_map)
|
||||
duals = Dict(k => JuMP.dual.(v) for (k, v) in constraint_map)
|
||||
|
||||
return (primal = primal, dual = duals, slack = slack)
|
||||
return (primal=primal, dual=duals, slack=slack)
|
||||
end
|
||||
|
||||
function solve(m::JuMP.Model, optimizer, warmstart = nothing)
|
||||
function solve(m::JuMP.Model, optimizer, warmstart=nothing)
|
||||
|
||||
JuMP.set_optimizer(m, optimizer)
|
||||
MOIU.attach_optimizer(m)
|
||||
|
||||
@ -47,7 +45,8 @@ function solve(m::JuMP.Model, optimizer, warmstart = nothing)
|
||||
return status, getwarmstart(m)
|
||||
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))
|
||||
|
||||
Base.flush(Base.stdout)
|
||||
@ -56,7 +55,7 @@ function solve(solverlog::String, m::JuMP.Model, optimizer, warmstart = nothing)
|
||||
redirect_stdout(logfile) do
|
||||
status, warmstart = solve(m, optimizer, warmstart)
|
||||
Base.Libc.flush_cstdio()
|
||||
return status, warmstart
|
||||
status, warmstart
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
205
src/sos_sdps.jl
205
src/sos_sdps.jl
@ -6,13 +6,15 @@ See also [sos_problem_primal](@ref).
|
||||
"""
|
||||
function sos_problem_dual(
|
||||
elt::StarAlgebras.AlgebraElement,
|
||||
order_unit::StarAlgebras.AlgebraElement = zero(elt);
|
||||
lower_bound = -Inf,
|
||||
order_unit::StarAlgebras.AlgebraElement=zero(elt);
|
||||
lower_bound=-Inf
|
||||
)
|
||||
@assert parent(elt) == parent(order_unit)
|
||||
algebra = parent(elt)
|
||||
moment_matrix = let m = algebra.mstructure
|
||||
[m[-i, j] for i in axes(m, 1), j in axes(m, 2)]
|
||||
mstructure = if StarAlgebras._istwisted(algebra.mstructure)
|
||||
algebra.mstructure
|
||||
else
|
||||
StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
|
||||
end
|
||||
|
||||
# 1 variable for every primal constraint
|
||||
@ -20,21 +22,60 @@ function sos_problem_dual(
|
||||
# Symmetrized:
|
||||
# 1 dual variable for every orbit of G acting on basis
|
||||
model = Model()
|
||||
JuMP.@variable(model, y[1:length(basis(algebra))])
|
||||
JuMP.@constraint(model, λ_dual, dot(order_unit, y) == 1)
|
||||
JuMP.@constraint(model, psd, y[moment_matrix] in PSDCone())
|
||||
@variable model y[1:length(basis(algebra))]
|
||||
@constraint model λ_dual dot(order_unit, y) == 1
|
||||
@constraint(model, psd, y[mstructure] in PSDCone())
|
||||
|
||||
if !isinf(lower_bound)
|
||||
throw("Not Implemented yet")
|
||||
JuMP.@variable(model, λ_ub_dual)
|
||||
JuMP.@objective(model, Min, dot(elt, y) + lower_bound * λ_ub_dual)
|
||||
@variable model λ_ub_dual
|
||||
@objective model Min dot(elt, y) + lower_bound * λ_ub_dual
|
||||
else
|
||||
JuMP.@objective(model, Min, dot(elt, y))
|
||||
@objective model Min dot(elt, y)
|
||||
end
|
||||
|
||||
return model
|
||||
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])
|
||||
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(
|
||||
elt::StarAlgebras.AlgebraElement,
|
||||
order_unit::StarAlgebras.AlgebraElement = zero(elt);
|
||||
upper_bound = Inf,
|
||||
augmented::Bool = iszero(StarAlgebras.aug(elt)) &&
|
||||
iszero(StarAlgebras.aug(order_unit)),
|
||||
order_unit::StarAlgebras.AlgebraElement=zero(elt);
|
||||
upper_bound=Inf,
|
||||
augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(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"
|
||||
end
|
||||
|
||||
A = constraints(parent(elt); augmented = augmented)
|
||||
A = constraints(parent(elt), augmented=augmented, twisted=true)
|
||||
|
||||
if !iszero(order_unit)
|
||||
λ = JuMP.@variable(model, λ)
|
||||
@ -95,11 +135,11 @@ end
|
||||
function invariant_constraint!(
|
||||
result::AbstractMatrix,
|
||||
basis::StarAlgebras.AbstractBasis,
|
||||
cnstrs::AbstractDict{K,<:ConstraintMatrix},
|
||||
cnstrs::AbstractDict{K,CM},
|
||||
invariant_vec::SparseVector,
|
||||
) where {K}
|
||||
) where {K,CM<:ConstraintMatrix}
|
||||
result .= zero(eltype(result))
|
||||
@inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
|
||||
for i in SparseArrays.nonzeroinds(invariant_vec)
|
||||
g = basis[i]
|
||||
A = cnstrs[g]
|
||||
for (idx, v) in nzpairs(A)
|
||||
@ -109,47 +149,25 @@ function invariant_constraint!(
|
||||
return result
|
||||
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)
|
||||
U = SymbolicWedderburn.image_basis(ds)
|
||||
return isapprox(U * U', I)
|
||||
end
|
||||
|
||||
function sos_problem_primal(
|
||||
sos_problem_primal(
|
||||
elt::StarAlgebras.AlgebraElement,
|
||||
wedderburn::WedderburnDecomposition;
|
||||
kwargs...,
|
||||
)
|
||||
return sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
|
||||
end
|
||||
kwargs...
|
||||
) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
|
||||
|
||||
function __fast_recursive_dot!(
|
||||
res::JuMP.AffExpr,
|
||||
Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
|
||||
Ms::AbstractArray{<:AbstractSparseMatrix},
|
||||
weights,
|
||||
Ms::AbstractArray{<:AbstractSparseMatrix};
|
||||
)
|
||||
@assert length(Ps) == length(Ms)
|
||||
|
||||
for (w, A, P) in zip(weights, Ms, Ps)
|
||||
for (A, P) in zip(Ms, Ps)
|
||||
iszero(Ms) && continue
|
||||
rows = rowvals(A)
|
||||
vals = nonzeros(A)
|
||||
@ -157,21 +175,13 @@ function __fast_recursive_dot!(
|
||||
for i in nzrange(A, cidx)
|
||||
ridx = rows[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
|
||||
return res
|
||||
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
|
||||
__show_itrs(n, total) = () -> [(Symbol("constraint"), "$n/$total")]
|
||||
|
||||
@ -179,38 +189,21 @@ function sos_problem_primal(
|
||||
elt::StarAlgebras.AlgebraElement,
|
||||
orderunit::StarAlgebras.AlgebraElement,
|
||||
wedderburn::WedderburnDecomposition;
|
||||
upper_bound = Inf,
|
||||
augmented = iszero(StarAlgebras.aug(elt)) &&
|
||||
iszero(StarAlgebras.aug(orderunit)),
|
||||
check_orthogonality = true,
|
||||
show_progress = false,
|
||||
upper_bound=Inf,
|
||||
augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
|
||||
check_orthogonality=true,
|
||||
show_progress=false
|
||||
)
|
||||
|
||||
@assert parent(elt) === parent(orderunit)
|
||||
if check_orthogonality
|
||||
if any(!isorth_projection, direct_summands(wedderburn))
|
||||
error(
|
||||
"Wedderburn decomposition contains a non-orthogonal projection",
|
||||
)
|
||||
error("Wedderburn decomposition contains a non-orthogonal projection")
|
||||
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)
|
||||
|
||||
# problem creation starts here
|
||||
model = JuMP.Model()
|
||||
if !feasibility_problem # add λ or not?
|
||||
λ = JuMP.@variable(model, λ)
|
||||
@ -224,52 +217,58 @@ function sos_problem_primal(
|
||||
end
|
||||
end
|
||||
|
||||
# semidefinite constraints as described by wedderburn
|
||||
Ps = map(direct_summands(wedderburn)) do ds
|
||||
P = map(direct_summands(wedderburn)) do ds
|
||||
dim = size(ds, 1)
|
||||
P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
|
||||
JuMP.@constraint(model, P in PSDCone())
|
||||
return P
|
||||
P
|
||||
end
|
||||
|
||||
begin # Ms are preallocated for the constraints loop
|
||||
begin # preallocating
|
||||
T = eltype(wedderburn)
|
||||
Ms = [spzeros.(T, size(p)...) for p in Ps]
|
||||
_eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
|
||||
Ms = [spzeros.(T, size(p)...) for p in P]
|
||||
M_orb = zeros(T, size(parent(elt).mstructure)...)
|
||||
end
|
||||
|
||||
X = StarAlgebras.coeffs(elt)
|
||||
U = StarAlgebras.coeffs(orderunit)
|
||||
X = convert(Vector{T}, StarAlgebras.coeffs(elt))
|
||||
U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
|
||||
|
||||
# defining constraints based on the multiplicative structure
|
||||
cnstrs = constraints(parent(elt); augmented = augmented)
|
||||
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))
|
||||
ProgressMeter.next!(prog; showvalues = __show_itrs(i, prog.n))
|
||||
augmented && i == id_one && continue
|
||||
# i == 500 && break
|
||||
ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n))
|
||||
|
||||
x = dot(X, 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
|
||||
JuMP.@constraint(model, x == _dot(Ps, Ms))
|
||||
JuMP.@constraint(
|
||||
model,
|
||||
x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
|
||||
)
|
||||
else
|
||||
JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms))
|
||||
JuMP.@constraint(
|
||||
model,
|
||||
x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
|
||||
)
|
||||
end
|
||||
end
|
||||
ProgressMeter.finish!(prog)
|
||||
return model, Ps
|
||||
return model, P
|
||||
end
|
||||
|
||||
@ -1,8 +1,9 @@
|
||||
@testset "1703.09680 Examples" begin
|
||||
|
||||
@testset "SL(2,Z)" begin
|
||||
N = 2
|
||||
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
|
||||
RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -14,15 +15,15 @@
|
||||
|
||||
status, certified, λ = check_positivity(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = ub,
|
||||
halfradius = 2,
|
||||
optimizer = scs_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 5_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
),
|
||||
unit,
|
||||
upper_bound=ub,
|
||||
halfradius=2,
|
||||
optimizer=scs_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=5_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
)
|
||||
)
|
||||
|
||||
@test status == JuMP.ALMOST_OPTIMAL
|
||||
@ -32,10 +33,8 @@
|
||||
|
||||
@testset "SL(3,F₅)" begin
|
||||
N = 3
|
||||
G = MatrixGroups.SpecialLinearGroup{N}(
|
||||
SymbolicWedderburn.Characters.FiniteFields.GF{5},
|
||||
)
|
||||
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
|
||||
G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{5})
|
||||
RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
|
||||
|
||||
Δ = let RG = RG, S = S
|
||||
RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -47,15 +46,15 @@
|
||||
|
||||
status, certified, λ = check_positivity(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = ub,
|
||||
halfradius = 2,
|
||||
optimizer = scs_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 5_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
),
|
||||
unit,
|
||||
upper_bound=ub,
|
||||
halfradius=2,
|
||||
optimizer=scs_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=5_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
)
|
||||
)
|
||||
|
||||
@test status == JuMP.OPTIMAL
|
||||
@ -63,15 +62,12 @@
|
||||
@test λ > 1
|
||||
|
||||
m = PropertyT.sos_problem_dual(elt, unit)
|
||||
PropertyT.solve(
|
||||
m,
|
||||
cosmo_optimizer(;
|
||||
eps = 1e-6,
|
||||
max_iters = 5_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
),
|
||||
)
|
||||
PropertyT.solve(m, cosmo_optimizer(
|
||||
eps=1e-6,
|
||||
max_iters=5_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
))
|
||||
|
||||
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
|
||||
@test JuMP.objective_value(m) ≈ 1.5 atol = 1e-2
|
||||
@ -80,7 +76,7 @@
|
||||
@testset "SAut(F₂)" begin
|
||||
N = 2
|
||||
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
|
||||
RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -92,46 +88,53 @@
|
||||
|
||||
status, certified, λ = check_positivity(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = ub,
|
||||
halfradius = 2,
|
||||
optimizer = scs_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 5_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
),
|
||||
unit,
|
||||
upper_bound=ub,
|
||||
halfradius=2,
|
||||
optimizer=scs_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=5_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
)
|
||||
)
|
||||
|
||||
@test status == JuMP.ALMOST_OPTIMAL
|
||||
@test λ < 0
|
||||
@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(
|
||||
sos_problem,
|
||||
cosmo_optimizer(;
|
||||
eps = 1e-7,
|
||||
max_iters = 10_000,
|
||||
accel = 0,
|
||||
alpha = 1.9,
|
||||
),
|
||||
cosmo_optimizer(
|
||||
eps=1e-7,
|
||||
max_iters=10_000,
|
||||
accel=0,
|
||||
alpha=1.9,
|
||||
)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
P = JuMP.value.(sos_problem[:P])
|
||||
Q = real.(sqrt(P))
|
||||
certified, λ_cert =
|
||||
PropertyT.certify_solution(elt, zero(elt), 0.0, Q; halfradius = 2)
|
||||
certified, λ_cert = PropertyT.certify_solution(
|
||||
elt,
|
||||
zero(elt),
|
||||
0.0,
|
||||
Q,
|
||||
halfradius=2,
|
||||
)
|
||||
@test !certified
|
||||
@test λ_cert < 0
|
||||
|
||||
end
|
||||
|
||||
@testset "SL(3,Z) has (T)" begin
|
||||
n = 3
|
||||
|
||||
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)
|
||||
|
||||
@ -142,18 +145,18 @@
|
||||
|
||||
opt_problem = PropertyT.sos_problem_primal(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = ub,
|
||||
augmented = false,
|
||||
unit,
|
||||
upper_bound=ub,
|
||||
augmented=false,
|
||||
)
|
||||
|
||||
status, _ = PropertyT.solve(
|
||||
opt_problem,
|
||||
cosmo_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 10_000,
|
||||
accel = 0,
|
||||
alpha = 1.5,
|
||||
cosmo_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=10_000,
|
||||
accel=0,
|
||||
alpha=1.5,
|
||||
),
|
||||
)
|
||||
|
||||
@ -167,13 +170,13 @@
|
||||
elt,
|
||||
unit,
|
||||
λ,
|
||||
Q;
|
||||
halfradius = 2,
|
||||
augmented = false,
|
||||
Q,
|
||||
halfradius=2,
|
||||
augmented=false,
|
||||
)
|
||||
|
||||
@test certified
|
||||
@test isapprox(λ_cert, λ, rtol = 1e-5)
|
||||
@test isapprox(λ_cert, λ, rtol=1e-5)
|
||||
end
|
||||
|
||||
@testset "augmented formulation" begin
|
||||
@ -183,18 +186,18 @@
|
||||
|
||||
opt_problem = PropertyT.sos_problem_primal(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = ub,
|
||||
augmented = true,
|
||||
unit,
|
||||
upper_bound=ub,
|
||||
augmented=true,
|
||||
)
|
||||
|
||||
status, _ = PropertyT.solve(
|
||||
opt_problem,
|
||||
scs_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 10_000,
|
||||
accel = -10,
|
||||
alpha = 1.5,
|
||||
scs_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=10_000,
|
||||
accel=-10,
|
||||
alpha=1.5,
|
||||
),
|
||||
)
|
||||
|
||||
@ -207,13 +210,13 @@
|
||||
elt,
|
||||
unit,
|
||||
λ,
|
||||
Q;
|
||||
halfradius = 2,
|
||||
augmented = true,
|
||||
Q,
|
||||
halfradius=2,
|
||||
augmented=true,
|
||||
)
|
||||
|
||||
@test certified
|
||||
@test isapprox(λ_cert, λ, rtol = 1e-5)
|
||||
@test isapprox(λ_cert, λ, rtol=1e-5)
|
||||
@test λ_cert > 2 // 10
|
||||
end
|
||||
end
|
||||
|
||||
@ -1,9 +1,9 @@
|
||||
@testset "1712.07167 Examples" begin
|
||||
|
||||
@testset "SAut(F₃)" begin
|
||||
N = 3
|
||||
G = SpecialAutomorphismGroup(FreeGroup(N))
|
||||
@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)
|
||||
|
||||
P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
|
||||
Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
|
||||
@ -15,7 +15,6 @@
|
||||
basis(RG),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
|
||||
)
|
||||
@info wd
|
||||
|
||||
Δ = let RG = RG, S = S
|
||||
RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -28,14 +27,14 @@
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
unit,
|
||||
wd;
|
||||
upper_bound = ub,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(;
|
||||
eps = 1e-7,
|
||||
max_iters = 10_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
wd,
|
||||
upper_bound=ub,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(
|
||||
eps=1e-7,
|
||||
max_iters=10_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
),
|
||||
)
|
||||
|
||||
@ -48,8 +47,7 @@
|
||||
n = 3
|
||||
|
||||
SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
|
||||
@info "running tests for" SL
|
||||
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)
|
||||
|
||||
@ -64,7 +62,6 @@
|
||||
basis(RSL),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
|
||||
)
|
||||
@info wd
|
||||
|
||||
elt = Δ^2
|
||||
unit = Δ
|
||||
@ -74,8 +71,8 @@
|
||||
elt,
|
||||
unit,
|
||||
wd,
|
||||
upper_bound = ub,
|
||||
augmented = false,
|
||||
upper_bound=ub,
|
||||
augmented=false,
|
||||
)
|
||||
|
||||
wdfl = SymbolicWedderburn.WedderburnDecomposition(
|
||||
@ -89,18 +86,18 @@
|
||||
model, varP = PropertyT.sos_problem_primal(
|
||||
elt,
|
||||
unit,
|
||||
wdfl;
|
||||
upper_bound = ub,
|
||||
augmented = false,
|
||||
wdfl,
|
||||
upper_bound=ub,
|
||||
augmented=false,
|
||||
)
|
||||
|
||||
status, warm = PropertyT.solve(
|
||||
model,
|
||||
cosmo_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 20_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
cosmo_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=20_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
),
|
||||
)
|
||||
|
||||
@ -108,34 +105,34 @@
|
||||
|
||||
status, _ = PropertyT.solve(
|
||||
model,
|
||||
scs_optimizer(;
|
||||
eps = 1e-10,
|
||||
max_iters = 100,
|
||||
accel = -20,
|
||||
alpha = 1.2,
|
||||
scs_optimizer(
|
||||
eps=1e-10,
|
||||
max_iters=100,
|
||||
accel=-20,
|
||||
alpha=1.2,
|
||||
),
|
||||
warm,
|
||||
warm
|
||||
)
|
||||
|
||||
@test status == JuMP.OPTIMAL
|
||||
|
||||
Q = @time let varP = varP
|
||||
Qs = map(varP) do P
|
||||
return real.(sqrt(JuMP.value.(P)))
|
||||
real.(sqrt(JuMP.value.(P)))
|
||||
end
|
||||
PropertyT.reconstruct(Qs, wdfl)
|
||||
end
|
||||
λ = 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(
|
||||
elt,
|
||||
unit,
|
||||
λ,
|
||||
Q;
|
||||
halfradius = 2,
|
||||
augmented = false,
|
||||
Q,
|
||||
halfradius=2,
|
||||
augmented=false,
|
||||
)
|
||||
|
||||
@test certified
|
||||
@ -158,23 +155,22 @@
|
||||
basis(RSL),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
|
||||
)
|
||||
@info wdfl
|
||||
|
||||
opt_problem, varP = PropertyT.sos_problem_primal(
|
||||
elt,
|
||||
unit,
|
||||
wdfl;
|
||||
upper_bound = ub,
|
||||
wdfl,
|
||||
upper_bound=ub,
|
||||
# augmented = true # since both elt and unit are augmented
|
||||
)
|
||||
|
||||
status, _ = PropertyT.solve(
|
||||
opt_problem,
|
||||
scs_optimizer(;
|
||||
eps = 1e-8,
|
||||
max_iters = 20_000,
|
||||
accel = 0,
|
||||
alpha = 1.9,
|
||||
scs_optimizer(
|
||||
eps=1e-8,
|
||||
max_iters=20_000,
|
||||
accel=0,
|
||||
alpha=1.9,
|
||||
),
|
||||
)
|
||||
|
||||
@ -182,7 +178,7 @@
|
||||
|
||||
Q = @time let varP = varP
|
||||
Qs = map(varP) do P
|
||||
return real.(sqrt(JuMP.value.(P)))
|
||||
real.(sqrt(JuMP.value.(P)))
|
||||
end
|
||||
PropertyT.reconstruct(Qs, wdfl)
|
||||
end
|
||||
@ -191,8 +187,8 @@
|
||||
elt,
|
||||
unit,
|
||||
JuMP.objective_value(opt_problem),
|
||||
Q;
|
||||
halfradius = 2,
|
||||
Q,
|
||||
halfradius=2,
|
||||
# augmented = true # since both elt and unit are augmented
|
||||
)
|
||||
|
||||
|
||||
@ -21,9 +21,8 @@ using SparseArrays
|
||||
@testset "Sq, Adj, Op in SL(4,Z)" begin
|
||||
N = 4
|
||||
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
|
||||
RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -40,7 +39,6 @@ using SparseArrays
|
||||
basis(RG),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
|
||||
)
|
||||
@info wd
|
||||
ivs = SymbolicWedderburn.invariant_vectors(wd)
|
||||
|
||||
sq, adj, op = PropertyT.SqAdjOp(RG, N)
|
||||
@ -84,7 +82,7 @@ using SparseArrays
|
||||
@testset "SAut(F₃)" begin
|
||||
n = 3
|
||||
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)
|
||||
|
||||
@test sq(one(G)) == 216
|
||||
@ -100,8 +98,7 @@ end
|
||||
n = 3
|
||||
|
||||
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)
|
||||
|
||||
Δ = RG(length(S)) - sum(RG(s) for s in S)
|
||||
|
||||
@ -116,7 +113,6 @@ end
|
||||
basis(RG),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
|
||||
)
|
||||
@info wd
|
||||
|
||||
sq, adj, op = PropertyT.SqAdjOp(RG, n)
|
||||
|
||||
@ -127,10 +123,10 @@ end
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
Δ,
|
||||
wd;
|
||||
upper_bound = UB,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
|
||||
wd,
|
||||
upper_bound=UB,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
@test certified
|
||||
@ -144,10 +140,10 @@ end
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
Δ,
|
||||
wd;
|
||||
upper_bound = UB,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
|
||||
wd,
|
||||
upper_bound=UB,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
@test certified
|
||||
@ -156,11 +152,10 @@ end
|
||||
m, _ = PropertyT.sos_problem_primal(elt, wd)
|
||||
PropertyT.solve(
|
||||
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
|
||||
(JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
|
||||
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL, JuMP.ITERATION_LIMIT)
|
||||
@test abs(JuMP.objective_value(m)) < 1e-3
|
||||
end
|
||||
|
||||
@ -173,10 +168,10 @@ end
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
Δ,
|
||||
wd;
|
||||
upper_bound = UB,
|
||||
halfradius = 2,
|
||||
optimizer = scs_optimizer(; accel = 50, alpha = 1.9),
|
||||
wd,
|
||||
upper_bound=UB,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
@test !certified
|
||||
@ -188,8 +183,7 @@ end
|
||||
n = 4
|
||||
|
||||
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)
|
||||
|
||||
Δ = RG(length(S)) - sum(RG(s) for s in S)
|
||||
|
||||
@ -204,7 +198,6 @@ end
|
||||
basis(RG),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
|
||||
)
|
||||
@info wd
|
||||
|
||||
sq, adj, op = PropertyT.SqAdjOp(RG, n)
|
||||
|
||||
@ -215,10 +208,10 @@ end
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
Δ,
|
||||
wd;
|
||||
upper_bound = UB,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
|
||||
wd,
|
||||
upper_bound=UB,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
@test certified
|
||||
@ -232,10 +225,10 @@ end
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
Δ,
|
||||
wd;
|
||||
upper_bound = UB,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
|
||||
wd,
|
||||
upper_bound=UB,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
@test certified
|
||||
@ -249,10 +242,10 @@ end
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
Δ,
|
||||
wd;
|
||||
upper_bound = UB,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
|
||||
wd,
|
||||
upper_bound=UB,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(accel=50, alpha=1.9)
|
||||
)
|
||||
@test status == JuMP.OPTIMAL
|
||||
@test !certified
|
||||
|
||||
@ -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
|
||||
@ -3,7 +3,7 @@ function test_action(basis, group, act)
|
||||
return @testset "action definition" begin
|
||||
@test all(basis) do b
|
||||
e = one(group)
|
||||
return action(act, e, b) == b
|
||||
action(act, e, b) == b
|
||||
end
|
||||
|
||||
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)
|
||||
x = action(act, h, action(act, g, a))
|
||||
y = action(act, g * h, a)
|
||||
return x == y
|
||||
x == y
|
||||
end
|
||||
|
||||
if act isa SymbolicWedderburn.ByPermutations
|
||||
@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
|
||||
@ -44,18 +44,21 @@ end
|
||||
## Testing
|
||||
|
||||
@testset "Actions on SL(3,ℤ)" begin
|
||||
|
||||
n = 3
|
||||
|
||||
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
|
||||
|
||||
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
|
||||
ΓpA = PropertyT.action_by_conjugation(SL, Γ)
|
||||
|
||||
test_action(basis(RSL), Γ, ΓpA)
|
||||
|
||||
@testset "mps is successful" begin
|
||||
|
||||
charsΓ =
|
||||
SymbolicWedderburn.Character{
|
||||
Rational{Int},
|
||||
@ -63,9 +66,9 @@ end
|
||||
|
||||
RΓ = SymbolicWedderburn._group_algebra(Γ)
|
||||
|
||||
@time mps, ranks =
|
||||
@time mps, simple =
|
||||
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
|
||||
@test all(isone, ranks)
|
||||
@test all(simple)
|
||||
end
|
||||
|
||||
@testset "Wedderburn decomposition" begin
|
||||
@ -74,17 +77,17 @@ end
|
||||
Γ,
|
||||
ΓpA,
|
||||
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 SymbolicWedderburn.size.(direct_summands(wd), 1) ==
|
||||
[40, 23, 18]
|
||||
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [40, 23, 18]
|
||||
@test all(issimple, direct_summands(wd))
|
||||
end
|
||||
end
|
||||
|
||||
@testset "Wreath action" begin
|
||||
|
||||
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
|
||||
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
|
||||
end
|
||||
@ -94,6 +97,7 @@ end
|
||||
test_action(basis(RSL), Γ, ΓpA)
|
||||
|
||||
@testset "mps is successful" begin
|
||||
|
||||
charsΓ =
|
||||
SymbolicWedderburn.Character{
|
||||
Rational{Int},
|
||||
@ -101,9 +105,9 @@ end
|
||||
|
||||
RΓ = SymbolicWedderburn._group_algebra(Γ)
|
||||
|
||||
@time mps, ranks =
|
||||
@time mps, simple =
|
||||
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
|
||||
@test all(isone, ranks)
|
||||
@test all(simple)
|
||||
end
|
||||
|
||||
@testset "Wedderburn decomposition" begin
|
||||
@ -112,30 +116,32 @@ end
|
||||
Γ,
|
||||
ΓpA,
|
||||
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 SymbolicWedderburn.size.(direct_summands(wd), 1) ==
|
||||
[14, 9, 6, 14, 12]
|
||||
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [14, 9, 6, 14, 12]
|
||||
@test all(issimple, direct_summands(wd))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@testset "Actions on SAut(F4)" begin
|
||||
|
||||
n = 4
|
||||
|
||||
SAutFn = SpecialAutomorphismGroup(FreeGroup(n))
|
||||
RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn; halfradius = 1)
|
||||
RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn, halfradius=1, twisted=true)
|
||||
|
||||
@testset "Permutation action" begin
|
||||
|
||||
Γ = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
|
||||
ΓpA = PropertyT.action_by_conjugation(SAutFn, Γ)
|
||||
|
||||
test_action(basis(RSAutFn), Γ, ΓpA)
|
||||
|
||||
@testset "mps is successful" begin
|
||||
|
||||
charsΓ =
|
||||
SymbolicWedderburn.Character{
|
||||
Rational{Int},
|
||||
@ -143,9 +149,9 @@ end
|
||||
|
||||
RΓ = SymbolicWedderburn._group_algebra(Γ)
|
||||
|
||||
@time mps, ranks =
|
||||
@time mps, simple =
|
||||
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
|
||||
@test all(isone, ranks)
|
||||
@test all(simple)
|
||||
end
|
||||
|
||||
@testset "Wedderburn decomposition" begin
|
||||
@ -154,17 +160,17 @@ end
|
||||
Γ,
|
||||
ΓpA,
|
||||
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 SymbolicWedderburn.size.(direct_summands(wd), 1) ==
|
||||
[4, 8, 5, 4]
|
||||
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [4, 8, 5, 4]
|
||||
@test all(issimple, direct_summands(wd))
|
||||
end
|
||||
end
|
||||
|
||||
@testset "Wreath action" begin
|
||||
|
||||
Γ = let P = PermGroup(perm"(1,2)", Perm(circshift(1:n, -1)))
|
||||
PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
|
||||
end
|
||||
@ -174,6 +180,7 @@ end
|
||||
test_action(basis(RSAutFn), Γ, ΓpA)
|
||||
|
||||
@testset "mps is successful" begin
|
||||
|
||||
charsΓ =
|
||||
SymbolicWedderburn.Character{
|
||||
Rational{Int},
|
||||
@ -181,9 +188,9 @@ end
|
||||
|
||||
RΓ = SymbolicWedderburn._group_algebra(Γ)
|
||||
|
||||
@time mps, ranks =
|
||||
@time mps, simple =
|
||||
SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
|
||||
@test all(isone, ranks)
|
||||
@test all(simple)
|
||||
end
|
||||
|
||||
@testset "Wedderburn decomposition" begin
|
||||
@ -192,12 +199,11 @@ end
|
||||
Γ,
|
||||
ΓpA,
|
||||
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 SymbolicWedderburn.size.(direct_summands(wd), 1) ==
|
||||
[1, 1, 2, 2, 1, 2, 2, 1]
|
||||
@test SymbolicWedderburn.size.(direct_summands(wd), 1) == [1, 1, 2, 2, 1, 2, 2, 1]
|
||||
@test all(issimple, direct_summands(wd))
|
||||
end
|
||||
end
|
||||
|
||||
@ -1,12 +1,6 @@
|
||||
function check_positivity(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = Inf,
|
||||
halfradius = 2,
|
||||
optimizer,
|
||||
)
|
||||
function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
|
||||
@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)
|
||||
P = JuMP.value.(sos_problem[:P])
|
||||
@ -15,23 +9,16 @@ function check_positivity(
|
||||
elt,
|
||||
unit,
|
||||
JuMP.objective_value(sos_problem),
|
||||
Q;
|
||||
halfradius = halfradius,
|
||||
Q,
|
||||
halfradius=halfradius,
|
||||
)
|
||||
return status, certified, λ_cert
|
||||
end
|
||||
|
||||
function check_positivity(
|
||||
elt,
|
||||
unit,
|
||||
wd;
|
||||
upper_bound = Inf,
|
||||
halfradius = 2,
|
||||
optimizer,
|
||||
)
|
||||
function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimizer)
|
||||
@assert aug(elt) == aug(unit) == 0
|
||||
@time sos_problem, Ps =
|
||||
PropertyT.sos_problem_primal(elt, unit, wd; upper_bound = upper_bound)
|
||||
PropertyT.sos_problem_primal(elt, unit, wd, upper_bound=upper_bound)
|
||||
|
||||
@time status, _ = PropertyT.solve(sos_problem, optimizer)
|
||||
|
||||
@ -42,7 +29,13 @@ function check_positivity(
|
||||
|
||||
λ = JuMP.value(sos_problem[:λ])
|
||||
|
||||
certified, λ_cert =
|
||||
PropertyT.certify_solution(elt, unit, λ, Q; halfradius = halfradius)
|
||||
certified, λ_cert = PropertyT.certify_solution(
|
||||
elt,
|
||||
unit,
|
||||
λ,
|
||||
Q,
|
||||
halfradius=halfradius
|
||||
)
|
||||
return status, certified, λ_cert
|
||||
end
|
||||
|
||||
|
||||
@ -1,17 +1,12 @@
|
||||
@testset "ConstraintMatrix" begin
|
||||
@test PropertyT.ConstraintMatrix{Float64}(
|
||||
[-1, 2, -1, 1, 4, 2, 6],
|
||||
3,
|
||||
2,
|
||||
π,
|
||||
) isa AbstractMatrix
|
||||
@test PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π) isa AbstractMatrix
|
||||
|
||||
cm = PropertyT.ConstraintMatrix{Float64}([-1, 2, -1, 1, 4, 2, 6], 3, 2, π)
|
||||
|
||||
@test cm == Float64[
|
||||
-π π
|
||||
2π 0
|
||||
0 π
|
||||
2π 0.0
|
||||
0.0 π
|
||||
]
|
||||
|
||||
@test collect(PropertyT.nzpairs(cm)) == [
|
||||
@ -23,7 +18,4 @@
|
||||
1 => -3.141592653589793
|
||||
1 => -3.141592653589793
|
||||
]
|
||||
|
||||
@test PropertyT.ConstraintMatrix{Float64}([-9:-1; 1:9], 3, 3, 1.0) ==
|
||||
zeros(3, 3)
|
||||
end
|
||||
|
||||
@ -1,9 +1,10 @@
|
||||
@testset "Adj via grading" begin
|
||||
|
||||
@testset "SL(n,Z) & Aut(F₄)" begin
|
||||
n = 4
|
||||
halfradius = 1
|
||||
SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
|
||||
RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = halfradius)
|
||||
RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=halfradius, twisted=true)
|
||||
|
||||
Δ = RSL(length(S)) - sum(RSL(s) for s in S)
|
||||
|
||||
@ -21,9 +22,10 @@
|
||||
@test PropertyT.Adj(Δs, :A₂) == adj
|
||||
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
|
||||
|
||||
|
||||
halfradius = 1
|
||||
G = SpecialAutomorphismGroup(FreeGroup(n))
|
||||
RG, S, sizes = PropertyT.group_algebra(G; halfradius = halfradius)
|
||||
RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true)
|
||||
|
||||
Δ = RG(length(S)) - sum(RG(s) for s in S)
|
||||
|
||||
@ -42,6 +44,7 @@
|
||||
@test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
|
||||
end
|
||||
|
||||
|
||||
@testset "Symplectic group" begin
|
||||
@testset "Sp2(ℤ)" begin
|
||||
genus = 2
|
||||
@ -49,8 +52,7 @@
|
||||
|
||||
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
|
||||
|
||||
RSpN, S_sp, sizes_sp =
|
||||
PropertyT.group_algebra(SpN; halfradius = halfradius)
|
||||
RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
|
||||
|
||||
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
|
||||
Δ = RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -71,8 +73,7 @@
|
||||
|
||||
SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
|
||||
|
||||
RSpN, S_sp, sizes_sp =
|
||||
PropertyT.group_algebra(SpN; halfradius = halfradius)
|
||||
RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
|
||||
|
||||
Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
|
||||
Δ = RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -85,14 +86,9 @@
|
||||
end
|
||||
|
||||
@testset "Adj numerics for genus=$genus" begin
|
||||
|
||||
all_subtypes = (
|
||||
:A₁,
|
||||
:C₁,
|
||||
Symbol("A₁×A₁"),
|
||||
Symbol("C₁×C₁"),
|
||||
Symbol("A₁×C₁"),
|
||||
:A₂,
|
||||
:C₂,
|
||||
:A₁, :C₁, Symbol("A₁×A₁"), Symbol("C₁×C₁"), Symbol("A₁×C₁"), :A₂, :C₂
|
||||
)
|
||||
|
||||
@test PropertyT.Adj(Δs, :A₂)[one(SpN)] == 384
|
||||
@ -105,8 +101,8 @@
|
||||
@test isinteger(PropertyT.Adj(Δs, subtype)[one(SpN)] / 16)
|
||||
end
|
||||
end
|
||||
@test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) ==
|
||||
Δ^2
|
||||
@test sum(PropertyT.Adj(Δs, subtype) for subtype in all_subtypes) == Δ^2
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@ -4,12 +4,12 @@ import JuMP
|
||||
import SCS
|
||||
|
||||
function scs_optimizer(;
|
||||
accel = 10,
|
||||
alpha = 1.5,
|
||||
eps = 1e-9,
|
||||
max_iters = 100_000,
|
||||
verbose = true,
|
||||
linear_solver = SCS.DirectSolver,
|
||||
accel=10,
|
||||
alpha=1.5,
|
||||
eps=1e-9,
|
||||
max_iters=100_000,
|
||||
verbose=true,
|
||||
linear_solver=SCS.DirectSolver
|
||||
)
|
||||
return JuMP.optimizer_with_attributes(
|
||||
SCS.Optimizer,
|
||||
@ -28,27 +28,27 @@ end
|
||||
import COSMO
|
||||
|
||||
function cosmo_optimizer(;
|
||||
accel = 15,
|
||||
alpha = 1.6,
|
||||
eps = 1e-9,
|
||||
max_iters = 100_000,
|
||||
verbose = true,
|
||||
verbose_timing = verbose,
|
||||
decompose = false,
|
||||
accel=15,
|
||||
alpha=1.6,
|
||||
eps=1e-9,
|
||||
max_iters=100_000,
|
||||
verbose=true,
|
||||
verbose_timing=verbose,
|
||||
decompose=false
|
||||
)
|
||||
return JuMP.optimizer_with_attributes(
|
||||
COSMO.Optimizer,
|
||||
"accelerator" =>
|
||||
COSMO.with_options(COSMO.AndersonAccelerator; mem = max(accel, 2)),
|
||||
"accelerator" => COSMO.with_options(
|
||||
COSMO.AndersonAccelerator,
|
||||
mem=max(accel, 2)
|
||||
),
|
||||
"alpha" => alpha,
|
||||
"decompose" => decompose,
|
||||
"eps_abs" => eps,
|
||||
"eps_rel" => eps,
|
||||
"eps_prim_inf" => eps,
|
||||
"eps_dual_inf" => eps,
|
||||
"max_iter" => max_iters,
|
||||
"verbose" => verbose,
|
||||
"verbose_timing" => verbose_timing,
|
||||
"check_termination" => 250,
|
||||
"check_termination" => 200,
|
||||
)
|
||||
end
|
||||
|
||||
@ -1,12 +1,11 @@
|
||||
@testset "Quick tests" begin
|
||||
|
||||
@testset "SL(2,F₇)" begin
|
||||
N = 2
|
||||
p = 7
|
||||
halfradius = 3
|
||||
G = MatrixGroups.SpecialLinearGroup{N}(
|
||||
SymbolicWedderburn.Characters.FiniteFields.GF{p},
|
||||
)
|
||||
RG, S, sizes = PropertyT.group_algebra(G; halfradius = 3)
|
||||
G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p})
|
||||
RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true)
|
||||
|
||||
Δ = let RG = RG, S = S
|
||||
RG(length(S)) - sum(RG(s) for s in S)
|
||||
@ -19,15 +18,15 @@
|
||||
@testset "standard formulation" begin
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
unit;
|
||||
upper_bound = ub,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(;
|
||||
eps = 1e-7,
|
||||
max_iters = 5_000,
|
||||
accel = 50,
|
||||
alpha = 1.95,
|
||||
),
|
||||
unit,
|
||||
upper_bound=ub,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(
|
||||
eps=1e-7,
|
||||
max_iters=5_000,
|
||||
accel=50,
|
||||
alpha=1.95,
|
||||
)
|
||||
)
|
||||
|
||||
@test status == JuMP.OPTIMAL
|
||||
@ -35,18 +34,14 @@
|
||||
@test λ_cert > 5857 // 10000
|
||||
|
||||
m = PropertyT.sos_problem_dual(elt, unit)
|
||||
PropertyT.solve(
|
||||
m,
|
||||
cosmo_optimizer(;
|
||||
eps = 1e-7,
|
||||
max_iters = 10_000,
|
||||
accel = 50,
|
||||
alpha = 1.95,
|
||||
),
|
||||
)
|
||||
PropertyT.solve(m, cosmo_optimizer(
|
||||
eps=1e-7,
|
||||
max_iters=10_000,
|
||||
accel=50,
|
||||
alpha=1.95,
|
||||
))
|
||||
|
||||
@test JuMP.termination_status(m) in
|
||||
(JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
|
||||
@test JuMP.termination_status(m) in (JuMP.ALMOST_OPTIMAL, JuMP.OPTIMAL)
|
||||
@test JuMP.objective_value(m) ≈ λ_cert atol = 1e-2
|
||||
end
|
||||
|
||||
@ -60,22 +55,20 @@
|
||||
Σ,
|
||||
act,
|
||||
basis(RG),
|
||||
StarAlgebras.Basis{UInt16}(
|
||||
@view basis(RG)[1:sizes[halfradius]]
|
||||
),
|
||||
StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[halfradius]]),
|
||||
)
|
||||
|
||||
status, certified, λ_cert = check_positivity(
|
||||
elt,
|
||||
unit,
|
||||
wd;
|
||||
upper_bound = ub,
|
||||
halfradius = 2,
|
||||
optimizer = cosmo_optimizer(;
|
||||
eps = 1e-7,
|
||||
max_iters = 10_000,
|
||||
accel = 50,
|
||||
alpha = 1.9,
|
||||
wd,
|
||||
upper_bound=ub,
|
||||
halfradius=2,
|
||||
optimizer=cosmo_optimizer(
|
||||
eps=1e-7,
|
||||
max_iters=10_000,
|
||||
accel=50,
|
||||
alpha=1.9,
|
||||
),
|
||||
)
|
||||
|
||||
|
||||
@ -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
|
||||
@ -24,8 +24,6 @@ if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI")
|
||||
include("1712.07167.jl")
|
||||
include("1812.03456.jl")
|
||||
|
||||
include("roots.jl")
|
||||
include("graded_adj.jl")
|
||||
include("Chevalley.jl")
|
||||
end
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user