diff --git a/Project.toml b/Project.toml
index d4c0302..a44307f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Vincent Du"]
 version = "0.1.0"
 
 [deps]
+CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
@@ -24,9 +25,11 @@ ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 SBMLToolkit = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 Satisfiability = "160ab843-0bc6-4ba4-9585-b7478b70f443"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
+CDDLib = "0.9.4"
 Catalyst = "14.4.1"
 DynamicPolynomials = "0.6.0"
 Graphs = "1.12.0"
@@ -40,20 +43,22 @@ ModelingToolkit = "9.46.1"
 MultivariatePolynomials = "0.5.7"
 Nemo = "0.47"
 Oscar = "1.1.1"
+Polyhedra = "0.7.8"
 PolynomialRoots = "1.0.0"
 ReactionNetworkImporters = "0.15.1"
 SBMLToolkit = "0.1.29"
 SafeTestsets = "0.1.0"
 Satisfiability = "0.2"
 SparseArrays = "1"
+StaticArrays = "1.9.8"
 SymbolicUtils = "3.7.2"
 julia = "1.10"
 
 [extras]
 SBMLImporter = "210efffb-c3c8-456d-a807-6f55560b12fe"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [targets]
diff --git a/src/CatalystNetworkAnalysis.jl b/src/CatalystNetworkAnalysis.jl
index e6c2f1f..866fc37 100644
--- a/src/CatalystNetworkAnalysis.jl
+++ b/src/CatalystNetworkAnalysis.jl
@@ -11,6 +11,9 @@ using SparseArrays
 using IterTools
 using StaticArrays 
 
+using Polyhedra
+import CDDLib
+
 import ModelingToolkit as MT
 
 const M::Float64 = 1E4
diff --git a/src/fluxbalanceanalysis.jl b/src/fluxbalanceanalysis.jl
deleted file mode 100644
index 6c2a977..0000000
--- a/src/fluxbalanceanalysis.jl
+++ /dev/null
@@ -1,19 +0,0 @@
-"""
-    elementaryfluxmodes(rn::ReactionSystem)
-
-    Given a reaction network, return the set of elementary flux modes of the reaction network. 
-"""
-function elementaryfluxmodes(rn::ReactionSystem)
-    S = netstoichmat(rn)
-    m, n = size(S)
-    hyperplanes = [Polyhedra.HyperPlane(S[i, :], 0) for i in 1:m]
-    halfspaces = [Polyhedra.HalfSpace(-I(n)[i, :], 0) for i in 1:n]
-    polycone = Polyhedra.polyhedron(hrep(hyperplanes, halfspaces))
-
-    Polyhedra.vrep(polycone)
-    Polyhedra.removevredundancy!(polycone)
-
-    EFMs = reduce(hcat, map(x->x.a, polycone.vrep.rays.rays))
-end
-
-
diff --git a/src/lp_utils.jl b/src/lp_utils.jl
index 6808b3f..182e968 100644
--- a/src/lp_utils.jl
+++ b/src/lp_utils.jl
@@ -113,3 +113,21 @@ function is_extreme_idxset(S::M, idxs::Vector{Int}) where {M <: AbstractMatrix}
     cone_mat_eq = @view cone_mat[rowidxs, :]
     return rank(cone_mat_eq) == n - 1
 end
+
+"""
+    elementaryfluxmodes(rn::ReactionSystem)
+
+    Given a reaction network, return the set of elementary flux modes of the reaction network. 
+"""
+function elementary_flux_modes(rn::ReactionSystem)
+    S = netstoichmat(rn)
+    m, n = size(S)
+    hyperplanes = [Polyhedra.HyperPlane(S[i, :], 0) for i in 1:m]
+    halfspaces = [Polyhedra.HalfSpace(-I(n)[i, :], 0) for i in 1:n]
+    polycone = Polyhedra.polyhedron(hrep(hyperplanes, halfspaces))
+
+    Polyhedra.vrep(polycone)
+    Polyhedra.removevredundancy!(polycone)
+
+    EFMs = reduce(hcat, map(x->x.a, polycone.vrep.rays.rays))
+end
diff --git a/src/persistence.jl b/src/persistence.jl
index a719783..6eb9b3a 100644
--- a/src/persistence.jl
+++ b/src/persistence.jl
@@ -4,7 +4,6 @@
 #   permanence in networks
 #   improve efficiency of siphon detection
 
