ILP solver

ilp.py

import sys
import numpy as np


class Simplex:

    def __init__(self):
        self.eps = 1e-9
        self.DEBUG = 0

    def parse_text(self, inputfile):
        text = open(inputfile).read().strip().split('\n')
        text = [l.split() for l in text]

        self.m, self.n = map(int, text[0])
        self.B = np.array(map(int, text[1]))
        self.N = np.array(map(int, text[2]))
        b = np.array(map(float, text[3]))
        A = np.array([map(float, l) for l in text[4:4 + self.m]])
        t = map(float, text[4 + self.m])
        z = t[0]
        c = np.array(t[1:])

        self.A = np.hstack((b[:, np.newaxis], A))
        self.c = np.hstack((z, c))

    def output(self):
        m, n, B, N, A, c = self.m, self.n, self.B, self.N, self.A, self.c
        sol = {}
        for i, v in enumerate(B):
            sol[v] = A[i][0]
        for j, v in enumerate(N):
            if c[j + 1] < 0:  # +1 fix gap between c and B
                sol[v] = 0.0
        print 'Solution:'
        print 'max z = {}'.format(c[0])
        print 'with'
        for v in sorted(sol.keys()):
            print 'x{} = {}'.format(v, sol[v])

    def inspect(self, title='DEBUG'):
        if not self.DEBUG:
            return
        print '-' * 30, title, '-' * 30
        print self.A
        print 'c', self.c
        print 'm:{} n:{}'.format(self.m, self.n)
        print 'B', self.B
        print 'N', self.N
        print

    def _pivot(self, j, i):
        B, N, A, c = self.B, self.N, self.A, self.c

        # adjust N, B
        N[j - 1], B[i] = B[i], N[j - 1]  # -1 fix gap between c and B

        # adjust A[i, :]
        pivot = A[i][j]
        A[i] /= -pivot
        A[i][j] = 1 / pivot

        # adjust c
        p = c[j]
        c += p * A[i]
        c[j] = p / pivot

        # adjust A
        for k in xrange(self.m):
            if k != i and A[k][j] != 0:
                p = A[k, j]
                A[k] += p * A[i]
                A[k, j] = p / pivot

    def finalize(self, init=False):
        while 1:
            if np.all(self.c[1:] < self.eps):  # finished
                return 'ok'

            # select column j (x0 first)
            if init and 0 in self.N:
                j = np.where(self.N == 0)[0][0]  # N.index(0)
            else:
                ok = np.where(self.c[1:] > self.eps)[0]
                if len(ok) == 1:  # fix error when len(ok)==1
                    j = ok[0]
                else:
                    order = np.argsort(self.N[ok])
                    j = ok[order][0]
                # fix gap by z
                j += 1

            # select row i
            col, b = self.A[:, j], self.A[:, 0]
            ok = np.where(col < -self.eps)[0]
            if len(ok) == 0:
                return 'unbounded'
            i = ok[np.argmax(b[ok] / col[ok])]

            # do the pivot
            self._pivot(j, i)
            self.inspect('pivoting')

    def need_init(self):
        return self.A[:, 0].min() < 0

    def initialize(self):
        # add x0 column (all 1) in A
        self.A = np.hstack([self.A, np.ones((self.m, 1))])

        # change opt object
        self.old_N = self.N
        self.old_c = self.c
        self.c = np.zeros(self.n + 2)
        self.c[-1] = -1

        # change record
        self.N = np.append(self.N, 0)

        # magic pivot
        j, i = self.n + 1, self.A[:, 0].argmin()
        self.inspect('aux')
        self._pivot(j, i)
        self.inspect('magic pivot')

    def de_initialize(self):
        # reconstruct c (lecture 69)
        idx0 = np.where(self.N == 0)[0][0]
        self.A = np.hstack((self.A[:, :idx0 + 1], self.A[:, idx0 + 2:]))
        self.N = np.hstack((self.N[:idx0], self.N[idx0 + 1:]))
        c = np.zeros(self.n + 1)
        for i, n in enumerate(self.old_N):
            if n not in self.B:
                idx = np.where(self.N == n)[0][0]
                c[idx + 1] += self.old_c[i + 1]
            else:
                idx = np.where(self.B == n)[0][0]
                c += self.old_c[i + 1] * self.A[idx, :]
        self.c = c

    def is_int(self, n):
        return abs(n - round(n)) <= self.eps

    def add_cutting_plane(self):
        B, N, A, c = self.B, self.N, self.A, self.c
        it = len(B) + len(N) + 1

        for i in xrange(len(B)):
            if not self.is_int(A[i, 0]):
                aa = -A[i, 1:]
                aa = aa - np.floor(aa)
                bb = A[i, 0]
                bb = - (bb - np.floor(bb))
                A = np.vstack([A, np.hstack([bb, aa])])
                B = np.append(B, it)
                it += 1
        self.A, self.B = A, B

    def dual(self):
        c = - np.hstack([self.c[0], self.A[:, 0]])
        A = - np.vstack([self.c[1:], self.A[:, 1:]]).transpose()
        self.m, self.n = self.n, self.m
        self.N, self.B = self.B, self.N
        self.c, self.A = c, A

    def solve(self):
        B, N, A, c = self.B, self.N, self.A, self.c

        self.inspect('first')
        if self.DEBUG:
            print 'need init:', self.need_init()

        if self.need_init():
            self.initialize()
            res = self.finalize()
            self.inspect('aux final')
            if self.c[0] > self.eps:
                return 'infeasible'
            self.de_initialize()
            self.inspect('origin')

        res = self.finalize()
        self.inspect('first relaxation')

        while not all(map(self.is_int, self.A[:, 0])):
            self.add_cutting_plane()
            self.inspect('cutting')
            self.dual()
            self.inspect('dual')
            res = self.finalize()
            self.inspect('dual final')
            if res == 'unbounded':
                return 'infeasible'
            self.dual()
            self.inspect('primal final')

        if self.DEBUG:
            self.output()
        return self.c[0]

if __name__ == '__main__':
    s = Simplex()
    # s.parse_text(sys.argv[1])
    # s.parse_text('unitTests/ilpTest1')
    s.parse_text('assignmentTests/part1.dict')
    print s.solve()