simplex.py

# python3

"""
Simplex algorithm.

nonbasics = [1, 2, 3]
basics = [4, 5, 6]

matrix = [
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 1, 1, 3, 0, 0, 0],
    [0, 2, 2, 5, 0, 0, 0],
    [0, 4, 1, 2, 0, 0, 0],
]

b = [0, 0, 0, 0, 30, 24, 36]
c = [0, 3, 1, 1, 0, 0, 0]

print(pivot(nonbasics, basics, matrix, b, c, 0, 6, 1))




maximize:     2x1 - x2
subject to:
              2x1 - x2 <= 2
              x1 - 5x2 <= -4
                x1, x2 >= 0
mmatrix=[
  [0, 2, -1],
  [0, 1, -5],
]

bb=[2, -4]
cc=[0, 2, -1]

nonbasics, basics, matrix, b, c, v = initialize(matrix=mmatrix, b=bb, c=cc, m=2)
pprint.pprint(simplex(nonbasics, basics, matrix, b, c, v))


Solve the following linear program using SIMPLEX:
maximize: x1 + 3x2
subject to
          x1 - x2  <= 8
         -x1 - x2  <= -3
         -x1 + 4x2 <= 2
          x1, x2   >= 0

mmatrix = [
    [0,  1, -1],
    [0, -1, -1],
    [0, -1,  4],
]

bb = [8, -3, 2]
cc = [0, 1, 3]

nonbasics, basics, matrix, b, c, v = initialize(matrix=mmatrix, b=bb, c=cc, m=2)

print("initial feasible slack form:")
pprint.pprint((nonbasics, basics, matrix, b, c, v))

print("solution:")

nonbasics = [2, 4]
basics = [3, 5, 1]

mmatrix = [
  [0,  0,    0,    0,  0,   0],
  [0,  0,    1.0,  0, -1.0, 0],
  [0,  0,    0,    0,  0,   0],
  [0,  2.0, -2.0,  0,  1.0, 0],
  [0, -1,   -1,    0,  0,   0],
  [0,  0.0,  5.0,  0, -1.0, 0],
]

bb = [3.0, 3.0, 0,   5.0, -3,   5.0]
cc = [0,   1,   2.0, 0,    1.0, 0]

vv = 3.0

pprint.pprint(simplex(nonbasics, basics, matrix, b, c, v))

##################################################
Solve the following linear program using SIMPLEX:
maximize: x1 + x2
subject to
          x1 + x2  <=  1
         -x1 - x2  <= -2

infeasible solution:
mmatrix = [
    [0,  1,  1],
    [0, -1, -1],
]

bb = [1, -2]
cc = [0, 1, 1]

mmatrix = [
    [0, -1, -1],
    [0,  1,  0],
    [0,  0,  1],
]

bb = [-1, 2, 2]
cc = [0, -1, 2]

nonbasics, basics, matrix, b, c, v = None, None, None, None, None, None

try:
    nonbasics, basics, matrix, b, c, v = initialize(matrix=mmatrix, b=bb, c=cc, m=2)
except InfeasibleException:
    print("No solution")  # infeasible
    exit(1)
except UnboundedException:
    print("Infinity")  # unbounded
    exit(1)

print("initial feasible slack form:")
pprint.pprint((nonbasics, basics, matrix, b, c, v))

print("solution:")
pprint.pprint(simplex(nonbasics, basics, matrix, b, c, v))


standard form:
we're given n real numbers c1, c2, ..., cn;
m real numbers b1, b2, ..., bm for aij for i = 1, 2, ..., m and j = 1, 2, ..., n.

def humanize(nonbasics, basics, matrix, b, c):
    for i in basics:
        print(f"x{i} = ", end='')

        first = True

        for j, k in enumerate(matrix[i]):
            if j in nonbasics:
                if first:
                    print(f"{k}*x{j}", end='')
                else:
                    print(f" + {k}*x{j}", end='')

                first = False

        print(f" + {b[i]}", end='')

        print()
"""

import heapq
from collections import deque


inf = 10 ** 9


class InfeasibleException(Exception):
    pass


class UnboundedException(Exception):
    pass


