Phase one basis solver working

Author: scheinerman

development/random-examples.jl

@@ -1,18 +1,19 @@
 using SimplexTableaux, LinearAlgebra, ProgressMeter
 
 """
-    random_example(n_vars::Int, n_cons::Int, modulus::Int=10)
+    random_example(n_vars::Int, n_cons::Int, lo::Int, hi::Int, canonical::Bool=true)::Tableau
 
 Create a random instance of a `Tableau` with `n_vars` variables, 
 `n_cons` constraints in which all data are chosen IID uniformly
-from `{1,2,...,modulus}`.
+from `{lo, ..., hi}`. 
 """
-function random_example(n_vars::Int, n_cons::Int, modulus::Int=10)::Tableau
-    f(x::Int) = mod1(x, modulus)
-    A = f.(rand(Int, n_cons, n_vars))
-    b = f.(rand(Int, n_cons))
-    c = f.(rand(Int, n_vars))
-    return Tableau(A, b, c)
+function random_example(
+    n_vars::Int, n_cons::Int, lo::Int, hi::Int, canonical::Bool=true
+)::Tableau
+    A = rand(lo:hi, n_cons, n_vars)
+    b = rand(lo:hi, n_cons)
+    c = rand(lo:hi, n_vars)
+    return Tableau(A, b, c, canonical)
 end
 
 """

src/Bases.jl

@@ -3,18 +3,52 @@
 
 Return a feasible basis for the LP in `T` or `nothing` if none exists.
 """
-function find_a_basis(T::Tableau)::Vector{Int}
-    r, c = size(T.A)
-    for B in combinations(1:c, r)
-        TT = set_basis!(T, B)
-        if is_feasible(TT)
-            return B
+function find_a_basis(T::Tableau)
+    A, b, c = get_Abc(T)
+    m = T.n_cons
+    n = T.n_vars
+
+    # ensure RHS is nonnegative
+    for i in 1:m
+        if b[i] < 0
+            b[i] = -b[i]
+            A[i, :] = -A[i, :]
         end
     end
-    @info "No basis found"
-    return zeros(Int, r)
+
+    # glue an identity matrix to RHS
+    A = hcat(A, eye(Int, m))
+
+    # create artifical c 
+    c = vcat(0*c, ones(Int, m))
+
+    TT = Tableau(A, b, c, false)
+
+    B = collect((n + 1):(n + m))
+    set_basis!(TT, B)
+
+    simplex_solve!(TT, false)
+    v = value(TT)
+    if v>0
+        @info "No basis found."
+        return 0*get_basis(TT)
+    end
+
+    return TT.B
 end
 
+# function find_a_basis(T::Tableau)::Vector{Int}
+#     r, c = size(T.A)
+#     for B in combinations(1:c, r)
+#         TT = set_basis!(T, B)
+#         if is_feasible(TT)
+#             return B
+#         end
+#     end
+#     @info "No basis found"
+#     return zeros(Int, r)
+# end
+
 """
     find_all_bases(T::Tableau)
 

src/LPsolve.jl

@@ -23,7 +23,7 @@ function lp_solve(T::Tableau, verbose::Bool=true)
     end
     if verbose
         obj_val = objective_value(MOD)
-        println("Optimum value = $obj_val\n")
+        println("Minimial objective value = $obj_val\n")
     end
 
     xval = JuMP.value.(x)

src/Solver.jl

@@ -120,7 +120,7 @@ function simplex_solve!(T::Tableau, verbose::Bool=true)
 
     if verbose
         println("Optimality reached. Pivot count = $pivot_count")
-        println("Value = $v = $(Float64(value(T)))")
+        println("Minimal value = $v = $(Float64(value(T)))")
     end
 
     return x