diff --git a/scripts/G₂_Adj.jl b/scripts/G₂_Adj.jl index efc6770..339d99c 100644 --- a/scripts/G₂_Adj.jl +++ b/scripts/G₂_Adj.jl @@ -1,29 +1,34 @@ using LinearAlgebra -BLAS.set_num_threads(1) -ENV["OMP_NUM_THREADS"] = 4 +BLAS.set_num_threads(8) -using MKL_jll -include(joinpath(@__DIR__, "../test/optimizers.jl")) +ENV["OMP_NUM_THREADS"] = 4 using Groups import Groups.MatrixGroups + +include(joinpath(@__DIR__, "../test/optimizers.jl")) using PropertyT -using SymbolicWedderburn -using SymbolicWedderburn.StarAlgebras -using PermutationGroups +using PropertyT.SymbolicWedderburn +using PropertyT.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 -const UPPER_BOUND = Inf - +@info "computing group algebra structure" 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 - -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 +@info wd 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 - @time status, warm = PropertyT.solve( - model, - scs_optimizer(; - linear_solver = SCS.MKLDirectSolver, - eps = 1e-10, - max_iters = 50_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 - -# Δ² - 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 -# ] +solve_in_loop( + model, + wd, + varP; + logdir = "./log/G2/r=$HALFRADIUS/Adj-InfΔ", + optimizer = scs_optimizer(; + eps = 1e-10, + max_iters = 50_000, + accel = 50, + alpha = 1.95, + ), + data = (elt = elt, unit = unit, halfradius = HALFRADIUS), +) diff --git a/scripts/G₂_has_T.jl b/scripts/G₂_has_T.jl new file mode 100644 index 0000000..e4e6c45 --- /dev/null +++ b/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 diff --git a/scripts/utils.jl b/scripts/utils.jl index 1357690..9a3f5a5 100644 --- a/scripts/utils.jl +++ b/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