simplex algorithm

simplex.py

from __future__ import division
try:
    from numpypy import *
except ImportError:
    from numpy import *
 
class Tableau:

    def __init__(self, obj):
        self.obj, self.rows, self.cons = obj, [], []

    def add_constraint(self, expression, value):
        self.rows.append(expression)
        self.cons.append(value)

    def display(self):
        print vstack((self.obj, self.tab))

    def _pivot(self, c, r):
        self.tab[r] /= self.tab[r][c]
        self.obj =  self.obj - self.obj[c] * self.tab[r]
        for i in xrange(len(self.rows)):
            if i!=r and self.tab[i][c]:
                self.tab[i] = self.tab[i] - self.tab[i][c] * self.tab[r]

    def solve(self):
        self.obj = array(self.obj + [0]*(len(self.rows)+1))
        self.tab = vstack(array(r, dtype=float) for r in self.rows)
        self.tab = hstack((self.tab, eye(len(self.rows))))
        self.tab = hstack((self.tab, array(self.cons, dtype=float).reshape(-1,1)))
        #self.display()
        
        while self.obj[:-1].min()<0:
            c = self.obj[:-1].argmin()
            r = (self.tab[:,-1] / self.tab[:,c]).argmin()
            self._pivot(c, r)
            #print '\ncolumn:{} row:{}'.format(c, r+1)
            #self.display()

        return vstack((self.obj, self.tab))


if __name__ == '__main__':
    """
    max z = 2x + 3y + 2z
    st
    2x + y + z <= 4
    x + 2y + z <= 7
    z          <= 5
    x,y,z >= 0
    """
 
    t = Tableau([-2,-3,-2])
    t.add_constraint([2, 1, 1], 4)
    t.add_constraint([1, 2, 1], 7)
    t.add_constraint([0, 0, 1], 5)
    print t.solve()