From 750ff5e123cf03b99d8dcc51a8a79eaaeaafd2db Mon Sep 17 00:00:00 2001
From: araujoms <maltusan@gmail.com>
Date: Sun, 29 Dec 2024 14:13:13 +0100
Subject: [PATCH] add multipartite tsirelson_bound

---
 docs/src/api.md  |  1 -
 src/games.jl     | 26 +++++++++++++++++++++++++-
 src/tsirelson.jl | 39 ++++++++++++++++++++-------------------
 test/nonlocal.jl |  4 +++-
 4 files changed, 48 insertions(+), 22 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index 9651c30..6672ca4 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -75,7 +75,6 @@ cglmp
 inn22
 local_bound
 tsirelson_bound
-tsirelson_bound_fc
 seesaw
 tensor_probability
 tensor_collinsgisin
diff --git a/src/games.jl b/src/games.jl
index 681ef31..b82bed7 100644
--- a/src/games.jl
+++ b/src/games.jl
@@ -26,7 +26,7 @@ chsh(d::Integer = 2) = chsh(Float64, d)
 export chsh
 
 """
-    cglmp([T=Float64,] d::Integer)
+    cglmp([T=Float64,] d::Integer = 3)
 
 CGLMP nonlocal game in full probability notation. If `T` is an integer type the game is unnormalized.
 
@@ -76,3 +76,27 @@ function inn22(::Type{T}, n) where {T}
 end
 inn22(n::Integer = 3) = inn22(Int, n)
 export inn22
+
+"""
+    gyni([T=Float64,] n::Integer)
+
+Guess your neighbour's input nonlocal game in full probability notation. If `T` is an integer type the game is unnormalized.
+
+Reference: Mafalda et al., [arXiv:1003.3844](https://arxiv.org/abs/1003.3844).
+"""
+function gyni(::Type{T}, n::Integer = 3) where {T}
+    G = zeros(T, 2 * ones(Int, 2 * n)...)
+
+    normalization = T <: Integer ? 1 : inv(T(2^(n - 1)))
+
+    for x in CartesianIndices(Tuple(2 * ones(Int, n)))
+        if sum(x.I) % 2 == 1
+            a = (x.I[2:n]..., x.I[1])
+            G[a..., x] = normalization
+        end
+    end
+
+    return G
+end
+gyni(n::Integer = 3) = gyni(Float64, n)
+export gyni
diff --git a/src/tsirelson.jl b/src/tsirelson.jl
index 3fc9fc8..e30f5b8 100644
--- a/src/tsirelson.jl
+++ b/src/tsirelson.jl
@@ -1,22 +1,23 @@
 """
-    tsirelson_bound(CG::Matrix, scenario::AbstractVecOrTuple, level::Integer)
+    tsirelson_bound(CG::Matrix, scenario::AbstractVecOrTuple, level)
 
-Upper bounds the Tsirelson bound of a bipartite Bell funcional game `CG`, written in Collins-Gisin notation.
-`scenario` is a vector detailing the number of inputs and outputs, in the order [oa, ob, ia, ib].
-`level` is an integer determining the level of the NPA hierarchy.
+Upper bounds the Tsirelson bound of a multipartite Bell funcional `CG`, written in Collins-Gisin notation.
+`scenario` is a tuple detailing the number of inputs and outputs, in the order (oa, ob, ..., ia, ib, ...).
+`level` is an integer or string determining the level of the NPA hierarchy.
 """
-function tsirelson_bound(CG::Matrix{T}, scenario::AbstractVecOrTuple{<:Integer}, level::Integer; solver = Hypatia.Optimizer{_solver_type(T)}) where {T <: Number}
-    oa, ob, ia, ib = scenario
-    A = QuantumNPA.projector(1, 1:(oa - 1), 1:ia)
-    B = QuantumNPA.projector(2, 1:(ob - 1), 1:ib)
+function tsirelson_bound(CG::Array{T, N}, scenario::AbstractVecOrTuple{<:Integer}, level; solver = Hypatia.Optimizer{_solver_type(T)}) where {T <: Number, N}
+    outs = scenario[1:N]
+    ins = 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
 
-    aind(a, x) = 1 + a + (x - 1) * (oa - 1)
-    bind(b, y) = 1 + b + (y - 1) * (ob - 1)
-
-    bell_functional = sum(CG[aind(a, x), bind(b, y)] * A[a, x] * B[b, y] for a in 1:(oa - 1), b in 1:(ob - 1), x in 1:ia, y in 1:ib)
-    bell_functional += sum(CG[aind(a, x), 1] * A[a, x] for a in 1:(oa - 1), x in 1:ia)
-    bell_functional += sum(CG[1, bind(b, y)] * B[b, y] for b in 1:(ob - 1), y in 1:ib)
-    bell_functional += CG[1, 1] * QuantumNPA.Id
+    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)
+        end
+    end
 
     Q = _npa(_solver_type(T), bell_functional, level; solver)
     return Q
@@ -24,12 +25,12 @@ end
 export tsirelson_bound
 
 """
-    tsirelson_bound_fc(FC::Matrix, level::Integer)
+    tsirelson_bound(FC::Matrix, level::Integer)
 
-Upper bounds the Tsirelson bound of a bipartite Bell funcional game `FC`, written in full correlation notation.
-`level` is an integer determining the level of the NPA hierarchy.
+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(FC::Matrix{T}, level::Integer; solver = Hypatia.Optimizer{_solver_type(T)}) where {T <: Number}
+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)
diff --git a/test/nonlocal.jl b/test/nonlocal.jl
index 200e73f..e579cf2 100644
--- a/test/nonlocal.jl
+++ b/test/nonlocal.jl
@@ -16,7 +16,9 @@
         τ 1  1;
         0 1 -1
     ]
-    @test tsirelson_bound_fc(tilted_chsh_fc, 3) ≈ 3.80128907501837942169727948014219026
+    @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
     for T in [Float64, Double64, Float128, BigFloat]
         @test eltype(chsh(T)) <: T
         @test eltype(cglmp(T)) <: T