rewrite scripts for G2

scripts/G₂_Adj.jl

@@ -1,29 +1,34 @@
 using LinearAlgebra
-BLAS.set_num_threads(1)
+BLAS.set_num_threads(8)
-ENV["OMP_NUM_THREADS"] = 4
 
-using MKL_jll
+ENV["OMP_NUM_THREADS"] = 4
-include(joinpath(@__DIR__, "../test/optimizers.jl"))
 
 using Groups
 import Groups.MatrixGroups
+
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
 using PropertyT
 
-using SymbolicWedderburn
+using PropertyT.SymbolicWedderburn
-using SymbolicWedderburn.StarAlgebras
+using PropertyT.PermutationGroups
-using PermutationGroups
+using PropertyT.StarAlgebras
 
-include(joinpath(@__DIR__, "G₂_gens.jl"))
+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"]
+
+include(joinpath(@__DIR__, "./G₂_gens.jl"))
 
 G, roots, Weyl = G₂_roots_weyl()
+@info "Running Adj² - λ·Δ sum of squares decomposition for G₂"
 
-const HALFRADIUS = 2
+@info "computing group algebra structure"
-const UPPER_BOUND = Inf
-
 RG, S, sizes = @time PropertyT.group_algebra(G, halfradius = HALFRADIUS)
 
-Δ = RG(length(S)) - sum(RG(s) for s in S)
+@info "computing WedderburnDecomposition"
-
 wd = let Σ = Weyl, RG = RG
     act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
         Dict(g => PermutationGroups.perm(g) for g in Σ),
@@ -38,57 +43,7 @@ wd = let Σ = Weyl, RG = RG
         semisimple = false,
     )
 end
-
+@info wd
-elt = Δ^2
-unit = Δ
-
-@time model, varP = PropertyT.sos_problem_primal(
-    elt,
-    unit,
-    wd;
-    upper_bound = UPPER_BOUND,
-    augmented = true,
-    show_progress = true,
-)
-warm = nothing
-
-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 begin
-        wd = wd
-        Ps = [JuMP.value.(P) for P in varP]
-        if any(any(isnan, P) for P in Ps)
-            throw("solver was probably interrupted, no valid solution available")
-        end
-        Qs = real.(sqrt.(Ps))
-        PropertyT.reconstruct(Qs, wd)
-    end
-    P = Q' * Q
-
-    @info "certifying the solution"
-    @time certified, λ = PropertyT.certify_solution(
-        elt,
-        unit,
-        JuMP.objective_value(model),
-        Q;
-        halfradius = HALFRADIUS,
-        augmented = true,
-    )
-end
-
-### grading below
 
 function desubscriptify(symbol::Symbol)
     digits = [
@@ -107,6 +62,7 @@ function PropertyT.grading(g::MatrixGroups.MatrixElt, roots = roots)
     return roots[id]
 end
 
+Δ = RG(length(S)) - sum(RG(s) for s in S)
 Δs = PropertyT.laplacians(
     RG,
     S,
@@ -114,7 +70,7 @@ end
 )
 
 elt = PropertyT.Adj(Δs)
-elt == Δ^2 - PropertyT.Sq(Δs)
+@assert elt == Δ^2 - PropertyT.Sq(Δs)
 unit = Δ
 
 @time model, varP = PropertyT.sos_problem_primal(
@@ -123,57 +79,21 @@ unit = Δ
     wd;
     upper_bound = UPPER_BOUND,
     augmented = true,
+    show_progress = true,
 )
 
 warm = nothing
 
-begin
+solve_in_loop(
-    @time status, warm = PropertyT.solve(
     model,
-        scs_optimizer(;
+    wd,
-            linear_solver = SCS.MKLDirectSolver,
+    varP;
+    logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ",
+    optimizer = scs_optimizer(;
         eps = 1e-10,
         max_iters = 50_000,
         accel = 50,
         alpha = 1.95,
     ),
-        warm,
+    data = (elt = elt, unit = unit, halfradius = HALFRADIUS),
-    )
+)
-
-    @info "reconstructing the solution"
-    Q = @time begin
-        wd = wd
-        Ps = [JuMP.value.(P) for P in varP]
-        if any(any(isnan, P) for P in Ps)
-            throw("solver was probably interrupted, no valid solution available")
-        end
-        Qs = real.(sqrt.(Ps))
-        PropertyT.reconstruct(Qs, wd)
-    end
-    P = Q' * Q
-
-    @info "certifying the solution"
-    @time certified, λ = PropertyT.certify_solution(
-        elt,
-        unit,
-        JuMP.objective_value(model),
-        Q;
-        halfradius = HALFRADIUS,
-        augmented = true,
-    )
-end
-
-# Δ² - 1 / 1 · Sq → -0.8818044647162608
-# Δ² - 2 / 3 · Sq → -0.1031738
-# Δ² - 1 / 2 · Sq → 0.228296213895906
-# Δ² - 1 / 3 · Sq → 0.520
-# Δ² - 0 / 1 · Sq → 0.9676851592000731
-# Sq → 0.333423
-
-# vals = [
-#     1.0 -0.8818
-#     2/3 -0.1032
-#     1/2  0.2282
-#     1/3  0.520
-#     0    0.9677
-# ]

scripts/G₂_has_T.jl

@@ -0,0 +1,97 @@
+using LinearAlgebra
+BLAS.set_num_threads(8)
+
+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"]
+
+include(joinpath(@__DIR__, "./G₂_gens.jl"))
+
+G, roots, Weyl = G₂_roots_weyl()
+@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 Σ = Weyl, RG = RG
+    act = PropertyT.AlphabetPermutation{eltype(Σ),Int64}(
+        Dict(g => PermutationGroups.perm(g) for g in Σ),
+    )
+
+    @time SymbolicWedderburn.WedderburnDecomposition(
+        Float64,
+        Σ,
+        act,
+        basis(RG),
+        StarAlgebras.Basis{UInt16}(@view basis(RG)[1:sizes[HALFRADIUS]]),
+        semisimple = false,
+    )
+end
+@info wd
+
+Δ = RG(length(S)) - sum(RG(s) for s in S)
+elt = Δ^2
+unit = Δ
+
+@time model, varP = PropertyT.sos_problem_primal(
+    elt,
+    unit,
+    wd;
+    upper_bound = UPPER_BOUND,
+    augmented = true,
+    show_progress = false,
+)
+warm = nothing
+status = JuMP.OPTIMIZE_NOT_CALLED
+
+while status ≠ JuMP.OPTIMAL
+    @time status, warm = PropertyT.solve(
+        model,
+        scs_optimizer(;
+            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/utils.jl

@@ -68,19 +68,19 @@ function solve_in_loop(model::JuMP.Model, args...; logdir, optimizer, data)
         end
 
         if flag == true && certified_λ ≥ 0
-            @info "Certification done with λ = $certified_λ"
+            @info "Certification done with λ = $certified_λ" certified_λ rel_change status
             return certified_λ
         else
             rel_change =
                 abs(certified_λ - old_lambda) /
                 (abs(certified_λ) + abs(old_lambda))
-            @info "Certification failed with λ = $λ" certified_λ rel_change
+            @info "Certification failed with λ = $λ" certified_λ rel_change status
         end
 
         old_lambda = certified_λ
 
         if rel_change < 1e-9
-            @info "No progress detected, breaking"
+            @info "No progress detected, breaking" certified_λ rel_change status
             break
         end
     end