diff --git a/CITATION.bib b/CITATION.bib
index fc4d008..1194c3e 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.1},
+	version = {v0.1.2},
 	year    = {2020},
 	month   = {6}
 }
diff --git a/Project.toml b/Project.toml
index 69efd13..aa10ed3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BellScenario"
 uuid = "4e58fcc1-4405-48af-b0f2-42df636d9190"
 authors = ["Brian Doolittle"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -9,6 +9,7 @@ Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 QBase = "e52e8ede-12bf-4731-8af7-b01f6064cb11"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
@@ -18,6 +19,7 @@ XPORTA = "8c143463-af6f-456f-8aed-72447cb569d2"
 Combinatorics = "1"
 Convex = "0.14"
 JSON = "0.21"
+Polyhedra = "0.6"
 QBase = "0.2"
 SCS = "0.7"
 XPORTA = "0.1"
diff --git a/docs/src/LocalPolytope/vertices.md b/docs/src/LocalPolytope/vertices.md
index e8b1dfd..b5a8d22 100644
--- a/docs/src/LocalPolytope/vertices.md
+++ b/docs/src/LocalPolytope/vertices.md
@@ -5,10 +5,17 @@ CurrentModule = LocalPolytope
 
 Vertices are extreme points of the local polytope.
 They correspond to deterministic strategies.
+The vertex representation of a convex polytope can be computed via the
+[Polyhedra.jl](https://github.com/JuliaPolyhedra/Polyhedra.jl) interface.
+
+```@docs
+vrep
+```
 
 ## Vertex Enumeration
 
-Enumerate the local polytope vertices for the specified Bell [`Scenario`](@ref).
+The vertices of each local polytope can be enumerated for the specified
+Bell [`Scenario`](@ref).
 
 ```@docs
 vertices
diff --git a/docs/src/assets/scenario_images/classical_chsh_scenario.png b/docs/src/assets/scenario_images/classical_chsh_scenario.png
new file mode 100644
index 0000000..c623413
Binary files /dev/null and b/docs/src/assets/scenario_images/classical_chsh_scenario.png differ
diff --git a/docs/src/user_guide.md b/docs/src/user_guide.md
index 7171baf..48190e2 100644
--- a/docs/src/user_guide.md
+++ b/docs/src/user_guide.md
@@ -1,5 +1,102 @@
 # User Guide
 
+## Quickstart
+
 ```julia
-using Pkg; Pkg.add("BellScenario")
+julia> using Pkg; Pkg.add("BellScenario")
+```
+
+```@example tutorial
+using BellScenario
+```
+
+## CHSH Scenario
+
+The CHSH scenario is a [`BipartiteNonSignaling`](@ref) scenario where Alice and Bob
+each have a black-box with binary inputs and outputs.
+
+![Classical CHSH Scenario](assets/scenario_images/classical_chsh_scenario.png)
+
+This scenario is significant because it is the simplest Bell scenario in which
+quantum nonlocality can be observed.
+We will use BellScenario.jl to compute the CH Bell inequality and optimize quantum
+measurements to violate the CH inequality.
+First, create a CHSH `Scenario` to specify the black-box arrangement in the figure
+above.
+
+```@example tutorial
+# (num_out_A, num_out_B, num_in_A, num_in_B)
+chsh_scenario = BipartiteNonSignaling(2,2,2,2)
+```
+
+Bell inequalities bound the set of local (classical) correlations.
+The local set is a convex polytope referred to as the *local polytope* and the
+facets of the local polytope are Bell inequalities.
+The standard method of computing Bell inequalities is to first compute the local
+polytope vertices, then apply a polytope transformation algorithm to compute the
+Bell inequalities.
+
+The BellScenario.jl package provides the [`LocalPolytope`](@ref) module to compute
+Bell inequalities.
+The first step is to commpute the vertex representation for the CHSH scenario.
+
+```@example tutorial
+chsh_polytope = LocalPolytope.vrep(chsh_scenario)
+```
+
+Then, the Bell inequalities can computed using the [`LocalPolytope.facets`](@ref)
+function.
+
+```@example tutorial
+chsh_facets = LocalPolytope.facets(chsh_polytope)
+```
+
+We'll take ``15^{th}`` facet as it represents the CH inequality
+
+```math
+- P_A(1|2) - P_B(1|1) + P(11|11) - P(11|12) + P(11|21) + P(11|22) \leq 0.
+```
+
+In fact, this inequality is equivalent to the more celebrated CHSH inequality.
+The difference is that the CH inequality is expressed in terms of probabilities
+whereas the CHSH inequality is expressed in terms of bipartite correlators.
+
+```@example tutorial
+ch_inequality = chsh_facets[15]
+```
+
+Now that we have computed a Bell inequality, we can find a quantum violation using
+the [`Nonlocality`](@ref) module.
+In this example, we will fix Alice's measurement and the quantum state shared
+between Alice and Bob.
+
+```@example tutorial
+# maximally entangled state
+ρ_AB = [1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1]/2
+
+# Alice's measurement bases
+Π_ax = [
+    [[1 0;0 0], [0 0;0 1]],      # Pauli Z basis
+    [[1 1;1 1]/2, [1 -1;-1 1]/2] # Pauli X basis
+]
+```
+
+Then, we convert the `ch_inequality` into a general representation of a [`BellGame`](@ref).
+
+```@example tutorial
+ch_game = convert(BellGame, ch_inequality, chsh_scenario)
+```
+
+Finally, we optimize Bob's measurement with respect to the fixed state and measurements.
+
+```@example tutorial
+opt_dict = Nonlocality.optimize_measurement(
+    chsh_scenario, ch_game, ρ_AB, A_POVMs=Π_ax
+)
+```
+
+We see that the inequality is violated for the optimized measurement and states.
+
+```@example tutorial
+opt_dict["violation"]
 ```
diff --git a/src/LocalPolytope/LocalPolytope.jl b/src/LocalPolytope/LocalPolytope.jl
index 1020e09..53a5279 100644
--- a/src/LocalPolytope/LocalPolytope.jl
+++ b/src/LocalPolytope/LocalPolytope.jl
@@ -45,21 +45,45 @@ however, it is important to note that a Bell violation can simply mean that more
 resources were used than anticipated.
 
 ### Module Exports:
-* [`vertices`](@ref) - Compute the set of extreme-points for the Local Polytope.
-* [`facets`](@ref) - Compute the linear inequalities which bound the Local Polytope.
+* [`vertices`](@ref) - Compute the set of extreme-points for the local polytope.
+* [`num_vertices`](@ref) - The number of vertices for the local polytope.
+* [`vrep`](@ref) - Construct a [`Polyhedron`](https://juliapolyhedra.github.io/Polyhedra.jl/stable/polyhedron/)
+    in the vertex representation.
+* [`facets`](@ref) - Compute the linear inequalities which bound the local polytope.
 * [`generator_vertex`](@ref) - Provide a canonical form for a vertex.
 * [`generator_facet`](@ref) - Provide a canonical form for a facet.
-* [`adjacency_decomposition`](@ref) - Efficiently compute the generating facets for the Local
-    Polytope using the adjacency decomposition technique.
+* [`adjacency_decomposition`](@ref) - Efficiently compute the generator facets for the local
+    polytope using the adjacency decomposition technique.
 """
 module LocalPolytope
 
 using LinearAlgebra, Combinatorics
-using XPORTA
+using XPORTA, Polyhedra
+
+import Polyhedra: vrep
 
 using ..BellScenario
 
 include("./vertices.jl")
+
+"""
+    vrep(scenario::Scenario; vertices_kwargs...) :: XPORTA.Polyhedron
+
+Constructs a `Polyhedron` using the vertex representation.
+See [Polyhedra.jl](https://github.com/JuliaPolyhedra/Polyhedra.jl)
+for more details.
+The `vertices_kwargs` keyword arguments are passed through to the `vertices`
+function for each `Scenario`.
+
+!!! note "Return Type"
+    This function differs from the Polyhedra.jl implementation in that it returns
+    a `Polyhedron` type rather than a `V-Representation`.
+    This is done to reduce the number of steps required to construct a polyhedron.
+"""
+function vrep(scenario::Scenario; vertices_kwargs...) :: XPORTA.Polyhedron
+    polyhedron(vrep(vertices(scenario; vertices_kwargs...)), XPORTA.Library())
+end
+
 include("./facets.jl")
 include("./generators.jl")
 include("./adjacency_decomposition.jl")
diff --git a/src/LocalPolytope/facets.jl b/src/LocalPolytope/facets.jl
index b9ddfb9..11c6dbb 100644
--- a/src/LocalPolytope/facets.jl
+++ b/src/LocalPolytope/facets.jl
@@ -55,3 +55,15 @@ function facets(
         "equalities" => map(row_id -> equalities[row_id,:], 1:size(equalities,1))
     )
 end
+
+"""
+    facets(poly::Polyhedron) :: Vector{Vector{Int64}}
+
+Returns the facets of the local polytope `poly`. If the facets have not already
+been computed, then they are transformed into the half-space representation from
+the vertex representation.
+"""
+function facets(poly::Polyhedron) :: Vector{Vector{Int64}}
+    inequalities = convert.(Int64, hrep(poly).inequalities)
+    map(row_id -> inequalities[row_id,:] , 1:size(inequalities,1))
+end
diff --git a/src/Nonlocality/optimize_measurements.jl b/src/Nonlocality/optimize_measurements.jl
index 0a12a1e..a23c23a 100644
--- a/src/Nonlocality/optimize_measurements.jl
+++ b/src/Nonlocality/optimize_measurements.jl
@@ -1,12 +1,12 @@
 export optimize_measurement
 
-@doc raw"""
+"""
 [`LocalSignaling`](@ref) scenario:
 
     optimize_measurement(
         scenario :: LocalSignaling,
         game :: BellGame,
-        ρ_states :: Vector{<:State}
+        ρ_states :: Vector{<:AbstractMatrix}
     )
 
 Finds the measurement that optimizes the score of the [`BellGame`](@ref) against
@@ -14,25 +14,26 @@ the set of quantum states `ρ_states`.
 The optimization is performed with the following semi-definite program:
 
 ```math
-\begin{aligned}
-&\max_{\{\Pi_y\}} \sum_{x,y} G_{x,y} \text{Tr}[\Pi_y \rho_x] \\
-&s.t. \quad \sum_y \Pi_y = \mathbb{I}, \quad \Pi_y \geq 0
-\end{aligned}
+\\begin{aligned}
+&\\max_{\\{\\Pi_y\\}} \\sum_{x,y} G_{x,y} \\text{Tr}[\\Pi_y \\rho_x] \\\\
+&s.t. \\quad \\sum_y \\Pi_y = \\mathbb{I}, \\quad \\Pi_y \\geq 0
+\\end{aligned}
 ```
 """
 function optimize_measurement(
     scenario::LocalSignaling,
     game::BellGame,
-    ρ_states::Vector{<:State},
+    ρ_states::Vector{<:AbstractMatrix};
+    atol=1e-7::Float64
 ) :: Dict
     if scenario.X != length(ρ_states)
         throw(DomainError(scenario, "expected length of `ρ_states` is $(scenario.X)), but got $(length(ρ_states)) instead"))
-    end
-
-    if size(ρ_states[1]) != (scenario.d,scenario.d)
+    elseif size(ρ_states[1]) != (scenario.d,scenario.d)
         throw(DomainError(ρ_states, "dimension of `ρ_states` is not $(scenario.d)"))
     end
 
+    states = map(ρ -> ρ isa State ? ρ : State(ρ, atol=atol), ρ_states)
+
     norm_game_vector = convert(Vector{Int64}, game)
     norm_bound = norm_game_vector[end]
     norm_game = reshape(norm_game_vector[1:(end-1)], (scenario.Y-1, scenario.X))
@@ -43,7 +44,7 @@ function optimize_measurement(
     constraints += (sum(map(Π_y -> imag(Π_y), Π_vars)) == zeros(Float64, scenario.d, scenario.d))
 
     # sum up the state contibutions for each row
-    H_y = map(row_id -> sum(norm_game[row_id,:] .* ρ_states), 1:scenario.Y-1)
+    H_y = map(row_id -> sum(norm_game[row_id,:] .* states), 1:scenario.Y-1)
 
     # add the objective
     objective = maximize(real(tr(sum(Π_vars[1:end-1] .* H_y))), constraints)
@@ -61,18 +62,18 @@ function optimize_measurement(
         "povm" => Π_opt,
         "game" => game,
         "scenario" => scenario,
-        "states" => ρ_states
+        "states" => states
     )
 end
 
-@doc raw"""
+"""
 [`BipartiteNonSignaling`](@ref) scenario:
 
     optimize_measurement(
         game :: BellGame,
         scenario :: BipartiteNonSignaling,
-        ρ_AB :: State;
-        A_POVMs :: Vector{<:POVM},
+        ρ_AB :: AbstractMatrix;
+        A_POVMs :: Vector{<:AbstractVector{<:AbstractMatrix}},
     ) :: Dict
 
 Find Bob's measurement which optimizes a `BellGame`'s score for the shared quantum
@@ -80,17 +81,17 @@ state `ρ_AB` and POVM measurement applied by Alice.
 The following semi-definite program optimizes the Bob's POVM:
 
 ```math
-\begin{aligned}
-    &\max_{\{\Pi_b^y\}} \sum_{a,b,x,y} G_{a,b,x,y}\text{Tr}[(\Pi_a^x \otimes \Pi_b^y)\rho_AB] \\
-    &s.t. \quad \sum_b \Pi_b^y = \mathbb{I},\quad \Pi_b^y \geq 0 \quad \forall\; y
-\end{aligned}
+\\begin{aligned}
+    &\\max_{\\{\\Pi_b^y\\}} \\sum_{a,b,x,y} G_{a,b,x,y}\\text{Tr}[(\\Pi_a^x \\otimes \\Pi_b^y)\\rho_{AB}] \\\\
+    &s.t. \\quad \\sum_b \\Pi_b^y = \\mathbb{I},\\quad \\Pi_b^y \\geq 0 \\quad \\forall\\; y
+\\end{aligned}
 ```
 
     optimize_measurement(
         game :: BellGame,
         scenario :: BipartiteNonSignaling,
-        ρ_AB :: State;
-        B_POVMs :: Vector{<:POVM},
+        ρ_AB :: AbstractMatrix;
+        B_POVMs :: Vector{<:AbstractVector{<:AbstractMatrix}},
     ) :: Dict
 
 Find Alice's measurement which optimizes a `BellGame`'s score for the shared quantum
@@ -98,23 +99,26 @@ state `ρ_{AB}` and POVM measurement applied by Bob.
 The following semi-definite program optimizes the Alice's POVM:
 
 ```math
-\begin{aligned}
-&\max_{\{\Pi_a^x\}} \sum_{a,b,x,y} G_{a,b,x,y}\text{Tr}[(\Pi_a^x \otimes \Pi_b^y)\rho_{AB}] \\
-&s.t. \quad \sum_a \Pi_a^x = \mathbb{I},\quad \Pi_a^x \geq 0 \quad \forall \;x
-\end{aligned}
+\\begin{aligned}
+&\\max_{\\{\\Pi_a^x\\}} \\sum_{a,b,x,y} G_{a,b,x,y}\\text{Tr}[(\\Pi_a^x \\otimes \\Pi_b^y)\\rho_{AB}] \\\\
+&s.t. \\quad \\sum_a \\Pi_a^x = \\mathbb{I},\\quad \\Pi_a^x \\geq 0 \\quad \\forall \\;x
+\\end{aligned}
 ```
 """
 function optimize_measurement(
     scenario::BipartiteNonSignaling,
     game::BellGame,
-    ρ_AB :: State;
-    A_POVMs=Vector{POVM}(undef,0) :: Vector{<:POVM},
-    B_POVMs=Vector{POVM}(undef,0) :: Vector{<:POVM}
+    ρ_AB :: AbstractMatrix;
+    A_POVMs=Vector{POVM}(undef,0) :: Vector{<:AbstractVector{<:AbstractMatrix}},
+    B_POVMs=Vector{POVM}(undef,0) :: Vector{<:AbstractVector{<:AbstractMatrix}},
+    atol=1e-7::Float64
 ) :: Dict
+    ρ = ρ_AB isa State ? ρ_AB : State(ρ_AB, atol=atol)
+
     if length(A_POVMs) > 0
-        return _optimize_measurement_B(scenario, game, ρ_AB, A_POVMs)
+        return _optimize_measurement_B(scenario, game, ρ, A_POVMs, atol=atol)
     elseif length(B_POVMs) > 0
-        return _optimize_measurement_A(scenario, game, ρ_AB, B_POVMs)
+        return _optimize_measurement_A(scenario, game, ρ, B_POVMs, atol=atol)
     else
         throw(DomainError((A_POVMs,B_POVMs), "either `A_POVMs` or `B_POVMs` must be specified."))
     end
@@ -126,22 +130,23 @@ end
 function _optimize_measurement_B(
     scenario::BipartiteNonSignaling,
     game::BellGame,
-    ρ_AB :: State,
-    A_POVMs :: Vector{<:POVM}
+    ρ_AB :: AbstractMatrix,
+    A_POVMs :: Vector{<:AbstractVector{<:AbstractMatrix}};
+    atol=1e-7::Float64
 ) :: Dict
     if !(scenario.X == length(A_POVMs))
         throw( DomainError(
             (scenario.X, " != length(A_POVMs)"),
             "A distinct POVM is not specified for input `X` of `BipartiteNonSignaling` scenario"
         ))
-    end
-
-    if !(all(i -> scenario.A == length(A_POVMs[i]), 1:scenario.X ))
+    elseif !(all(i -> scenario.A == length(A_POVMs[i]), 1:scenario.X ))
         throw( DomainError(
             A_POVMs, "Each POVM must have $(scenario.A) outputs."
         ))
     end
 
+    Π_Ax = map(Π -> Π isa POVM ? Π : POVM(Π, atol=atol), A_POVMs)
+
     ρ_dim = size(ρ_AB,1)
     Π_A_dim = size(A_POVMs[1][1],1)
 
@@ -156,7 +161,7 @@ function _optimize_measurement_B(
     # Objective function
     problem = maximize(real(
         sum(x -> sum(y -> sum(a -> sum(b ->
-            game[(a-1)*scenario.B+b,(x-1)*scenario.Y+y]*tr(kron(A_POVMs[x][a],B_POVMs[y][b])*ρ_AB.M),
+            game[(a-1)*scenario.B+b,(x-1)*scenario.Y+y]*tr(kron(Π_Ax[x][a],B_POVMs[y][b])*ρ_AB.M),
         1:scenario.B), 1:scenario.A), 1:scenario.Y), 1:scenario.X)
     ))
 
@@ -181,7 +186,7 @@ function _optimize_measurement_B(
         "game" => game,
         "scenario" => scenario,
         "state" => ρ_AB,
-        "A_POVMs" => A_POVMs,
+        "A_POVMs" => Π_Ax,
         "B_POVMs" => Π_B_opt
     )
 end
@@ -192,22 +197,23 @@ end
 function _optimize_measurement_A(
     scenario::BipartiteNonSignaling,
     game::BellGame,
-    ρ_AB :: State,
-    B_POVMs :: Vector{<:POVM}
+    ρ_AB :: AbstractMatrix,
+    B_POVMs :: Vector{<:AbstractVector{<:AbstractMatrix}};
+    atol=1e-7::Float64
 ) :: Dict
     if !(scenario.Y == length(B_POVMs))
         throw( DomainError(
             (scenario.Y, " != length(B_POVMs)"),
             "A distinct POVM is not specified for input `Y` of `BipartiteNonSignaling` scenario"
         ))
-    end
-
-    if !(all(i -> scenario.B == length(B_POVMs[i]), 1:scenario.Y ))
+    elseif !(all(i -> scenario.B == length(B_POVMs[i]), 1:scenario.Y ))
         throw( DomainError(
             B_POVMs, "Each POVM must have $(scenario.B) outputs."
         ))
     end
 
+    Π_By = map(Π -> Π isa POVM ? Π : POVM(Π, atol=atol), B_POVMs)
+
     ρ_dim = size(ρ_AB,1)
     Π_B_dim = size(B_POVMs[1][1],1)
 
@@ -222,7 +228,7 @@ function _optimize_measurement_A(
     # objective function
     problem = maximize(real(
         sum(x -> sum(y -> sum(a -> sum(b ->
-            game[(a-1)*scenario.B+b,(x-1)*scenario.Y+y]*tr(kron(A_POVMs[x][a],B_POVMs[y][b])*ρ_AB.M),
+            game[(a-1)*scenario.B+b,(x-1)*scenario.Y+y]*tr(kron(A_POVMs[x][a],Π_By[y][b])*ρ_AB.M),
         1:scenario.B), 1:scenario.A), 1:scenario.Y), 1:scenario.X)
     ))
 
@@ -248,14 +254,14 @@ function _optimize_measurement_A(
         "scenario" => scenario,
         "state" => ρ_AB,
         "A_POVMs" => Π_A_opt,
-        "B_POVMs" => B_POVMs
+        "B_POVMs" => Π_By
     )
 end
 
 # """
 # Helper function for converting POVM optimization variables to a POVM
 # """
-function _opt_vars_to_povm(povm::Vector{Matrix{Complex{Float64}}}; atol=1) :: POVM
+function _opt_vars_to_povm(povm::Vector{<:AbstractMatrix}; atol=1) :: POVM
     Π = is_povm(povm,atol=atol) ? POVM(povm,atol=atol) : throw(DomainError(povm, "not a povm"))
 
     Π
diff --git a/src/quantum_strategies.jl b/src/quantum_strategies.jl
index 3c02828..10bbc0f 100644
--- a/src/quantum_strategies.jl
+++ b/src/quantum_strategies.jl
@@ -4,41 +4,116 @@ export quantum_strategy
 Constructs a strategy matrix given quantum states and measurements. The supported
 scenarios include:
 
-`BlackBox` scenarios
+[`BlackBox`](@ref) scenarios:
 
     quantum_strategy(
-        Π :: POVM,
-        ρ_states :: Vector{<:State}
+        Π :: AbstractVector{<:AbstractMatrix},
+        ρ_states :: Vector{<:State};
+        atol=1e-7::Float64
     ) :: Strategy
 
-`LocalSignaling` scenarios
+For a quantum system the conditional proabilities are constructed as
+
+```math
+    P(y|x) = \\text{Tr}[\\Pi_y\\rho_x].
+```
+"""
+function quantum_strategy(
+    Π :: AbstractVector{<:AbstractMatrix},
+    ρ_states :: Vector{<:AbstractMatrix};
+    atol=1e-7::Float64
+) :: Strategy
+    povm = Π isa POVM ? Π : POVM(Π, atol=atol)
+    states = map(ρ -> ρ isa State ? ρ : State(ρ, atol=atol), ρ_states)
+
+    conditionals = measure(povm, states)
+    Strategy(conditionals.distribution, atol=atol)
+end
+
+"""
+[`LocalSignaling`](@ref) scenarios:
 
     quantum_strategy(
-        Π :: POVM,
-        ρ_states :: Vector{<:State},
-        scenario :: LocalSignaling
+        Π :: AbstractVector{<:AbstractMatrix},
+        ρ_states :: Vector{<:AbstractMatrix},
+        scenario :: LocalSignaling;
+        atol=1e-7::Float64
     ) :: Strategy
 
+For quantum systems the conditional probabilities are construct as
+
+```math
+    P(y|x) = \\text{Tr}[\\Pi_y\\rho_x].
+```
+
 A `DomainError` is thrown if the provided states and measurements are not compatible
 with the specified scenario.
 """
 function quantum_strategy(
-    Π :: POVM,
-    ρ_states :: Vector{<:State}
+    Π :: AbstractVector{<:AbstractMatrix},
+    ρ_states :: Vector{<:AbstractMatrix},
+    scenario :: LocalSignaling;
+    atol=1e-7::Float64
 ) :: Strategy
-    conditionals = measure(Π, ρ_states)
-    scenario = BlackBox(length(ρ_states), length(Π))
-    Strategy(conditionals, scenario)
+    if (size(Π[1]) != (scenario.d, scenario.d)) || (size(ρ_states[1]) != (scenario.d, scenario.d))
+        throw(DomainError((Π, ρ_states), "POVM or States are not dimension `d=$(scenario.d)`."))
+    end
+
+    povm = Π isa POVM ? Π : POVM(Π, atol=atol)
+    states = map(ρ -> ρ isa State ? ρ : State(ρ, atol=atol), ρ_states)
+
+    conditionals = measure(povm, states)
+    Strategy(conditionals.distribution, scenario, atol=atol)
 end
 
+"""
+[`BipartiteNonSignaling`](@ref) scenarios:
+
+    quantum_strategy(
+            ρ_AB :: AbstractMatrix,
+            Π_Ax :: Vector{<:AbstractVector{<:AbstractMatrix}},
+            Π_By :: Vector{<:AbstractVector{<:AbstractMatrix}},
+            scenario :: BipartiteNonSignaling;
+            atol::Float64 = 1e-7
+    )
+
+Constructs a strategy matrix in the generalized representation for the quantum
+system with conditional probabilities.
+
+```math
+    P(ab|xy) = \\text{Tr}[(\\Pi_a^x\\otimes\\Pi_b^y)\\rho_{AB}]
+```
+
+A `DomainError` is thrown if
+* The length of each POVM does not match the scenarios number of outputs
+* The number of each party's POVMS doesn't match the the number of inputs.
+"""
 function quantum_strategy(
-    Π :: POVM,
-    ρ_states :: Vector{<:State},
-    scenario :: LocalSignaling
+        ρ_AB :: AbstractMatrix,
+        Π_Ax :: Vector{<:AbstractVector{<:AbstractMatrix}},
+        Π_By :: Vector{<:AbstractVector{<:AbstractMatrix}},
+        scenario :: BipartiteNonSignaling;
+        atol::Float64 = 1e-7
 ) :: Strategy
-    if (size(Π[1]) != (scenario.d, scenario.d)) || (size(ρ_states[1]) != (scenario.d, scenario.d))
-        throw(DomainError((Π, ρ_states), "POVM or States are not dimension `d=$(scenario.d)`."))
+    if scenario.X != length(Π_Ax)
+        throw(DomainError(Π_Ax, "`length(Π_Ax) == $(scenario.X)` is not satisfied."))
+    elseif scenario.Y != length(Π_By)
+        throw(DomainError(Π_By, "`length(Π_By) == $(scenario.Y)` is not satisfied."))
+    elseif all(Π -> length(Π) != scenario.A , Π_Ax)
+        throw(DomainError(Π_Ax, "Each POVM in Π_Ax must have $(scenario.A) elements"))
+    elseif all(Π -> length(Π) != scenario.B, Π_By)
+        throw(DomainError(Π_By, "Each POVM in Π_By must have $(scenario.B) elements"))
     end
-    conditionals = measure(Π, ρ_states)
-    Strategy(conditionals, scenario)
+
+    povm_Ax = map(Π -> Π isa POVM ? Π : POVM(Π, atol=atol), Π_Ax)
+    povm_By = map(Π -> Π isa POVM ? Π : POVM(Π, atol=atol), Π_By)
+    state = ρ_AB isa State ? ρ_AB : State(ρ_AB, atol=atol)
+
+    strat_mat = zeros(Float64, (scenario.A,scenario.B,scenario.X,scenario.Y))
+    for a in 1:scenario.A, b in 1:scenario.B, x in 1:scenario.X, y in 1:scenario.Y
+        Π_AB = kron(povm_Ax[x][a].M, povm_By[y][b].M)
+        strat_mat[b,a,y,x] = tr(Π_AB * state)
+    end
+
+    Strategy(reshape(strat_mat, strategy_dims(scenario)), scenario, atol=atol)
 end
diff --git a/src/scenarios.jl b/src/scenarios.jl
index b12c195..960b8a4 100644
--- a/src/scenarios.jl
+++ b/src/scenarios.jl
@@ -95,7 +95,6 @@ struct BipartiteNonSignaling <: Scenario
         ) ? new(A, B, X, Y) : throw(
             DomainError((A, B, X, Y), "inputs must be ≥ 1")
         )
-
 end
 
 """
diff --git a/test/Project.toml b/test/Project.toml
index 1f3215b..911d5f0 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,7 +1,15 @@
 [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/unit/LocalPolytope.jl b/test/unit/LocalPolytope.jl
index 18cfe1a..2a3519b 100644
--- a/test/unit/LocalPolytope.jl
+++ b/test/unit/LocalPolytope.jl
@@ -1,4 +1,4 @@
-using Test, LinearAlgebra
+using Test, LinearAlgebra, Polyhedra, XPORTA
 
 @testset "test/unit/LocalPolytope.jl" begin
 
@@ -6,14 +6,38 @@ using BellScenario
 
 @testset "running all tests for module LocalPolytope.jl" begin
 
-  println("running LocalPolytope.jl unit tests:")
-  for test in readdir("./test/unit/LocalPolytope/")
-      # run only julia files in test directory
-      if occursin(r"^.*\.jl$", test)
-          println("./unit/LocalPolytope/$test")
-          @time include("./LocalPolytope/$test")
-      end
-  end
+    @testset "LocalPolytope.vrep()" begin
+        @testset "no arguments" begin
+            scenario = LocalSignaling(3,3,2)
+            local_poly = vrep(scenario)
+
+            @test local_poly isa XPORTA.Polyhedron
+            @test npoints(local_poly) == 21
+
+            vertices = map(p -> convert.(Int64,p), points(local_poly))
+            @test LocalPolytope.vertices(scenario) == vertices
+        end
+
+        @testset "passing arguments" begin
+            scenario = LocalSignaling(3,3,2)
+            local_poly = vrep(scenario, rank_d_only=true)
+
+            @test local_poly isa XPORTA.Polyhedron
+            @test npoints(local_poly) == 18
+
+            vertices = map(p -> convert.(Int64,p), points(local_poly))
+            @test LocalPolytope.vertices(scenario, rank_d_only=true) == vertices
+        end
+    end
+
+    println("running LocalPolytope.jl unit tests:")
+    for test in readdir("./test/unit/LocalPolytope/")
+    # run only julia files in test directory
+    if occursin(r"^.*\.jl$", test)
+        println("./unit/LocalPolytope/$test")
+        @time include("./LocalPolytope/$test")
+    end
+    end
 
 end
 
diff --git a/test/unit/LocalPolytope/facets.jl b/test/unit/LocalPolytope/facets.jl
index 1e55302..c0472e7 100644
--- a/test/unit/LocalPolytope/facets.jl
+++ b/test/unit/LocalPolytope/facets.jl
@@ -1,4 +1,4 @@
-using Test
+using Test, Polyhedra
 
 @testset "./src/LocalPolytope/facets.jl" begin
 
@@ -47,4 +47,21 @@ using BellScenario
     end
 end
 
+@testset "facets(::Polyhedron)" begin
+    scenario = LocalSignaling(3,3,2)
+    local_poly = LocalPolytope.vrep(scenario)
+
+    @test !hrepiscomputed(local_poly)
+
+    facets = LocalPolytope.facets(local_poly)
+    @test length(facets) == 15
+
+    @test hrepiscomputed(local_poly)
+
+    match_vertices = LocalPolytope.vertices(scenario)
+    match_facets = LocalPolytope.facets(match_vertices)["facets"]
+
+    @test facets == match_facets
+end
+
 end
diff --git a/test/unit/Nonlocality/optimize_measurements.jl b/test/unit/Nonlocality/optimize_measurements.jl
index 71e6591..ee326eb 100644
--- a/test/unit/Nonlocality/optimize_measurements.jl
+++ b/test/unit/Nonlocality/optimize_measurements.jl
@@ -5,6 +5,17 @@ using Test, QBase
 using BellScenario
 
 @testset "Nonlocality.optimize_measurement(LocalSignaling)" begin
+    @testset "non State types" begin
+        scenario = LocalSignaling(3,3,2)
+        game = BellGame([1 0 0;0 1 0;0 0 1], 2)
+        ρ_states = [[1 0;0 0],[0 0;0 1],[0 0;0 1]]
+
+        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]]
+    end
+
     @testset "trine states" begin
         scenario = LocalSignaling(3,3,2)
         game = BellGame([1 0 0;0 1 0;0 0 1], 2)
@@ -59,10 +70,38 @@ using BellScenario
 
         scenario = LocalSignaling(4,4,3)
         @test_throws DomainError  Nonlocality.optimize_measurement(scenario, game, states)
+
+        scenario = LocalSignaling(3,3,2)
+        game = BellGame([1 0 0;0 1 0;0 0 1],2)
+        states = [[1 0;0 1],[0 0;0 1],[0.5 0;0 0.5]]
+        @test_throws DomainError Nonlocality.optimize_measurement(scenario, game, states)
     end
 end
 
 @testset "Nonlocality.optimize_measurement(::BipartiteNonSignaling)" begin
+    @testset "matrix states and povms" begin
+        scenario = BipartiteNonSignaling(2,2,2,2)
+        game = BellGame([0 1 1 0;0 0 0 1;0 0 0 1;1 0 0 1],2)
+        ρ_AB = [1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1]/2
+        POVMs = [
+            [[1 0;0 0],[0 0;0 1]],
+            [[0.5 0.5;0.5 0.5],[0.5 -0.5;-0.5 0.5]]
+        ]
+
+        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 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 opt_dictB["B_POVMs"] isa Vector{<:POVM}
+        @test opt_dictB["state"] isa State
+    end
+
     @testset "CHSH inequality" begin
         scenario = BipartiteNonSignaling(2,2,2,2)
         game = BellGame([0 1 1 0;0 0 0 1;0 0 0 1;1 0 0 1],2)
@@ -127,6 +166,16 @@ end
 
         @test_throws DomainError Nonlocality.optimize_measurement(scenarioA, game, ρ_AB, A_POVMs=POVMs)
         @test_throws DomainError Nonlocality.optimize_measurement(scenarioB, game, ρ_AB, B_POVMs=POVMs)
+
+        @test_throws DomainError Nonlocality.optimize_measurement(
+            scenario, game, [1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1],A_POVMs=POVMs
+        )
+        @test_throws DomainError Nonlocality.optimize_measurement(
+            scenario, game, ρ_AB, A_POVMs = [[[1 0;0 1],[1 0;0 1]],[[1 0;0 0],[1 0;0 0]]]
+        )
+        @test_throws DomainError Nonlocality.optimize_measurement(
+            scenario, game, ρ_AB, B_POVMs = [[[1 0;0 1],[1 0;0 1]],[[1 0;0 0],[1 0;0 0]]]
+        )
     end
 end
 
diff --git a/test/unit/quantum_strategies.jl b/test/unit/quantum_strategies.jl
index 1320f17..55ed50e 100644
--- a/test/unit/quantum_strategies.jl
+++ b/test/unit/quantum_strategies.jl
@@ -5,33 +5,125 @@ using Test, QBase
 using BellScenario
 
 @testset "quantum_strategy()" begin
-    q_strat = quantum_strategy(
-        POVM([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]]),
-        State.([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]])
-    )
-
-    @test q_strat isa Strategy
-    @test q_strat.scenario == BlackBox(3,3)
-    @test q_strat == [1 0 0;0 1 0;0 0 1]
+    @testset "QBase types" begin
+        q_strat = quantum_strategy(
+            POVM([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]]),
+            State.([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]])
+        )
+
+        @test q_strat isa Strategy
+        @test q_strat.scenario == BlackBox(3,3)
+        @test q_strat == [1 0 0;0 1 0;0 0 1]
+    end
+
+    @testset "matrix types" begin
+        q_strat = quantum_strategy(
+            [[1 0;0 0],[0 0;0 1]],
+            [[1 0;0 0],[0 0;0 1]]
+        )
+
+        @test q_strat isa Strategy
+        @test q_strat.scenario == BlackBox(2,2)
+        @test q_strat == [1 0;0 1]
+    end
 end
 
 @testset "quantum_strategy(LocalSignaling)" begin
-    scenario = LocalSignaling(3,3,2)
-    q_strat = quantum_strategy(
-        mirror_symmetric_qubit_3povm(π/3),
-        trine_qubit_states(),
-        scenario
-    )
-
-    @test q_strat isa Strategy
-    @test q_strat ≈ [2/3 1/6 1/6;1/6 2/3 1/6;1/6 1/6 2/3]
-    @test q_strat.scenario isa LocalSignaling
-
-    @test_throws DomainError quantum_strategy(
-        POVM([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]]),
-        State.([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]]),
-        scenario
-    )
+    @testset "QBase types" begin
+        scenario = LocalSignaling(3,3,2)
+        q_strat = quantum_strategy(
+            mirror_symmetric_qubit_3povm(π/3),
+            trine_qubit_states(),
+            scenario
+        )
+
+        @test q_strat isa Strategy
+        @test q_strat ≈ [2/3 1/6 1/6;1/6 2/3 1/6;1/6 1/6 2/3]
+        @test q_strat.scenario isa LocalSignaling
+    end
+
+    @testset "Matrix types" begin
+        scenario = LocalSignaling(2,2,2)
+
+        q_strat = quantum_strategy(
+            [[1 0;0 0],[0 0;0 1]],
+            [[1 0;0 0],[0 0;0 1]],
+            scenario
+        )
+
+        @test q_strat isa Strategy
+        @test q_strat.scenario == LocalSignaling(2,2,2)
+        @test q_strat == [1 0;0 1]
+    end
+
+    @testset "DomainErrors" begin
+        scenario = LocalSignaling(3,3,2)
+        @test_throws DomainError quantum_strategy(
+            POVM([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]]),
+            State.([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]]),
+            scenario
+        )
+    end
+end
+
+@testset "quantum_strategy(BipartiteNonSignaling)" begin
+    @testset "QBase types ch violation" begin
+        scenario = BipartiteNonSignaling(2,2,2,2)
+        ρ_AB = State([1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1]/2)
+        Π_Ax = [
+            POVM([[1. 0;0 0],[0 0;0 1.]]),
+            POVM([[1 1;1 1]/2,[1 -1;-1 1]/2])
+        ]
+        Π_By = [
+            POVM([bloch_qubit_state(π/4),bloch_qubit_state(-3π/4)]),
+            POVM([bloch_qubit_state(3π/4),bloch_qubit_state(-π/4)])
+        ]
+
+        q_strat = quantum_strategy(ρ_AB, Π_Ax, Π_By, scenario)
+        @test q_strat isa Strategy
+        @test q_strat.scenario isa BipartiteNonSignaling
+
+        # generalized rep of CH inequality
+        g = BellGame([1 0 1 1;1 1 0 0;0 1 1 0;1 1 1 0], 3)
+        @test sum(g .* q_strat) ≈ 3.207 atol=2e-4
+    end
+
+    @testset "matrix types classical states" begin
+        scenario = BipartiteNonSignaling(2,2,2,2)
+        ρ_AB = kron([1 0;0 0],[0 0;0 1])
+        Π_Ax = [
+            [[1 0;0 0],[0 0;0 1]],
+            [[0 0;0 1],[1 0;0 0]]
+        ]
+        Π_By = [
+            [[1 0;0 0],[0 0;0 1]],
+            [[0 0;0 1],[1 0;0 0]]
+        ]
+
+        q_strat = quantum_strategy(ρ_AB, Π_Ax, Π_By, scenario)
+        @test q_strat isa Strategy
+        @test q_strat == [0 1 0 0;1 0 0 0;0 0 0 1;0 0 1 0]
+        @test q_strat.scenario == scenario
+    end
+
+    @testset "Domain Errors" begin
+        ρ_AB = State(kron([1 0;0 0],[0 0;0 1]))
+        Π = POVM([[1 0;0 0],[0 0;0 1]])
+        Π_Ax = [Π, Π]
+        Π_By = [Π, Π]
+
+        scenario = BipartiteNonSignaling(3,2,2,2)
+        @test_throws DomainError quantum_strategy(ρ_AB, Π_Ax, Π_By, scenario)
+
+        scenario = BipartiteNonSignaling(2,3,2,2)
+        @test_throws DomainError quantum_strategy(ρ_AB, Π_Ax, Π_By, scenario)
+
+        scenario = BipartiteNonSignaling(2,2,3,2)
+        @test_throws DomainError quantum_strategy(ρ_AB, Π_Ax, Π_By, scenario)
+
+        scenario = BipartiteNonSignaling(2,2,2,3)
+        @test_throws DomainError quantum_strategy(ρ_AB, Π_Ax, Π_By, scenario)
+    end
 end
 
 end