Merge pull request #20 from ChitambarLab/v0-1-3

V0.1.3
ChitambarLab · Dec 6, 2023 · 8561043 · 8561043 · bdoolittle · Dec 6, 2023
2 parents 74932c1 + cb28d32
commit 8561043
Show file tree

Hide file tree

Showing 17 changed files with 201 additions and 59 deletions.
diff --git a/.github/workflows/deploy_docs.yml b/.github/workflows/deploy_docs.yml
@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.4'
+          version: '1.8'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

diff --git a/.gitignore b/.gitignore
@@ -12,3 +12,4 @@ docs/Manifest.toml
 
 porta_tmp
 porta.log
+Highs.log
diff --git a/CITATION.bib b/CITATION.bib
@@ -3,7 +3,7 @@ @misc{BellScenario.jl
 	title   = {BellScenario.jl},
 	howpublished =  {\url{https://github.com/ChitambarLab/BellScenario.jl}},
 	url     = {https://github.com/ChitambarLab/BellScenario.jl},
-	version = {v0.1.2},
+	version = {v0.1.3},
 	year    = {2020},
 	month   = {12}
 }
diff --git a/Project.toml b/Project.toml
@@ -1,13 +1,15 @@
 name = "BellScenario"
 uuid = "4e58fcc1-4405-48af-b0f2-42df636d9190"
 authors = ["Brian Doolittle"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 QBase = "e52e8ede-12bf-4731-8af7-b01f6064cb11"
@@ -17,10 +19,23 @@ XPORTA = "8c143463-af6f-456f-8aed-72447cb569d2"
 
 [compat]
 Combinatorics = "1"
-Convex = "0.14"
+Convex = "0.15"
+HiGHS = "1"
 JSON = "0.21"
-Polyhedra = "0.6"
+JuMP = "1"
+Polyhedra = "0.7"
 QBase = "0.2"
-SCS = "0.7"
-XPORTA = "0.1"
-julia = "1.4"
+SCS = "2"
+XPORTA = "0.1.3"
+julia = "1"
+
+[extras]
+Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
+QBase = "e52e8ede-12bf-4731-8af7-b01f6064cb11"
+SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+XPORTA = "8c143463-af6f-456f-8aed-72447cb569d2"
+
+[targets]
+test = ["Test", "SafeTestsets", "Polyhedra", "XPORTA", "QBase", "SCS"]
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,6 +1,8 @@
 [deps]
+BellScenario = "4e58fcc1-4405-48af-b0f2-42df636d9190"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 QBase = "e52e8ede-12bf-4731-8af7-b01f6064cb11"
+XPORTA = "8c143463-af6f-456f-8aed-72447cb569d2"
 
 [compat]
-Documenter = "0.25.5"
+Documenter = "1"
diff --git a/docs/make.jl b/docs/make.jl
@@ -1,3 +1,5 @@
+push!(LOAD_PATH,"../src/")
+
 using Documenter, BellScenario
 
 DocMeta.setdocmeta!(BellScenario, :DocTestSetup, :(using BellScenario); recursive=true)
@@ -34,6 +36,7 @@ makedocs(;
     repo="https://github.com/ChitambarLab/BellScenario.jl/blob/{commit}{path}#L{line}",
     sitename="BellScenario.jl",
     authors="Brian Doolittle",
+    warnonly = true,
 )
 
 deploydocs(;

diff --git a/docs/src/LocalPolytope/facets.md b/docs/src/LocalPolytope/facets.md
@@ -6,3 +6,9 @@ CurrentModule = LocalPolytope
 ```@docs
 facets
 ```
+
+## Linear Programming of Facets
+
+```@docs
+linear_nonclassicality_witness
+```
diff --git a/src/LocalPolytope/LocalPolytope.jl b/src/LocalPolytope/LocalPolytope.jl
@@ -58,7 +58,8 @@ resources were used than anticipated.
 module LocalPolytope
 
 using LinearAlgebra, Combinatorics
-using XPORTA, Polyhedra
+using JuMP, HiGHS  # Linear Programming
+using XPORTA, Polyhedra  # Convex Polytope
 
 import Polyhedra: vrep
 
@@ -88,6 +89,7 @@ include("./facets.jl")
 include("./generators.jl")
 include("./adjacency_decomposition.jl")
 include("./utils.jl")
+include("./linear_nonclassicality_witnesses.jl")
 
 # legacy code
 include("../Legacy/LocalPolytope.jl")

diff --git a/src/LocalPolytope/linear_nonclassicality_witnesses.jl b/src/LocalPolytope/linear_nonclassicality_witnesses.jl
@@ -0,0 +1,76 @@
+"""
+    linear_nonclassicality_witness(
+        vertices :: Vector{Vector{T}} where T <: Real,
+        test_behavior :: Vector{<:Real};
+        verbose=false :: Bool
+    ) :: Vector{Float64}
+
+Obtains a linear inequality ``(\\mathbf{G}^\\star, \\beta^\\star)`` bounding the convex hull of the
+set of `vertices` (``\\mathcal{V}``) and is violated by the `test_behavior` ``(\\mathbf{P}\\notin \\text{Conv}(\\mathcal{V}))``.
+This task is achieved using the  linear program described by Brunner et al. in Eq. 19 of
+[https://arxiv.org/abs/1303.2849](https://arxiv.org/abs/1303.2849). The linear program is
+
+
+```math
+\\begin{align}
+    \\min_{(\\mathbf{G}, \\beta)} \\quad & \\langle \\mathbf{G},\\mathbf{P}\\rangle - \\beta & \\notag\\\\
+    \\text{s.t.} \\quad  & \\langle \\mathbf{G}, \\mathbf{V} \\rangle - \\beta \\leq 0 & \\quad \\forall \\;\\; \\mathbf{V} \\in \\mathcal{V}  \\notag\\\\
+    & \\langle \\mathbf{G}, \\mathbf{P} \\rangle \\leq 1 & \\notag\\\\
+\\end{align}
+```
+
+The solution to the linear program is a linear nonclassicality witness ``(\\mathbf{G}^\\star, \\beta^\\star)`` becuase the
+constraint ``\\lange\\mathbf{G}^\\star, \\mathbf{V} \\rangle - \\beta^\\star \\leq 0`` ensures that no behavior in the polytope
+``\\text{Conv}(\\mathcal{V})`` can violate the inequality. Provided that ``\\mathbf{P} \\notin \\text{Conv}(\\mathcal{V})``
+technique the output linear inequality witnesses the nonclassicality of the `test_behavior`.
+
+The optimized facet inequality ``(\\mathbf{G}^\\star, \\beta^\\star)`` is returned as a vector ``(G^\\star_{0,0}, \\dots, G^\\star_{Y,X}, -\\beta^\\star)``
+where ``G^\\star_{y,x}`` are the elements of ``\\mathbf{G}^\\star``.
+
+!!! note "Supporting Software"
+    The linear programming is performed using [HiGHS](https://highs.dev/)
+    solver via the [`JuMP`](https://jump.dev/JuMP.jl/stable/)
+    interface. Please refer to the source code for more details.
+
+!!! note "Converting Output into Bell Game"
+    The linear programming software outputs numerical values that have numerical error. Moreover, the linear inequality is
+    scaled such that the classical bound is zero and the `test_behavior` score is one. In order to convert the output
+    inequality into a `BellGame`, care must be taken to obtain the correct scaling factor to ensure that elements are integers.
+
+!!! note "Classical Test Behavior"
+    If the `test_behavior` ``\\mathbf{P}`` is classical, meaning it satisfies ``\\mathbf{P}\\in\\text{Conv}(\\mathcal{V})``,
+    then the zero vector is returned as the optimal solution. Note that if all elements of ``\\mathbf{G}^\\star``
+    satisfy ``G^\\star_{y,x}=0``, then all behaviors ``\\mathbf{P} \\in \\text{Conv}(\\mathcal{V})`` are trivially optimal
+    as ``\\langle\\mathbf{G}^{\\star}, \\mathbf{P} \\rangle - \\beta^\\star \\leq 0``.
+"""
+function linear_nonclassicality_witness(
+    vertices :: Vector{Vector{T}} where T <: Real,
+    test_behavior :: Vector{<:Real};
+    verbose=false :: Bool
+) :: Vector{Float64}
+    dim_v = length(vertices[1]) + 1
+
+    # initializing model
+    model = Model(HiGHS.Optimizer)
+    set_attribute(model, MOI.Silent(), !verbose)
+
+    # adding variable for inequality
+    @variable(model, s[1:dim_v])
+
+    # adding constraints to model
+    for v in vertices
+        va = [v..., -1]
+        @constraint(model, sum(s.*va) <= 0)
+    end
+
+    @constraint(model, c, sum(s.*[test_behavior..., -1]) <= 1)
+
+    # defining the optimization objective
+    @objective(model, Max, sum(s.*[test_behavior..., -1]))
+
+    # optimizing
+    optimize!(model)
+
+    # return optimized linear inequality in vector form
+    return value.(s)
+end
diff --git a/src/Nonlocality/optimize_measurements.jl b/src/Nonlocality/optimize_measurements.jl
@@ -49,8 +49,7 @@ function optimize_measurement(
     # add the objective
     objective = maximize(real(tr(sum(Π_vars[1:end-1] .* H_y))), constraints)
 
-    # optimize model
-    solve!(objective, SCS.Optimizer(verbose=0))
+    solve!(objective, SCS.Optimizer())
 
     # parse/return results
     score = objective.optval
@@ -172,7 +171,7 @@ function _optimize_measurement_B(
     end
 
     # optimize model
-    solve!(problem, SCS.Optimizer(verbose=0))
+    solve!(problem, SCS.Optimizer())
 
     # parse/return results
     score = problem.optval
@@ -239,7 +238,7 @@ function _optimize_measurement_A(
     end
 
     # optimize model
-    solve!(problem, SCS.Optimizer(verbose=0))
+    solve!(problem, SCS.Optimizer())
 
     # parse/return results
     score = problem.optval

diff --git a/src/combinatorics.jl b/src/combinatorics.jl
@@ -31,8 +31,8 @@ Each partition is a vector containing a  set of `k` vectors designating each gro
 E.g.
 
 ```jldoctest
-julia> stirling2_partitions( 4, 2 )
-7-element Array{Array{Array{Int64,1},1},1}:
+julia> stirling2_partitions(4, 2)
+7-element Vector{Vector{Vector{Int64}}}:
  [[1, 2, 3], [4]]
  [[3], [1, 2, 4]]
  [[1, 2], [3, 4]]
@@ -96,8 +96,8 @@ that the column id is grouped into the corresponding row.
 E.g.
 
 ```jldoctest
-julia> stirling2_matrices( 4, 2 )
-7-element Array{Array{Bool,2},1}:
+julia> stirling2_matrices(4, 2)
+7-element Vector{Matrix{Bool}}:
  [1 1 1 0; 0 0 0 1]
  [0 0 1 0; 1 1 0 1]
  [1 1 0 0; 0 0 1 1]
@@ -135,8 +135,8 @@ Generates the set of square permutation matrices of dimension `dim`.
 E.g.
 
 ```jldoctest
-julia> permutation_matrices( 3 )
-6-element Array{Array{Bool,2},1}:
+julia> permutation_matrices(3)
+6-element Vector{Matrix{Bool}}:
  [1 0 0; 0 1 0; 0 0 1]
  [1 0 0; 0 0 1; 0 1 0]
  [0 1 0; 1 0 0; 0 0 1]
@@ -160,7 +160,7 @@ E.g.
 
 ```jldoctest
 julia> n_choose_k_matrices( 4, 2 )
-6-element Array{Array{Bool,2},1}:
+6-element Vector{Matrix{Bool}}:
  [1 0; 0 1; 0 0; 0 0]
  [1 0; 0 0; 0 1; 0 0]
  [1 0; 0 0; 0 0; 0 1]

diff --git a/test/Project.toml b/test/Project.toml
diff --git a/test/integration/LocalPolytope/adjacency_decomposition.jl b/test/integration/LocalPolytope/adjacency_decomposition.jl
@@ -39,15 +39,15 @@ end
         skip = [BellGame([1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],1)]
         dict = LocalPolytope.adjacency_decomposition(vertices, BG, scenario, skip_games=skip, dir=test_dir)
 
-        # positvity included, but is skipped in computation
-        @test collect(keys(dict)) == [
+        # non-negativity included, but is skipped in computation
+        @test issetequal(collect((keys(dict))), [
             [1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],
             [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1],
             [1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],
             [2 0 0 0;1 1 1 0;0 2 0 0;0 0 1 1],
             [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0],
             [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1],
-        ]
+        ])
 
         @test dict[[1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0]]["skipped"]
 
@@ -63,14 +63,14 @@ end
 
         dict = LocalPolytope.adjacency_decomposition(vertices, BG, scenario, dir=test_dir)
 
-        @test collect(keys(dict)) == [
+        @test issetequal(collect(keys(dict)), [
             [1 0 0 0;1 0 0 0;0 1 0 0;0 0 1 0],
             [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1],
             [1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0],
             [2 0 0 0;1 1 1 0;0 2 0 0;0 0 1 1],
             [2 0 0 0;1 1 1 0;0 2 0 0;0 0 2 0],
             [1 1 0 0;1 0 1 0;0 1 1 0;0 0 0 1]
-        ]
+        ])
 
         @test dict[[1 0 0 0;1 0 0 0;1 0 0 0;0 0 0 0]]["skipped"] == false
     end

diff --git a/test/unit/LocalPolytope.jl b/test/unit/LocalPolytope.jl
@@ -9,7 +9,7 @@ using BellScenario
     @testset "LocalPolytope.vrep()" begin
         @testset "no arguments" begin
             scenario = LocalSignaling(3,3,2)
-            local_poly = vrep(scenario)
+            local_poly = LocalPolytope.vrep(scenario)
 
             @test local_poly isa XPORTA.Polyhedron
             @test npoints(local_poly) == 21
@@ -20,7 +20,7 @@ using BellScenario
 
         @testset "passing arguments" begin
             scenario = LocalSignaling(3,3,2)
-            local_poly = vrep(scenario, rank_d_only=true)
+            local_poly = LocalPolytope.vrep(scenario, rank_d_only=true)
 
             @test local_poly isa XPORTA.Polyhedron
             @test npoints(local_poly) == 18

diff --git a/test/unit/LocalPolytope/facets.jl b/test/unit/LocalPolytope/facets.jl
@@ -1,4 +1,4 @@
-using Test, Polyhedra
+using Test, XPORTA
 
 @testset "./src/LocalPolytope/facets.jl" begin
 

diff --git a/test/unit/LocalPolytope/linear_nonclassicality_witnesses.jl b/test/unit/LocalPolytope/linear_nonclassicality_witnesses.jl
@@ -0,0 +1,53 @@
+using Test
+
+@testset "./src/LocalPolytope/linear_nonclassicality_witnesses.jl" begin
+
+using BellScenario
+
+@testset "linear_nonclassicality_witness" begin
+    @testset "CHSH Inequality Correlation Vertices" begin
+
+        chsh_correlation_vertices = [
+            [-1, -1, -1, -1,  1,  1,  1,  1],
+            [-1, -1, -1,  1,  1, -1,  1, -1],
+            [-1, -1,  1, -1, -1,  1, -1,  1],
+            [-1, -1,  1,  1, -1, -1, -1, -1],
+            [-1,  1, -1, -1,  1,  1, -1, -1],
+            [-1,  1, -1,  1,  1, -1, -1,  1],
+            [-1,  1,  1, -1, -1,  1,  1, -1],
+            [-1,  1,  1,  1, -1, -1,  1,  1],
+            [ 1, -1, -1, -1, -1, -1,  1,  1],
+            [ 1, -1, -1,  1, -1,  1,  1, -1],
+            [ 1, -1,  1, -1,  1, -1, -1,  1],
+            [ 1, -1,  1,  1,  1,  1, -1, -1],
+            [ 1,  1, -1, -1, -1, -1, -1, -1],
+            [ 1,  1, -1,  1, -1,  1, -1,  1],
+            [ 1,  1,  1, -1,  1, -1,  1, -1],
+            [ 1,  1,  1,  1,  1,  1,  1,  1],
+        ]
+
+        pr_box_test_behavior = [0, 0, 0, 0, 1, 1, 1, -1]
+
+        facet_vec = LocalPolytope.linear_nonclassicality_witness(chsh_correlation_vertices, pr_box_test_behavior[:])
+
+        @test facet_vec ≈ [0, 0, 0, 0, 0.5, 0.5, 0.5, -0.5, 1]  # [A0, A1, B0, B1, AB00, AB01, AB10, AB11]
+    end
+
+    @testset "CH Inequality Probability Vertices" begin
+        chsh_scenario = BipartiteNonSignaling(2, 2, 2, 2)
+
+        chsh_vertices = LocalPolytope.vertices(chsh_scenario, "non-signaling")
+
+        pr_box_test_behavior = [1, 1, 1, 1, 1, 1, 1, 0] / 2  # nonlocal test behavior
+        local_test_behavior = [1/2, 1/2, 1/2, 1/2, 1/4, 1/4, 1/4, 1/4]  # white noise behavior
+
+        chsh_facet_vec = LocalPolytope.linear_nonclassicality_witness(chsh_vertices, pr_box_test_behavior[:])
+        local_game_vec = LocalPolytope.linear_nonclassicality_witness(chsh_vertices, local_test_behavior[:])
+
+        # -2*PA(0|0) - 2*PB(0|0) + 2*PAB(00|00) + 2*PAB(01|00) + 2*PAB(10|00) - 2*PAB(11|00) <= 0
+        @test chsh_facet_vec ≈ [-2, 0, -2, 0, 2, 2, 2, -2, 0]
+        @test local_game_vec ≈ [0, 0, 0, 0, 0, 0, 0, 0, 0]  # zero vector returned as optimal solution.
+    end
+end
+
+end