"""Part 8: the simplex method, from scratch, in exact fractions. Solves minimise c.x subject to A_i . x (<=, >= or =) b_i for every row i, x >= 0. Every number is a Fraction, so there is no rounding anywhere and "optimal" means optimal. The method is the textbook one, a dense tableau in two phases: - Phase I adds an artificial variable to every row that has no obvious starting variable and minimises their sum. If that sum can't reach 0 the problem is infeasible; otherwise the basis it ends with is a corner of the feasible region to start from. - Phase II minimises the real objective from that corner, one pivot at a time. Pivots follow Bland's rule: the entering variable is the lowest-numbered one that improves the objective, and ties for the leaving row go to the lowest-numbered basic variable. That rule can never cycle, so the method always stops (Bland, 1977). Alongside the answer it returns the dual prices y, one per row, read off the final tableau. They satisfy b.y = c.x and A^T y <= c, which is a certificate anyone can check that no better x exists. """ from fractions import Fraction from math import gcd LE, GE, EQ = '<=', '>=', '=' class Result: def __init__(self, status, x=None, value=None, y=None, pivots=None, tableau=None): self.status = status # 'optimal', 'infeasible' or 'unbounded' self.tableau = tableau # the final tableau, for gomory_cuts self.x = x # values of the original variables self.value = value # c.x self.y = y # dual prices, one per row self.pivots = pivots # [(phase, entering, leaving, objective after, x after), ...] def solve(c, rows, trace=False): """c: costs, one per variable. rows: [(coefficients {j: a_ij}, op, b_i), ...].""" n = len(c) c = [Fraction(v) for v in c] # Every row with b >= 0, flipping it if need be; remember the flips for the duals. norm = [] for coef, op, b in rows: coef, b = {j: Fraction(a) for j, a in coef.items() if a}, Fraction(b) flip = b < 0 if flip: coef, b = {j: -a for j, a in coef.items()}, -b op = {LE: GE, GE: LE, EQ: EQ}[op] norm.append((coef, op, b, flip)) m = len(norm) # Columns: the n originals, one slack (<=) or surplus (>=) per inequality, one artificial for # every >= and = row. unit[i] is the column that starts as row i's basic variable; its column # in the final tableau gives row i's dual price. cols = n extra, basis, unit, artificial = [], [], [], set() for i, (coef, op, b, flip) in enumerate(norm): if op == LE: extra.append((cols, i, 1)) basis.append(cols) unit.append(cols) cols += 1 else: if op == GE: extra.append((cols, i, -1)) cols += 1 extra.append((cols, i, 1)) artificial.add(cols) basis.append(cols) unit.append(cols) cols += 1 T = [[Fraction(0)] * (cols + 1) for _ in range(m)] for i, (coef, op, b, flip) in enumerate(norm): for j, a in coef.items(): T[i][j] = a T[i][cols] = b for j, i, s in extra: T[i][j] = Fraction(s) pivots = [] if trace else None def run(cost, phase, allowed): """Minimise cost.x from the current basis; returns False if unbounded.""" while True: # reduced costs r_j = c_j - c_B . column j r = [cost[j] - sum(cost[basis[i]] * T[i][j] for i in range(m)) for j in range(cols)] enter = next((j for j in range(cols) if allowed(j) and j not in basis and r[j] < 0), None) if enter is None: return True best, leave = None, None for i in range(m): if T[i][enter] > 0: ratio = T[i][cols] / T[i][enter] if best is None or ratio < best or (ratio == best and basis[i] < basis[leave]): best, leave = ratio, i if leave is None: return False out = basis[leave] pivot(leave, enter) if trace: value = sum(cost[basis[i]] * T[i][cols] for i in range(m)) now = [Fraction(0)] * n for i in range(m): if basis[i] < n: now[basis[i]] = T[i][cols] pivots.append((phase, enter, out, value, now)) def pivot(i, j): p = T[i][j] T[i] = [v / p for v in T[i]] for k in range(m): if k != i and T[k][j]: f = T[k][j] T[k] = [a - f * b for a, b in zip(T[k], T[i])] basis[i] = j # Phase I: minimise the sum of the artificial variables. if artificial: cost1 = [Fraction(1) if j in artificial else Fraction(0) for j in range(cols)] run(cost1, 1, lambda j: True) if sum(T[i][cols] for i in range(m) if basis[i] in artificial) > 0: return Result('infeasible', pivots=pivots) # Drive any artificial still basic (at 0) out of the basis. If its row has nothing else in # it, the row was redundant: it stays, and no later pivot can change it. for i in range(m): if basis[i] in artificial: j = next((j for j in range(cols) if j not in artificial and T[i][j] != 0), None) if j is not None: pivot(i, j) # Phase II: the real objective; artificial columns may not re-enter. cost2 = c + [Fraction(0)] * (cols - n) if not run(cost2, 2, lambda j: j not in artificial): return Result('unbounded', pivots=pivots) x = [Fraction(0)] * cols for i in range(m): x[basis[i]] = T[i][cols] # Dual prices: y_i = c_B B^-1 e_i, and column unit[i] of the tableau is B^-1 e_i. y = [] for i, (coef, op, b, flip) in enumerate(norm): yi = sum(cost2[basis[k]] * T[k][unit[i]] for k in range(m)) y.append(-yi if flip else yi) value = sum(ci * xi for ci, xi in zip(c, x[:n])) # What each column stands for, in the original variables: (coefficients, constant). A slack is # b - a.x, a surplus a.x - b; artificial columns are 0 at the optimum and stand for nothing. stands = [({j: Fraction(1)}, Fraction(0)) for j in range(n)] + [None] * (cols - n) for j, i, s in extra: if j not in artificial: coef, op, b, flip = norm[i] stands[j] = ({k: -a * s for k, a in coef.items()}, b * s) tableau = {'T': T, 'basis': basis, 'cols': cols, 'stands': stands} return Result('optimal', x[:n], value, y, pivots, tableau) def gomory_cuts(result): """Gomory's fractional cuts from an optimal tableau (Gomory, 1958). Valid when every variable must be a whole number and every row has whole-number coefficients, so that every slack is a whole number too. Take a row whose basic variable has a fractional value: x_B + sum_k a_k y_k = b over the nonbasic columns y_k >= 0. Every whole-number answer satisfies sum_k frac(a_k) y_k >= frac(b), and the current corner, where the y_k are 0, doesn't. Returns each such cut as a row in the original variables, scaled to whole numbers so that its own slack is a whole number as well and later cuts stay valid.""" T, basis, cols, stands = (result.tableau[k] for k in ('T', 'basis', 'cols', 'stands')) rows = [] for i, row in enumerate(T): f0 = row[cols] - (row[cols].numerator // row[cols].denominator) if f0 == 0 or stands[basis[i]] is None: continue coef, const = {}, Fraction(0) for k in range(cols): if k in basis or stands[k] is None: continue f = row[k] - (row[k].numerator // row[k].denominator) if f: for j, a in stands[k][0].items(): coef[j] = coef.get(j, 0) + f * a const += f * stands[k][1] coef = {j: a for j, a in coef.items() if a} rhs = f0 - const scale = 1 for q in list(coef.values()) + [rhs]: scale = scale * q.denominator // gcd(scale, q.denominator) rows.append(({j: int(a * scale) for j, a in coef.items()}, GE, int(rhs * scale))) return rows def check_certificate(c, rows, x, y): """True when x is feasible, y is dual feasible with the right signs, and c.x = b.y: then x is optimal, whatever program produced it. Works in exact fractions.""" n = len(c) if any(v < 0 for v in x): return False for (coef, op, b), yi in zip(rows, y): ax = sum(Fraction(a) * x[j] for j, a in coef.items()) if (op == LE and ax > b) or (op == GE and ax < b) or (op == EQ and ax != b): return False if (op == GE and yi < 0) or (op == LE and yi > 0): return False for j in range(n): if sum(Fraction(coef.get(j, 0)) * yi for (coef, op, b), yi in zip(rows, y)) > c[j]: return False return sum(Fraction(ci) * xi for ci, xi in zip(c, x)) == sum(Fraction(b) * yi for (_, _, b), yi in zip(rows, y))