add script for SpN_Adj.jl

scripts/SpN_Adj.jl

@@ -0,0 +1,80 @@
+using LinearAlgebra
+BLAS.set_num_threads(8)
+
+ENV["OMP_NUM_THREADS"] = 1
+
+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)
+
+RG, S, sizes =
+    @time PropertyT.group_algebra(G, halfradius=HALFRADIUS, twisted=true)
+
+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)
+    # Σ = P
+    act = PropertyT.action_by_conjugation(G, Σ)
+    @info "Computing WedderburnDecomposition"
+
+    wdfl = @time SymbolicWedderburn.WedderburnDecomposition(
+        Float64,
+        Σ,
+        act,
+        basis(RG),
+        StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
+    )
+end
+
+Δ = RG(length(S)) - sum(RG(s) for s in S)
+Δs = PropertyT.laplacians(
+    RG,
+    S,
+    x -> (gx = PropertyT.grading(x); Set([gx, -gx])),
+)
+
+# elt = Δ^2
+elt = PropertyT.Adj(Δs, :C₂)
+unit = Δ
+
+@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/r=$HALFRADIUS/Sp($N,Z)/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)
+)
+

scripts/argparse.jl

@@ -0,0 +1,19 @@
+using ArgParse
+
+args_settings = ArgParseSettings()
+@add_arg_table! args_settings begin
+    "-N"
+    help = "the degree/genus/etc. parameter for a group"
+    arg_type = Int
+    default = 3
+    "--halfradius", "-R"
+    help = "the halfradius on which perform the sum of squares decomposition"
+    arg_type = Int
+    default = 2
+    "--upper_bound", "-u"
+    help = "set upper bound for the optimization problem to speed-up the convergence"
+    arg_type = Float64
+    default = Inf
+end
+
+parsed_args = parse_args(ARGS, args_settings)

scripts/utils.jl

@@ -0,0 +1,85 @@
+using Dates
+using Serialization
+using Logging
+
+import JuMP
+
+function get_solution(model)
+    λ = JuMP.value(model[:λ])
+    Q = real.(sqrt(JuMP.value.(model[:P])))
+    solution = Dict(:λ => λ, :Q => Q)
+    return solution
+end
+
+function get_solution(model, wd, varP; logdir)
+    λ = 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
+
+function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
+    @info "logging to $logdir"
+    status = JuMP.UNKNOWN_RESULT_STATUS
+    warm = try
+        solution = deserialize(joinpath(logdir, "solution.sjl"))
+        warm = solution[:warm]
+        @info "trying to warm-start model with λ=$(solution[:λ])..."
+        warm
+    catch
+        nothing
+    end
+    old_lambda = 0.0
+    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
+            # 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
+
+            serialize(joinpath(logdir, "solution_$date.sjl"), solution)
+            serialize(joinpath(logdir, "solution.sjl"), solution)
+
+            flag, λ_cert = open(log_file, append=true) do io
+                with_logger(SimpleLogger(io)) do
+                    PropertyT.certify_solution(
+                        data.elt,
+                        data.unit,
+                        solution[:λ],
+                        solution[:Q],
+                        halfradius=data.halfradius,
+                    )
+                end
+            end
+
+            solution[:λ], flag, λ_cert
+        end
+
+        if flag == true && certified_λ ≥ 0
+            @info "Certification done with λ = $certified_λ"
+            return certified_λ
+        else
+            rel_change = abs(certified_λ - old_lambda) / (abs(certified_λ) + abs(old_lambda))
+            @info "Certification failed with λ = $λ" certified_λ rel_change
+        end
+
+        old_lambda = certified_λ
+
+        if rel_change < 1e-9
+            @info "No progress detected, breaking"
+            break
+        end
+    end
+
+    return status == JuMP.OPTIMAL ? old_lambda : NaN
+end