lots of re-formatting

.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

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]),
+        (sizes[halfradius], sizes[halfradius]);
-        precompute=false,
+        precompute = false,
     )
     return RG, S, sizes
 end

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(α))

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 =
-    should_we_augment = !augmented && StarAlgebras.aug(elt) == StarAlgebras.aug(orderunit) == 0
+        !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

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} =
+function ConstraintMatrix(
-    ConstraintMatrix{T}(nzeros, n, m, val)
+    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

src/gradings.jl

@@ -1,7 +1,8 @@
 ## something about roots
 
-Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N} =
+function Roots.Root(e::MatrixGroups.ElementaryMatrix{N}) where {N}
-    Roots.𝕖(N, e.i) - Roots.𝕖(N, e.j)
+    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.classify_sub_root_system(roots, first(α), first(β)) == subtype
-                roots,
+        );
-                first(α),
+        init = zero(first(values(rootsystem))),
-                first(β),
+    )
-            ) == subtype
+end
-        ),
+
-        init=zero(first(values(rootsystem))),
+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

src/reconstruct.jl

@@ -7,35 +7,31 @@ end
 
 function reconstruct(
     Ms::AbstractVector{<:AbstractMatrix},
-    wbdec::WedderburnDecomposition;
+    wbdec::WedderburnDecomposition,
-    atol=eps(real(eltype(wbdec))) * 10__outer_dim(wbdec)
 )
     n = __outer_dim(wbdec)
     res = sum(zip(Ms, SymbolicWedderburn.direct_summands(wbdec))) do (M, ds)
         res = similar(M, n, n)
-        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;
+    hom,
-    atol=eps(real(eltype(ds))) * 10max(size(res)...)
 )
-    res .= zero(eltype(res))
     U = SymbolicWedderburn.image_basis(ds)
     d = SymbolicWedderburn.degree(ds)
-    tmp = (U' * M * U) .* d
+    Θπ = (U' * M * U) .* d
 
-    res = average!(res, tmp, G, hom)
+    res .= zero(eltype(res))
-    if eltype(res) <: AbstractFloat
+    Θπ = average!(res, Θπ, G, hom)
-        __droptol!(res, atol) # TODO: is this really necessary?!
+    return Θπ
-    end
-    return res
 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

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} =
+function Plane(α::Root, β::Root)
-    Plane(α, β, [a * α + b * β for a in -3:3 for b in -3:3])
+    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
 

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])
-    for (ct, idx) in pairs(constraint_map)
+            JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
-        JuMP.set_start_value.(idx, warmstart.slack[ct])
+        end
-        JuMP.set_dual_start_value.(idx, warmstart.dual[ct])
+    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
 

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);
+    order_unit::StarAlgebras.AlgebraElement = zero(elt);
-    lower_bound=-Inf
+    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);
+    order_unit::StarAlgebras.AlgebraElement = zero(elt);
-    upper_bound=Inf,
+    upper_bound = Inf,
-    augmented::Bool=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(order_unit))
+    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...
+    kwargs...,
-) = sos_problem_primal(elt, zero(elt), wedderburn; 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,
+    upper_bound = Inf,
-    augmented=iszero(StarAlgebras.aug(elt)) && iszero(StarAlgebras.aug(orderunit)),
+    augmented = iszero(StarAlgebras.aug(elt)) &&
-    check_orthogonality=true,
+        iszero(StarAlgebras.aug(orderunit)),
-    show_progress=false
+    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)),
+        length(invariant_vectors(wedderburn));
-        dt=1,
+        dt = 1,
-        desc="Adding constraints... ",
+        desc = "Adding constraints: ",
-        enabled=show_progress,
+        enabled = show_progress,
-        barlen=60,
+        barlen = 60,
-        showspeed=true
+        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)

test/optimizers.jl

@@ -4,12 +4,12 @@ import JuMP
 import SCS
 
 function scs_optimizer(;
-    accel=10,
+    accel = 10,
-    alpha=1.5,
+    alpha = 1.5,
-    eps=1e-9,
+    eps = 1e-9,
-    max_iters=100_000,
+    max_iters = 100_000,
-    verbose=true,
+    verbose = true,
-    linear_solver=SCS.DirectSolver
+    linear_solver = SCS.DirectSolver,
 )
     return JuMP.optimizer_with_attributes(
         SCS.Optimizer,
@@ -28,20 +28,18 @@ end
 import COSMO
 
 function cosmo_optimizer(;
-    accel=15,
+    accel = 15,
-    alpha=1.6,
+    alpha = 1.6,
-    eps=1e-9,
+    eps = 1e-9,
-    max_iters=100_000,
+    max_iters = 100_000,
-    verbose=true,
+    verbose = true,
-    verbose_timing=verbose,
+    verbose_timing = verbose,
-    decompose=false
+    decompose = false,
 )
     return JuMP.optimizer_with_attributes(
         COSMO.Optimizer,
-        "accelerator" => COSMO.with_options(
+        "accelerator" =>
-            COSMO.AndersonAccelerator,
+            COSMO.with_options(COSMO.AndersonAccelerator; mem = max(accel, 2)),
-            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