add/update scripts for SLNZ/SpNZ

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

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

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