-
 """
     ispersistent(rs::ReactionSystem)
 
@@ -94,9 +93,9 @@ function minimalsiphons_smt(rs::ReactionSystem)
     return removesupersets(siphons)
 end
 
-function removesupersets(indexsets::AbstractArray{<:Array})
+function removesupersets(indexsets)
     indexsets = sort(indexsets, by=x->length(x))
-    minimalsets = Array[]
+    minimalsets = Array{Int64}[]
 
     for s in indexsets
         if !any(ms->issubset(ms, s), minimalsets)
@@ -148,5 +147,5 @@ function iscritical(siphon::Vector, S::Matrix)
 
     # If there is a non-negative vector in the nullspace of S_red', then there is a positive conservation law with a support that is the subset of the siphon, and the siphon is not critical 
     conslaws_r = conservationlaws(S_r)
-    !haspositivesolution(copy(conslaws_r'), nonneg = true)
+    !has_positive_solution(copy(conslaws_r'), nonneg = true)
 end
diff --git a/src/translated.jl b/src/translated.jl
index a977030..0c69590 100644
--- a/src/translated.jl
+++ b/src/translated.jl
@@ -7,8 +7,8 @@
 struct Translation{T<:Int} 
     rn::ReactionSystem
     Y_T::Matrix{T} # The new complex-composition matrix.
-    δ_K::T
-    δ_E::T
+    δ_K::T # Kinetic deficiency 
+    δ_E::T # Effective deficiency
 end
 
 """
@@ -35,10 +35,12 @@ function WRDZ_translation(rn::ReactionSystem)
 
         # Update reactant complexes
         @label step_3
-        rxs = findall(==(1), rr_adj[i, :])
+        rxs = findall(==(1), @view rr_adj[i, :])
         for j in rxs
             p, s = (dst(e[i]), src(e[j]))
-            (@view Y_T[:, j]) .= Y_K[:, p] - Y_K[:, s] + Y_T[:, i]
+            for k in 1:size(Y_T, 1)
+                Y_T[k, j] = Y_K[k, p] - Y_K[k, s] + Y_T[k, i]
+            end
             P3 = setdiff(push!(P3, j), P2)
         end
 
@@ -113,7 +115,7 @@ function rrgraph(rn::ReactionSystem)
     end
 
     optimize!(model)
-    is_solved_and_feasible(model) ? return edge : return
+    is_solved_and_feasible(model) && return edge 
 end
 
 # Generate partitions of the reaction-to-reaction graph
diff --git a/src/utils.jl b/src/utils.jl
index 4119e40..2e2c7eb 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -18,7 +18,7 @@ end
 """
     removespec(rn, spec)
 
-    Return the net stoichiometric matrix and incidence matrix obtained when removing a species at index spec from the reaction network. 
+    Return the net stoichiometric matrix and incidence matrix obtained when removing a species at index spec from the reaction network.
 """
 function removespec(rn::ReactionSystem, spec::Int64) 
     s, r = size(netstoichmat(rn))
@@ -39,9 +39,8 @@ end
 """
     transitiveclosure(arr, relation)
 
-    Given an iterable array and a (symmetric, reflexive) relation (a function that takes two objects, and returns true if they belong to the relation), return the partitions of the set under the transitive closure of the given relation. 
+    Given an iterable array and a (symmetric, reflexive) relation (a function that takes two objects, and returns true if they belong to the relation), return the partitions of the set under the transitive closure of the given relation.
 """
-
 function transitiveclosure(arr, relation) 
     adjmat = SparseArrays.spzeros(Bool, length(arr), length(arr))
     for idx in CartesianIndices(adjmat)
diff --git a/test/lp_utils.jl b/test/lp_utils.jl
new file mode 100644
index 0000000..4ddf54e
--- /dev/null
+++ b/test/lp_utils.jl
@@ -0,0 +1,29 @@
+# Testing of the linear programming utilities in CatalystNetworkAnalysis
+
+import CatalystNetworkAnalysis as C
+# Test that matrices with positive element in their image space are identified.
+let
+    a = rand(1:10, 10, 10)
+    @test has_positive_solution(a)
+    
+    b = -rand(1:10, 10, 10)
+    @test has_positive_solution(b)
+
+    a = [1 -1; -1 1]
+    @test !has_positive_solution(b)
+
+    b = [-1 1 1 1;
+         1 -1 1 1;
+         1 1 -1 1
+         1 1 1 -1]
+    @test !has_positive_solution(b)
+    @test has_positive_solution(b, nonneg=true)
+end
+
+# Testing that extreme rays are properly identified.
+let
+end
+
+# Testing elementary flux modes are found.
+let
+end