From a596fd78f50b002734b7bcab770129385c4d4768 Mon Sep 17 00:00:00 2001
From: Marek Kaluba <kalmar@amu.edu.pl>
Date: Tue, 27 Feb 2024 18:14:59 +0100
Subject: [PATCH] add scripts/PRA_has_T.jl

---
 scripts/PRA_has_T.jl | 162 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 162 insertions(+)
 create mode 100644 scripts/PRA_has_T.jl

diff --git a/scripts/PRA_has_T.jl b/scripts/PRA_has_T.jl
new file mode 100644
index 0000000..68a41e6
--- /dev/null
+++ b/scripts/PRA_has_T.jl
@@ -0,0 +1,162 @@
+using LinearAlgebra
+BLAS.set_num_threads(4)
+ENV["OMP_NUM_THREADS"] = 4
+include(joinpath(@__DIR__, "../test/optimizers.jl"))
+using SCS_MKL_jll
+
+using Groups
+import Groups.MatrixGroups
+
+using PropertyT
+
+import PropertyT.SW as SW
+using PropertyT.PG
+using PropertyT.SA
+
+include(joinpath(@__DIR__, "argparse.jl"))
+
+const N = parsed_args["N"]
+const HALFRADIUS = parsed_args["halfradius"]
+const UPPER_BOUND = parsed_args["upper_bound"]
+
+# fixes/hacks
+import Groups.KnuthBendix
+KnuthBendix.ordering(o::KnuthBendix.WordOrdering) = o
+function KnuthBendix.rewrite!(
+    u::KnuthBendix.AbstractWord,
+    w::KnuthBendix.AbstractWord,
+    o::KnuthBendix.WordOrdering,
+)
+    return KnuthBendix.rewrite!(u, w, KnuthBendix.alphabet(o))
+end
+
+struct Letter{T} <: Groups.GSymbol # letter of an Alphabet
+    elt::T
+end
+
+Base.show(io::IO, tt::Letter) = show(io, tt.elt)
+Base.inv(tt::Letter) = Letter(inv(tt.elt))
+Base.:(==)(tt::Letter, ss::Letter) = tt.elt == ss.elt
+Base.hash(tt::Letter, h::UInt) = hash(tt.elt, hash(Letter, h))
+
+Base.Base.@propagate_inbounds function Groups.evaluate!(
+    v::Tuple{Vararg{T,N}},
+    tt::Letter,
+    tmp = one(first(v)),
+) where {T,N}
+    return Groups.evaluate!(v, tt.elt, tmp)
+end
+
+function PropertyT._conj(tt::Letter, g)
+    G = parent(tt.elt)
+    A = alphabet(G)
+
+    w = [A[PropertyT._conj(A[l], g)] for l in word(tt.elt)]
+    return Letter(G(w))
+end
+
+G = let G = SpecialAutomorphismGroup(FreeGroup(N + 1))
+    A = alphabet(G)
+    lambdas = [Groups.λ(1, i) for i in 2:N+1]
+    append!(lambdas, [Groups.λ(i, 1) for i in 2:N+1])
+    rhos = [Groups.ϱ(1, i) for i in 2:N+1]
+    append!(rhos, [Groups.ϱ(i, 1) for i in 2:N+1])
+
+    _alph = eltype(G)[]
+
+    for i in 2:N+1
+        for j in 2:N+1
+            i == j && continue
+            g = G([A[Groups.ϱ(1, i)], A[Groups.ϱ(j, 1)]])
+            h = G([A[Groups.λ(1, i)], A[Groups.λ(j, 1)]])
+            push!(_alph, g, h)
+        end
+    end
+
+    alph = Letter.(_alph)
+    AutomorphismGroup(
+        FreeGroup(N + 1),
+        alph,
+        KnuthBendix.LenLex(Groups.Alphabet(alph)),
+        Groups.domain(one(G)),
+    )
+end
+# @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"(2,3)", Perm([1; 1 .+ circshift(1:N, -1)]))
+    Σ = Groups.Constructions.WreathProduct(PermGroup(perm"(1,2)"), P)
+    act = PropertyT.action_by_conjugation(G, P)
+
+    wdfl = @time SW.WedderburnDecomposition(
+        Float64,
+        P,
+        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,
+    show_progress = true,
+)
+
+let status = JuMP.OPTIMIZE_NOT_CALLED, warm = warm, eps = 1e-10
+    certified, λ = false, 0.0
+    while status ≠ JuMP.OPTIMAL
+        @time status, warm = PropertyT.solve(
+            model,
+            scs_optimizer(;
+                linear_solver = SCS.MKLDirectSolver,
+                eps = eps,
+                max_iters = N * 10_000,
+                accel = 50,
+                alpha = 1.95,
+            ),
+            warm,
+        )
+
+        @info "reconstructing the solution"
+        Q = @time let wd = wd, Ps = [JuMP.value.(P) for P in varP], eps = 1e-10
+            PropertyT.__droptol!.(Ps, 100eps)
+            Qs = real.(sqrt.(Ps))
+            PropertyT.__droptol!.(Qs, eps)
+
+            PropertyT.reconstruct(Qs, wd)
+        end
+
+        @info "certifying the solution"
+        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
+end