diff --git a/.github/workflows/deploy_docs.yml b/.github/workflows/deploy_docs.yml
index 3204579..a790487 100644
--- 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
index cf79245..5fe9211 100644
--- 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
index 3c70bcf..950e571 100644
--- 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
index aa10ed3..08326b4 100644
--- 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
index e8f10b1..ec88457 100644
--- 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
index 4260ef4..db16e17 100644
--- 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
index 77c86d8..5747933 100644
--- 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
index 53a5279..1d9f275 100644
--- 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
new file mode 100644
index 0000000..5b999be
--- /dev/null
+++ 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
index a23c23a..4ae35f1 100644
--- 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
index 55f2cce..1c3a155 100644
--- 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
deleted file mode 100644
index 911d5f0..0000000
--- a/test/Project.toml
+++ /dev/null
@@ -1,15 +0,0 @@
-[deps]
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
-QBase = "e52e8ede-12bf-4731-8af7-b01f6064cb11"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-XPORTA = "8c143463-af6f-456f-8aed-72447cb569d2"
-
-[compat]
-JSON = "0.21"
-Polyhedra = "0.6"
-QBase = "0.2"
-XPORTA = "0.1"
-julia = "1.4"
diff --git a/test/integration/LocalPolytope/adjacency_decomposition.jl b/test/integration/LocalPolytope/adjacency_decomposition.jl
index 57519eb..d0b3ce6 100644
--- 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
index 2a3519b..b67a082 100644
--- 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
index c0472e7..afeb4fa 100644
--- 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
new file mode 100644
index 0000000..8350840
--- /dev/null
+++ 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
\ No newline at end of file
diff --git a/test/unit/Nonlocality/optimize_measurements.jl b/test/unit/Nonlocality/optimize_measurements.jl
index ee326eb..ad14fb8 100644
--- a/test/unit/Nonlocality/optimize_measurements.jl
+++ b/test/unit/Nonlocality/optimize_measurements.jl
@@ -12,8 +12,8 @@ using BellScenario
 
         dict = Nonlocality.optimize_measurement(scenario, game, ρ_states)
 
-        @test isapprox(dict["violation"], 0.0, atol=1e-6)
-        @test dict["povm"] ≈ [[1 0;0 0],[0 0;0 0.5],[0 0;0 0.5]]
+        @test isapprox(dict["violation"], 0.0, atol=1e-4)
+        @test isapprox(dict["povm"], [[1 0;0 0],[0 0;0 0.5],[0 0;0 0.5]], atol=1e-4)
     end
 
     @testset "trine states" begin
@@ -23,10 +23,10 @@ using BellScenario
 
         dict = Nonlocality.optimize_measurement(scenario, game, ρ_states)
 
-        @test isapprox(dict["violation"], 0.0, atol=1e-6)
-        @test all(isapprox.(dict["povm"][1], 2/3*ρ_states[1], atol=1e-6))
-        @test all(isapprox.(dict["povm"][2], 2/3*ρ_states[2], atol=1e-6))
-        @test all(isapprox.(dict["povm"][3], 2/3*ρ_states[3], atol=1e-6))
+        @test isapprox(dict["violation"], 0.0, atol=1e-5)
+        @test all(isapprox.(dict["povm"][1], 2/3*ρ_states[1], atol=1e-5))
+        @test all(isapprox.(dict["povm"][2], 2/3*ρ_states[2], atol=1e-5))
+        @test all(isapprox.(dict["povm"][3], 2/3*ρ_states[3], atol=1e-5))
 
         @test dict["states"] == ρ_states
         @test dict["game"] == game
@@ -40,7 +40,7 @@ using BellScenario
 
         dict = Nonlocality.optimize_measurement(scenario, game, ρ_states)
 
-        @test isapprox(dict["violation"], 0.0, atol=1e-6)
+        @test isapprox(dict["violation"], 0.0, atol=1e-5)
         @test all(isapprox.(dict["povm"][1], 1/2*ρ_states[1], atol=1e-3))
         @test all(isapprox.(dict["povm"][2], 1/2*ρ_states[2], atol=1e-3))
         @test all(isapprox.(dict["povm"][3], 1/2*ρ_states[3], atol=1e-3))
@@ -91,13 +91,13 @@ end
         opt_dictA = Nonlocality.optimize_measurement(scenario, game, ρ_AB, A_POVMs=POVMs)
         opt_dictB = Nonlocality.optimize_measurement(scenario, game, ρ_AB, B_POVMs=POVMs)
 
-        @test opt_dictA["score"] ≈ 2.20710279
-        @test opt_dictA["violation"] ≈ 0.207102796
+        @test isapprox(opt_dictA["score"], 2.20710279, atol=1e-5)
+        @test isapprox(opt_dictA["violation"], 0.207102796, atol=1e-5)
         @test opt_dictA["A_POVMs"] isa Vector{<:POVM}
         @test opt_dictA["state"] isa State
 
-        @test opt_dictB["score"] ≈ 2.20710279
-        @test opt_dictB["violation"] ≈ 0.207102796
+        @test isapprox(opt_dictB["score"], 2.20710279, atol=1e-5)
+        @test isapprox(opt_dictB["violation"], 0.207102796, atol=1e-5)
         @test opt_dictB["B_POVMs"] isa Vector{<:POVM}
         @test opt_dictB["state"] isa State
     end
@@ -120,8 +120,8 @@ end
         @test opt_dictA["scenario"] == scenario
         @test opt_dictA["game"] == game
         @test opt_dictA["game"].β == 2
-        @test opt_dictA["score"] ≈ 2.20710279
-        @test opt_dictA["violation"] ≈ 0.207102796
+        @test isapprox(opt_dictA["score"], 2.20710279, atol=1e-5)
+        @test isapprox(opt_dictA["violation"], 0.207102796, atol=1e-5)
         @test opt_dictA["state"] == ρ_AB
         @test opt_dictA["A_POVMs"] == POVMs
         @test isapprox(
@@ -137,8 +137,8 @@ end
         @test opt_dictB["scenario"] == scenario
         @test opt_dictB["game"] == game
         @test opt_dictB["game"].β == 2
-        @test opt_dictB["score"] ≈ 2.2071028
-        @test opt_dictB["violation"] ≈ 0.207102796
+        @test isapprox(opt_dictB["score"], 2.2071028, atol=1e-5)
+        @test isapprox(opt_dictB["violation"], 0.207102796, atol=1e-5)
         @test opt_dictB["state"] == ρ_AB
         @test opt_dictB["B_POVMs"] == POVMs
         @test isapprox(