Linear & Integer Programming

"""Dense-tableau simplex, transcribed literally.

Solves   minimize c.x   subject to   A x <= b,  x >= 0,  with b >= 0.

b >= 0 is what lets us skip phase one: the slack variables form a ready-made
feasible starting basis (x = 0). Each block below is one step of the algorithm.
"""

TOL = 1e-9
MAX_PIVOTS = 10_000


def simplex(c, A, b):
    m = len(A)            # number of constraints
    n = len(c)            # number of decision variables
    width = n + m + 1     # decision columns, slack columns, right-hand side

    # Tableau rows are [ A | I | b ]; the last row is the cost row [ c | 0 | 0 ].
    tableau = []
    for i in range(m):
        row = list(A[i])
        row += [1.0 if k == i else 0.0 for k in range(m)]
        row.append(b[i])
        tableau.append(row)
    tableau.append(list(c) + [0.0] * m + [0.0])

    basis = [n + i for i in range(m)]     # slack i starts out basic

    for _ in range(MAX_PIVOTS):
        # 1. Pricing: entering column = the most negative reduced cost.
        enter = -1
        most_negative = -TOL
        for j in range(n + m):
            if tableau[m][j] < most_negative:
                most_negative = tableau[m][j]
                enter = j
        if enter == -1:
            break                         # no improving edge: optimal

        # 2. Ratio test: how far can we move before a basic variable hits 0?
        leave = -1
        best_ratio = float("inf")
        for i in range(m):
            if tableau[i][enter] > TOL:
                ratio = tableau[i][width - 1] / tableau[i][enter]
                if ratio < best_ratio:
                    best_ratio = ratio
                    leave = i
        if leave == -1:
            raise ValueError("unbounded: the objective falls without limit")

        # 3. Pivot: Gauss-Jordan elimination on (leave, enter).
        pivot = tableau[leave][enter]
        for j in range(width):
            tableau[leave][j] /= pivot
        for i in range(m + 1):
            if i == leave:
                continue
            factor = tableau[i][enter]
            for j in range(width):
                tableau[i][j] -= factor * tableau[leave][j]
        basis[leave] = enter
    else:
        raise RuntimeError("pivot limit reached; the problem may be cycling")

    x = [0.0] * n
    for i in range(m):
        if basis[i] < n:
            x[basis[i]] = tableau[i][width - 1]
    objective = -tableau[m][width - 1]
    # Reduced cost of slack i is -y_i, and y_i = d(objective)/d(b_i).
    shadow_prices = [-tableau[m][n + i] for i in range(m)]
    return x, objective, shadow_prices


if __name__ == "__main__":
    # minimize -3 x1 - 5 x2  (i.e. maximize 3 x1 + 5 x2)
    #   x1 <= 4,   2 x2 <= 12,   3 x1 + 2 x2 <= 18
    x, z, y = simplex([-3.0, -5.0], [[1, 0], [0, 2], [3, 2]], [4.0, 12.0, 18.0])
    print(x, z, y)    # [2.0, 6.0] -36.0 [0.0, -1.5, -1.0]

"""Linear programming - a typed, validated two-phase simplex API.

Infeasible and unbounded problems are legitimate answers, not crashes, so each
has its own exception type the caller can catch without guessing from a message.
"""

from dataclasses import dataclass

import numpy as np
from numpy.typing import NDArray

Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]

TOL = 1e-9
MAX_PIVOTS = 50_000


class LinearProgramError(Exception):
    """Base class for every outcome the solver reports on purpose."""


class InfeasibleError(LinearProgramError):
    """No x satisfies every constraint."""


class UnboundedError(LinearProgramError):
    """The objective falls without limit along a feasible ray."""


class PivotLimitError(LinearProgramError):
    """Cycling, or an instance harder than MAX_PIVOTS allows."""


@dataclass(frozen=True)
class LinearProgram:
    """minimize c.x  subject to  A x <= b,  x >= 0.  Validated on construction."""

    c: Vector
    A: Matrix
    b: Vector

    def __post_init__(self) -> None:
        if self.A.ndim != 2 or self.c.ndim != 1 or self.b.ndim != 1:
            raise ValueError("expected A (m, n), c (n,), b (m,)")
        rows, cols = self.A.shape
        if rows == 0 or cols == 0:
            raise ValueError("cannot solve an empty program")
        if self.c.shape != (cols,) or self.b.shape != (rows,):
            raise ValueError(
                f"shape mismatch: A={self.A.shape}, c={self.c.shape}, b={self.b.shape}"
            )
        for name, values in (("c", self.c), ("A", self.A), ("b", self.b)):
            if not np.isfinite(values).all():
                raise ValueError(f"{name} contains NaN or infinity")


@dataclass(frozen=True)
class Solution:
    x: Vector
    objective: float
    shadow_prices: Vector      # d(objective)/d(b_i), valid while the basis holds


def _pivot(tableau: Matrix, row: int, col: int) -> None:
    tableau[row] /= tableau[row, col]
    factors = tableau[:, col].copy()
    factors[row] = 0.0
    tableau -= np.outer(factors, tableau[row])      # one rank-1 update


