An extremely naive implementation of the simplex algorithm.

simplex.py

from collections import namedtuple
import numpy as np

# This is an implementation of the version of the simplex algorithm we are
# expected to execute by hand in LPMS. There is just one important difference:
# because python and numpy use zero-based indexing, this does too. So for
# example you need to subtract 1 from each element of your basic and nonbasic
# sets.

LPState = namedtuple("LPState", ["basic", "nonbasic"])
LPSolution = namedtuple("LPSolution", ["basic_variables", "basic_variable_values", "objective"])

class LP:
    def __init__(self, costs, A, b):
        """
        Constructs a linear programming problem of the form

            maximize f = (costs)^T . x
            subject to A . x <= b
                       x >= 0

        with x being the vector of decision variables.
        """
        self._check_dims(costs, A, b)
        self.costs = costs
        self.A = A
        self.b = b

        # Computed properties
        self.num_constraints = A.shape[0]
        self.num_variables = A.shape[1]
        self.Abar = np.concatenate([A, np.identity(self.num_constraints)], axis=1)
        self.full_costs = np.concatenate([self.costs, np.repeat(0, self.num_constraints)])

    def _check_dims(self, costs, A, b):
        if len(costs.shape) != 1:
            raise ValueError("Expected costs to be a vector")
        if len(A.shape) != 2:
            raise ValueError("Expected A to be a matrix")
        if len(b.shape) != 1:
            raise ValueError("Expected b to be a vector")

        if A.shape[1] != costs.shape[0]:
            raise ValueError("Size mismatch between A and costs")

        if A.shape[0] != b.shape[0]:
            raise ValueError("Size mismatch between A and b")

    def __repr__(self):
        def show_constraint(j):
            return "  " + " + ".join([
                repr(coeff) + "*x_" + repr(i) for (i, coeff) in enumerate(self.A[j])
                ]) + " <= " + repr(self.b[j])

        objective_fn = " + ".join([
            repr(c) + "*x_" + repr(i) for (i, c) in enumerate(self.costs)
            ])
        constraints = [show_constraint(i) for i in range(0, self.num_constraints)]
        return "\n".join([
            "maximize f = " + objective_fn,
            "subject to:"] +
            constraints +
            ["  x_i >= 0"]
            )

    def initial_state(self):
        # TODO: ensure this is feasible
        m = self.num_constraints
        n = self.num_variables
        return LPState(
                basic=range(n, n+m),
                nonbasic=range(0, n)
                )

    def dual(self):
        return LP(-self.b,
                -np.transpose(self.A),
                -self.costs
                )


eg2 = LP(
    # costs
    np.array([1, 2, 3]),
    # A
    np.array([[1, 0, 2],
              [0, 1, 2]]),
    # b
    np.array([3, 2])
    )

eg2_state = LPState(
        basic=[0,2],
        nonbasic=[1,3,4]
        )

assignment2 = LP(
    # costs
    np.array([1, 10]),
    # A
    np.array([[1, 20],
              [0, 1 ]]),
    # b
    np.array([100, 1])
    )

def iterate(lp, state, verbose=True):
    """Returns an LPSolution if the given state comprises a feasible optimal
    basis, a new LPState if it is feasible but non-optimal, and throws an error
    if the state is non-feasible or the problem is unbounded."""
    def info(msg):
        if verbose:
            print(msg)

    B = lp.Abar[:,state.basic]
    N = lp.Abar[:,state.nonbasic]

    info("B = {}".format(B))
    info("N = {}".format(N))

    bhat = np.linalg.solve(B, lp.b)

    info("bhat = {}".format(bhat))

    if not is_feasible(bhat):
        raise ValueError("unfeasible basis; bhat = {}".format(bhat))

    basic_costs = lp.full_costs[state.basic]
    nonbasic_costs = lp.full_costs[state.nonbasic]

    info("basic_costs = {}".format(basic_costs))
    info("nonbasic_costs = {}".format(nonbasic_costs))

    pi = np.linalg.solve(B.T, basic_costs)

    info("pi = {}".format(pi))

    reduced_costs = nonbasic_costs - np.matmul(N.T, pi)

    info("reduced_costs = {}".format(reduced_costs))

    if is_optimal(reduced_costs):
        info("Optimal!")
        return LPSolution(
                basic_variables=state.basic,
                basic_variable_values=bhat,
                objective=np.dot(basic_costs, bhat)
                )

    q = np.argmax(reduced_costs)
    qp = state.nonbasic[q]
    # q-th column of N
    a_q = N[:,q]

    info("q = {}".format(q))
    info("qp = {}".format(qp))
    info("a_q = {}".format(a_q))

    a_q_hat = np.linalg.solve(B, a_q)
    info("a_q_hat = {}".format(a_q_hat))

    if is_unbounded(a_q_hat):
        raise ValueError("Unbounded!")

    possible_ps = np.where(a_q_hat > 0, bhat / a_q_hat, np.infty)

    info("possible_ps = {}".format(possible_ps))

    p = np.argmin(possible_ps)
    pp = state.basic[p]

    info("p = {}".format(p))
    info("pp = {}".format(pp))

    new_nonbasic = list(state.nonbasic)
    new_nonbasic[q] = pp

    new_basic = list(state.basic)
    new_basic[p] = qp

    return LPState(
            nonbasic=new_nonbasic,
            basic=new_basic
            )

def is_feasible(bhat):
    return (bhat > 0).all()

def is_optimal(reduced_costs):
    return (reduced_costs <= 0).all()

def is_unbounded(a_q_hat):
    return (a_q_hat <= 0).all()

def main():
    "This runs through 1.a.ii of May 2016 paper."
    may2016_q1 = LP(
        # costs
        np.array([3, 1, 2]),
        # A
        np.array([[1, -2, 4],
                  [2, -3, -1]]),
        # b
        np.array([4, 12])
    )
    initial_state = LPState(
        basic=[0,4],
        nonbasic=[1,2,3]
        )

    new_state = iterate(may2016_q1, initial_state)
    iterate(may2016_q1, new_state)

if __name__ == '__main__':
    main()