Fix nasty inplace bug and factorise out _bfw_init

Author: sebastiendesignolle

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
-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

src/bell_frank_wolfe.jl

@@ -67,47 +67,14 @@ function bell_frank_wolfe(
         kwargs...,
     ) where {T <: Number, N}
     Random.seed!(seed)
-    if !prob
-        LMO = BellCorrelationsLMO
-        DS = BellCorrelationsDS
-        m = collect(size(p))
-        if o === nothing
-            o = zeros(T, size(p))
-            o[end] = marg
-        end
-        reynolds = reynolds_permutedims
-        build_deflate_inflate = build_deflate_inflate_permutedims
-    else
-        LMO = BellProbabilitiesLMO
-        DS = BellProbabilitiesDS
-        m = collect(size(p)[(N ÷ 2 + 1):end])
-        if o === nothing
-            o = ones(T, size(p)) / prod(size(p)[1:(N ÷ 2)])
-        end
-        reynolds = reynolds_permutelastdims
-        build_deflate_inflate = build_deflate_inflate_permutelastdims
-    end
-    # symmetry detection
-    if sym === nothing
-        if all(diff(m) .== 0) && p ≈ reynolds(p) && (v0 == 1 || o ≈ reynolds(o))
-            deflate, inflate = build_deflate_inflate(p)
-            sym = true
-        else
-            sym = false
-        end
-    end
+    LMO, DS, m, o, sym, deflate, inflate = _bfw_init(p, v0, prob, marg, o, sym, deflate, inflate, verbose)
     if verbose > 0
         println("Visibility: ", v0)
     end
     # choosing the point on the line between o and p according to the visibility v0
-    ro = deflate(o)
-    rp = deflate(p)
+    ro = deflate(copy(o))
+    rp = deflate(copy(p))
     vp = v0 * rp + (one(T) - v0) * ro
-    if verbose > 1
-        println("   #Inputs: ", all(diff(m) .== 0) ? m[end] - (marg && !prob) : m .- (marg && !prob))
-        println(" Symmetric: ", sym)
-        println(" Dimension: ", length(vp))
-    end
     # create the LMO
     if sym
         lmo = FrankWolfe.SubspaceLMO(LMO(p, vp; mode, nb, marg), deflate, inflate)

src/nonlocality_threshold.jl

@@ -24,49 +24,17 @@ function nonlocality_threshold(
         sym = nothing,
         deflate = identity,
         inflate = identity,
+        verbose = 0,
         kwargs...,
     ) where {T <: Number, N}
-    if !prob
-        m = collect(size(p))
-        if o === nothing
-            o = zeros(T, size(p))
-            o[end] = marg
-        end
-        reynolds = reynolds_permutedims
-        build_deflate_inflate = build_deflate_inflate_permutedims
-    else
-        m = collect(size(p)[(N ÷ 2 + 1):end])
-        if o === nothing
-            o = ones(T, size(p)) / prod(size(p)[1:(N ÷ 2)])
-        end
-        reynolds = reynolds_permutelastdims
-        build_deflate_inflate = build_deflate_inflate_permutelastdims
-    end
-    # symmetry detection
-    if sym === nothing
-        if all(diff(m) .== 0) && p ≈ reynolds(p) && (v0 == 1 || o ≈ reynolds(o))
-            deflate, inflate = build_deflate_inflate(p)
-            sym = true
-        else
-            sym = false
-        end
-    end
+    _, _, _, o, sym, deflate, inflate = _bfw_init(p, 0, prob, marg, o, sym, deflate, inflate, verbose)
     lower_bound = zero(T)
     upper_bound = one(T)
     local_model = nothing
     bell_inequality = nothing
     while upper_bound - lower_bound > 10.0^(-precision)
-        x, ds, primal, dual_gap, as, M, β = bell_frank_wolfe(
-            p;
-            v0,
-            epsilon,
-            marg,
-            o,
-            sym,
-            deflate,
-            inflate,
-            kwargs...,
-        )
+        res = bell_frank_wolfe(p; v0, epsilon, marg, o, sym, deflate, inflate, kwargs...)
+        x, ds, primal, dual_gap, as, M, β = res
         if primal > 10epsilon && dual_gap > 10epsilon
             @warn "Please increase nb or max_iteration"
         end

src/utils.jl

@@ -1120,3 +1120,41 @@ function move_marg(FC::AbstractArray{T, N}, sense = -1) where {T, N}
     return circshift(FC, ntuple(i -> sense, Val(N)))
 end
 export move_marg
+
+function _bfw_init(p::Array{T, N}, v0, prob, marg, o, sym, deflate, inflate, verbose) where {T <: Number, N}
+    if !prob
+        LMO = BellCorrelationsLMO
+        DS = BellCorrelationsDS
+        m = collect(size(p))
+        if o === nothing
+            o = zeros(T, size(p))
+            o[end] = marg
+        end
+        reynolds = reynolds_permutedims
+        build_deflate_inflate = build_deflate_inflate_permutedims
+    else
+        LMO = BellProbabilitiesLMO
+        DS = BellProbabilitiesDS
+        m = collect(size(p)[(N ÷ 2 + 1):end])
+        if o === nothing
+            o = ones(T, size(p)) / prod(size(p)[1:(N ÷ 2)])
+        end
+        reynolds = reynolds_permutelastdims
+        build_deflate_inflate = build_deflate_inflate_permutelastdims
+    end
+    # symmetry detection
+    if sym === nothing
+        if all(diff(m) .== 0) && p ≈ reynolds(p) && (v0 == 1 || o ≈ reynolds(o))
+            deflate, inflate = build_deflate_inflate(p)
+            sym = true
+        else
+            sym = false
+        end
+    end
+    if verbose > 1
+        println("   #Inputs: ", all(diff(m) .== 0) ? m[end] - (marg && !prob) : m .- (marg && !prob))
+        println(" Symmetric: ", sym)
+        println(" Dimension: ", length(deflate(p)))
+    end
+    return LMO, DS, m, o, sym, deflate, inflate
+end