def _run_simplex(tableau: Matrix, basis: NDArray[np.int_], eligible: int) -> None:
    """Pivot until optimal. Only columns below 'eligible' may enter the basis."""
    rows = tableau.shape[0] - 1
    for _ in range(MAX_PIVOTS):
        reduced = tableau[rows, :eligible]
        enter = int(np.argmin(reduced))
        if reduced[enter] >= -TOL:
            return
        column = tableau[:rows, enter]
        positive = column > TOL
        if not positive.any():
            raise UnboundedError("objective decreases without limit")
        ratios = np.full(rows, np.inf)
        ratios[positive] = tableau[:rows, -1][positive] / column[positive]
        leave = int(np.argmin(ratios))
        _pivot(tableau, leave, enter)
        basis[leave] = enter
    raise PivotLimitError(f"no optimum within {MAX_PIVOTS} pivots")


def solve(lp: LinearProgram) -> Solution:
    """Solve exactly. Raises InfeasibleError or UnboundedError, never returns garbage."""
    m, n = lp.A.shape

    # A row with b < 0 is negated so its right-hand side is >= 0; its slack then
    # enters with coefficient -1 and cannot start basic, so it gets an artificial.
    sign = np.where(lp.b < 0.0, -1.0, 1.0)
    art_rows = np.flatnonzero(sign < 0.0)
    k = art_rows.size
    art_cols = n + m + np.arange(k)

    tableau = np.zeros((m + 1, n + m + k + 1))
    tableau[:m, :n] = sign[:, None] * lp.A
    tableau[:m, n : n + m] = np.diag(sign)
    tableau[:m, -1] = sign * lp.b
    tableau[art_rows, art_cols] = 1.0
    basis = np.arange(n, n + m)
    basis[art_rows] = art_cols

    if k > 0:
        # Phase I: minimize the sum of the artificials. Zero means feasible.
        tableau[m, art_cols] = 1.0
        for row in art_rows:
            tableau[m] -= tableau[row]
        _run_simplex(tableau, basis, eligible=n + m + k)
        if -tableau[m, -1] > TOL:
            raise InfeasibleError("no point satisfies every constraint")
        for row in np.flatnonzero(basis >= n + m):      # drive artificials out
            real = np.flatnonzero(np.abs(tableau[row, : n + m]) > TOL)
            if real.size > 0:
                _pivot(tableau, int(row), int(real[0]))
                basis[row] = real[0]

    # Phase II: the true cost, priced out against the current basis.
    tableau[m] = 0.0
    tableau[m, :n] = lp.c
    for row, col in enumerate(basis):
        if col < n:
            tableau[m] -= lp.c[col] * tableau[row]
    _run_simplex(tableau, basis, eligible=n + m)

    x = np.zeros(n)
    for row, col in enumerate(basis):
        if col < n:
            x[col] = tableau[row, -1]
    return Solution(x=x, objective=float(-tableau[m, -1]), shadow_prices=-tableau[m, n : n + m])

"""Revised simplex with factorisation reuse - then the production MILP path."""

import numpy as np
import scipy.linalg as la
from numpy.typing import NDArray
from scipy import sparse
from scipy.optimize import Bounds, LinearConstraint, milp

Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]


def revised_simplex(
    A: Matrix, b: Vector, c: Vector, basis: NDArray[np.int_],
    tol: float = 1e-9, max_pivots: int = 50_000,
) -> tuple[Vector, float, Vector]:
    """minimize c.x  s.t.  A x = b, x >= 0, from a feasible starting basis.

    No tableau and no explicit inverse. Each iteration factors the m x m basis
    ONCE and reuses that factorisation for all three solves it needs: the primal
    values, the duals, and the entering column.
    """
    m, total = A.shape
    basis = basis.copy()
    ratios = np.empty(m)                                   # reused every iteration

    for _ in range(max_pivots):
        factor = la.lu_factor(A[:, basis])                 # one O(m^3) factorisation...
        x_basic = la.lu_solve(factor, b)                   # ...B x_B = b
        duals = la.lu_solve(factor, c[basis], trans=1)     # ...B^T y = c_B
        reduced = c - A.T @ duals                          # price EVERY column in one matvec
        reduced[basis] = 0.0

        enter = int(np.argmin(reduced))
        if reduced[enter] >= -tol:
            x = np.zeros(total)
            x[basis] = x_basic
            return x, float(c @ x), duals

        direction = la.lu_solve(factor, A[:, enter])       # ...B d = a_enter
        positive = direction > tol
        if not positive.any():
            raise ValueError("unbounded")
        ratios.fill(np.inf)
        np.divide(x_basic, direction, out=ratios, where=positive)
        basis[int(np.argmin(ratios))] = enter
    raise RuntimeError("pivot limit reached")


def solve_milp(
    c: Vector, A_ub: Matrix, b_ub: Vector, is_integer: NDArray[np.bool_],
    time_limit_s: float, mip_rel_gap: float,
) -> Vector:
    """Production path once any variable must be integer: HiGHS via SciPy.

    Presolve, cutting planes and primal heuristics are the part a hand-rolled
    simplex does not have, and they dominate MILP solve time.
    """
    result = milp(
        c=c,
        constraints=LinearConstraint(sparse.csr_array(A_ub), ub=b_ub),
        integrality=is_integer.astype(int),
        bounds=Bounds(lb=0.0),
        options={"time_limit": time_limit_s, "mip_rel_gap": mip_rel_gap},
    )
    if result.status == 2:
        raise ValueError("infeasible")
    if result.status == 3:
        raise ValueError("unbounded")
    if result.x is None:
        raise TimeoutError(f"no incumbent within {time_limit_s}s: {result.message}")
    return result.x