Simplex Algorithm implementation in Python (3.x)

simplex.py

import numpy as np


def simplex(c, A, b, basis):
  z = 0
  tableau = None
  foundOptimal = False
  numRound = 1
  # we'll be using tableau because it's easier and makes more sense.

  # using the basis, check if A is in canonical form.
  # i.e: the basis variables used create the identity matrix.
  if isCanonicalForm(A,basis):
    # convert c, A, and b to tableau form.
    tableau = componmentsToTableau(-c,A,b)
  
  while not foundOptimal:
    print("Round:", numRound)
    numRound += 1
    # to get the entering variable, look @ the negative entries in
    # the top row (which is c). Use Blands rule to get the entering
    # variable.
    enteringVariable = getEnteringVariable(tableau)
    # now we know the column to go to, we need to look into what
    # row to pivot from.
    pivotLocation = findMin(A, enteringVariable, b)
    if pivotLocation == -1:
      print("The solution is unbounded..")
      return -1

    print("Pivot Location [regards to A]: (" + str(enteringVariable+1)+","+str(pivotLocation+1)+")")
    # now we can pivot.
    tableau = pivot(tableau,pivotLocation+2,enteringVariable+2)
    print("new Tableau:\n",tableau, "\n")
    # extract the components of the equation
    (A,b,c,z) = tableauToComponments(tableau)
    # find the leaving column. We can tell which one is the
    # leaving column by looking at the columns that make up A
    # and removing ones that don't make the basis.
    leavingVariable, basis = findLeavingColumn(A,basis)

    print("Entering, Leaving Columns (for Bases):", enteringVariable+1, leavingVariable)
    # at this stage we've just finished a round of simplex.
    # check if C is all positive. if they are, then we've found the optimal
    # solution.
    if min(c) >= 0: foundOptimal = True
    
    print("---------------------------------Best: ", z)

  
  print()
  print("Found the (optimal?) solution:")
  print("z = ", z)
  print("Where x = ", produceX(A,basis,b))

  return False

def produceX(A,basis,b):
  # use the basis to produce the columns in A
  # iterate through the basis (it starts from 1 so -1 to them)
  basis = [x-1 for x in basis]

  # find how many items go in x.
  numberOfXValues = A.shape[1]
  result = []
  # iterate through the indexes
  for x in range(numberOfXValues):
    if x in basis:
      pos = A[:,x].T.tolist()[0].index(1)
      # add b's value to x
      result.append(b.tolist()[0][pos])
    else:
      # fill it with a zero.
      result.append(0.0)

  return result

def getEnteringVariable(t):
  # repull c.
  c = t[0].tolist()[0][1:-1]
  # using bland's rule (Choose the entering variable k with
  #  the SMALLEST INDEX) and the negative entries, get the
  #  first one that pops up.
  for i in range(len(c)):
    if c[i] < 0:
      return i
  # otherwise we can't find a solution.
  return -1

def isCanonicalForm(A,basis):
  """
  Check if the basis variables actually create
  the identity matrix.
  """
  # get the basic variables X_b
  X_b = []
  for b in basis:
    X_b.append(np.asarray(np.transpose(A[:,b-1])).tolist()[0])
  # convert X_b into a numpy matrix
  X_b = np.matrix(X_b)
  
  # get the size of the basis to create an identity matrix.
  dim = X_b.shape[0]
  I = np.eye(dim)

  # check if they're the same.
  if np.array_equal(X_b, I):
    return True
  else:
    return False

def findMin(A, ent, b):
  # convert those subsections into lists
  # so we can traverse them easier.
  a_col = A[:,ent].T.tolist()[0]
  b_cop = b.copy().tolist()[0]

  negatives = []
  contenders = []
  for i in range(len(a)):
    division = b_cop[i]/a_col[i]
    # if it's a negative value then we mess up the number
    # (we want the smallest positive)
    if division < 0:
      division = 2036854775807
      negatives.append(division)
    contenders.append(division)
  # if they're all negative then our problem is unbounded.
  if len(negatives) == len(contenders):
    return -1
  # now we want the index of the smallest positive number.
  index = contenders.index(min(contenders))
  return index

