1
0
mirror of https://github.com/kalmarek/PropertyT.jl.git synced 2024-11-19 15:25:29 +01:00
PropertyT.jl/test/1812.03456.jl

246 lines
7.9 KiB
Julia

@testset "Sq, Adj, Op" begin
function isconstant_on_orbit(v, orb)
isempty(orb) && return true
k = v[first(orb)]
return all(v[o] == k for o in orb)
end
@testset "unit tests" begin
for N in [3,4]
M = MatrixAlgebra(zz, N)
@test PropertyT.EltaryMat(M, 1, 2) isa MatAlgElem
e12 = PropertyT.EltaryMat(M, 1, 2)
@test e12[1,2] == 1
@test inv(e12)[1,2] == -1
S = PropertyT.generating_set(M)
@test e12 S
@test length(PropertyT.generating_set(M)) == 2N*(N-1)
@test all(inv(s) S for s in S)
A = SAut(FreeGroup(N))
@test length(PropertyT.generating_set(A)) == 4N*(N-1)
S = PropertyT.generating_set(A)
@test all(inv(s) S for s in S)
end
@test PropertyT.isopposite(perm"(1,2,3)(4)", perm"(1,4,2)")
@test PropertyT.isadjacent(perm"(1,2,3)", perm"(1,2)(3)")
@test !PropertyT.isopposite(perm"(1,2,3)", perm"(1,2)(3)")
@test !PropertyT.isadjacent(perm"(1,4)", perm"(2,3)(4)")
@test isconstant_on_orbit([1,1,1,2,2], [2,3])
@test !isconstant_on_orbit([1,1,1,2,2], [2,3,4])
end
@testset "Sq, Adj, Op" begin
N = 4
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
Δ = PropertyT.Laplacian(S, 2)
RG = parent(Δ)
autS = WreathProduct(PermGroup(2), PermGroup(N))
orbits = PropertyT.orbit_decomposition(autS, RG.basis)
@test PropertyT.Sq(RG) isa GroupRingElem
sq = PropertyT.Sq(RG)
@test all(isconstant_on_orbit(sq, orb) for orb in orbits)
@test PropertyT.Adj(RG) isa GroupRingElem
adj = PropertyT.Adj(RG)
@test all(isconstant_on_orbit(adj, orb) for orb in orbits)
@test PropertyT.Op(RG) isa GroupRingElem
op = PropertyT.Op(RG)
@test all(isconstant_on_orbit(op, orb) for orb in orbits)
sq, adj, op = PropertyT.SqAdjOp(RG, N)
@test sq == PropertyT.Sq(RG)
@test adj == PropertyT.Adj(RG)
@test op == PropertyT.Op(RG)
e = one(M)
g = PropertyT.EltaryMat(M, 1,2)
h = PropertyT.EltaryMat(M, 1,3)
k = PropertyT.EltaryMat(M, 3,4)
edges = N*(N-1)÷2
@test sq[e] == 20*edges
@test sq[g] == sq[h] == -8
@test sq[g^2] == sq[h^2] == 1
@test sq[g*h] == sq[h*g] == 0
# @test adj[e] == ...
@test adj[g] == adj[h] # == ...
@test adj[g^2] == adj[h^2] == 0
@test adj[g*h] == adj[h*g] # == ...
# @test op[e] == ...
@test op[g] == op[h] # == ...
@test op[g^2] == op[h^2] == 0
@test op[g*h] == op[h*g] == 0
@test op[g*k] == op[k*g] # == ...
@test op[h*k] == op[k*h] == 0
end
end
@testset "1812.03456 examples" begin
function SOS_residual(x::GroupRingElem, Q::Matrix)
RG = parent(x)
@time sos = PropertyT.compute_SOS(RG, Q);
return x - sos
end
function check_positivity(elt, Δ, orbit_data, upper_bound, warm=nothing; with_solver=with_SCS(20_000, accel=10))
SDP_problem, varP = PropertyT.SOS_problem(elt, Δ, orbit_data; upper_bound=upper_bound)
status, warm = PropertyT.solve(SDP_problem, with_solver, warm);
Base.Libc.flush_cstdio()
@info "Optimization status:" status
λ = value(SDP_problem[])
Ps = [value.(P) for P in varP]
Qs = real.(sqrt.(Ps));
Q = PropertyT.reconstruct(Qs, orbit_data);
b = SOS_residual(elt - λ*Δ, Q)
return b, λ, warm
end
@testset "SL(3,Z)" begin
N = 3
halfradius = 2
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
Δ = PropertyT.Laplacian(S, halfradius)
RG = parent(Δ)
orbit_data = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
orbit_data = PropertyT.decimate(orbit_data);
@testset "Sq₃ is SOS" begin
elt = PropertyT.Sq(RG)
UB = 0.05 # 0.105?
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < 2λ/100
end
@testset "Adj₃ is SOS" begin
elt = PropertyT.Adj(RG)
UB = 0.1 # 0.157?
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ
@test 2^2*norm(residual, 1) < λ/100
end
@testset "Op₃ is empty, so can not be certified" begin
elt = PropertyT.Op(RG)
UB = Inf
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) > λ
end
end
@testset "SL(4,Z)" begin
N = 4
halfradius = 2
M = MatrixAlgebra(zz, N)
S = PropertyT.generating_set(M)
Δ = PropertyT.Laplacian(S, halfradius)
RG = parent(Δ)
orbit_data = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
orbit_data = PropertyT.decimate(orbit_data);
@testset "Sq₄ is SOS" begin
elt = PropertyT.Sq(RG)
UB = 0.2 # 0.3172
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < λ/100
end
@testset "Adj₄ is SOS" begin
elt = PropertyT.Adj(RG)
UB = 0.3 # 0.5459?
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < λ/100
end
@testset "we can't cerify that Op₄ SOS" begin
elt = PropertyT.Op(RG)
UB = 2.0
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB,
with_solver=with_SCS(20_000, accel=10, eps=2e-10))
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) > λ # i.e. we can't certify positivity
end
@testset "Adj₄ + Op₄ is SOS" begin
elt = PropertyT.Adj(RG) + PropertyT.Op(RG)
UB = 0.6 # 0.82005
residual, λ, _ = check_positivity(elt, Δ, orbit_data, UB)
Base.Libc.flush_cstdio()
@info "obtained λ and residual" λ norm(residual, 1)
@test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
@test 2^2*norm(residual, 1) < λ/100
end
end
# @testset "Adj₄ + 100 Op₄ ∈ ISAut(F₄) is SOS" begin
# N = 4
# halfradius = 2
# M = SAut(FreeGroup(N))
# S = PropertyT.generating_set(M)
# Δ = PropertyT.Laplacian(S, halfradius)
# RG = parent(Δ)
# orbit_data = PropertyT.OrbitData(RG, WreathProduct(PermGroup(2), PermGroup(N)))
# orbit_data = PropertyT.decimate(orbit_data);
#
# @time elt = PropertyT.Adj(RG) + 100PropertyT.Op(RG)
# UB = 0.05
#
# warm = nothing
#
# residual, λ, warm = check_positivity(elt, Δ, orbit_data, UB, warm, with_solver=with_SCS(20_000, accel=10))
# @info "obtained λ and residual" λ norm(residual, 1)
#
# @test 2^2*norm(residual, 1) < λ # i.e. we can certify positivity
# @test 2^2*norm(residual, 1) < λ/100
# end
end