From 4e43811ea3385a3d7abaab8ca4796bd03fe50dec Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Fri, 17 Mar 2023 15:53:36 +0100
Subject: [PATCH 01/27] fix a nasty bug with negatives in ConstraintMatrix

---
 src/constraint_matrix.jl     | 13 +++++++++----
 test/constratint_matrices.jl | 14 +++++++++++---
 2 files changed, 20 insertions(+), 7 deletions(-)

diff --git a/src/constraint_matrix.jl b/src/constraint_matrix.jl
index 8fbd36a..5734837 100644
--- a/src/constraint_matrix.jl
+++ b/src/constraint_matrix.jl
@@ -50,10 +50,15 @@ ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} =
 
 Base.size(cm::ConstraintMatrix) = cm.size
 
-__get_positive(cm::ConstraintMatrix, idx::Integer) =
-    convert(eltype(cm), cm.val * length(searchsorted(cm.pos, idx)))
-__get_negative(cm::ConstraintMatrix, idx::Integer) =
-    convert(eltype(cm), cm.val * length(searchsorted(cm.neg, idx)))
+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
 
 Base.@propagate_inbounds function Base.getindex(
     cm::ConstraintMatrix,
diff --git a/test/constratint_matrices.jl b/test/constratint_matrices.jl
index afdee61..1d04876 100644
--- a/test/constratint_matrices.jl
+++ b/test/constratint_matrices.jl
@@ -1,12 +1,17 @@
 @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
-        0.0 π
+        2π 0
+        0 π
     ]
 
     @test collect(PropertyT.nzpairs(cm)) == [
@@ -18,4 +23,7 @@
         1 => -3.141592653589793
         1 => -3.141592653589793
     ]
+
+    @test PropertyT.ConstraintMatrix{Float64}([-9:-1; 1:9], 3, 3, 1.0) ==
+          zeros(3, 3)
 end

From 1fb324b49ae492b42c828acd479c552ab0f06d23 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 20:26:42 +0100
Subject: [PATCH 02/27] update constraints to StarAlgebras-0.2

---
 Project.toml             |  4 +--
 src/PropertyT.jl         |  4 +--
 src/certify.jl           | 29 ++++++++++-----------
 src/constraint_matrix.jl | 47 ++++++++++++++++++++++++++++++++++
 src/sos_sdps.jl          | 55 +++++-----------------------------------
 5 files changed, 72 insertions(+), 67 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8255f05..75f97da 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,9 +19,9 @@ Groups = "0.7"
 IntervalArithmetic = "0.20"
 JuMP = "1.3"
 ProgressMeter = "1.7"
-SCS = "1.1.0"
+SCS = "1.1"
 StaticArrays = "1"
-SymbolicWedderburn = "0.3.2"
+SymbolicWedderburn = "0.3.4"
 julia = "1.6"
 
 [extras]
diff --git a/src/PropertyT.jl b/src/PropertyT.jl
index 6b97d5d..16430ee 100644
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@@ -26,13 +26,13 @@ include("gradings.jl")
 
 include("actions/actions.jl")
 
-function group_algebra(G::Groups.Group, S=gens(G); halfradius::Integer, twisted::Bool)
+function group_algebra(G::Groups.Group, S = gens(G); halfradius::Integer)
     S = union!(S, inv.(S))
     @info "generating wl-metric ball of radius $(2halfradius)"
     @time E, sizes = Groups.wlmetric_ball(S, radius=2halfradius)
     @info "sizes = $(sizes)"
     @info "computing the *-algebra structure for G"
