diff --git a/Project.toml b/Project.toml
index aa10ed3..dc8146e 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,18 @@ 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"
+SCS = "2"
 XPORTA = "0.1"
-julia = "1.4"
+julia = "1"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]
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..df73a69
--- /dev/null
+++ b/src/LocalPolytope/linear_nonclassicality_witnesses.jl
@@ -0,0 +1,43 @@
+"""
+    linear_nonclassicality_witness(
+        vertices :: Vector{Vector{T}} where T <: Real,
+        test_behavior :: Vector{<:Real};
+        verbose=false :: Bool
+    ) :: Vector{Float64}
+
+Implements a linear program that obtains a linear inequality that witnesses the nonclassciality of
+the provided `test_behavior` (see https://arxiv.org/abs/1303.2849 Eq. 19). A complete set of
+enumerated vertices are required as input. The optimized facet inequality is returned in vector
+form with the classical bound being appended to the end.
+"""
+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/test/unit/LocalPolytope/linear_nonclassicality_witnesses.jl b/test/unit/LocalPolytope/linear_nonclassicality_witnesses.jl
new file mode 100644
index 0000000..7794402
--- /dev/null
+++ b/test/unit/LocalPolytope/linear_nonclassicality_witnesses.jl
@@ -0,0 +1,50 @@
+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
+
+        facet_vec = LocalPolytope.linear_nonclassicality_witness(chsh_vertices, pr_box_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 facet_vec ≈ [-2, 0, -2, 0, 2, 2, 2, -2, 0]
+    end
+end
+
+end
\ No newline at end of file