L P

# solves *bounded* LPs of the form:
        # max cx
        # sub to: Ax <= b
from sympy import *
from itertools import combinations
        # enumerates all the vertices of {x | Ax <= b}
def enumeratevertices(A, b):
    m, n = A.rows, A.cols     
    for rowlist in combinations(range(m), n):
        Ap = A.extract(rowlist, range(n))
        bp = b.extract(rowlist, [0])
        if Ap.det() != 0:
            xp = Ap.LUsolve(bp)
            d = A * xp - b
            feasible = True
            for i in range(m):
                if d[i] > 0:
                    feasible = False
            if feasible:
                yield xp       
        # finds the optimum using vertex enumeration
def findoptimum(A, b, c):
    m, n = A.rows, A.cols
    bestvalue, bestvertex = None, None
    for vertex in enumeratevertices(A, b):
        if bestvalue is None or (vertex.T*c)[0] > bestvalue:
            bestvalue = (vertex.T * c)[0]
            bestvertex = vertex
    return bestvertex
def solve(A, b, c):
    x = findoptimum(A, b, c)
    if not x:
        print 'LP is infeasible'
    else:
        print 'Vertex', x.T, 'is optimal'
        print 'Optimal value is', c.T*x
if __name__ == '__main__':
    A = Matrix([[-10,  -6, -9, -10],
                [  8,  -6, -5,  -5],
                [ -7,  -1, -9,   3],
                [ -1,  -4,  5,  10],
                [  1,   2,  0,  10],
                [  2,  -9,  3,  -8],
                [ -8,  -1, -8,   1],
                [  7, -10,  4,  -4],
                [-10,   2,  5,   8],
                [ -7,   9,  4,  -4],
                [ -1,   0,  0,   0],
                [  0,  -1,  0,   0],
                [  0,   0, -1,   0],
                [  0,   0,  0,  -1]])
    b = Matrix([9, 7, 3, 4, 8, 0, 3, 2, 4, 8, 0, 0, 0, 0])
    c = Matrix([2, -2, -3, 8])
    solve(A, b, c)