From a5f5a4ea35539240e73b5cba4b75c80c6c4201f5 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Sun, 19 Mar 2023 20:33:28 +0100 Subject: [PATCH] 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