-    @time RG = StarAlgebras.StarAlgebra{twisted}(
+    @time RG = StarAlgebras.StarAlgebra(
         G,
         StarAlgebras.Basis{UInt32}(E),
         (sizes[halfradius], sizes[halfradius]),
diff --git a/src/certify.jl b/src/certify.jl
index cdebcc9..2079791 100644
--- a/src/certify.jl
+++ b/src/certify.jl
@@ -8,26 +8,25 @@ end
 function __sos_via_sqr!(
     res::StarAlgebras.AlgebraElement,
     P::AbstractMatrix;
-    augmented::Bool
+    augmented::Bool,
+    id = (b = basis(parent(res)); b[one(first(b))]),
 )
-    StarAlgebras.zero!(res)
     A = parent(res)
-    b = basis(A)
-    @assert size(A.mstructure) == size(P)
-    e = b[one(b[1])]
+    mstr = A.mstructure
+    @assert size(mstr) == size(P)
 
-    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)
+    StarAlgebras.zero!(res)
+    for j in axes(mstr, 2)
+        for i in axes(mstr, 1)
             p = P[i, j]
-            xy = b[A.mstructure[i, j]]
-            # either result += P[x,y]*(x*y)
-            res[xy] += p
+            x_star_y = mstr[-i, j]
+            res[x_star_y] += p
+            # either result += P[x,y]*(x'*y)
             if augmented
-                # or result += P[x,y]*(1-x)*(1-y) == P[x,y]*(2-x-y+xy)
-                y = b[j]
-                res[e] += p
-                res[x] -= p
+                # 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
                 res[y] -= p
             end
         end
diff --git a/src/constraint_matrix.jl b/src/constraint_matrix.jl
index 5734837..73f3238 100644
--- a/src/constraint_matrix.jl
+++ b/src/constraint_matrix.jl
@@ -132,3 +132,50 @@ 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
diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index 2f6faea..04a8916 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -11,10 +11,8 @@ function sos_problem_dual(
 )
     @assert parent(elt) == parent(order_unit)
     algebra = parent(elt)
-    mstructure = if StarAlgebras._istwisted(algebra.mstructure)
-        algebra.mstructure
-    else
-        StarAlgebras.MTable{true}(basis(algebra), table_size=size(algebra.mstructure))
+    moment_matrix = let m = algebra.mstructure
+        [m[-i, j] for i in axes(m, 1) for j in axes(m, 2)]
     end
 
     # 1 variable for every primal constraint
@@ -24,7 +22,7 @@ function sos_problem_dual(
     model = Model()
     @variable model y[1:length(basis(algebra))]
     @constraint model λ_dual dot(order_unit, y) == 1
-    @constraint(model, psd, y[mstructure] in PSDCone())
+    @constraint(model, psd, y[moment_matrix] in PSDCone())
 
     if !isinf(lower_bound)
         throw("Not Implemented yet")
@@ -37,45 +35,6 @@ function sos_problem_dual(
     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`.
@@ -111,7 +70,7 @@ function sos_problem_primal(
         @warn "Setting `upper_bound` together with zero `order_unit` has no effect"
     end
 
-    A = constraints(parent(elt), augmented=augmented, twisted=true)
+    A = constraints(parent(elt); augmented = augmented)
 
     if !iszero(order_unit)
         λ = JuMP.@variable(model, λ)
@@ -135,9 +94,9 @@ end
 function invariant_constraint!(
     result::AbstractMatrix,
     basis::StarAlgebras.AbstractBasis,
-    cnstrs::AbstractDict{K,CM},
+    cnstrs::AbstractDict{K,<:ConstraintMatrix},
     invariant_vec::SparseVector,
-) where {K,CM<:ConstraintMatrix}
+) where {K}
     result .= zero(eltype(result))
     for i in SparseArrays.nonzeroinds(invariant_vec)
         g = basis[i]
@@ -234,7 +193,7 @@ function sos_problem_primal(
     U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
 
     # defining constraints based on the multiplicative structure
-    cnstrs = constraints(parent(elt), augmented=augmented, twisted=true)
+    cnstrs = constraints(parent(elt); augmented = augmented)
 
     prog = ProgressMeter.Progress(
         length(invariant_vectors(wedderburn)),

From b5fa1ac0ef1c5dea5bd66ede94ed18dafca7d18e Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 20:36:02 +0100
Subject: [PATCH 03/27] use Cartan matrix to classify root-subsystems

---
 src/roots.jl | 85 +++++++++++++++++++++++++++++++++-------------------
 1 file changed, 54 insertions(+), 31 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index 3c5d6f3..f8f809e 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -71,39 +71,55 @@ end
 𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
 𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
 
+reflection(α::Root, β::Root) = β - Int(2dot(α, β) / dot(α, α)) * α
+function cartan(α, β)
+    return [
+        length(reflection(a, b) - b) / length(a) for a in (α, β), b in (α, β)
+    ]
+end
+
 """
     classify_root_system(α, β)
 Return the symbol of smallest system generated by roots `α` and `β`.
 
-The classification is based only on roots length and
-proportionality/orthogonality.
+The classification is based only on roots length,
+proportionality/orthogonality and Cartan matrix.
 """
-function classify_root_system(α::AbstractRoot, β::AbstractRoot)
-    lα, lβ = length(α), length(β)
+function classify_root_system(
+    α::AbstractRoot,
+    β::AbstractRoot,
+    long::Tuple{Bool,Bool},
+)
     if isproportional(α, β)
-        if lα ≈ lβ ≈ √2
-            return :A₁
-        elseif lα ≈ lβ ≈ 2.0
+        if all(long)
             return :C₁
+        elseif all(.!long) # both short
+            return :A₁
         else
+            @error "Proportional roots of different length"
             error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
         end
     elseif isorthogonal(α, β)
-        if lα ≈ lβ ≈ √2
-            return Symbol("A₁×A₁")
-        elseif lα ≈ lβ ≈ 2.0
+        if all(long)
             return Symbol("C₁×C₁")
-        elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
+        elseif all(.!long) # both short
+            return Symbol("A₁×A₁")
+        elseif any(long)
             return Symbol("A₁×C₁")
-        else
-            error("Unknown root system ⟨α, β⟩:\n α = $α\n β = $β")
         end
     else # ⟨α, β⟩ is 2-dimensional, but they're not orthogonal
-        if lα ≈ lβ ≈ √2
+        a, b, c, d = abs.(cartan(α, β))
+        @assert a == d == 2
+        b, c = b < c ? (b, c) : (c, b)
+        if b == c == 1
             return :A₂
-        elseif (lα ≈ 2.0 && lβ ≈ √2) || (lα ≈ √2 && lβ ≈ 2)
+        elseif b == 1 && c == 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
@@ -130,12 +146,17 @@ function Base.in(r::R, plane::Plane{R}) where {R}
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
+function _islong(α::Root, Ω)
+    lα = length(α)
+    return any(r -> lα - length(r) > eps(lα), Ω)
+end
+
 function classify_sub_root_system(
     Ω::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(Ω, β)
 
@@ -146,28 +167,30 @@ function classify_sub_root_system(
     l = length(subsystem)
     if l == 1
         x = first(subsystem)
-        return classify_root_system(x, x)
+        long = _islong(x, Ω)
+        return classify_root_system(x, -x, (long, long))
     elseif l == 2
-        return classify_root_system(subsystem...)
+        x, y = subsystem
+        return classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω)))
     elseif l == 3
-        a = classify_root_system(subsystem[1], subsystem[2])
-        b = classify_root_system(subsystem[2], subsystem[3])
-        c = classify_root_system(subsystem[1], subsystem[3])
+        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))
 
-        if a == b == c # it's only A₂
+        if :A₂ == a == b == c # it's only A₂
             return a
         end
 
-        C = (:C₂, Symbol("C₁×C₁"))
-        if (a ∈ C && b ∈ C && c ∈ C) && (:C₂ ∈ (a, b, c))
-            return :C₂
-        end
+        throw("Unknown subroot system! $((x,y,z))")
     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
+        subtypes = [
+            classify_root_system(x, y, (_islong(x, Ω), _islong(y, Ω))) for
+            x in subsystem for y in subsystem if x ≠ y
+        ]
+        if :C₂ in subtypes
+            return :C₂
         end
     end
     @error "Unknown root subsystem generated by" α β

From a1de0ecc85cf8e0a12079eed4949f6b039c600b2 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 20:37:46 +0100
Subject: [PATCH 04/27] improve __fast_recursive_dot a bit

---
 src/sos_sdps.jl | 25 ++++++++++++++-----------
 1 file changed, 14 insertions(+), 11 deletions(-)

diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index 04a8916..e256bd0 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -122,11 +122,12 @@ sos_problem_primal(
 function __fast_recursive_dot!(
     res::JuMP.AffExpr,
     Ps::AbstractArray{<:AbstractMatrix{<:JuMP.VariableRef}},
-    Ms::AbstractArray{<:AbstractSparseMatrix};
+    Ms::AbstractArray{<:AbstractSparseMatrix},
+    weights,
 )
     @assert length(Ps) == length(Ms)
 
-    for (A, P) in zip(Ms, Ps)
+    for (w, A, P) in zip(weights, Ms, Ps)
         iszero(Ms) && continue
         rows = rowvals(A)
         vals = nonzeros(A)
@@ -134,13 +135,21 @@ function __fast_recursive_dot!(
             for i in nzrange(A, cidx)
                 ridx = rows[i]
                 val = vals[i]
-                JuMP.add_to_expression!(res, P[ridx, cidx], val)
+                JuMP.add_to_expression!(res, P[ridx, cidx], w * 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")]
 
@@ -217,15 +226,9 @@ function sos_problem_primal(
         # @info [nnz(m) / length(m) for m in Ms]
 
         if feasibility_problem
-            JuMP.@constraint(
-                model,
-                x == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
-            )
+            JuMP.@constraint(model, x == _dot(P, Ms))
         else
-            JuMP.@constraint(
-                model,
-                x - λ * u == __fast_recursive_dot!(JuMP.AffExpr(), P, Ms)
-            )
+            JuMP.@constraint(model, x - λ * u == _dot(P, Ms))
         end
     end
     ProgressMeter.finish!(prog)

From a5f5a4ea35539240e73b5cba4b75c80c6c4201f5 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 20:33:28 +0100
Subject: [PATCH 05/27] lots of re-formatting

---
 .JuliaFormatter.toml     |  9 ++++++++
 src/PropertyT.jl         |  6 ++---
 src/actions/actions.jl   |  3 +--
 src/certify.jl           | 38 +++++++++++++++++++------------
 src/constraint_matrix.jl | 19 ++++++++++++----
 src/gradings.jl          | 29 ++++++++++++++++--------
 src/reconstruct.jl       | 32 +++++++++++++-------------
 src/roots.jl             | 18 +++++++--------
 src/solve.jl             | 23 ++++++++++---------
 src/sos_sdps.jl          | 49 ++++++++++++++++++++++------------------
 test/optimizers.jl       | 34 +++++++++++++---------------
 11 files changed, 150 insertions(+), 110 deletions(-)
 create mode 100644 .JuliaFormatter.toml

diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..f196f10
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1,9 @@
+# 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
diff --git a/src/PropertyT.jl b/src/PropertyT.jl
index 16430ee..d4a664a 100644
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@@ -29,14 +29,14 @@ include("actions/actions.jl")
 function group_algebra(G::Groups.Group, S = gens(G); halfradius::Integer)
     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(
         G,
         StarAlgebras.Basis{UInt32}(E),
-        (sizes[halfradius], sizes[halfradius]),
-        precompute=false,
+        (sizes[halfradius], sizes[halfradius]);
+        precompute = false,
     )
     return RG, S, sizes
 end
diff --git a/src/actions/actions.jl b/src/actions/actions.jl
index 9cc6872..8205e8b 100644
--- a/src/actions/actions.jl
+++ b/src/actions/actions.jl
@@ -1,5 +1,4 @@
 import SymbolicWedderburn.action
-StarAlgebras.star(g::Groups.GroupElement) = inv(g)
 
 include("alphabet_permutation.jl")
 
@@ -10,7 +9,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(α))
diff --git a/src/certify.jl b/src/certify.jl
index 2079791..43648aa 100644
--- a/src/certify.jl
+++ b/src/certify.jl
@@ -35,7 +35,11 @@ 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²)
@@ -43,10 +47,14 @@ function __sos_via_cnstr!(res::StarAlgebras.AlgebraElement, Q²::AbstractMatrix,
     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
 
@@ -80,13 +88,12 @@ 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(
@@ -95,17 +102,18 @@ 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
@@ -122,16 +130,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 && inf(λ_certified) > 0.0, inf(λ_certified)
+    return check && inf(λ_certified) > 0.0, λ_certified
 end
diff --git a/src/constraint_matrix.jl b/src/constraint_matrix.jl
index 73f3238..1e95aa1 100644
--- a/src/constraint_matrix.jl
+++ b/src/constraint_matrix.jl
@@ -28,7 +28,12 @@ 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
 
@@ -45,8 +50,14 @@ struct ConstraintMatrix{T,I} <: AbstractMatrix{T}
     end
 end
 
-ConstraintMatrix(nzeros::AbstractArray{<:Integer}, n, m, val::T) where {T} =
-    ConstraintMatrix{T}(nzeros, n, m, val)
+function ConstraintMatrix(
+    nzeros::AbstractArray{<:Integer},
+    n,
+    m,
+    val::T,
+) where {T}
+    return ConstraintMatrix{T}(nzeros, n, m, val)
+end
 
 Base.size(cm::ConstraintMatrix) = cm.size
 
@@ -80,7 +91,7 @@ Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
 Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
 
 # TODO: iterate over (idx=>val) pairs combining vals
-function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int}=(1, 1))
+function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int} = (1, 1))
     k = iterate(itr.m.pos, state[1])
     isnothing(k) && return iterate(itr, state[2])
     idx, st = k
diff --git a/src/gradings.jl b/src/gradings.jl
index 8f34efa..b51de1a 100644
--- a/src/gradings.jl
+++ b/src/gradings.jl
@@ -1,7 +1,8 @@
 ## something about roots
 
-Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
-    Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
+function Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N}
+    return Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
+end
 
 function Roots.Root(s::MatrixGroups.ElementarySymplectic{N}) where {N}
     if s.symbol === :A
@@ -51,20 +52,28 @@ 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))),
+            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))),
     )
 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
diff --git a/src/reconstruct.jl b/src/reconstruct.jl
index 75451ee..3971d0a 100644
--- a/src/reconstruct.jl
+++ b/src/reconstruct.jl
@@ -7,35 +7,31 @@ end
 
 function reconstruct(
     Ms::AbstractVector{<:AbstractMatrix},
-    wbdec::WedderburnDecomposition;
-    atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
+    wbdec::WedderburnDecomposition,
 )
     n = __outer_dim(wbdec)
     res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
         res = similar(M, n, n)
-        reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom, atol=atol)
+        res = _reconstruct!(res, M, ds, __group_of(wbdec), wbdec.hom)
+        return res
     end
     return res
 end
 
-function reconstruct!(
+function _reconstruct!(
     res::AbstractMatrix,
     M::AbstractMatrix,
     ds::SymbolicWedderburn.DirectSummand,
     G,
-    hom;
-    atol=eps(real(eltype(ds))) * 10max(size(res)...)
+    hom,
 )
-    res .= zero(eltype(res))
     U = SymbolicWedderburn.image_basis(ds)
     d = SymbolicWedderburn.degree(ds)
-    tmp = (U' * M * U) .* d
+    Θπ = (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
+    res .= zero(eltype(res))
+    Θπ = average!(res, Θπ, G, hom)
+    return Θπ
 end
 
 function __droptol!(M::AbstractMatrix, tol)
@@ -52,14 +48,18 @@ function average!(
     res::AbstractMatrix,
     M::AbstractMatrix,
     G::Groups.Group,
-    hom::SymbolicWedderburn.InducedActionHomomorphism{<:SymbolicWedderburn.ByPermutations}
+    hom::SymbolicWedderburn.InducedActionHomomorphism{
+        <:SymbolicWedderburn.ByPermutations,
+    },
 )
     @assert size(M) == size(res)
     for g in G
-        gext = SymbolicWedderburn.induce(hom, g)
+        p = SymbolicWedderburn.induce(hom, g)
         Threads.@threads for c in axes(res, 2)
+            # note: p is a permutation,
+            # so accesses below are guaranteed to be disjoint
             for r in axes(res, 1)
-                res[r, c] += M[r^gext, c^gext]
+                res[r^p, c^p] += M[r, c]
             end
         end
     end
diff --git a/src/roots.jl b/src/roots.jl
index f8f809e..acbcdb2 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -30,7 +30,7 @@ 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.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))
@@ -38,13 +38,13 @@ cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
 function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
     N == M || return false
     val = abs(cos_angle(α, β))
-    return isapprox(val, one(val), atol=eps(one(val)))
+    return isapprox(val, one(val); atol = eps(one(val)))
 end
 
 function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
     N == M || return false
     val = cos_angle(α, β)
-    return isapprox(val, zero(val), atol=eps(one(val)))
+    return isapprox(val, zero(val); atol = eps(one(val)))
 end
 
 function _positive_direction(α::Root{N}) where {N}
@@ -53,19 +53,18 @@ function _positive_direction(α::Root{N}) where {N}
 end
 
 function positive(roots::AbstractVector{<:Root{N}}) where {N}
-    # return those roots for which dot(α, Root([½, ¼, …])) > 0.0
     pd = _positive_direction(first(roots))
     return filter(α -> dot(α, pd) > 0.0, roots)
 end
 
 function Base.show(io::IO, r::Root)
-    print(io, "Root$(r.coord)")
+    return print(io, "Root$(r.coord)")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
     lngth² = sum(x -> x^2, r.coord)
     l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
-    print(io, "Root in ℝ^$N of length $l\n", r.coord)
+    return print(io, "Root in ℝ^$N of length $l\n", r.coord)
 end
 
 𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
@@ -139,10 +138,11 @@ struct Plane{R<:Root}
     vectors::Vector{R}
 end
 
-Plane(α::R, β::R) where {R<:Root} =
-    Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
+function Plane(α::Root, β::Root)
+    return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
+end
 
-function Base.in(r::R, plane::Plane{R}) where {R}
+function Base.in(r::Root, plane::Plane)
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
diff --git a/src/solve.jl b/src/solve.jl
index e7df555..c872617 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -7,12 +7,15 @@ 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)
 
-    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])
+        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
     end
     return model
 end
@@ -28,11 +31,10 @@ 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)
 
@@ -45,8 +47,7 @@ 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)
@@ -55,7 +56,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()
-            status, warmstart
+            return status, warmstart
         end
     end
 
diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index e256bd0..eb86dcf 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -6,8 +6,8 @@ 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)
@@ -55,9 +55,10 @@ 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)
 
@@ -113,11 +114,13 @@ function isorth_projection(ds::SymbolicWedderburn.DirectSummand)
     return isapprox(U * U', I)
 end
 
-sos_problem_primal(
+function sos_problem_primal(
     elt::StarAlgebras.AlgebraElement,
     wedderburn::WedderburnDecomposition;
-    kwargs...
-) = sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
+    kwargs...,
+)
+    return sos_problem_primal(elt, zero(elt), wedderburn; kwargs...)
+end
 
 function __fast_recursive_dot!(
     res::JuMP.AffExpr,
@@ -157,16 +160,18 @@ 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
 
@@ -189,7 +194,7 @@ function sos_problem_primal(
         dim = size(ds, 1)
         P = JuMP.@variable(model, [1:dim, 1:dim], Symmetric)
         JuMP.@constraint(model, P in PSDCone())
-        P
+        return P
     end
 
     begin # preallocating
@@ -205,16 +210,16 @@ function sos_problem_primal(
     cnstrs = constraints(parent(elt); augmented = augmented)
 
     prog = ProgressMeter.Progress(
-        length(invariant_vectors(wedderburn)),
-        dt=1,
-        desc="Adding constraints... ",
-        enabled=show_progress,
-        barlen=60,
-        showspeed=true
+        length(invariant_vectors(wedderburn));
+        dt = 1,
+        desc = "Adding constraints: ",
+        enabled = show_progress,
+        barlen = 60,
+        showspeed = true,
     )
 
     for (i, iv) in enumerate(invariant_vectors(wedderburn))
-        ProgressMeter.next!(prog, showvalues=__show_itrs(i, prog.n))
+        ProgressMeter.next!(prog; showvalues = __show_itrs(i, prog.n))
 
         x = dot(X, iv)
         u = dot(U, iv)
diff --git a/test/optimizers.jl b/test/optimizers.jl
index f415de7..5832ca9 100644
--- a/test/optimizers.jl
+++ b/test/optimizers.jl
@@ -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,20 +28,18 @@ 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,
@@ -49,6 +47,6 @@ function cosmo_optimizer(;
         "max_iter" => max_iters,
         "verbose" => verbose,
         "verbose_timing" => verbose_timing,
-        "check_termination" => 200,
+        "check_termination" => 250,
     )
 end

From bacd170504f4be48635ce71390a77a6f03991b75 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 20:39:43 +0100
Subject: [PATCH 06/27] preserve trace when diagonalizing M_orb

---
 src/sos_sdps.jl | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index eb86dcf..4eecd09 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -225,8 +225,14 @@ function sos_problem_primal(
         u = dot(U, 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)...))
+
+        Ms = SymbolicWedderburn.diagonalize!(
+            Ms,
+            M_orb,
+            wedderburn;
+            trace_preserving = true,
+        )
+        # SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
 
         # @info [nnz(m) / length(m) for m in Ms]
 

From fb3b51fd6e81f563a01d17f8632dcffd5a1977b8 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 23:26:44 +0100
Subject: [PATCH 07/27] fix sos_problem_dual

---
 src/sos_sdps.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index 4eecd09..a91b240 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -12,7 +12,7 @@ function sos_problem_dual(
     @assert parent(elt) == parent(order_unit)
     algebra = parent(elt)
     moment_matrix = let m = algebra.mstructure
-        [m[-i, j] for i in axes(m, 1) for j in axes(m, 2)]
+        [m[-i, j] for i in axes(m, 1), j in axes(m, 2)]
     end
 
     # 1 variable for every primal constraint
@@ -20,16 +20,16 @@ function sos_problem_dual(
     # Symmetrized:
     # 1 dual variable for every orbit of G acting on basis
     model = Model()
-    @variable model y[1:length(basis(algebra))]
-    @constraint model λ_dual dot(order_unit, y) == 1
-    @constraint(model, psd, y[moment_matrix] in PSDCone())
+    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())
 
     if !isinf(lower_bound)
         throw("Not Implemented yet")
-        @variable model λ_ub_dual
-        @objective model Min dot(elt, y) + lower_bound * λ_ub_dual
+        JuMP.@variable(model, λ_ub_dual)
+        JuMP.@objective(model, Min, dot(elt, y) + lower_bound * λ_ub_dual)
     else
-        @objective model Min dot(elt, y)
+        JuMP.@objective(model, Min, dot(elt, y))
     end
 
     return model

From 914b06807041ee42d9c2e6fd5618d5b44855b0bc Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Sun, 19 Mar 2023 23:28:36 +0100
Subject: [PATCH 08/27] update and reformat tests

---
 test/1703.09680.jl       | 155 ++++++++++++++++++------------------
 test/1712.07167.jl       |  88 ++++++++++----------
 test/1812.03456.jl       |  67 +++++++++-------
 test/Chevalley.jl        | 168 +++++++++++++++++++++++++++++++++++++++
 test/actions.jl          |  56 ++++++-------
 test/check_positivity.jl |  35 ++++----
 test/graded_adj.jl       |  26 +++---
 test/quick_tests.jl      |  63 ++++++++-------
 test/runtests.jl         |   1 +
 9 files changed, 424 insertions(+), 235 deletions(-)
 create mode 100644 test/Chevalley.jl

diff --git a/test/1703.09680.jl b/test/1703.09680.jl
index 6ed0f79..e801ab9 100644
--- a/test/1703.09680.jl
+++ b/test/1703.09680.jl
@@ -1,9 +1,8 @@
 @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, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -15,15 +14,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
@@ -33,8 +32,10 @@
 
     @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, twisted=true)
+        G = MatrixGroups.SpecialLinearGroup{N}(
+            SymbolicWedderburn.Characters.FiniteFields.GF{5},
+        )
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -46,15 +47,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
@@ -62,12 +63,15 @@
         @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
@@ -76,7 +80,7 @@
     @testset "SAut(F₂)" begin
         N = 2
         G = SpecialAutomorphismGroup(FreeGroup(N))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -88,53 +92,46 @@
 
         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, twisted=true)
+        RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
 
         Δ = RSL(length(S)) - sum(RSL(s) for s in S)
 
@@ -145,18 +142,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,
                 ),
             )
 
@@ -170,13 +167,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
@@ -186,18 +183,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,
                 ),
             )
 
@@ -210,13 +207,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
diff --git a/test/1712.07167.jl b/test/1712.07167.jl
index be479d9..abc19b8 100644
--- a/test/1712.07167.jl
+++ b/test/1712.07167.jl
@@ -1,9 +1,9 @@
 @testset "1712.07167 Examples" begin
-
     @testset "SAut(F₃)" begin
         N = 3
         G = SpecialAutomorphismGroup(FreeGroup(N))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        @info "running tests for" G
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1)))
         Σ = PropertyT.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
@@ -15,6 +15,7 @@
             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)
@@ -27,14 +28,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,
             ),
         )
 
@@ -47,7 +48,8 @@
         n = 3
 
         SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
+        @info "running tests for" SL
+        RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
 
         Δ = RSL(length(S)) - sum(RSL(s) for s in S)
 
@@ -62,6 +64,7 @@
                 basis(RSL),
                 StarAlgebras.Basis{UInt16}(@view basis(RSL)[1:sizes[2]]),
             )
+            @info wd
 
             elt = Δ^2
             unit = Δ
@@ -71,8 +74,8 @@
                 elt,
                 unit,
                 wd,
-                upper_bound=ub,
-                augmented=false,
+                upper_bound = ub,
+                augmented = false,
             )
 
             wdfl = SymbolicWedderburn.WedderburnDecomposition(
@@ -86,18 +89,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,
                 ),
             )
 
@@ -105,34 +108,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
-                    real.(sqrt(JuMP.value.(P)))
+                    return 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
@@ -155,22 +158,23 @@
                 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,
                 ),
             )
 
@@ -178,7 +182,7 @@
 
             Q = @time let varP = varP
                 Qs = map(varP) do P
-                    real.(sqrt(JuMP.value.(P)))
+                    return real.(sqrt(JuMP.value.(P)))
                 end
                 PropertyT.reconstruct(Qs, wdfl)
             end
@@ -187,8 +191,8 @@
                 elt,
                 unit,
                 JuMP.objective_value(opt_problem),
-                Q,
-                halfradius=2,
+                Q;
+                halfradius = 2,
                 # augmented = true # since both elt and unit are augmented
             )
 
diff --git a/test/1812.03456.jl b/test/1812.03456.jl
index b0d9482..46802b9 100644
--- a/test/1812.03456.jl
+++ b/test/1812.03456.jl
@@ -21,8 +21,9 @@ 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, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -39,6 +40,7 @@ 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)
@@ -82,7 +84,7 @@ using SparseArrays
     @testset "SAut(F₃)" begin
         n = 3
         G = SpecialAutomorphismGroup(FreeGroup(n))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
         sq, adj, op = PropertyT.SqAdjOp(RG, n)
 
         @test sq(one(G)) == 216
@@ -98,7 +100,8 @@ end
         n = 3
 
         G = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        @info "running tests for" G
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = RG(length(S)) - sum(RG(s) for s in S)
 
@@ -113,6 +116,7 @@ end
             basis(RG),
             StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
         )
+        @info wd
 
         sq, adj, op = PropertyT.SqAdjOp(RG, n)
 
@@ -123,10 +127,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
@@ -140,10 +144,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
@@ -152,10 +156,11 @@ 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
 
@@ -168,10 +173,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
@@ -183,7 +188,8 @@ end
         n = 4
 
         G = MatrixGroups.SpecialLinearGroup{n}(Int8)
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=2, twisted=true)
+        @info "running tests for" G
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 2)
 
         Δ = RG(length(S)) - sum(RG(s) for s in S)
 
@@ -198,6 +204,7 @@ end
             basis(RG),
             StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[2]]),
         )
+        @info wd
 
         sq, adj, op = PropertyT.SqAdjOp(RG, n)
 
@@ -208,10 +215,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
@@ -225,10 +232,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
@@ -242,10 +249,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
diff --git a/test/Chevalley.jl b/test/Chevalley.jl
new file mode 100644
index 0000000..4bb5862
--- /dev/null
+++ b/test/Chevalley.jl
@@ -0,0 +1,168 @@
+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.0, 1.0, 0.0, 0.0)
+                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, 0.0, 0.0)
+                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 / 2 .* (1.0, 1.0, 1.0, 1.0)
+                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(length(a), sqrt(2))
+        b = F4[6]
+        @test isapprox(length(b), sqrt(2))
+        c = a + b
+        @test isapprox(length(c), 2.0)
+        @test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
+
+        long = F4[findfirst(r -> length(r) == sqrt(2), F4)]
+        short = F4[findfirst(r -> 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.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+                    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 / 2 .* ntuple(i -> 1.0, 8)
+                    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 -> 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
diff --git a/test/actions.jl b/test/actions.jl
index af7d123..08e6a03 100644
--- a/test/actions.jl
+++ b/test/actions.jl
@@ -3,7 +3,7 @@ function test_action(basis, group, act)
     return @testset "action definition" begin
         @test all(basis) do b
             e = one(group)
-            action(act, e, b) == b
+            return 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)
-            x == y
+            return x == y
         end
 
         if act isa SymbolicWedderburn.ByPermutations
             @test all(basis) do b
-                action(act, g, b) ∈ basis && action(act, h, b) ∈ basis
+                return action(act, g, b) ∈ basis && action(act, h, b) ∈ basis
             end
         end
     end
@@ -44,21 +44,18 @@ end
 ## Testing
 
 @testset "Actions on SL(3,ℤ)" begin
-
     n = 3
 
     SL = MatrixGroups.SpecialLinearGroup{n}(Int8)
-    RSL, S, sizes = PropertyT.group_algebra(SL, halfradius=2, twisted=true)
+    RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = 2)
 
     @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},
@@ -66,9 +63,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -77,17 +74,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
@@ -97,7 +94,6 @@ end
         test_action(basis(RSL), Γ, ΓpA)
 
         @testset "mps is successful" begin
-
             charsΓ =
                 SymbolicWedderburn.Character{
                     Rational{Int},
@@ -105,9 +101,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -116,32 +112,30 @@ 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, twisted=true)
+    RSAutFn, S, sizes = PropertyT.group_algebra(SAutFn; halfradius = 1)
 
     @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},
@@ -149,9 +143,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -160,17 +154,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
@@ -180,7 +174,6 @@ end
         test_action(basis(RSAutFn), Γ, ΓpA)
 
         @testset "mps is successful" begin
-
             charsΓ =
                 SymbolicWedderburn.Character{
                     Rational{Int},
@@ -188,9 +181,9 @@ end
 
             RΓ = SymbolicWedderburn._group_algebra(Γ)
 
-            @time mps, simple =
+            @time mps, ranks =
                 SymbolicWedderburn.minimal_projection_system(charsΓ, RΓ)
-            @test all(simple)
+            @test all(isone, ranks)
         end
 
         @testset "Wedderburn decomposition" begin
@@ -199,11 +192,12 @@ 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
diff --git a/test/check_positivity.jl b/test/check_positivity.jl
index 974077b..d57e80f 100644
--- a/test/check_positivity.jl
+++ b/test/check_positivity.jl
@@ -1,6 +1,12 @@
-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])
@@ -9,16 +15,23 @@ function check_positivity(elt, unit; upper_bound=Inf, halfradius=2, optimizer)
         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)
 
@@ -29,13 +42,7 @@ function check_positivity(elt, unit, wd; upper_bound=Inf, halfradius=2, optimize
 
     λ = 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
-
diff --git a/test/graded_adj.jl b/test/graded_adj.jl
index bf5a804..447b58c 100644
--- a/test/graded_adj.jl
+++ b/test/graded_adj.jl
@@ -1,10 +1,9 @@
 @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, twisted=true)
+        RSL, S, sizes = PropertyT.group_algebra(SL; halfradius = halfradius)
 
         Δ = RSL(length(S)) - sum(RSL(s) for s in S)
 
@@ -22,10 +21,9 @@
         @test PropertyT.Adj(Δs, :A₂) == adj
         @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
 
-
         halfradius = 1
         G = SpecialAutomorphismGroup(FreeGroup(n))
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=halfradius, twisted=true)
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = halfradius)
 
         Δ = RG(length(S)) - sum(RG(s) for s in S)
 
@@ -44,7 +42,6 @@
         @test PropertyT.Adj(Δs, Symbol("A₁×A₁")) == op
     end
 
-
     @testset "Symplectic group" begin
         @testset "Sp2(ℤ)" begin
             genus = 2
@@ -52,7 +49,8 @@
 
             SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
 
-            RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
+            RSpN, S_sp, sizes_sp =
+                PropertyT.group_algebra(SpN; halfradius = halfradius)
 
             Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
                 Δ = RG(length(S)) - sum(RG(s) for s in S)
@@ -73,7 +71,8 @@
 
         SpN = MatrixGroups.SymplecticGroup{2genus}(Int8)
 
-        RSpN, S_sp, sizes_sp = PropertyT.group_algebra(SpN, halfradius=halfradius, twisted=true)
+        RSpN, S_sp, sizes_sp =
+            PropertyT.group_algebra(SpN; halfradius = halfradius)
 
         Δ, Δs = let RG = RSpN, S = S_sp, ψ = identity
             Δ = RG(length(S)) - sum(RG(s) for s in S)
@@ -86,9 +85,14 @@
         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
@@ -101,8 +105,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
-
diff --git a/test/quick_tests.jl b/test/quick_tests.jl
index 9a38b5c..9478286 100644
--- a/test/quick_tests.jl
+++ b/test/quick_tests.jl
@@ -1,11 +1,12 @@
 @testset "Quick tests" begin
-
     @testset "SL(2,F₇)" begin
         N = 2
         p = 7
         halfradius = 3
-        G = MatrixGroups.SpecialLinearGroup{N}(SymbolicWedderburn.Characters.FiniteFields.GF{p})
-        RG, S, sizes = PropertyT.group_algebra(G, halfradius=3, twisted=true)
+        G = MatrixGroups.SpecialLinearGroup{N}(
+            SymbolicWedderburn.Characters.FiniteFields.GF{p},
+        )
+        RG, S, sizes = PropertyT.group_algebra(G; halfradius = 3)
 
         Δ = let RG = RG, S = S
             RG(length(S)) - sum(RG(s) for s in S)
@@ -18,15 +19,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
@@ -34,14 +35,18 @@
             @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
 
@@ -55,20 +60,22 @@
                 Σ,
                 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,
                 ),
             )
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 7003c1e..5e6a9be 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -25,5 +25,6 @@ if haskey(ENV, "FULL_TEST") || haskey(ENV, "CI")
         include("1812.03456.jl")
 
         include("graded_adj.jl")
+        include("Chevalley.jl")
     end
 end

From 92d4ba0c207ab71ba63250f043ddca6e74c36474 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Mon, 20 Mar 2023 01:40:40 +0100
Subject: [PATCH 09/27] update script for SpN_Adj.jl

---
 scripts/SpN_Adj.jl | 30 +++++++++++++++---------------
 scripts/utils.jl   | 24 ++++++++++++++----------
 2 files changed, 29 insertions(+), 25 deletions(-)

diff --git a/scripts/SpN_Adj.jl b/scripts/SpN_Adj.jl
index c5922de..70ff7fe 100644
--- a/scripts/SpN_Adj.jl
+++ b/scripts/SpN_Adj.jl
@@ -24,8 +24,7 @@ const GENUS = 2N
 
 G = MatrixGroups.SymplecticGroup{GENUS}(Int8)
 
-RG, S, sizes =
-    @time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true)
+RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
 
 wd = let RG = RG, N = N
     G = StarAlgebras.object(RG)
@@ -42,6 +41,8 @@ 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)
@@ -58,23 +59,22 @@ 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₂-InfΔ",
-    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),
 )
-
diff --git a/scripts/utils.jl b/scripts/utils.jl
index df891b3..1357690 100644
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@@ -11,12 +11,13 @@ function get_solution(model)
     return solution
 end
 
-function get_solution(model, wd, varP; logdir)
+function get_solution(model, wd, varP)
     λ = JuMP.value(model[:λ])
 
     Qs = [real.(sqrt(JuMP.value.(P))) for P in varP]
     Q = PropertyT.reconstruct(Qs, wd)
     solution = Dict(:λ => λ, :Q => Q)
+
     return solution
 end
 
@@ -42,22 +43,23 @@ 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)
-
-            solution = get_solution(model, args...; logdir=logdir)
-            solution[:warm] = warm
+            status, warm =
+                @time PropertyT.solve(log_file, model, optimizer, warm)
 
+            solution = get_solution(model, args...)
             serialize(joinpath(logdir, "solution_$date.sjl"), solution)
             serialize(joinpath(logdir, "solution.sjl"), solution)
 
-            flag, λ_cert = open(log_file, append=true) do io
+            solution[:warm] = warm
+
+            flag, λ_cert = open(log_file; append = true) do io
                 with_logger(SimpleLogger(io)) do
-                    PropertyT.certify_solution(
+                    return PropertyT.certify_solution(
                         data.elt,
                         data.unit,
                         solution[:λ],
-                        solution[:Q],
-                        halfradius=data.halfradius,
+                        solution[:Q];
+                        halfradius = data.halfradius,
                     )
                 end
             end
@@ -69,7 +71,9 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
             @info "Certification done with λ = $certified_λ"
             return certified_λ
         else
-            rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda))
+            rel_change =
+                abs(certified_λ - old_lambda) /
+                (abs(certified_λ) + abs(old_lambda))
             @info "Certification failed with λ = $λ" certified_λ rel_change
         end
 

From a14c6d2669d826b960c415be50beafdc7ac8970c Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Mon, 20 Mar 2023 01:40:59 +0100
Subject: [PATCH 10/27] =?UTF-8?q?add=20G=E2=82=82=20script?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 scripts/G₂_Adj.jl  | 179 +++++++++++++++++++++++++
 scripts/G₂_gens.jl | 308 +++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 487 insertions(+)
 create mode 100644 scripts/G₂_Adj.jl
 create mode 100644 scripts/G₂_gens.jl

diff --git a/scripts/G₂_Adj.jl b/scripts/G₂_Adj.jl
new file mode 100644
index 0000000..efc6770
--- /dev/null
+++ b/scripts/G₂_Adj.jl
@@ -0,0 +1,179 @@
+using LinearAlgebra
+BLAS.set_num_threads(1)
+ENV["OMP_NUM_THREADS"] = 4
+
+using MKL_jll
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
+
+using Groups
+import Groups.MatrixGroups
+using PropertyT
+
+using SymbolicWedderburn
+using SymbolicWedderburn.StarAlgebras
+using PermutationGroups
+
+include(joinpath(@__DIR__, "G₂_gens.jl"))
+
+G, roots, Weyl = G₂_roots_weyl()
+
+const HALFRADIUS = 2
+const UPPER_BOUND = Inf
+
+RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
+
+Δ = RG(length(S)) - sum(RG(s) for s in S)
+
+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
+
+elt = Δ^2
+unit = Δ
+
+@time model, varP = PropertyT.sos_problem_primal(
+    elt,
+    unit,
+    wd;
+    upper_bound = UPPER_BOUND,
+    augmented = true,
+    show_progress = true,
+)
+warm = nothing
+
+begin
+    @time status, warm = PropertyT.solve(
+        model,
+        scs_optimizer(;
+            linear_solver = SCS.MKLDirectSolver,
+            eps = 1e-10,
+            max_iters = 20_000,
+            accel = 50,
+            alpha = 1.95,
+        ),
+        warm,
+    )
+
+    @info "reconstructing the solution"
+    Q = @time begin
+        wd = wd
+        Ps = [JuMP.value.(P) for P in varP]
+        if any(any(isnan, P) for P in Ps)
+            throw("solver was probably interrupted, no valid solution available")
+        end
+        Qs = real.(sqrt.(Ps))
+        PropertyT.reconstruct(Qs, wd)
+    end
+    P = Q' * Q
+
+    @info "certifying the solution"
+    @time certified, λ = PropertyT.certify_solution(
+        elt,
+        unit,
+        JuMP.objective_value(model),
+        Q;
+        halfradius = HALFRADIUS,
+        augmented = true,
+    )
+end
+
+### grading below
+
+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
+
+Δs = PropertyT.laplacians(
+    RG,
+    S,
+    x -> (gx = PropertyT.grading(x); Set([gx, -gx])),
+)
+
+elt = PropertyT.Adj(Δs)
+elt == Δ^2 - PropertyT.Sq(Δs)
+unit = Δ
+
+@time model, varP = PropertyT.sos_problem_primal(
+    elt,
+    unit,
+    wd;
+    upper_bound = UPPER_BOUND,
+    augmented = true,
+)
+
+warm = nothing
+
+begin
+    @time status, warm = PropertyT.solve(
+        model,
+        scs_optimizer(;
+            linear_solver = SCS.MKLDirectSolver,
+            eps = 1e-10,
+            max_iters = 50_000,
+            accel = 50,
+            alpha = 1.95,
+        ),
+        warm,
+    )
+
+    @info "reconstructing the solution"
+    Q = @time begin
+        wd = wd
+        Ps = [JuMP.value.(P) for P in varP]
+        if any(any(isnan, P) for P in Ps)
+            throw("solver was probably interrupted, no valid solution available")
+        end
+        Qs = real.(sqrt.(Ps))
+        PropertyT.reconstruct(Qs, wd)
+    end
+    P = Q' * Q
+
+    @info "certifying the solution"
+    @time certified, λ = PropertyT.certify_solution(
+        elt,
+        unit,
+        JuMP.objective_value(model),
+        Q;
+        halfradius = HALFRADIUS,
+        augmented = true,
+    )
+end
+
+# Δ² - 1 / 1 · Sq → -0.8818044647162608
+# Δ² - 2 / 3 · Sq → -0.1031738
+# Δ² - 1 / 2 · Sq → 0.228296213895906
+# Δ² - 1 / 3 · Sq → 0.520
+# Δ² - 0 / 1 · Sq → 0.9676851592000731
+# Sq → 0.333423
+
+# vals = [
+#     1.0 -0.8818
+#     2/3 -0.1032
+#     1/2  0.2282
+#     1/3  0.520
+#     0    0.9677
+# ]
diff --git a/scripts/G₂_gens.jl b/scripts/G₂_gens.jl
new file mode 100644
index 0000000..48677d5
--- /dev/null
+++ b/scripts/G₂_gens.jl
@@ -0,0 +1,308 @@
+#= 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

From 15286c0c4ae9213c686da3fa42f78114b360ae23 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Mon, 20 Mar 2023 01:42:41 +0100
Subject: [PATCH 11/27] add Manifest with a branch of SymbolicWedderburn

---
 Manifest.toml | 636 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 636 insertions(+)
 create mode 100644 Manifest.toml

diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000..458ec6b
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,636 @@
+# This file is machine-generated - editing it directly is not advised
+
+julia_version = "1.8.5"
+manifest_format = "2.0"
+project_hash = "5452e4eb8836d03e65c342a444826b894c3e14dc"
+
+[[deps.AbstractAlgebra]]
+deps = ["GroupsCore", "InteractiveUtils", "LinearAlgebra", "MacroTools", "Markdown", "Random", "RandomExtensions", "SparseArrays", "Test"]
+git-tree-sha1 = "29e65c331f97db9189ef00a4c7aed8127c2fd2d4"
+uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
+version = "0.27.10"
+
+[[deps.Accessors]]
+deps = ["Compat", "CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Requires", "StaticArrays", "Test"]
+git-tree-sha1 = "beabc31fa319f9de4d16372bff31b4801e43d32c"
+uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+version = "0.1.28"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.6.1"
+
+[[deps.ArgCheck]]
+git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
+uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
+version = "2.3.0"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.1"
+
+[[deps.Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[deps.BangBang]]
+deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
+git-tree-sha1 = "7fe6d92c4f281cf4ca6f2fba0ce7b299742da7ca"
+uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
+version = "0.3.37"
+
+[[deps.Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[deps.Baselet]]
+git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
+uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
+version = "0.1.1"
+
+[[deps.BenchmarkTools]]
+deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "1.3.2"
+
+[[deps.Bzip2_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
+uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
+version = "1.0.8+0"
+
+[[deps.CRlibm]]
+deps = ["CRlibm_jll"]
+git-tree-sha1 = "32abd86e3c2025db5172aa182b982debed519834"
+uuid = "96374032-68de-5a5b-8d9e-752f78720389"
+version = "1.0.1"
+
+[[deps.CRlibm_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "e329286945d0cfc04456972ea732551869af1cfc"
+uuid = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
+version = "1.0.1+0"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.15.7"
+
+[[deps.ChangesOfVariables]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
+git-tree-sha1 = "485193efd2176b88e6622a39a246f8c5b600e74e"
+uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+version = "0.1.6"
+
+[[deps.CodecBzip2]]
+deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
+git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
+uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
+version = "0.7.2"
+
+[[deps.CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.1"
+
+[[deps.CommonSubexpressions]]
+deps = ["MacroTools", "Test"]
+git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.3.0"
+
+[[deps.Compat]]
+deps = ["Dates", "LinearAlgebra", "UUIDs"]
+git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.6.1"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.0.1+0"
+
+[[deps.CompositionsBase]]
+git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
+uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
+version = "0.1.1"
+
+[[deps.ConstructionBase]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "89a9db8d28102b094992472d333674bd1a83ce2a"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "1.5.1"
+
+[[deps.Cyclotomics]]
+deps = ["LRUCache", "Memoize", "Primes", "SparseArrays", "Test"]
+git-tree-sha1 = "dc2e5fd64c188399434e83fa5c10c6fa4eff962a"
+uuid = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
+version = "0.3.2"
+
+[[deps.DataAPI]]
+git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.14.0"
+
+[[deps.DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.13"
+
+[[deps.DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
+
+[[deps.Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[deps.DefineSingletons]]
+git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
+uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
+version = "0.1.2"
+
+[[deps.DiffResults]]
+deps = ["StaticArraysCore"]
+git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "1.1.0"
+
+[[deps.DiffRules]]
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "a4ad7ef19d2cdc2eff57abbbe68032b1cd0bd8f8"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "1.13.0"
+
+[[deps.Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[deps.DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.9.3"
+
+[[deps.Downloads]]
+deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
+
+[[deps.ErrorfreeArithmetic]]
+git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
+uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
+version = "0.5.2"
+
+[[deps.ExternalDocstrings]]
+git-tree-sha1 = "1224740fc4d07c989949e1c1b508ebd49a65a5f6"
+uuid = "e189563c-0753-4f5e-ad5c-be4293c83fb4"
+version = "0.1.1"
+
+[[deps.FastRounding]]
+deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
+git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
+uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f"
+version = "0.3.1"
+
+[[deps.FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
+
+[[deps.Folds]]
+deps = ["Accessors", "BangBang", "Baselet", "DefineSingletons", "Distributed", "ExternalDocstrings", "InitialValues", "MicroCollections", "Referenceables", "Requires", "Test", "ThreadedScans", "Transducers"]
+git-tree-sha1 = "638109532de382a1f99b1aae1ca8b5d08515d85a"
+uuid = "41a02a25-b8f0-4f67-bc48-60067656b558"
+version = "0.2.8"
+
+[[deps.ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
+git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.35"
+
+[[deps.Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
+[[deps.Groups]]
+deps = ["Folds", "GroupsCore", "KnuthBendix", "LinearAlgebra", "Logging", "OrderedCollections", "PermutationGroups", "StaticArrays"]
+git-tree-sha1 = "5dbf642ee0048e6ad5f0bda11af17e40b8e8dd2f"
+uuid = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
+version = "0.7.5"
+
+[[deps.GroupsCore]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "9e1a5e9f3b81ad6a5c613d181664a0efc6fe6dd7"
+uuid = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
+version = "0.4.0"
+
+[[deps.InitialValues]]
+git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3"
+uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c"
+version = "0.3.1"
+
+[[deps.IntegerMathUtils]]
+git-tree-sha1 = "f366daebdfb079fd1fe4e3d560f99a0c892e15bc"
+uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
+version = "0.1.0"
+
+[[deps.InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[deps.IntervalArithmetic]]
+deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "RoundingEmulator", "SetRounding", "StaticArrays"]
+git-tree-sha1 = "c1c88395d09366dae431556bcb598ad08fa1392b"
+uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+version = "0.20.8"
+
+[[deps.InverseFunctions]]
+deps = ["Test"]
+git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.8"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.2.2"
+
+[[deps.IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
+
+[[deps.JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[deps.JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.3"
+
+[[deps.JuMP]]
+deps = ["LinearAlgebra", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays"]
+git-tree-sha1 = "611b9f12f02c587d860c813743e6cec6264e94d8"
+uuid = "4076af6c-e467-56ae-b986-b466b2749572"
+version = "1.9.0"
+
+[[deps.KnuthBendix]]
+deps = ["MacroTools", "ProgressMeter"]
+git-tree-sha1 = "a5ef62c9f2b4461b246a610b7402f5f2dbffbfaa"
+uuid = "c2604015-7b3d-4a30-8a26-9074551ec60a"
+version = "0.4.0"
+
+[[deps.LRUCache]]
+git-tree-sha1 = "d862633ef6097461037a00a13f709a62ae4bdfdd"
+uuid = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
+version = "1.4.0"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.3"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "7.84.0+0"
+
+[[deps.LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.10.2+0"
+
+[[deps.Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LogExpFunctions]]
+deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "0a1b7c2863e44523180fdb3146534e265a91870b"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.23"
+
+[[deps.Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[deps.MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.10"
+
+[[deps.Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[deps.MathOptInterface]]
+deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
+git-tree-sha1 = "f219b62e601c2f2e8adb7b6c48db8a9caf381c82"
+uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+version = "1.13.1"
+
+[[deps.MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+version = "2.28.0+0"
+
+[[deps.Memoize]]
+deps = ["MacroTools"]
+git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
+uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
+version = "0.4.4"
+
+[[deps.MicroCollections]]
+deps = ["BangBang", "InitialValues", "Setfield"]
+git-tree-sha1 = "629afd7d10dbc6935ec59b32daeb33bc4460a42e"
+uuid = "128add7d-3638-4c79-886c-908ea0c25c34"
+version = "0.1.4"
+
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[deps.MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+version = "2022.2.1"
+
+[[deps.MutableArithmetics]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "3295d296288ab1a0a2528feb424b854418acff57"
+uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+version = "1.2.3"
+
+[[deps.NaNMath]]
+deps = ["OpenLibm_jll"]
+git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "1.0.2"
+
+[[deps.NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+version = "1.2.0"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.20+0"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+version = "0.8.1+0"
+
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[deps.Parsers]]
+deps = ["Dates", "SnoopPrecompile"]
+git-tree-sha1 = "478ac6c952fddd4399e71d4779797c538d0ff2bf"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "2.5.8"
+
+[[deps.PermutationGroups]]
+deps = ["AbstractAlgebra", "GroupsCore", "Markdown", "Random"]
+git-tree-sha1 = "1bfba1a836e2c085270d3f5657803e0902e20564"
+uuid = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
+version = "0.3.3"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.8.0"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
+
+[[deps.Primes]]
+deps = ["IntegerMathUtils"]
+git-tree-sha1 = "311a2aa90a64076ea0fac2ad7492e914e6feeb81"
+uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+version = "0.5.3"
+
+[[deps.Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[deps.Profile]]
+deps = ["Printf"]
+uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+
+[[deps.ProgressMeter]]
+deps = ["Distributed", "Printf"]
+git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9"
+uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
+version = "1.7.2"
+
+[[deps.REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[deps.Random]]
+deps = ["SHA", "Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[deps.RandomExtensions]]
+deps = ["Random", "SparseArrays"]
+git-tree-sha1 = "062986376ce6d394b23d5d90f01d81426113a3c9"
+uuid = "fb686558-2515-59ef-acaa-46db3789a887"
+version = "0.4.3"
+
+[[deps.RecipesBase]]
+deps = ["SnoopPrecompile"]
+git-tree-sha1 = "261dddd3b862bd2c940cf6ca4d1c8fe593e457c8"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "1.3.3"
+
+[[deps.Referenceables]]
+deps = ["Adapt"]
+git-tree-sha1 = "e681d3bfa49cd46c3c161505caddf20f0e62aaa9"
+uuid = "42d2dcc6-99eb-4e98-b66c-637b7d73030e"
+version = "0.1.2"
+
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.0"
+
+[[deps.RoundingEmulator]]
+git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
+uuid = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
+version = "0.2.1"
+
+[[deps.SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+version = "0.7.0"
+
+[[deps.Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[deps.SetRounding]]
+git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
+uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
+version = "0.2.1"
+
+[[deps.Setfield]]
+deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
+git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
+uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+version = "1.1.1"
+
+[[deps.SnoopPrecompile]]
+deps = ["Preferences"]
+git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
+uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+version = "1.0.3"
+
+[[deps.Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[deps.SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.2.0"
+
+[[deps.SplittablesBase]]
+deps = ["Setfield", "Test"]
+git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5"
+uuid = "171d559e-b47b-412a-8079-5efa626c420e"
+version = "0.1.15"
+
+[[deps.StarAlgebras]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "265b89a5dfb38fe94ad48b997a253b2393fce6f1"
+uuid = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
+version = "0.2.0"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
+git-tree-sha1 = "6aa098ef1012364f2ede6b17bf358c7f1fbe90d4"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.5.17"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.0"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[deps.SymbolicWedderburn]]
+deps = ["Cyclotomics", "GroupsCore", "LinearAlgebra", "PermutationGroups", "Primes", "SparseArrays", "StarAlgebras"]
+git-tree-sha1 = "3a098aae261c4c57d6ee41dab62973065f0cb68b"
+repo-rev = "enh/towards_simple_projections"
+repo-url = "https://github.com/kalmarek/SymbolicWedderburn.jl.git"
+uuid = "858aa9a9-4c7c-4c62-b466-2421203962a2"
+version = "0.3.4"
+
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
+
+[[deps.TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[deps.Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
+git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.10.1"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.1"
+
+[[deps.Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[deps.ThreadedScans]]
+deps = ["ArgCheck"]
+git-tree-sha1 = "ca1ba3000289eacba571aaa4efcefb642e7a1de6"
+uuid = "24d252fe-5d94-4a69-83ea-56a14333d47a"
+version = "0.1.0"
+
+[[deps.TranscodingStreams]]
+deps = ["Random", "Test"]
+git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.9.11"
+
+[[deps.Transducers]]
+deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
+git-tree-sha1 = "c42fa452a60f022e9e087823b47e5a5f8adc53d5"
+uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
+version = "0.4.75"
+
+[[deps.UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[deps.Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[deps.Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+version = "1.2.12+3"
+
+[[deps.ZygoteRules]]
+deps = ["MacroTools"]
+git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
+uuid = "700de1a5-db45-46bc-99cf-38207098b444"
+version = "0.2.2"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.1.1+0"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.48.0+0"
+
+[[deps.p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
+version = "17.4.0+0"

From f9f852439f8a113e429336a0e73f0ae15824ed1b Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Mon, 20 Mar 2023 02:29:19 +0100
Subject: [PATCH 12/27] rewrite scripts for G2

---
 scripts/G₂_Adj.jl   | 150 ++++++++++--------------------------------
 scripts/G₂_has_T.jl |  97 +++++++++++++++++++++++++++
 scripts/utils.jl      |   6 +-
 3 files changed, 135 insertions(+), 118 deletions(-)
 create mode 100644 scripts/G₂_has_T.jl

diff --git a/scripts/G₂_Adj.jl b/scripts/G₂_Adj.jl
index efc6770..339d99c 100644
--- a/scripts/G₂_Adj.jl
+++ b/scripts/G₂_Adj.jl
@@ -1,29 +1,34 @@
 using LinearAlgebra
-BLAS.set_num_threads(1)
-ENV["OMP_NUM_THREADS"] = 4
+BLAS.set_num_threads(8)
 
-using MKL_jll
-include(joinpath(@__DIR__, "../test/optimizers.jl"))
+ENV["OMP_NUM_THREADS"] = 4
 
 using Groups
 import Groups.MatrixGroups
+
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
 using PropertyT
 
-using SymbolicWedderburn
-using SymbolicWedderburn.StarAlgebras
-using PermutationGroups
+using PropertyT.SymbolicWedderburn
+using PropertyT.PermutationGroups
+using PropertyT.StarAlgebras
 
-include(joinpath(@__DIR__, "G₂_gens.jl"))
+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₂"
 
-const HALFRADIUS = 2
-const UPPER_BOUND = Inf
-
+@info "computing group algebra structure"
 RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
 
-Δ = RG(length(S)) - sum(RG(s) for s in S)
-
+@info "computing WedderburnDecomposition"
 wd = let Σ = Weyl, RG = RG
     act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
         Dict(g => PermutationGroups.perm(g) for g in Σ),
@@ -38,57 +43,7 @@ wd = let Σ = Weyl, RG = RG
         semisimple = false,
     )
 end
-
-elt = Δ^2
-unit = Δ
-
-@time model, varP = PropertyT.sos_problem_primal(
-    elt,
-    unit,
-    wd;
-    upper_bound = UPPER_BOUND,
-    augmented = true,
-    show_progress = true,
-)
-warm = nothing
-
-begin
-    @time status, warm = PropertyT.solve(
-        model,
-        scs_optimizer(;
-            linear_solver = SCS.MKLDirectSolver,
-            eps = 1e-10,
-            max_iters = 20_000,
-            accel = 50,
-            alpha = 1.95,
-        ),
-        warm,
-    )
-
-    @info "reconstructing the solution"
-    Q = @time begin
-        wd = wd
-        Ps = [JuMP.value.(P) for P in varP]
-        if any(any(isnan, P) for P in Ps)
-            throw("solver was probably interrupted, no valid solution available")
-        end
-        Qs = real.(sqrt.(Ps))
-        PropertyT.reconstruct(Qs, wd)
-    end
-    P = Q' * Q
-
-    @info "certifying the solution"
-    @time certified, λ = PropertyT.certify_solution(
-        elt,
-        unit,
-        JuMP.objective_value(model),
-        Q;
-        halfradius = HALFRADIUS,
-        augmented = true,
-    )
-end
-
-### grading below
+@info wd
 
 function desubscriptify(symbol::Symbol)
     digits = [
@@ -107,6 +62,7 @@ function PropertyT.grading(g::MatrixGroups.MatrixElt, roots = roots)
     return roots[id]
 end
 
+Δ = RG(length(S)) - sum(RG(s) for s in S)
 Δs = PropertyT.laplacians(
     RG,
     S,
@@ -114,7 +70,7 @@ end
 )
 
 elt = PropertyT.Adj(Δs)
-elt == Δ^2 - PropertyT.Sq(Δs)
+@assert elt == Δ^2 - PropertyT.Sq(Δs)
 unit = Δ
 
 @time model, varP = PropertyT.sos_problem_primal(
@@ -123,57 +79,21 @@ unit = Δ
     wd;
     upper_bound = UPPER_BOUND,
     augmented = true,
+    show_progress = true,
 )
 
 warm = nothing
 
-begin
-    @time status, warm = PropertyT.solve(
-        model,
-        scs_optimizer(;
-            linear_solver = SCS.MKLDirectSolver,
-            eps = 1e-10,
-            max_iters = 50_000,
-            accel = 50,
-            alpha = 1.95,
-        ),
-        warm,
-    )
-
-    @info "reconstructing the solution"
-    Q = @time begin
-        wd = wd
-        Ps = [JuMP.value.(P) for P in varP]
-        if any(any(isnan, P) for P in Ps)
-            throw("solver was probably interrupted, no valid solution available")
-        end
-        Qs = real.(sqrt.(Ps))
-        PropertyT.reconstruct(Qs, wd)
-    end
-    P = Q' * Q
-
-    @info "certifying the solution"
-    @time certified, λ = PropertyT.certify_solution(
-        elt,
-        unit,
-        JuMP.objective_value(model),
-        Q;
-        halfradius = HALFRADIUS,
-        augmented = true,
-    )
-end
-
-# Δ² - 1 / 1 · Sq → -0.8818044647162608
-# Δ² - 2 / 3 · Sq → -0.1031738
-# Δ² - 1 / 2 · Sq → 0.228296213895906
-# Δ² - 1 / 3 · Sq → 0.520
-# Δ² - 0 / 1 · Sq → 0.9676851592000731
-# Sq → 0.333423
-
-# vals = [
-#     1.0 -0.8818
-#     2/3 -0.1032
-#     1/2  0.2282
-#     1/3  0.520
-#     0    0.9677
-# ]
+solve_in_loop(
+    model,
+    wd,
+    varP;
+    logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
+    optimizer = scs_optimizer(;
+        eps = 1e-10,
+        max_iters = 50_000,
+        accel = 50,
+        alpha = 1.95,
+    ),
+    data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
+)
diff --git a/scripts/G₂_has_T.jl b/scripts/G₂_has_T.jl
new file mode 100644
index 0000000..e4e6c45
--- /dev/null
+++ b/scripts/G₂_has_T.jl
@@ -0,0 +1,97 @@
+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
+
+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
diff --git a/scripts/utils.jl b/scripts/utils.jl
index 1357690..9a3f5a5 100644
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@@ -68,19 +68,19 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
         end
 
         if flag == true && certified_λ ≥ 0
-            @info "Certification done with λ = $certified_λ"
+            @info "Certification done with λ = $certified_λ" certified_λ rel_change status
             return certified_λ
         else
             rel_change =
                 abs(certified_λ - old_lambda) /
                 (abs(certified_λ) + abs(old_lambda))
-            @info "Certification failed with λ = $λ" certified_λ rel_change
+            @info "Certification failed with λ = $λ" certified_λ rel_change status
         end
 
         old_lambda = certified_λ
 
         if rel_change < 1e-9
-            @info "No progress detected, breaking"
+            @info "No progress detected, breaking" certified_λ rel_change status
             break
         end
     end

From 1a43a1b1bef1e5c69fe89cda7227052e9610899b Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 4 Apr 2023 19:58:51 +0200
Subject: [PATCH 13/27] in reconstruct: average the sum, not sum the averages!

---
 src/reconstruct.jl | 27 +++++++++++++--------------
 1 file changed, 13 insertions(+), 14 deletions(-)

diff --git a/src/reconstruct.jl b/src/reconstruct.jl
index 3971d0a..ce55619 100644
--- a/src/reconstruct.jl
+++ b/src/reconstruct.jl
@@ -12,9 +12,10 @@ function reconstruct(
     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, __group_of(wbdec), wbdec.hom)
+        res = _reconstruct!(res, M, ds)
         return res
     end
+    res = average!(zero(res), res, __group_of(wbdec), wbdec.hom)
     return res
 end
 
@@ -22,16 +23,14 @@ function _reconstruct!(
     res::AbstractMatrix,
     M::AbstractMatrix,
     ds::SymbolicWedderburn.DirectSummand,
-    G,
-    hom,
 )
-    U = SymbolicWedderburn.image_basis(ds)
-    d = SymbolicWedderburn.degree(ds)
-    Θπ = (U' * M * U) .* d
-
     res .= zero(eltype(res))
-    Θπ = average!(res, Θπ, G, hom)
-    return Θπ
+    if !iszero(M)
+        U = SymbolicWedderburn.image_basis(ds)
+        d = SymbolicWedderburn.degree(ds)
+        res = (U' * M * U) .* d
+    end
+    return res
 end
 
 function __droptol!(M::AbstractMatrix, tol)
@@ -52,18 +51,18 @@ function average!(
         <: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)
         Threads.@threads for c in axes(res, 2)
-            # note: p is a permutation,
-            # so accesses below are guaranteed to be disjoint
             for r in axes(res, 1)
-                res[r^p, c^p] += M[r, c]
+                if !iszero(M[r, c])
+                    res[r^p, c^p] += M[r, c] / o
+                end
             end
         end
     end
-    o = Groups.order(Int, G)
-    res ./= o
     return res
 end

From 132802feeb500d1636444f2accadb605c3a9e2e2 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 4 Apr 2023 23:14:56 +0200
Subject: [PATCH 14/27] use sparse matrices for invariant constraint

---
 src/sos_sdps.jl | 41 ++++++++++++++++++++++++++++++-----------
 1 file changed, 30 insertions(+), 11 deletions(-)

diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index a91b240..2a8880c 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -99,7 +99,7 @@ function invariant_constraint!(
     invariant_vec::SparseVector,
 ) where {K}
     result .= zero(eltype(result))
-    for i in SparseArrays.nonzeroinds(invariant_vec)
+    @inbounds for i in SparseArrays.nonzeroinds(invariant_vec)
         g = basis[i]
         A = cnstrs[g]
         for (idx, v) in nzpairs(A)
@@ -109,6 +109,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)
@@ -190,7 +209,8 @@ function sos_problem_primal(
         end
     end
 
-    P = map(direct_summands(wedderburn)) do ds
+    # semidefinite constraints as described by wedderburn
+    Ps = 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())
@@ -224,24 +244,23 @@ function sos_problem_primal(
         x = dot(X, iv)
         u = dot(U, iv)
 
-        M_orb = invariant_constraint!(M_orb, basis(parent(elt)), cnstrs, iv)
+        spM_orb = invariant_constraint(basis(parent(elt)), cnstrs, iv)
 
         Ms = SymbolicWedderburn.diagonalize!(
             Ms,
-            M_orb,
+            spM_orb,
             wedderburn;
             trace_preserving = true,
         )
-        # SparseArrays.droptol!.(Ms, 10 * eps(T) * max(size(M_orb)...))
-
-        # @info [nnz(m) / length(m) for m in Ms]
-
+        for M in Ms
+            SparseArrays.droptol!(M, _eps)
+        end
         if feasibility_problem
-            JuMP.@constraint(model, x == _dot(P, Ms))
+            JuMP.@constraint(model, x == _dot(Ps, Ms))
         else
-            JuMP.@constraint(model, x - λ * u == _dot(P, Ms))
+            JuMP.@constraint(model, x - λ * u == _dot(Ps, Ms))
         end
     end
     ProgressMeter.finish!(prog)
-    return model, P
+    return model, Ps
 end

From 3f2be20152432639028bbe3db3279dda27908215 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 4 Apr 2023 23:15:40 +0200
Subject: [PATCH 15/27] make nzpairs(::ConstraintMatrix) type stable

---
 src/constraint_matrix.jl | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/constraint_matrix.jl b/src/constraint_matrix.jl
index 1e95aa1..997a037 100644
--- a/src/constraint_matrix.jl
+++ b/src/constraint_matrix.jl
@@ -83,26 +83,29 @@ Base.@propagate_inbounds function Base.getindex(
     return pos - neg
 end
 
-struct NZPairsIter{T}
-    m::ConstraintMatrix{T}
+struct NZPairsIter{T,I}
+    m::ConstraintMatrix{T,I}
 end
 
 Base.eltype(::Type{NZPairsIter{T}}) where {T} = Pair{Int,T}
 Base.IteratorSize(::Type{<:NZPairsIter}) = Base.SizeUnknown()
 
 # TODO: iterate over (idx=>val) pairs combining vals
-function Base.iterate(itr::NZPairsIter, state::Tuple{Int,Int} = (1, 1))
+function Base.iterate(
+    itr::NZPairsIter,
+    state::Tuple{Int,Nothing} = (1, nothing),
+)
     k = iterate(itr.m.pos, state[1])
-    isnothing(k) && return iterate(itr, state[2])
+    isnothing(k) && return iterate(itr, (nothing, 1))
     idx, st = k
-    return idx => itr.m.val, (st, 1)
+    return idx => itr.m.val, (st, nothing)
 end
 
-function Base.iterate(itr::NZPairsIter, state::Int)
-    k = iterate(itr.m.neg, state[1])
+function Base.iterate(itr::NZPairsIter, state::Tuple{Nothing,Int})
+    k = iterate(itr.m.neg, state[2])
     isnothing(k) && return nothing
     idx, st = k
-    return idx => -itr.m.val, st
+    return idx => -itr.m.val, (nothing, st)
 end
 
 """

From 150b5c2cba28ed53444b47fc8f45362781c6314c Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 4 Apr 2023 23:48:59 +0200
Subject: [PATCH 16/27] skip the identity constraint if augmented

---
 src/sos_sdps.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index 2a8880c..11bf8b7 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -194,6 +194,10 @@ function sos_problem_primal(
         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
     feasibility_problem = iszero(orderunit)
 
     model = JuMP.Model()
@@ -240,6 +244,8 @@ function sos_problem_primal(
 
     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
 
         x = dot(X, iv)
         u = dot(U, iv)

From 5b4a7f68049e40ff2ef9b1cea114f531ebed19ce Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 4 Apr 2023 23:50:48 +0200
Subject: [PATCH 17/27] reshuffle sos_sdps for clarity

---
 src/sos_sdps.jl | 31 +++++++++++++++++--------------
 1 file changed, 17 insertions(+), 14 deletions(-)

diff --git a/src/sos_sdps.jl b/src/sos_sdps.jl
index 11bf8b7..7dfb306 100644
--- a/src/sos_sdps.jl
+++ b/src/sos_sdps.jl
@@ -198,8 +198,19 @@ function sos_problem_primal(
         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, λ)
@@ -221,27 +232,19 @@ function sos_problem_primal(
         return P
     end
 
-    begin # preallocating
+    begin # Ms are preallocated for the constraints loop
         T = eltype(wedderburn)
-        Ms = [spzeros.(T, size(p)...) for p in P]
-        M_orb = zeros(T, size(parent(elt).mstructure)...)
+        Ms = [spzeros.(T, size(p)...) for p in Ps]
+        _eps = 10 * eps(T) * max(size(parent(elt).mstructure)...)
     end
 
-    X = convert(Vector{T}, StarAlgebras.coeffs(elt))
-    U = convert(Vector{T}, StarAlgebras.coeffs(orderunit))
+    X = StarAlgebras.coeffs(elt)
+    U = StarAlgebras.coeffs(orderunit)
 
     # defining constraints based on the multiplicative structure
     cnstrs = constraints(parent(elt); augmented = augmented)
 
-    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

From 005ffc29cb8278648d85015b46743695c177ed03 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Wed, 5 Apr 2023 17:50:12 +0200
Subject: [PATCH 18/27] cosmo fails at high precision with zero condition

---
 test/1812.03456.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/1812.03456.jl b/test/1812.03456.jl
index 46802b9..383e01e 100644
--- a/test/1812.03456.jl
+++ b/test/1812.03456.jl
@@ -176,7 +176,7 @@ end
                 wd;
                 upper_bound = UB,
                 halfradius = 2,
-                optimizer = cosmo_optimizer(; accel = 50, alpha = 1.9),
+                optimizer = scs_optimizer(; accel = 50, alpha = 1.9),
             )
             @test status == JuMP.OPTIMAL
             @test !certified

From f0986982ced8a72c4f7ae207167a24b7c5535332 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Thu, 6 Apr 2023 11:39:54 +0200
Subject: [PATCH 19/27] reorganize Roots module

---
 src/roots.jl      | 133 ++++++++++++++++++++++++++--------------------
 test/Chevalley.jl |  22 ++++----
 2 files changed, 87 insertions(+), 68 deletions(-)

diff --git a/src/roots.jl b/src/roots.jl
index acbcdb2..d1fdb3d 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -7,73 +7,48 @@ export Root, isproportional, isorthogonal, ~, ⟂
 
 abstract type AbstractRoot{N,T} end
 
-struct Root{N,T} <: AbstractRoot{N,T}
-    coord::SVector{N,T}
-end
-
-Root(a) = Root(SVector(a...))
-
-function Base.:(==)(r::Root{N}, s::Root{M}) where {M,N}
-    M == N || return false
-    r.coord == s.coord || return false
-    return true
-end
-
-Base.hash(r::Root, h::UInt) = hash(r.coord, hash(Root, h))
-
-Base.:+(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord + s.coord)
-Base.:-(r::Root{N,T}, s::Root{N,T}) where {N,T} = Root{N,T}(r.coord - s.coord)
-Base.:-(r::Root{N}) where {N} = Root(-r.coord)
-
-Base.:*(a::Number, r::Root) = Root(a * r.coord)
-Base.:*(r::Root, a::Number) = a * r
-
-Base.length(r::AbstractRoot) = norm(r, 2)
-
-LinearAlgebra.norm(r::Root, p::Real = 2) = norm(r.coord, p)
-LinearAlgebra.dot(r::Root, s::Root) = dot(r.coord, s.coord)
+ℓ₂length(r::AbstractRoot) = norm(r, 2)
+ambient_dim(r::AbstractRoot) = length(r)
+Base.:*(r::AbstractRoot, a::Number) = a * r
 
 cos_angle(a, b) = dot(a, b) / (norm(a) * norm(b))
 
-function isproportional(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
-    N == M || return false
+function isproportional(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || return false
     val = abs(cos_angle(α, β))
     return isapprox(val, one(val); atol = eps(one(val)))
 end
 
-function isorthogonal(α::AbstractRoot{N}, β::AbstractRoot{M}) where {N,M}
-    N == M || return false
+function isorthogonal(α::AbstractRoot, β::AbstractRoot)
+    ambient_dim(α) == ambient_dim(β) || return false
     val = cos_angle(α, β)
     return isapprox(val, zero(val); atol = eps(one(val)))
 end
 
-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}
+function positive(roots::AbstractVector{<:AbstractRoot})
+    isempty(roots) && return empty(roots)
     pd = _positive_direction(first(roots))
     return filter(α -> dot(α, pd) > 0.0, roots)
 end
 
-function Base.show(io::IO, r::Root)
-    return print(io, "Root$(r.coord)")
+function Base.show(io::IO, r::AbstractRoot)
+    return print(io, "Root $(r.coord)")
 end
 
-function Base.show(io::IO, ::MIME"text/plain", r::Root{N}) where {N}
-    lngth² = sum(x -> x^2, r.coord)
-    l = isinteger(sqrt(lngth²)) ? "$(sqrt(lngth²))" : "√$(lngth²)"
+function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
+    l₂l = ℓ₂length(r)
+    l = isinteger(l₂l) ? "$(l₂l)" : "√$(l₂l^2)"
     return print(io, "Root in ℝ^$N of length $l\n", r.coord)
 end
 
-𝕖(N, i) = Root(ntuple(k -> k == i ? 1 : 0, N))
-𝕆(N, ::Type{T}) where {T} = Root(ntuple(_ -> zero(T), N))
-
-reflection(α::Root, β::Root) = β - Int(2dot(α, β) / dot(α, α)) * α
-function cartan(α, β)
+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 (α, β)
+        ℓ₂length(reflection(a, b) - b) / ℓ₂length(a) for a in (α, β),
+        b in (α, β)
     ]
 end
 
@@ -124,7 +99,10 @@ function classify_root_system(
     end
 end
 
-function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
+function proportional_root_from_system(
+    Ω::AbstractVector{<:AbstractRoot},
+    α::AbstractRoot,
+)
     k = findfirst(v -> isproportional(α, v), Ω)
     if isnothing(k)
         error("Line L_α not contained in root system Ω:\n α = $α\n Ω = $Ω")
@@ -132,31 +110,31 @@ function proportional_root_from_system(Ω::AbstractVector{<:Root}, α::Root)
     return Ω[k]
 end
 
-struct Plane{R<:Root}
+struct Plane{R<:AbstractRoot}
     v1::R
     v2::R
     vectors::Vector{R}
 end
 
-function Plane(α::Root, β::Root)
+function Plane(α::AbstractRoot, β::AbstractRoot)
     return Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
 end
 
-function Base.in(r::Root, plane::Plane)
+function Base.in(r::AbstractRoot, plane::Plane)
     return any(isproportional(r, v) for v in plane.vectors)
 end
 
-function _islong(α::Root, Ω)
-    lα = length(α)
-    return any(r -> lα - length(r) > eps(lα), Ω)
+function _islong(α::AbstractRoot, Ω)
+    lα = ℓ₂length(α)
+    return any(r -> lα - ℓ₂length(r) > eps(lα), Ω)
 end
 
 function classify_sub_root_system(
-    Ω::AbstractVector{<:Root{N}},
-    α::Root{N},
-    β::Root{N},
+    Ω::AbstractVector{<:AbstractRoot{N}},
+    α::AbstractRoot{N},
+    β::AbstractRoot{N},
 ) where {N}
-    @assert 1 ≤ length(unique(length, Ω)) ≤ 2
+    @assert 1 ≤ length(unique(ℓ₂length, Ω)) ≤ 2
     v = proportional_root_from_system(Ω, α)
     w = proportional_root_from_system(Ω, β)
 
@@ -197,4 +175,45 @@ function classify_sub_root_system(
     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
diff --git a/test/Chevalley.jl b/test/Chevalley.jl
index 4bb5862..e4fa863 100644
--- a/test/Chevalley.jl
+++ b/test/Chevalley.jl
@@ -22,7 +22,7 @@ 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.0, 1.0, 0.0, 0.0)
+            long = let x = (1, 1, 0, 0) .// 1
                 PropertyT.Roots.Root.(
                     union(
                         (x^g for g in Σ),
@@ -32,14 +32,14 @@ end
                 )
             end
 
-            short = let x = (1.0, 0.0, 0.0, 0.0)
+            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 / 2 .* (1.0, 1.0, 1.0, 1.0)
+            halfs = let x = (1, 1, 1, 1) .// 2
                 PropertyT.Roots.Root.(union(x .* sgn for sgn in signs))
             end
 
@@ -49,15 +49,15 @@ end
         @test length(F4) == 48
 
         a = F4[1]
-        @test isapprox(length(a), sqrt(2))
+        @test isapprox(PropertyT.Roots.ℓ₂length(a), sqrt(2))
         b = F4[6]
-        @test isapprox(length(b), sqrt(2))
+        @test isapprox(PropertyT.Roots.ℓ₂length(b), sqrt(2))
         c = a + b
-        @test isapprox(length(c), 2.0)
+        @test isapprox(PropertyT.Roots.ℓ₂length(c), 2.0)
         @test PropertyT.Roots.classify_root_system(b, c, (false, true)) == :C₂
 
-        long = F4[findfirst(r -> length(r) == sqrt(2), F4)]
-        short = F4[findfirst(r -> length(r) == 1.0, F4)]
+        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₁")])
 
@@ -94,7 +94,7 @@ end
                     perm"(1,2,3,4,5,6,7,8)",
                     perm"(1,2)",
                 )
-                long = let x = (1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+                long = let x = (1, 1, 0, 0, 0, 0, 0, 0) .// 1
                     PropertyT.Roots.Root.(
                         union(
                             (x^g for g in Σ),
@@ -108,7 +108,7 @@ end
                     p for p in Iterators.product(fill([-1, +1], 8)...) if
                     iseven(count(==(-1), p))
                 )
-                halfs = let x = 1 / 2 .* ntuple(i -> 1.0, 8)
+                halfs = let x = (1, 1, 1, 1, 1, 1, 1, 1) .// 2
                     rts = unique(PropertyT.Roots.Root(x .* sgn) for sgn in signs)
                 end
 
@@ -119,7 +119,7 @@ end
 
         @testset "E8" begin
             @test length(E8) == 240
-            @test all(r -> length(r) ≈ sqrt(2), E8)
+            @test all(r -> PropertyT.Roots.ℓ₂length(r) ≈ sqrt(2), E8)
 
             let Ω = E8, α = first(Ω)
                 counts = countmap([

From 17274f895f433a5f5aa86927ee53874d38827570 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Thu, 6 Apr 2023 11:40:49 +0200
Subject: [PATCH 20/27] set eps_*_inf for COSMO solver

---
 test/optimizers.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/optimizers.jl b/test/optimizers.jl
index 5832ca9..01b45db 100644
--- a/test/optimizers.jl
+++ b/test/optimizers.jl
@@ -44,6 +44,8 @@ function cosmo_optimizer(;
         "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,

From 58f0ccb141d816440ec7a55f57462ee890c24184 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Thu, 6 Apr 2023 13:17:52 +0200
Subject: [PATCH 21/27] use IntervalMatrices
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

IntervalMatrices use Rump algorithm to matrix multiplication
This brings time to Qint'*Qint down to ~40s which is
5-8 × slower than Q'*Q (for size n=2^13).

The naive version is ~100 × slower than Q'*Q even for n = 2^10.
---
 Manifest.toml    | 174 +++++++----------------------------------------
 Project.toml     |   2 +
 src/PropertyT.jl |   1 -
 src/certify.jl   |  13 ++--
 4 files changed, 34 insertions(+), 156 deletions(-)

diff --git a/Manifest.toml b/Manifest.toml
index 458ec6b..09a2e96 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "5452e4eb8836d03e65c342a444826b894c3e14dc"
+project_hash = "46bad61f58c1a6b455db53aa91fae7351cb279ab"
 
 [[deps.AbstractAlgebra]]
 deps = ["GroupsCore", "InteractiveUtils", "LinearAlgebra", "MacroTools", "Markdown", "Random", "RandomExtensions", "SparseArrays", "Test"]
@@ -10,23 +10,6 @@ git-tree-sha1 = "29e65c331f97db9189ef00a4c7aed8127c2fd2d4"
 uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
 version = "0.27.10"
 
-[[deps.Accessors]]
-deps = ["Compat", "CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Requires", "StaticArrays", "Test"]
-git-tree-sha1 = "beabc31fa319f9de4d16372bff31b4801e43d32c"
-uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.28"
-
-[[deps.Adapt]]
-deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef"
-uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.6.1"
-
-[[deps.ArgCheck]]
-git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
-uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
-version = "2.3.0"
-
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 version = "1.1.1"
@@ -34,20 +17,9 @@ version = "1.1.1"
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
-[[deps.BangBang]]
-deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
-git-tree-sha1 = "7fe6d92c4f281cf4ca6f2fba0ce7b299742da7ca"
-uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
-version = "0.3.37"
-
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
-[[deps.Baselet]]
-git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
-uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
-version = "0.1.1"
-
 [[deps.BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
 git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8"
@@ -113,48 +85,22 @@ deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 version = "1.0.1+0"
 
-[[deps.CompositionsBase]]
-git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
-uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
-version = "0.1.1"
-
-[[deps.ConstructionBase]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "89a9db8d28102b094992472d333674bd1a83ce2a"
-uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
-version = "1.5.1"
-
 [[deps.Cyclotomics]]
 deps = ["LRUCache", "Memoize", "Primes", "SparseArrays", "Test"]
 git-tree-sha1 = "dc2e5fd64c188399434e83fa5c10c6fa4eff962a"
 uuid = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 version = "0.3.2"
 
-[[deps.DataAPI]]
-git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630"
-uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.14.0"
-
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
 git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 version = "0.18.13"
 
-[[deps.DataValueInterfaces]]
-git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
-uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
-version = "1.0.0"
-
 [[deps.Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
-[[deps.DefineSingletons]]
-git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
-uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
-version = "0.1.2"
-
 [[deps.DiffResults]]
 deps = ["StaticArraysCore"]
 git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
@@ -187,11 +133,6 @@ git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
 uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
 version = "0.5.2"
 
-[[deps.ExternalDocstrings]]
-git-tree-sha1 = "1224740fc4d07c989949e1c1b508ebd49a65a5f6"
-uuid = "e189563c-0753-4f5e-ad5c-be4293c83fb4"
-version = "0.1.1"
-
 [[deps.FastRounding]]
 deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
 git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
@@ -201,27 +142,17 @@ version = "0.3.1"
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
-[[deps.Folds]]
-deps = ["Accessors", "BangBang", "Baselet", "DefineSingletons", "Distributed", "ExternalDocstrings", "InitialValues", "MicroCollections", "Referenceables", "Requires", "Test", "ThreadedScans", "Transducers"]
-git-tree-sha1 = "638109532de382a1f99b1aae1ca8b5d08515d85a"
-uuid = "41a02a25-b8f0-4f67-bc48-60067656b558"
-version = "0.2.8"
-
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
 git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
 version = "0.10.35"
 
-[[deps.Future]]
-deps = ["Random"]
-uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
-
 [[deps.Groups]]
-deps = ["Folds", "GroupsCore", "KnuthBendix", "LinearAlgebra", "Logging", "OrderedCollections", "PermutationGroups", "StaticArrays"]
-git-tree-sha1 = "5dbf642ee0048e6ad5f0bda11af17e40b8e8dd2f"
+deps = ["GroupsCore", "KnuthBendix", "LinearAlgebra", "Logging", "OrderedCollections", "PermutationGroups", "Random", "StaticArrays"]
+git-tree-sha1 = "47967c88672431b6d3624d1cb200abdc15aa0b47"
 uuid = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
-version = "0.7.5"
+version = "0.7.6"
 
 [[deps.GroupsCore]]
 deps = ["Markdown", "Random"]
@@ -229,11 +160,6 @@ git-tree-sha1 = "9e1a5e9f3b81ad6a5c613d181664a0efc6fe6dd7"
 uuid = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 version = "0.4.0"
 
-[[deps.InitialValues]]
-git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3"
-uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c"
-version = "0.3.1"
-
 [[deps.IntegerMathUtils]]
 git-tree-sha1 = "f366daebdfb079fd1fe4e3d560f99a0c892e15bc"
 uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
@@ -249,6 +175,12 @@ git-tree-sha1 = "c1c88395d09366dae431556bcb598ad08fa1392b"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 version = "0.20.8"
 
+[[deps.IntervalMatrices]]
+deps = ["IntervalArithmetic", "LinearAlgebra", "Random", "Reexport", "SparseArrays"]
+git-tree-sha1 = "b5a981918e38fb40aeea659ec6511187fa488c14"
+uuid = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
+version = "0.8.3"
+
 [[deps.InverseFunctions]]
 deps = ["Test"]
 git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f"
@@ -260,11 +192,6 @@ git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
 version = "0.2.2"
 
-[[deps.IteratorInterfaceExtensions]]
-git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
-uuid = "82899510-4779-5014-852e-03e436cf321d"
-version = "1.0.0"
-
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
 git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
@@ -279,9 +206,9 @@ version = "0.21.3"
 
 [[deps.JuMP]]
 deps = ["LinearAlgebra", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays"]
-git-tree-sha1 = "611b9f12f02c587d860c813743e6cec6264e94d8"
+git-tree-sha1 = "4ec0e68fecbbe1b78db2ddf1ac573963ed5adebc"
 uuid = "4076af6c-e467-56ae-b986-b466b2749572"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.KnuthBendix]]
 deps = ["MacroTools", "ProgressMeter"]
@@ -341,9 +268,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[deps.MathOptInterface]]
 deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
-git-tree-sha1 = "f219b62e601c2f2e8adb7b6c48db8a9caf381c82"
+git-tree-sha1 = "88551cb56065a7bf9276af1d1970d2bd0d4b6ed4"
 uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-version = "1.13.1"
+version = "1.14.0"
 
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -356,12 +283,6 @@ git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
 uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
 version = "0.4.4"
 
-[[deps.MicroCollections]]
-deps = ["BangBang", "InitialValues", "Setfield"]
-git-tree-sha1 = "629afd7d10dbc6935ec59b32daeb33bc4460a42e"
-uuid = "128add7d-3638-4c79-886c-908ea0c25c34"
-version = "0.1.4"
-
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
@@ -402,9 +323,9 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.4.1"
+version = "1.6.0"
 
 [[deps.Parsers]]
 deps = ["Dates", "SnoopPrecompile"]
@@ -469,17 +390,10 @@ git-tree-sha1 = "261dddd3b862bd2c940cf6ca4d1c8fe593e457c8"
 uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 version = "1.3.3"
 
-[[deps.Referenceables]]
-deps = ["Adapt"]
-git-tree-sha1 = "e681d3bfa49cd46c3c161505caddf20f0e62aaa9"
-uuid = "42d2dcc6-99eb-4e98-b66c-637b7d73030e"
-version = "0.1.2"
-
-[[deps.Requires]]
-deps = ["UUIDs"]
-git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
-uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.3.0"
+[[deps.Reexport]]
+git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.2.2"
 
 [[deps.RoundingEmulator]]
 git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
@@ -498,12 +412,6 @@ git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
 uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
 version = "0.2.1"
 
-[[deps.Setfield]]
-deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
-git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
-uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
-version = "1.1.1"
-
 [[deps.SnoopPrecompile]]
 deps = ["Preferences"]
 git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
@@ -523,12 +431,6 @@ git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
 version = "2.2.0"
 
-[[deps.SplittablesBase]]
-deps = ["Setfield", "Test"]
-git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5"
-uuid = "171d559e-b47b-412a-8079-5efa626c420e"
-version = "0.1.15"
-
 [[deps.StarAlgebras]]
 deps = ["LinearAlgebra", "SparseArrays"]
 git-tree-sha1 = "265b89a5dfb38fe94ad48b997a253b2393fce6f1"
@@ -537,9 +439,9 @@ version = "0.2.0"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
-git-tree-sha1 = "6aa098ef1012364f2ede6b17bf358c7f1fbe90d4"
+git-tree-sha1 = "b8d897fe7fa688e93aef573711cb207c08c9e11e"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.5.17"
+version = "1.5.19"
 
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
@@ -552,7 +454,7 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.SymbolicWedderburn]]
 deps = ["Cyclotomics", "GroupsCore", "LinearAlgebra", "PermutationGroups", "Primes", "SparseArrays", "StarAlgebras"]
-git-tree-sha1 = "3a098aae261c4c57d6ee41dab62973065f0cb68b"
+git-tree-sha1 = "7afdc5576b8fc49032c1f765c1514e5369fb77fe"
 repo-rev = "enh/towards_simple_projections"
 repo-url = "https://github.com/kalmarek/SymbolicWedderburn.jl.git"
 uuid = "858aa9a9-4c7c-4c62-b466-2421203962a2"
@@ -563,18 +465,6 @@ deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 version = "1.0.0"
 
-[[deps.TableTraits]]
-deps = ["IteratorInterfaceExtensions"]
-git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
-uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
-version = "1.0.1"
-
-[[deps.Tables]]
-deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
-git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec"
-uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.10.1"
-
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
@@ -584,24 +474,12 @@ version = "1.10.1"
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[[deps.ThreadedScans]]
-deps = ["ArgCheck"]
-git-tree-sha1 = "ca1ba3000289eacba571aaa4efcefb642e7a1de6"
-uuid = "24d252fe-5d94-4a69-83ea-56a14333d47a"
-version = "0.1.0"
-
 [[deps.TranscodingStreams]]
 deps = ["Random", "Test"]
 git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.9.11"
 
-[[deps.Transducers]]
-deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"]
-git-tree-sha1 = "c42fa452a60f022e9e087823b47e5a5f8adc53d5"
-uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999"
-version = "0.4.75"
-
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
@@ -614,12 +492,6 @@ deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
 version = "1.2.12+3"
 
-[[deps.ZygoteRules]]
-deps = ["MacroTools"]
-git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
-uuid = "700de1a5-db45-46bc-99cf-38207098b444"
-version = "0.2.2"
-
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
diff --git a/Project.toml b/Project.toml
index 75f97da..c13977e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ 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"
@@ -17,6 +18,7 @@ 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"
diff --git a/src/PropertyT.jl b/src/PropertyT.jl
index d4a664a..60b0aef 100644
--- a/src/PropertyT.jl
+++ b/src/PropertyT.jl
@@ -3,7 +3,6 @@ module PropertyT
 using LinearAlgebra
 using SparseArrays
 
-using IntervalArithmetic
 using JuMP
 
 using Groups
diff --git a/src/certify.jl b/src/certify.jl
index 43648aa..b93c843 100644
--- a/src/certify.jl
+++ b/src/certify.jl
@@ -1,3 +1,6 @@
+import IntervalArithmetic
+import IntervalMatrices
+
 function augment_columns!(Q::AbstractMatrix)
     for c in eachcol(Q)
         c .-= sum(c) ./ length(c)
@@ -63,7 +66,7 @@ function sufficient_λ(residual::StarAlgebras.AlgebraElement, λ; halfradius)
     suff_λ = λ - 2.0^(2ceil(log2(halfradius))) * L1_norm
 
     eq_sign = let T = eltype(residual)
-        if T <: Interval
+        if T <: IntervalArithmetic.Interval
             "∈"
         elseif T <: Union{Rational,Integer}
             "="
@@ -119,8 +122,10 @@ function certify_solution(
         return false, λ_flpoint
     end
 
-    λ_int = @interval(λ)
-    Q_int = [@interval(q) for q in Q]
+    λ_int = IntervalArithmetic.@interval(λ)
+    Q_int = IntervalMatrices.IntervalMatrix([
+        IntervalArithmetic.@interval(q) for q in Q
+    ])
 
     check, sos_int = @time if should_we_augment
         @info("Projecting columns of Q to the augmentation ideal...")
@@ -141,5 +146,5 @@ function certify_solution(
     λ_certified =
         sufficient_λ(elt, orderunit, λ_int, sos_int; halfradius = halfradius)
 
-    return check && inf(λ_certified) > 0.0, λ_certified
+    return check && IntervalArithmetic.inf(λ_certified) > 0.0, λ_certified
 end

From 56aed884161210688aa2ad6d5570a612a10ba73f Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Thu, 6 Apr 2023 16:16:27 +0200
Subject: [PATCH 22/27] =?UTF-8?q?use=20MKLDirect=20for=20G=E2=82=82?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 scripts/G₂_Adj.jl | 71 ++++++++++++++++++++++++++++++++++++---------
 1 file changed, 57 insertions(+), 14 deletions(-)

diff --git a/scripts/G₂_Adj.jl b/scripts/G₂_Adj.jl
index 339d99c..e9acfc9 100644
--- a/scripts/G₂_Adj.jl
+++ b/scripts/G₂_Adj.jl
@@ -1,6 +1,6 @@
 using LinearAlgebra
 BLAS.set_num_threads(8)
-
+using MKL_jll
 ENV["OMP_NUM_THREADS"] = 4
 
 using Groups
@@ -84,16 +84,59 @@ unit = Δ
 
 warm = nothing
 
-solve_in_loop(
-    model,
-    wd,
-    varP;
-    logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
-    optimizer = scs_optimizer(;
-        eps = 1e-10,
-        max_iters = 50_000,
-        accel = 50,
-        alpha = 1.95,
-    ),
-    data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
-)
+let status = JuMP.OPTIMIZE_NOT_CALLED, warm = warm, eps = 1e-9
+    certified, λ = false, 0.0
+    while status ≠ JuMP.OPTIMAL
+        @time status, warm = PropertyT.solve(
+            model,
+            scs_optimizer(;
+                linear_solver = SCS.MKLDirectSolver,
+                eps = eps,
+                max_iters = 100_000,
+                accel = 50,
+                alpha = 1.95,
+            ),
+            warm,
+        )
+
+        @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
+
+        @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
+end
+
+# solve_in_loop(
+#     model,
+#     wd,
+#     varP;
+#     logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
+#     optimizer = scs_optimizer(;
+#         eps = 1e-10,
+#         max_iters = 50_000,
+#         accel = 50,
+#         alpha = 1.95,
+#     ),
+#     data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
+# )

From 263af398eb99b08bac90f4e9c276df8bfa707f98 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 11 Apr 2023 10:31:54 +0200
Subject: [PATCH 23/27] =?UTF-8?q?G=E2=82=82=5FAdj:=20bring=20back=20solve?=
 =?UTF-8?q?=5Fin=5Floop?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 scripts/G₂_Adj.jl | 72 +++++++++------------------------------------
 1 file changed, 14 insertions(+), 58 deletions(-)

diff --git a/scripts/G₂_Adj.jl b/scripts/G₂_Adj.jl
index e9acfc9..8a499d0 100644
--- a/scripts/G₂_Adj.jl
+++ b/scripts/G₂_Adj.jl
@@ -82,61 +82,17 @@ unit = Δ
     show_progress = true,
 )
 
-warm = nothing
-
-let status = JuMP.OPTIMIZE_NOT_CALLED, warm = warm, eps = 1e-9
-    certified, λ = false, 0.0
-    while status ≠ JuMP.OPTIMAL
-        @time status, warm = PropertyT.solve(
-            model,
-            scs_optimizer(;
-                linear_solver = SCS.MKLDirectSolver,
-                eps = eps,
-                max_iters = 100_000,
-                accel = 50,
-                alpha = 1.95,
-            ),
-            warm,
-        )
-
-        @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
-
-        @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
-end
-
-# solve_in_loop(
-#     model,
-#     wd,
-#     varP;
-#     logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
-#     optimizer = scs_optimizer(;
-#         eps = 1e-10,
-#         max_iters = 50_000,
-#         accel = 50,
-#         alpha = 1.95,
-#     ),
-#     data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
-# )
+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),
+)

From 9afbccda1566d6dd209f1ebfd32d230795643eda Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 11 Apr 2023 10:32:24 +0200
Subject: [PATCH 24/27] drop zeros when reconstructing the solution

---
 scripts/utils.jl | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/scripts/utils.jl b/scripts/utils.jl
index 9a3f5a5..840eb80 100644
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@@ -11,11 +11,17 @@ function get_solution(model)
     return solution
 end
 
-function get_solution(model, wd, varP)
+function get_solution(model, wd, varP, eps = 1e-10)
     λ = JuMP.value(model[:λ])
 
-    Qs = [real.(sqrt(JuMP.value.(P))) for P in varP]
-    Q = PropertyT.reconstruct(Qs, wd)
+    @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
+
     solution = Dict(:λ => λ, :Q => Q)
 
     return solution

From ed832d69fb3fe75fc4f12324736e587a3bf216c0 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Thu, 13 Apr 2023 01:00:22 +0200
Subject: [PATCH 25/27] fix for solve_in_loop script

---
 scripts/utils.jl | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/scripts/utils.jl b/scripts/utils.jl
index 840eb80..5dd6c80 100644
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@@ -81,14 +81,13 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
                 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
         end
 
         old_lambda = certified_λ
-
-        if rel_change < 1e-9
-            @info "No progress detected, breaking" certified_λ rel_change status
-            break
-        end
     end
 
     return status == JuMP.OPTIMAL ? old_lambda : NaN

From 34e19768d4200e30396913a63d1101c8e6ab51c5 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Mon, 8 May 2023 17:49:32 +0200
Subject: [PATCH 26/27] rm Manifest and update scripts

---
 Manifest.toml         | 508 ------------------------------------------
 scripts/G₂_has_T.jl |   1 +
 scripts/SpN_Adj.jl    |   2 +-
 scripts/utils.jl      |   2 +-
 4 files changed, 3 insertions(+), 510 deletions(-)
 delete mode 100644 Manifest.toml

diff --git a/Manifest.toml b/Manifest.toml
deleted file mode 100644
index 09a2e96..0000000
--- a/Manifest.toml
+++ /dev/null
@@ -1,508 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
-julia_version = "1.8.5"
-manifest_format = "2.0"
-project_hash = "46bad61f58c1a6b455db53aa91fae7351cb279ab"
-
-[[deps.AbstractAlgebra]]
-deps = ["GroupsCore", "InteractiveUtils", "LinearAlgebra", "MacroTools", "Markdown", "Random", "RandomExtensions", "SparseArrays", "Test"]
-git-tree-sha1 = "29e65c331f97db9189ef00a4c7aed8127c2fd2d4"
-uuid = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
-version = "0.27.10"
-
-[[deps.ArgTools]]
-uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
-version = "1.1.1"
-
-[[deps.Artifacts]]
-uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
-
-[[deps.Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[deps.BenchmarkTools]]
-deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8"
-uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.3.2"
-
-[[deps.Bzip2_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
-uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
-version = "1.0.8+0"
-
-[[deps.CRlibm]]
-deps = ["CRlibm_jll"]
-git-tree-sha1 = "32abd86e3c2025db5172aa182b982debed519834"
-uuid = "96374032-68de-5a5b-8d9e-752f78720389"
-version = "1.0.1"
-
-[[deps.CRlibm_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "e329286945d0cfc04456972ea732551869af1cfc"
-uuid = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
-version = "1.0.1+0"
-
-[[deps.ChainRulesCore]]
-deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9"
-uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.15.7"
-
-[[deps.ChangesOfVariables]]
-deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
-git-tree-sha1 = "485193efd2176b88e6622a39a246f8c5b600e74e"
-uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
-version = "0.1.6"
-
-[[deps.CodecBzip2]]
-deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
-git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
-uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
-version = "0.7.2"
-
-[[deps.CodecZlib]]
-deps = ["TranscodingStreams", "Zlib_jll"]
-git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83"
-uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.7.1"
-
-[[deps.CommonSubexpressions]]
-deps = ["MacroTools", "Test"]
-git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
-uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-version = "0.3.0"
-
-[[deps.Compat]]
-deps = ["Dates", "LinearAlgebra", "UUIDs"]
-git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.6.1"
-
-[[deps.CompilerSupportLibraries_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.1+0"
-
-[[deps.Cyclotomics]]
-deps = ["LRUCache", "Memoize", "Primes", "SparseArrays", "Test"]
-git-tree-sha1 = "dc2e5fd64c188399434e83fa5c10c6fa4eff962a"
-uuid = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
-version = "0.3.2"
-
-[[deps.DataStructures]]
-deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
-uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.13"
-
-[[deps.Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[deps.DiffResults]]
-deps = ["StaticArraysCore"]
-git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
-uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-version = "1.1.0"
-
-[[deps.DiffRules]]
-deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "a4ad7ef19d2cdc2eff57abbbe68032b1cd0bd8f8"
-uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.13.0"
-
-[[deps.Distributed]]
-deps = ["Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[deps.DocStringExtensions]]
-deps = ["LibGit2"]
-git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
-uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.9.3"
-
-[[deps.Downloads]]
-deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
-uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
-version = "1.6.0"
-
-[[deps.ErrorfreeArithmetic]]
-git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
-uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
-version = "0.5.2"
-
-[[deps.FastRounding]]
-deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
-git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
-uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f"
-version = "0.3.1"
-
-[[deps.FileWatching]]
-uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
-
-[[deps.ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "00e252f4d706b3d55a8863432e742bf5717b498d"
-uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.35"
-
-[[deps.Groups]]
-deps = ["GroupsCore", "KnuthBendix", "LinearAlgebra", "Logging", "OrderedCollections", "PermutationGroups", "Random", "StaticArrays"]
-git-tree-sha1 = "47967c88672431b6d3624d1cb200abdc15aa0b47"
-uuid = "5d8bd718-bd84-11e8-3b40-ad14f4a32557"
-version = "0.7.6"
-
-[[deps.GroupsCore]]
-deps = ["Markdown", "Random"]
-git-tree-sha1 = "9e1a5e9f3b81ad6a5c613d181664a0efc6fe6dd7"
-uuid = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
-version = "0.4.0"
-
-[[deps.IntegerMathUtils]]
-git-tree-sha1 = "f366daebdfb079fd1fe4e3d560f99a0c892e15bc"
-uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
-version = "0.1.0"
-
-[[deps.InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[deps.IntervalArithmetic]]
-deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "RoundingEmulator", "SetRounding", "StaticArrays"]
-git-tree-sha1 = "c1c88395d09366dae431556bcb598ad08fa1392b"
-uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-version = "0.20.8"
-
-[[deps.IntervalMatrices]]
-deps = ["IntervalArithmetic", "LinearAlgebra", "Random", "Reexport", "SparseArrays"]
-git-tree-sha1 = "b5a981918e38fb40aeea659ec6511187fa488c14"
-uuid = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
-version = "0.8.3"
-
-[[deps.InverseFunctions]]
-deps = ["Test"]
-git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f"
-uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.8"
-
-[[deps.IrrationalConstants]]
-git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
-uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
-version = "0.2.2"
-
-[[deps.JLLWrappers]]
-deps = ["Preferences"]
-git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
-uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.4.1"
-
-[[deps.JSON]]
-deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
-uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.3"
-
-[[deps.JuMP]]
-deps = ["LinearAlgebra", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays"]
-git-tree-sha1 = "4ec0e68fecbbe1b78db2ddf1ac573963ed5adebc"
-uuid = "4076af6c-e467-56ae-b986-b466b2749572"
-version = "1.10.0"
-
-[[deps.KnuthBendix]]
-deps = ["MacroTools", "ProgressMeter"]
-git-tree-sha1 = "a5ef62c9f2b4461b246a610b7402f5f2dbffbfaa"
-uuid = "c2604015-7b3d-4a30-8a26-9074551ec60a"
-version = "0.4.0"
-
-[[deps.LRUCache]]
-git-tree-sha1 = "d862633ef6097461037a00a13f709a62ae4bdfdd"
-uuid = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
-version = "1.4.0"
-
-[[deps.LibCURL]]
-deps = ["LibCURL_jll", "MozillaCACerts_jll"]
-uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.3"
-
-[[deps.LibCURL_jll]]
-deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
-uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "7.84.0+0"
-
-[[deps.LibGit2]]
-deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-
-[[deps.LibSSH2_jll]]
-deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
-uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.10.2+0"
-
-[[deps.Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[deps.LinearAlgebra]]
-deps = ["Libdl", "libblastrampoline_jll"]
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-
-[[deps.LogExpFunctions]]
-deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "0a1b7c2863e44523180fdb3146534e265a91870b"
-uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.23"
-
-[[deps.Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[deps.MacroTools]]
-deps = ["Markdown", "Random"]
-git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2"
-uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.10"
-
-[[deps.Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[deps.MathOptInterface]]
-deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "Printf", "SnoopPrecompile", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
-git-tree-sha1 = "88551cb56065a7bf9276af1d1970d2bd0d4b6ed4"
-uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-version = "1.14.0"
-
-[[deps.MbedTLS_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.0+0"
-
-[[deps.Memoize]]
-deps = ["MacroTools"]
-git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
-uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
-version = "0.4.4"
-
-[[deps.Mmap]]
-uuid = "a63ad114-7e13-5084-954f-fe012c677804"
-
-[[deps.MozillaCACerts_jll]]
-uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2022.2.1"
-
-[[deps.MutableArithmetics]]
-deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "3295d296288ab1a0a2528feb424b854418acff57"
-uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
-version = "1.2.3"
-
-[[deps.NaNMath]]
-deps = ["OpenLibm_jll"]
-git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
-uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "1.0.2"
-
-[[deps.NetworkOptions]]
-uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
-version = "1.2.0"
-
-[[deps.OpenBLAS_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
-uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.20+0"
-
-[[deps.OpenLibm_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+0"
-
-[[deps.OpenSpecFun_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
-uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
-version = "0.5.5+0"
-
-[[deps.OrderedCollections]]
-git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282"
-uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.0"
-
-[[deps.Parsers]]
-deps = ["Dates", "SnoopPrecompile"]
-git-tree-sha1 = "478ac6c952fddd4399e71d4779797c538d0ff2bf"
-uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.5.8"
-
-[[deps.PermutationGroups]]
-deps = ["AbstractAlgebra", "GroupsCore", "Markdown", "Random"]
-git-tree-sha1 = "1bfba1a836e2c085270d3f5657803e0902e20564"
-uuid = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
-version = "0.3.3"
-
-[[deps.Pkg]]
-deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.8.0"
-
-[[deps.Preferences]]
-deps = ["TOML"]
-git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
-uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.3.0"
-
-[[deps.Primes]]
-deps = ["IntegerMathUtils"]
-git-tree-sha1 = "311a2aa90a64076ea0fac2ad7492e914e6feeb81"
-uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
-version = "0.5.3"
-
-[[deps.Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[deps.Profile]]
-deps = ["Printf"]
-uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
-
-[[deps.ProgressMeter]]
-deps = ["Distributed", "Printf"]
-git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9"
-uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
-version = "1.7.2"
-
-[[deps.REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
-uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[deps.Random]]
-deps = ["SHA", "Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[deps.RandomExtensions]]
-deps = ["Random", "SparseArrays"]
-git-tree-sha1 = "062986376ce6d394b23d5d90f01d81426113a3c9"
-uuid = "fb686558-2515-59ef-acaa-46db3789a887"
-version = "0.4.3"
-
-[[deps.RecipesBase]]
-deps = ["SnoopPrecompile"]
-git-tree-sha1 = "261dddd3b862bd2c940cf6ca4d1c8fe593e457c8"
-uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-version = "1.3.3"
-
-[[deps.Reexport]]
-git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
-uuid = "189a3867-3050-52da-a836-e630ba90ab69"
-version = "1.2.2"
-
-[[deps.RoundingEmulator]]
-git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
-uuid = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
-version = "0.2.1"
-
-[[deps.SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-version = "0.7.0"
-
-[[deps.Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[deps.SetRounding]]
-git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
-uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
-version = "0.2.1"
-
-[[deps.SnoopPrecompile]]
-deps = ["Preferences"]
-git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
-uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
-version = "1.0.3"
-
-[[deps.Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[deps.SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[deps.SpecialFunctions]]
-deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.2.0"
-
-[[deps.StarAlgebras]]
-deps = ["LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "265b89a5dfb38fe94ad48b997a253b2393fce6f1"
-uuid = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
-version = "0.2.0"
-
-[[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"]
-git-tree-sha1 = "b8d897fe7fa688e93aef573711cb207c08c9e11e"
-uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.5.19"
-
-[[deps.StaticArraysCore]]
-git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a"
-uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.0"
-
-[[deps.Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[deps.SymbolicWedderburn]]
-deps = ["Cyclotomics", "GroupsCore", "LinearAlgebra", "PermutationGroups", "Primes", "SparseArrays", "StarAlgebras"]
-git-tree-sha1 = "7afdc5576b8fc49032c1f765c1514e5369fb77fe"
-repo-rev = "enh/towards_simple_projections"
-repo-url = "https://github.com/kalmarek/SymbolicWedderburn.jl.git"
-uuid = "858aa9a9-4c7c-4c62-b466-2421203962a2"
-version = "0.3.4"
-
-[[deps.TOML]]
-deps = ["Dates"]
-uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
-version = "1.0.0"
-
-[[deps.Tar]]
-deps = ["ArgTools", "SHA"]
-uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
-version = "1.10.1"
-
-[[deps.Test]]
-deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[deps.TranscodingStreams]]
-deps = ["Random", "Test"]
-git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f"
-uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.9.11"
-
-[[deps.UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
-[[deps.Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
-
-[[deps.Zlib_jll]]
-deps = ["Libdl"]
-uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.12+3"
-
-[[deps.libblastrampoline_jll]]
-deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
-uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.1.1+0"
-
-[[deps.nghttp2_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.48.0+0"
-
-[[deps.p7zip_jll]]
-deps = ["Artifacts", "Libdl"]
-uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+0"
diff --git a/scripts/G₂_has_T.jl b/scripts/G₂_has_T.jl
index e4e6c45..cde9048 100644
--- a/scripts/G₂_has_T.jl
+++ b/scripts/G₂_has_T.jl
@@ -59,6 +59,7 @@ unit = Δ
 )
 warm = nothing
 status = JuMP.OPTIMIZE_NOT_CALLED
+certified, λ = false, nothing
 
 while status ≠ JuMP.OPTIMAL
     @time status, warm = PropertyT.solve(
diff --git a/scripts/SpN_Adj.jl b/scripts/SpN_Adj.jl
index 70ff7fe..35b9a6b 100644
--- a/scripts/SpN_Adj.jl
+++ b/scripts/SpN_Adj.jl
@@ -69,7 +69,7 @@ solve_in_loop(
     model,
     wd,
     varP;
-    logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-InfΔ",
+    logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Adj_C₂-$(UPPER_BOUND)Δ",
     optimizer = cosmo_optimizer(;
         eps = 1e-10,
         max_iters = 20_000,
diff --git a/scripts/utils.jl b/scripts/utils.jl
index 5dd6c80..4c4b394 100644
--- a/scripts/utils.jl
+++ b/scripts/utils.jl
@@ -74,7 +74,7 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
         end
 
         if flag == true && certified_λ ≥ 0
-            @info "Certification done with λ = $certified_λ" certified_λ rel_change status
+            @info "Certification done with λ = $certified_λ" certified_λ status
             return certified_λ
         else
             rel_change =

From d40a0fe1170e5b0743809bda3fe827cc5251011a Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 9 May 2023 01:02:39 +0200
Subject: [PATCH 27/27] fix error with printing of roots

---
 src/roots.jl     |  6 +++---
 test/roots.jl    | 10 ++++++++++
 test/runtests.jl |  1 +
 3 files changed, 14 insertions(+), 3 deletions(-)
 create mode 100644 test/roots.jl

diff --git a/src/roots.jl b/src/roots.jl
index d1fdb3d..92a915c 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -5,7 +5,7 @@ using LinearAlgebra
 
 export Root, isproportional, isorthogonal, ~, ⟂
 
-abstract type AbstractRoot{N,T} end
+abstract type AbstractRoot{N,T} end # <: AbstractVector{T} ?
 
 ℓ₂length(r::AbstractRoot) = norm(r, 2)
 ambient_dim(r::AbstractRoot) = length(r)
@@ -37,8 +37,8 @@ end
 
 function Base.show(io::IO, ::MIME"text/plain", r::AbstractRoot)
     l₂l = ℓ₂length(r)
-    l = isinteger(l₂l) ? "$(l₂l)" : "√$(l₂l^2)"
-    return print(io, "Root in ℝ^$N of length $l\n", r.coord)
+    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)
 end
 
 function reflection(α::AbstractRoot, β::AbstractRoot)
diff --git a/test/roots.jl b/test/roots.jl
new file mode 100644
index 0000000..7cb32d0
--- /dev/null
+++ b/test/roots.jl
@@ -0,0 +1,10 @@
+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
diff --git a/test/runtests.jl b/test/runtests.jl
index 5e6a9be..a306e55 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,6 +24,7 @@ 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