Add local_model and adapt nonlocality_threshold examples

Author: sebastiendesignolle

examples/CHSH.jl

@@ -30,8 +30,8 @@ println()
 println("Lower bound")
 println(lower_bound) # 0.7071
 println("Local model")
-display(local_model.x)
-println(local_model.x == sum(local_model.weights[i] * local_model.atoms[i] for i in 1:length(local_model))) # true
+display(local_model)
+# display(sum(weight * atom for (weight, atom) in local_model))
 
 println()
 

examples/GHZ.jl

@@ -19,7 +19,7 @@ measurements = [[(σ[1] - v[i, 1] * σ[2] - v[i, 2] * σ[3] - v[i, 3] * σ[4]) /
 p = tensor_correlation(rho, measurements, N; marg = false) # the marginals vanish in this cas
 
 # Frank-Wolfe
-@time lower_bound, upper_bound, local_model, bell_inequality = nonlocality_threshold(p; sym = false)
+lower_bound, upper_bound, local_model, bell_inequality = nonlocality_threshold(p)
 
 println("Correlation tensor")
 display(p[:, :, 1]) # only printing part of the tensor
@@ -29,8 +29,8 @@ println()
 println("Lower bound")
 println(lower_bound) # 0.49315
 println("Local model")
-display(local_model.x[:, :, 1])
-println(local_model.x == sum(local_model.weights[i] * local_model.atoms[i] for i in 1:length(local_model))) # true
+# display(local_model)
+display(sum(weight * atom for (weight, atom) in local_model)[:, :, 1])
 
 println()
 

src/fw_methods.jl

@@ -627,7 +627,7 @@ end
 function _muladd_memory_mode(as, d::AbstractArray{T}, a, v) where {T <: Number}
     idx_a = _unsafe_find_atom(as, a)
     idx_v = _unsafe_find_atom(as, v)
-    d[1] = typemax(T)
+    @inbounds d[1] = typemax(T)
     if idx_v > idx_a
         @inbounds d[2] = ((as.dots_x[idx_a] + as.dots_b[idx_a]) - (as.dots_x[idx_v] + as.dots_b[idx_v])) / (as.dots_A[idx_a][idx_a] + as.dots_A[idx_v][idx_v] - 2as.dots_A[idx_v][idx_a])
     else

src/nonlocality_threshold.jl

@@ -4,37 +4,41 @@
 Compute the nonlocality threshold of the probability/correlation tensor `p`.
 
 Returns:
- - `lower_bound`: a lower bound on the nonlocality threshold under the measurements provided in input,
+ - `lower_bound`: a (exact up to analyticity step) lower bound on the nonlocality threshold of `p`,
  - `upper_bound`: a (heuristic) upper bound on the nonlocality threshold of `p`
- - `local_model`: a decomposition of the tensor `p` with visibility `lower_bound`,
+ - `local_model`: a decomposition of the tensor `p` with visibility `lower_bound` (up to a distance `2√epsilon`),
  - `bell_inequality`: a (heuristic) Bell inequality corresponding to `upper_bound`.
 
 Optional arguments:
+ - ``
  - `precision`: number of digits of `lower_bound`, 4 by default,
  - for the other optional arguments, see `bell_frank_wolfe`.
 """
 function nonlocality_threshold(
         p::Array{T, N},
         lower_bound = zero(T),
         upper_bound = one(T);
+        upper::Bool = true,
         precision = 4,
         prob::Bool = false,
         marg::Bool = false,
-        v0 = one(T),
         epsilon = Base.rtoldefault(T),
         o = nothing,
         sym = nothing,
         deflate = identity,
         inflate = identity,
         verbose = 0,
         active_set = nothing,
-        shortcut = 10,
+        shortcut = 4,
         kwargs...,
     ) where {T <: Number, N}
+    expand_permutedims = sym === nothing
     _, _, _, o, sym, deflate, inflate = _bfw_init(p, 0, prob, marg, o, sym, deflate, inflate, verbose > 0)
-    local_model = nothing
+    expand_permutedims &= sym
+    v0 = upper ? upper_bound : lower_bound
+    ass = nothing
     bell_inequality = nothing
-    while log10(upper_bound - lower_bound) > -precision
+    while round(log10(upper_bound - lower_bound); digits = 4) > -precision
         res = bell_frank_wolfe(p; v0, epsilon, prob, marg, o, sym, deflate, inflate, verbose, verbose_init = false, active_set, shortcut, mode_last = -1, kwargs...)
         x, ds, primal, dual_gap, active_set, M, β = res
         if dual_gap ≥ primal && primal > 10epsilon && dual_gap > 10epsilon
@@ -55,13 +59,13 @@ function nonlocality_threshold(
             end
         else
             lower_bound = v0
-            local_model = ActiveSetStorage(active_set) # in-place mess otherwise
+            ass = ActiveSetStorage(active_set)
             if upper_bound < lower_bound
                 upper_bound = round(v0 + 2 * 10.0^(-precision); digits = precision)
             end
             v0 = (lower_bound + upper_bound) / 2
         end
     end
-    return lower_bound, upper_bound, load_active_set(local_model, T; deflate), bell_inequality
+    return lower_bound, upper_bound, local_model(ass; deflate, expand_permutedims), bell_inequality
 end
 export nonlocality_threshold

src/types.jl

@@ -577,8 +577,28 @@ function load_active_set(
     end
     weights = T2.(ass.weights)
     weights /= sum(weights)
-    res = FrankWolfe.ActiveSetQuadraticProductCaching([(weights[i], deflate(atoms[i])) for i in eachindex(ass.weights)], I, deflate(p))
-    return res
+    return FrankWolfe.ActiveSetQuadraticProductCaching([(weights[i], deflate(atoms[i])) for i in eachindex(ass.weights)], I, deflate(p))
+end
+
+function local_model(
+        ass::ActiveSetStorage{T, N, HasMarginals};
+        marg = HasMarginals,
+        deflate = identity,
+        expand_permutedims = false,
+    ) where {T <: Number, N, HasMarginals}
+    as = load_active_set(ass, T; marg, deflate = expand_permutedims ? identity : deflate)
+    if expand_permutedims
+        weights = repeat(as.weights / factorial(N); inner = factorial(N))
+        atoms = Vector{BellCorrelationsDS{T, N, marg}}(undef, factorial(N) * length(as))
+        i = 0
+        for atom in as.atoms, per in permutations(1:N)
+            i += 1
+            atoms[i] = BellCorrelationsDS(atom.ax[per], atom.lmo)
+        end
+        return collect(zip(weights, atoms))
+    else
+        return collect(zip(as.weights, as.atoms))
+    end
 end
 
 # for multi-outcome scenarios
@@ -630,6 +650,25 @@ function load_active_set(
     end
     weights = T2.(ass.weights)
     weights /= sum(weights)
-    res = FrankWolfe.ActiveSetQuadraticProductCaching([(weights[i], deflate(atoms[i])) for i in eachindex(ass.weights)], I, deflate(p))
-    return res
+    return FrankWolfe.ActiveSetQuadraticProductCaching([(weights[i], deflate(atoms[i])) for i in eachindex(ass.weights)], I, deflate(p))
+end
+
+function local_model(
+        ass::ActiveSetStorageMulti{T, N};
+        deflate = identity,
+        expand_permutedims = false,
+    ) where {T <: Number, N}
+    as = load_active_set(ass, T; deflate = expand_permutedims ? identity : deflate)
+    if expand_permutedims
+        weights = repeat(as.weights / factorial(N); inner = factorial(N))
+        atoms = Vector{BellProbabilitiesDS{T, N, marg}}(undef, factorial(N) * length(as))
+        i = 0
+        for atom in as.atoms, per in permutations(1:N)
+            i += 1
+            atoms[i] = BellCorrelationsDS(atom.ax[per], atom.lmo)
+        end
+        return collect(zip(weights, atoms))
+    else
+        return collect(zip(as.weights, as.atoms))
+    end
 end