Nearly there!

Author: scheinerman

Project.toml

@@ -4,9 +4,13 @@ authors = ["Ed Scheinerman <ers@jhu.edu>"]
 version = "0.0.0"
 
 [deps]
+ChooseOptimizer = "858a232f-1959-5553-8cfc-91e1fd5304e2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
+ChooseOptimizer = "0.3"
 DataFrames = "1"
+JuMP = "1"
 LinearAlgebra = "1"

src/LPsolve.jl

@@ -0,0 +1,23 @@
+"""
+    lp_solve(T::Tableau)
+
+Use a standard linear programming solver [HiGHS by defaut] to find the optimal solution to the
+LP in `T`. Returns the values of the variables. 
+"""
+function lp_solve(T::Tableau)
+    MOD = Model(get_solver())
+    @variable(MOD, x[1:(T.n_vars)] >= 0)
+    @objective(MOD, Max, sum(T.c[i] * x[i] for i in 1:(T.n_vars)))
+    for i in 1:(T.n_cons)
+        @constraint(MOD, sum(T.A[i, j] * x[j] for j in 1:(T.n_vars)) <= T.b[i])
+    end
+
+    optimize!(MOD)
+    status = Int(termination_status(MOD))
+
+    if status ≠ 1
+        error("Linear program is either infeasible or unbounded. Status = $status")
+    end
+
+    return value.(x)
+end

src/Pivoting.jl

@@ -65,7 +65,7 @@ end
 
 """
     pivot(T::Tableau, i::Int, j::Int)
-    pivot!(T::Tableau)
+    pivot(T::Tableau)
 
 Non-modifying version of `pivot!`
 """

src/SimpleTableaux.jl

@@ -1,14 +1,25 @@
 module SimpleTableaux
 
+using ChooseOptimizer
 using DataFrames
+using JuMP
 using LinearAlgebra
 
 import DataFrames: DataFrame
 import Base: show
 
 TabEntry = Rational{BigInt}
 
-export Tableau, find_pivot, find_pivot_column, find_pivot_row, pivot, pivot!, restore
+export Tableau,
+    find_pivot,
+    find_pivot_column,
+    find_pivot_row,
+    lp_solve,
+    pivot,
+    pivot!,
+    pivot_solve,
+    pivot_solve!,
+    restore
 
 """
 Tableau(A::Matrix, b::Vector, c::Vector)
@@ -40,9 +51,11 @@ end
     restore(T::Tableau)
 
 Create a new `Tableau` based on the original data used to create `T`.
+
+Recommended usage: `T = restore(T)` to reset `T` to its original state. 
 """
 function restore(T::Tableau)
-    return Tableau(T.A, T.b, T.C)
+    return Tableau(T.A, T.b, T.c)
 end
 
 include("Exact.jl")
@@ -54,5 +67,7 @@ function show(io::IO, T::Tableau)
 end
 
 include("Pivoting.jl")
+include("Solver.jl")
+include("LPsolve.jl")
 
 end # module SimpleTableaux

src/Solver.jl

@@ -0,0 +1,73 @@
+"""
+    pivot_solve!(T::Tableau, verbose::Bool=true)
+
+Same as `pivot_solve` but this version modifies `T`.
+"""
+function pivot_solve!(T::Tableau, verbose::Bool=true)
+    count = 0
+    while true
+        if verbose
+            println("Iteration $count")
+            println(T)
+            println()
+            count += 1
+        end
+        i, j = find_pivot(T)
+        if i == 0 || j == 0
+            break
+        end
+        pivot!(T, i, j)
+    end # end while
+
+    if verbose
+        val = Exact(T.M[end,end])
+        println("Optimum value = $val")
+    end
+
+    return _get_x(T)
+end # end pivot_solve
+
+"""
+    pivot_solve(T::Tableau, verbose::Bool=true)
+
+Solve the linear program in `T` by pivoting. 
+With `verbose` set to `true` all steps of the solution are shown.
+
+Returns the optimal solution to the linear program. 
+"""
+function pivot_solve(T::Tableau, verbose::Bool=true)
+    TT = deepcopy(T)
+    return pivot_solve!(TT, verbose)
+end
+
+"""
+    _is_basis_vector(x::Vector)::Bool
+
+Determine if `x` is vector with a single entry equal to 1 and the rest equal to 0.
+"""
+function _is_basis_vector(x::Vector)::Bool
+    x = sort(x)
+    n = length(x)
+    y = zeros(Int, n)
+    y[end] = 1
+
+    return x == y
+end
+
+"""
+    _get_x(T::Tableau)
+
+Extract the values of the solution to `T` after finish pivoting.
+"""
+function _get_x(T::Tableau)
+    x = zeros(TabEntry, T.n_vars)
+
+    for i in 1:(T.n_vars)
+        col = T.M[1:(T.n_cons), i]
+        if _is_basis_vector(col)
+            j = findfirst(col .== 1)
+            x[i] = T.M[j, end]
+        end
+    end
+    return x
+end

test/examples.jl

@@ -23,3 +23,19 @@ function example2()
     c = [3; 4; 1]
     return Tableau(A, b, c)
 end
+
+function example3()
+    A = [1 1 3; 2 2 5; 4 1 2]
+    b = [30; 24; 36]
+    c = [3; 1; 1]
+    return Tableau(A, b, c)
+end
+
+function unbounded_example()
+    A = [-1 0; 0 -1]
+    b = [-1; -1]
+    c = [1; 1]
+    return Tableau(A, b, c)
+end
+
+nothing

test/runtests.jl

@@ -3,11 +3,12 @@ using SimpleTableaux
 using LinearAlgebra
 using DataFrames
 
-@testset "Dumb start" begin
+@testset "Simple start" begin
     A = [1 2 3 4; 5 6 7 8]
     b = [1; 3//2]
     c = [5; 6; 9; 2]
     T = Tableau(A, b, c)
-    pivot!(T)
-    @test true
+    x = pivot_solve!(T, false)
+
+    @test dot(c, x) == T.M[end, end]
 end