From 2a84d04edc55c3b2601d672acdc2c2966da38281 Mon Sep 17 00:00:00 2001
From: kalmarek <kalmar@amu.edu.pl>
Date: Tue, 16 Apr 2019 17:05:34 +0200
Subject: [PATCH] test the correctness of hpc sos computation

---
 src/checksolution.jl    | 12 +++++---
 test/SOS_correctness.jl | 61 +++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl        |  7 +++--
 3 files changed, 73 insertions(+), 7 deletions(-)
 create mode 100644 test/SOS_correctness.jl

diff --git a/src/checksolution.jl b/src/checksolution.jl
index e99fca2..c102191 100644
--- a/src/checksolution.jl
+++ b/src/checksolution.jl
@@ -37,14 +37,18 @@ function compute_SOS(RG::GroupRing, Q::AbstractMatrix{<:Real})
     return GroupRingElem(result, RG)
 end
 
-function compute_SOS_square(RG::GroupRing, Q::AbstractMatrix{<:Real})
-    result = zeros(eltype(Q), maximum(RG.pm));
+function compute_SOS_square(pm::AbstractMatrix{<:Integer}, Q::AbstractMatrix{<:Real})
+    result = zeros(eltype(Q), maximum(pm));
 
     for i in 1:size(Q,2)
-        GroupRings.fmac!(result, view(Q,:,i), view(Q,:,i), RG.pm)
+        GroupRings.fmac!(result, view(Q,:,i), view(Q,:,i), pm)
     end
 
-    return GroupRingElem(result, RG)
+    return result
+end
+
+function compute_SOS_square(RG::GroupRing, Q::AbstractMatrix{<:Real})
+    return GroupRingElem(compute_SOS_square(RG.pm, Q), RG)
 end
 
 function augIdproj(Q::AbstractMatrix{T}) where {T<:Real}
diff --git a/test/SOS_correctness.jl b/test/SOS_correctness.jl
new file mode 100644
index 0000000..f991ee9
--- /dev/null
+++ b/test/SOS_correctness.jl
@@ -0,0 +1,61 @@
+using PropertyT.GroupRings
+
+@testset "Correctness of HPC SOS computation" begin
+
+    function prepare(G_name, λ, S_size)
+        pm = load("$G_name/delta.jld", "pm")
+        P = load("$G_name/$λ/solution.jld", "P")
+        @time Q = real(sqrt(P))
+
+        Δ_coeff = SparseVector(maximum(pm), collect(1:1+S_size), [S_size; ((-1.0) for i in 1:S_size)...])
+
+        Δ²_coeff = GroupRings.GRmul!(spzeros(length(Δ_coeff)), Δ_coeff, Δ_coeff, pm)
+
+        eoi = Δ²_coeff - λ*Δ_coeff
+
+        Q = PropertyT.augIdproj(Q)
+
+        return eoi, pm, Q
+    end
+
+    #########################################################
+    NAME = "SL(3,Z)"
+    eoi, pm, Q = prepare(NAME, 0.1, 3*2*2)
+
+    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
+    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
+
+    @test norm(sos_sqr - sos_hpc, 1) < 3e-12
+    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
+
+    #########################################################
+    NAME = "oSL(3,Z)"
+    eoi, pm, Q = prepare(NAME, 0.27, 3*2*2)
+
+    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
+    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
+
+    @test norm(sos_sqr - sos_hpc, 1) < 3e-12
+    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
+
+    #########################################################
+    NAME = "oSL(4,Z)"
+    eoi, pm, Q = prepare(NAME, 1.3, 4*3*2)
+
+    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
+    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
+
+    @test norm(sos_sqr - sos_hpc, 1) < 2e-10
+    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
+
+    #########################################################
+    NAME = "oSAut(F3)"
+    eoi, pm, Q = prepare(NAME, 0.15, 4*3*2*2)
+
+    @time sos_sqr = PropertyT.compute_SOS_square(pm, Q)
+    @time sos_hpc = PropertyT.compute_SOS(pm, Q)
+
+    @test norm(sos_sqr - sos_hpc, 1) < 6e-11
+    @info "$NAME:\nDifference in l₁-norm between square and hpc sos decompositions:" norm(eoi-sos_sqr,1) norm(eoi-sos_hpc,1) norm(sos_sqr - sos_hpc, 1)
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index ed689b2..773fe0a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,10 +15,11 @@ function Groups.gens(M::MatSpace)
     return S
 end
 
-solver(iters; accel=1) = 
-    with_optimizer(SCS.Optimizer, 
+solver(iters; accel=1) =
+    with_optimizer(SCS.Optimizer,
     linear_solver=SCS.Direct, max_iters=iters,
     acceleration_lookback=accel, eps=1e-10, warm_start=true)
-    
+
 include("1703.09680.jl")
 include("1712.07167.jl")
+include("SOS_correctness.jl")