# Let's give an algorithm to solve an LP by direct enumeration.
using Combinatorics # for combinations
using LinearAlgebra

# Check if a given point is a basic feasible point.
function basic_feasible_point(A::Matrix,b::Vector,set::Vector)
    m,n = size(A)
    @assert length(set) == m "need more indices to define a BFP"
    binds = set # basic variable indices
    ninds = setdiff(1:size(A,1),binds) # non-basic
    B = A[:,binds]
    N = A[:,ninds]
    
    if rank(B) != m
        return (:Infeasible, 0)
    end
    
    xb = B\b
    x = zeros(eltype(xb),n)
    x[binds] = xb
    x[ninds] .= zero(eltype(xb))    
    
    if any(xb .< 0)
        return (:Infeasible, x)
    else
        return (:Feasible, x)
    end
end

function exhaustive_linprog(c::Vector, A::Matrix, b::Vector)
    m,n = size(A)
    # add some lines of code to track optimal
    for inds in combinations(1:n,m) # for all combinations
        bfp = basic_feasible_point(A,b,inds)
        if bfp[1] == :Feasible
            # compute the objective
            x = bfp[2]
            obj = (c'*x)[1]
            @show (obj, inds, x)
        end
    end
end

A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

A = [A1 Matrix{Float64}(I,3,3)] # form the problem with slacks.
exhaustive_linprog([-1.0,-2.0,0.0,0.0,0.0],A,b)

(obj, inds, x) = (-13.0, [1, 2, 3], [3.0, 5.0, 3.0, 0.0, 0.0])
(obj, inds, x) = (-9.0, [1, 2, 5], [1.0, 4.0, 0.0, 0.0, 2.0])
(obj, inds, x) = (-3.0, [1, 3, 4], [3.0, 0.0, 8.0, 10.0, 0.0])
(obj, inds, x) = (-4.0, [2, 4, 5], [0.0, 2.0, 0.0, 3.0, 3.0])
(obj, inds, x) = (0.0, [3, 4, 5], [0.0, 0.0, 2.0, 7.0, 3.0])