From f4936dd50a14413d5d83b479fe7a2af69eb02812 Mon Sep 17 00:00:00 2001 From: Marek Kaluba Date: Tue, 30 May 2023 16:35:55 +0200 Subject: [PATCH] add/update scripts for SLNZ/SpNZ --- scripts/SLNZ_has_T.jl | 94 +++++++++++++++++++++++++++++++++++++++++++ scripts/SpNZ_has_T.jl | 83 ++++++++++++++++++++++++++++++++++++++ scripts/utils.jl | 21 +++++----- 3 files changed, 189 insertions(+), 9 deletions(-) create mode 100644 scripts/SLNZ_has_T.jl create mode 100644 scripts/SpNZ_has_T.jl diff --git a/scripts/SLNZ_has_T.jl b/scripts/SLNZ_has_T.jl new file mode 100644 index 0000000..0c36f4e --- /dev/null +++ b/scripts/SLNZ_has_T.jl @@ -0,0 +1,94 @@ +using LinearAlgebra +using MKL_jll +BLAS.set_num_threads(4) + +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")) + +const N = parsed_args["N"] +const HALFRADIUS = parsed_args["halfradius"] +const UPPER_BOUND = parsed_args["upper_bound"] + +G = MatrixGroups.SpecialLinearGroup{N}(Int8) +@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 RG = RG, N = N + G = StarAlgebras.object(RG) + P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) + Σ = Groups.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + act = PropertyT.action_by_conjugation(G, Σ) + + wdfl = @time SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]), + ) +end +@info wd + +Δ = RG(length(S)) - sum(RG(s) for s in S) +elt = Δ^2 +unit = Δ +warm = nothing + +@info "defining optimization problem" +@time model, varP = PropertyT.sos_problem_primal( + elt, + unit, + wd; + upper_bound = UPPER_BOUND, + augmented = true, +) +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 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/SpNZ_has_T.jl b/scripts/SpNZ_has_T.jl new file mode 100644 index 0000000..7a40d27 --- /dev/null +++ b/scripts/SpNZ_has_T.jl @@ -0,0 +1,83 @@ +using LinearAlgebra +using MKL_jll +BLAS.set_num_threads(4) + +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"] + +const GENUS = 2N + +G = MatrixGroups.SymplecticGroup{GENUS}(Int8) +@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 RG = RG, N = N + G = StarAlgebras.object(RG) + P = PermGroup(perm"(1,2)", Perm(circshift(1:N, -1))) + Σ = Groups.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P) + act = PropertyT.action_by_conjugation(G, Σ) + + wdfl = @time SymbolicWedderburn.WedderburnDecomposition( + Float64, + Σ, + act, + basis(RG), + StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]), + ) +end +@info wd + +Δ = RG(length(S)) - sum(RG(s) for s in S) +elt = Δ^2 +unit = Δ + +@info "defining optimization problem" +@time model, varP = PropertyT.sos_problem_primal( + elt, + unit, + wd; + upper_bound = UPPER_BOUND, + augmented = true, + show_progress = true, +) + +solve_in_loop( + model, + wd, + varP; + logdir = "./log/Sp($N,Z)/r=$HALFRADIUS/Δ²-$(UPPER_BOUND)Δ", + optimizer = scs_optimizer(; + linear_solver = SCS.MKLDirectSolver, + eps = 1e-10, + max_iters = 50_000, + accel = 50, + alpha = 1.95, + ), + data = (elt = elt, unit = unit, halfradius = HALFRADIUS), +) + +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 4c4b394..65fbf57 100644 --- a/scripts/utils.jl +++ b/scripts/utils.jl @@ -39,13 +39,14 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data) nothing end old_lambda = 0.0 + certified = false while status != JuMP.OPTIMAL date = now() log_file = joinpath(logdir, "solver_$date.log") @info "Current logfile is $log_file." isdir(dirname(log_file)) || mkpath(dirname(log_file)) - λ, flag, certified_λ = let + certified, certified_λ = let # logstream = current_logger().logger.stream # v = @ccall setvbuf(logstream.handle::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, 1::Cint, 0::Cint)::Cint # @warn v @@ -58,7 +59,7 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data) solution[:warm] = warm - flag, λ_cert = open(log_file; append = true) do io + certified, λ_cert = open(log_file; append = true) do io with_logger(SimpleLogger(io)) do return PropertyT.certify_solution( data.elt, @@ -70,25 +71,27 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data) end end - solution[:λ], flag, λ_cert + certified, λ_cert end - if flag == true && certified_λ ≥ 0 + if certified == true @info "Certification done with λ = $certified_λ" certified_λ status - return certified_λ + end + + if status == JuMP.OPTIMAL + return certified, certified_λ else rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda)) - @info "Certification failed with λ = $λ" certified_λ rel_change status + @info "Relatie improement for λ" rel_change if rel_change < 1e-9 @info "No progress detected, breaking" certified_λ rel_change status - break + return certified, certified_λ end end - old_lambda = certified_λ end - return status == JuMP.OPTIMAL ? old_lambda : NaN + return certified, old_lambda end