added multipartite tsirelson_bound for full correlation

src/games.jl

@@ -8,15 +8,11 @@ Reference: Buhrman and Massar, [arXiv:quant-ph/0409066](https://arxiv.org/abs/qu
 function chsh(::Type{T}, d::Integer = 2) where {T}
     G = zeros(T, d, d, d, d)
 
-    if T <: Integer
-        element = 1
-    else
-        element = inv(T(d^2))
-    end
+    normalization = T <: Integer ? 1 : inv(T(d^2))
 
     for a in 0:(d - 1), b in 0:(d - 1), x in 0:(d - 1), y in 0:(d - 1)
         if mod(a + b + x * y, d) == 0
-            G[a + 1, b + 1, x + 1, y + 1] = element
+            G[a + 1, b + 1, x + 1, y + 1] = normalization
         end
     end
 
@@ -35,11 +31,7 @@ References: [arXiv:quant-ph/0106024](https://arxiv.org/abs/quant-ph/0106024) for
 function cglmp(::Type{T}, d::Integer = 3) where {T}
     G = zeros(T, d, d, 2, 2)
 
-    if T <: Integer
-        normalization = 1
-    else
-        normalization = inv(T(4 * (d - 1)))
-    end
+    normalization = T <: Integer ? 1 : inv(T(4 * (d - 1)))
 
     for a in 0:(d - 1), b in 0:(d - 1), x in 0:1, y in 0:1, k in 0:(d - 2)
         if mod(a - b, d) == mod((-1)^mod(x + y, 2) * k + x * y, d)

src/tsirelson.jl

@@ -8,42 +8,38 @@ Upper bounds the Tsirelson bound of a multipartite Bell funcional `CG`, written
 function tsirelson_bound(CG::Array{T, N}, scenario::AbstractVecOrTuple{<:Integer}, level; solver = Hypatia.Optimizer{_solver_type(T)}) where {T <: Number, N}
     outs = Tuple(scenario[1:N])
     ins = Tuple(scenario[(N + 1):(2 * N)])
-    Π = [[QuantumNPA.projector(n, 1:(outs[n] - 1), 1:ins[n]) [QuantumNPA.Id for _ in 1:(outs[n] - 1)]] for n in 1:N]
-    cgindex(a, x) = (x .!= (ins .+ 1)) .* (a .+ (x .- 1) .* (outs .- 1)) .+ 1
+    Π = [[[QuantumNPA.Id for _ in 1:(outs[n] - 1)] QuantumNPA.projector(n, 1:(outs[n] - 1), 1:ins[n])] for n in 1:N]
+    cgindex(a, x) = (x .!= 1) .* (a .+ (x .- 2) .* (outs .- 1)) .+ 1
 
     bell_functional = QuantumNPA.Polynomial()
     for x in CartesianIndices(ins .+ 1)
-        cgiterators = map((a, b, c) -> a == b ? (1:1) : (1:c), x.I, ins .+ 1, outs .- 1)
-        for a in Iterators.product(cgiterators...)
-            bell_functional += CG[cgindex(a, x.I)...] * prod(Π[n][a[n], x.I[n]] for n in 1:N)
+        for a in CartesianIndices((x.I .!= 1) .* (outs .- 2) .+ 1)
+            bell_functional += CG[cgindex(a.I, x.I)...] * prod(Π[n][a.I[n], x.I[n]] for n in 1:N)
         end
     end
 
     Q = _npa(_solver_type(T), bell_functional, level; solver)
     return Q
 end
 export tsirelson_bound
-
 """
     tsirelson_bound(FC::Matrix, level::Integer)
 
 Upper bounds the Tsirelson bound of a bipartite Bell funcional `FC`, written in full correlation notation.
 `level` is an integer or string determining the level of the NPA hierarchy.
 """
-function tsirelson_bound(FC::Matrix{T}, level; solver = Hypatia.Optimizer{_solver_type(T)}) where {T <: Number}
-    ia, ib = size(FC) .- 1
-    A = QuantumNPA.dichotomic(1, 1:ia)
-    B = QuantumNPA.dichotomic(2, 1:ib)
+function tsirelson_bound(FC::Array{T, N}, level; solver = Hypatia.Optimizer{_solver_type(T)}) where {T <: Number, N}
+    ins = size(FC) .- 1
+    O = [[QuantumNPA.Id; QuantumNPA.dichotomic(n, 1:ins[n])] for n in 1:N]
 
-    bell_functional = sum(FC[x + 1, y + 1] * A[x] * B[y] for x in 1:ia, y in 1:ib)
-    bell_functional += sum(FC[x + 1, 1] * A[x] for x in 1:ia)
-    bell_functional += sum(FC[1, y + 1] * B[y] for y in 1:ib)
-    bell_functional += FC[1, 1] * QuantumNPA.Id
+    bell_functional = QuantumNPA.Polynomial()
+    for x in CartesianIndices(ins .+ 1)
+        bell_functional += FC[x] * prod(O[n][x[n]] for n in 1:N)
+    end
 
     Q = _npa(_solver_type(T), bell_functional, level; solver)
     return Q
 end
-export tsirelson_bound_fc
 
 function _npa(::Type{T}, functional, level; solver) where {T <: AbstractFloat}
     model = JuMP.GenericModel{T}()

test/nonlocal.jl

@@ -9,7 +9,7 @@
     cglmp_cg = tensor_collinsgisin(cglmp())
     @test seesaw(cglmp_cg, (3, 3, 2, 2), 3)[1] ≈ (15 + sqrt(33)) / 24
     @test seesaw(inn22(), (2, 2, 3, 3), 2)[1] ≈ 1.25
-    @test tsirelson_bound(cglmp_cg, (3, 3, 2, 2), 2) ≈ (15 + sqrt(33)) / 24 rtol = 1.0e-7
+    @test tsirelson_bound(cglmp_cg, (3, 3, 2, 2), "1 + A B") ≈ (15 + sqrt(33)) / 24 rtol = 1.0e-7
     τ = Double64(9) / 10
     tilted_chsh_fc = [
         0 τ  0;
@@ -19,6 +19,17 @@
     @test tsirelson_bound(tilted_chsh_fc, 3) ≈ 3.80128907501837942169727948014219026
     gyni_cg = tensor_collinsgisin(gyni())
     @test tsirelson_bound(gyni_cg, 2*ones(Int,6), 3) ≈ 0.25 rtol = 1e-6 #for some reason CI gives a different result
+    Śliwa18 = [0   0   0
+               1   1   0
+               1   1   0;;;
+               0  -2   0
+               1   0   1
+               1   0  -1;;;
+               0   0   2
+               0   1  -1
+               0  -1  -1]
+    @test tsirelson_bound(Śliwa18, 2) ≈ 2*(7-sqrt(17)) rtol = 1e-7
+
     for T in [Float64, Double64, Float128, BigFloat]
         @test eltype(chsh(T)) <: T
         @test eltype(cglmp(T)) <: T