def componmentsToTableau(c,A,b,z=None):
  """
  Converts the matrix into an augmented matrix
  based on the linear system z-c^(T)x = z'
  """
  tableau = None

  # create top
  top = [1.0]+np.transpose(c).tolist()[0]+[0.0]
  top = np.matrix(top)

  # create zeros
  lhs_zero = np.array([np.zeros(A.shape[0])]).T
  # transpose b
  b_t = b.T
  # initiate bottom variable
  bottom = np.asarray(A.copy())
  # concatenate those zeros and bottom
  bottom = np.hstack((lhs_zero, bottom))
  # concatenate bottom and b
  bottom = np.hstack((bottom,b_t))
  # now we concatenate top and bottom
  tableau = np.vstack((top, bottom))
  return tableau

def tableauToComponments(tableau):
  """
  Converts the tableau to regular components.
  """
  # let's create A
  A = tableau[1:]
  A = np.delete(A, 0, 1)
  A = np.delete(A, A.shape[1]-1, 1)
  # creating b
  b = tableau[:,-1][1:].T
  # creating c
  c = tableau[0].tolist()[0][1:-1]
  # creating z
  z = tableau[0,-1]
  return (A,b,c,z)

def findLeavingColumn(A, basis):
  """
  We'll use bland's other rule here
  i.e we get the smallest index first; or rather
  the first thing we find that isn't to scratch will suffice.
  """

  # iterate through the columns and see if we can
  # find the identity matrix. 
  new_basis = []
  for i in range(A.shape[1]):
    # get the column
    column = A[:,i].T.tolist()[0]
    count_zero = 0
    int_is_one = False
    for j in column:
      if j == 0: count_zero += 1
      else:
        if j == 1: int_is_one = True
    if int_is_one and count_zero == len(column)-1:
      new_basis.append(i+1)

  # create placeholder variable for our leaving
  # column
  odd_one_out = -1
  for i in basis:
    # this oneliner finds the first leaving value.
    if i not in new_basis: odd_one_out = i; break

  return (odd_one_out, new_basis)



def pivot(a, i, j):
  """
  Remember with elementary row operations you can do the following:
  1. Multiply a row with a constant.
  2. Add a row * constant to another row.
  3. Swap the positions of two rows.

  Pivoting on an entry T_(i,j) != 0 means transforming the matrix T into 
  a matrix T' by performing the following operations:
  - multiply the ith row of T by 1/T_(i,j) (to make T_(i,j) == 1)
  - For each k != i, add the ith row * constant to the kth row. (to 
    make T'_(k,j) = 0)
  """
  # make sure we simplify the code s.t i and j refer to np indexes.
  i, j = i - 1, j - 1
  # make a copy of the matrix we're going to manipulate
  T = a.copy()
  # get the pivot
  pivot = T[i,j]
  # 1. multiply the ith row of T by 1/T_(i,j) (to make T_(i,j) == 1)
  divider = 1/pivot
  T[i] = np.multiply(T[i], divider)


  # 2. For each k != i, add the ith row * -constant to the kth row. 
  #    (to make T'_(k,j) = 0)
  for k in range(T.shape[0]):
    if k != i:
      # find the negative constant
      constant = -T[k,j]
      multiplier = np.multiply(T[i],constant)
      T[k] = np.add(T[k], multiplier)
  return T

if __name__ == "__main__":
  # Bounded Example:
  c = np.matrix([[2],[3],[0],[0],[0]], dtype=float)
  a = np.matrix([
    [1,1,1,0,0],
    [2,1,0,1,0],
    [-1,1,0,0,1]
  ], dtype=float)
  b = np.matrix([6,10,4])
  basis = [3,4,5]
  # -----------

  # Unbounded Example:
  # c = np.matrix([[-1],[3],[0],[0],[1]],dtype=float)
  # a = np.matrix([
  #   [-2,4,1,0,1],
  #   [-3,7,0,1,1]
  # ], dtype=float)
  # b = np.matrix([1,3])
  # basis = [3,4]

  simplex(c, a, b, basis)

I've verified the program to detect unbounded examples; that being said I've only tested it on two LPs. I'll update it in the event that other LPs break it.