def initialize(matrix, b, c, m):
    n = len(b)
    k = b.index(min(b))

    if b[k] >= 0:
        for mmm in matrix:
            mmm += [0 for _ in range(n)]

        for i in range(n+m+1):
            matrix.append([0 for _ in range(len(c) + n)])

        items = deque(matrix)
        items.rotate(m+1)
        matrix = list(items)

        c += [0 for _ in range(n)]
        b = [0] * (m+1) + b

        return [i for i in range(1, m+1)], [i for i in range(m+1, n+m+1)], matrix, b, c, 0

    aux_matrix = [[0 for _ in range(n+m+1)] for _ in range(n+m+1)]

    for i, row in enumerate(matrix):
        for j, x in enumerate(row):
            aux_matrix[i+1+m][j] = x

        aux_matrix[i+1+m][0] = -1

    aux_b, aux_c = b.copy(), c.copy()

    aux_b = [0] * (m+1) + aux_b
    aux_c = [0] * (m+n+1)
    aux_c[0] = -1

    aux_nonbasics = [i for i in range(0, m+1)]
    aux_basics = [i for i in range(m+1, n+m+1)]
    aux_v = 0

    # leaving variable l, it must be basic.
    l = m + 1 + k

    aux_v = pivot(aux_nonbasics, aux_basics, aux_matrix, aux_b, aux_c, aux_v, l, 0)

    # use simplex to find an optimal solution to aux
    sol = simplex(aux_nonbasics, aux_basics, aux_matrix, aux_b, aux_c, aux_v)

    # if the optimal solution to aux sets x0 to 0
    if sol[0] == 0:
        # if x0 is basic perform one (degenerate) pivot to make it nonbasic
        # from the final slack form of aux, remove x0 from the constraints.
        if 0 in aux_basics:
            # entering variable can be any non-basic variable that is not 0.
            e = next(e for e in aux_nonbasics if c[e] > 0)
            nonbasic, basics, matrix, b, aux_c, v = pivot(aux_nonbasics, aux_basics, aux_matrix, aux_b, aux_c, aux_v, 0, e)

        # restore the original objective function of L
        aux_c = c
        for _ in range(len(b)):
            aux_c.append(0)

        aux_nonbasics.remove(0)

        # replace each basic variable in this objective function
        # by the right-hand side of its associated constraint.
        # example:
        #   original obj func: x1 + 3x2
        #   constraints: 3 - x2 + x4 = x1
        #   restored obj func => 3 - x2 + x4 + 3x2 = 3 - 2x2 + x4
        for i, j in enumerate(aux_c.copy()):
            if j == 0:
                continue

            if i in aux_basics:
                for k in range(1, len(aux_c)):
                    aux_c[k] = round(aux_c[k] - j * aux_matrix[i][k], 20)

                aux_c[i] = 0
                aux_v += round(j * aux_b[i], 10)

        for m in aux_matrix:
            m[0] = 0

        return aux_nonbasics, aux_basics, aux_matrix, aux_b, aux_c, aux_v
    else:
        raise InfeasibleException


def simplex(nonbasics, basics, matrix, b, c, v):
    while 0 < len(nonbasics):
        e = nonbasics[0]

        for i in nonbasics:
            if c[e] < c[i]:
                e = i

        if c[e] <= 0:
            break

        delta = []
        for i in basics:
            if matrix[i][e] > 0:
                heapq.heappush(delta, (round(b[i] / matrix[i][e], 9), i))
            else:
                heapq.heappush(delta, (inf, i))

        value, l = heapq.heappop(delta)

        if value == inf:
            raise UnboundedException
        else:
            v = pivot(nonbasics, basics, matrix, b, c, v, l, e)

    ans = [0] * (len(c) + 1)

    for i in range(0, len(c) + 1):
        if i in basics:
            ans[i] = b[i]
        else:
            ans[i] = 0

    return ans


def pivot(nonbasics, basics, matrix, b, c, v, l, e):
    b[e] = round(b[l] / matrix[l][e], 20)

    for j in nonbasics:
        if j == e:
            continue

        matrix[e][j] = round(matrix[l][j] / matrix[l][e], 20)

    matrix[e][l] = round(1 / matrix[l][e], 20)

    # compute the coefficients of the remaining constraints.
    for i in basics:
        if i == l:
            continue

        # example:
        # x3 = -3x0 + 1
        # x0 = 2 + x1
        # e = 0; i = 3
        # x3 = 1    -     -3       *  2
        b[i] = b[i] - matrix[i][e] * b[e]

        for j in nonbasics:
            if j == e:
                continue

            matrix[i][j] = round(matrix[i][j] - matrix[i][e] * matrix[e][j], 20)

        matrix[i][l] = round(-matrix[i][e] * matrix[e][l], 20)
        matrix[i][e] = 0

    # compute objective function.
    v = v + c[e] * b[e]

    for j in nonbasics:
        if j == e:
            continue

        c[j] = round(c[j] - c[e] * matrix[e][j], 20)

    c[l] = round(-c[e] * matrix[e][l], 20)
    c[e] = 0

    # compute new sets of basic and non-basic variables.
    nonbasics.remove(e)
    nonbasics.append(l)

    basics.remove(l)
    basics.append(e)

    return v


if __name__ == "__main__":
    nonbasics, basics, v = None, None, None

    def int_inputs():
        return list(map(int, input().split(" ")))

    n, m = int_inputs()

    matrix = [[0] for i in range(n)]

    for i in range(n):
        matrix[i] += int_inputs()

    b = int_inputs()
    c = [0] + int_inputs()

    result = None

    try:
        nonbasics, basics, matrix, b, c, v = initialize(matrix=matrix, b=b, c=c, m=m)
        result = simplex(nonbasics, basics, matrix, b, c, v)
    except InfeasibleException:
        print("No solution")  # infeasible
        exit(1)
    except UnboundedException:
        print("Infinity")  # unbounded
        exit(1)

    print("Bounded solution")
    for i in range(m):
        print(result[i+1], end=" ")