From 30b098f6bca59753588d904f8fa889b57b93fea8 Mon Sep 17 00:00:00 2001 From: Tomas Roun Date: Sat, 17 Feb 2024 22:15:09 +0100 Subject: [PATCH] Refactor code --- .gitignore | 1 + README.md | 12 +- pivotal/__init__.py | 3 +- pivotal/api.py | 391 ++------------------------------- pivotal/errors.py | 10 - pivotal/expressions.py | 255 +++++++++++++++++++++ pivotal/simplex.py | 488 +++++++++++++++++++++++++++-------------- pyproject.toml | 40 ++-- tests/test_api.py | 122 +++++------ tests/test_simplex.py | 194 ++++++++-------- 10 files changed, 798 insertions(+), 718 deletions(-) create mode 100644 pivotal/expressions.py diff --git a/.gitignore b/.gitignore index 41fd999..25fe812 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ __pycache__ *.egg-info +dist .pytest_cache/ diff --git a/README.md b/README.md index a21d5ef..034cf91 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,7 @@ sum(X) ``` Note that variables are considered equal if they have the same name, so -for examples this expression: +for example this expression: ```python Variable("x") + 2 + Variable("x") @@ -72,7 +72,7 @@ Variable("x") + 2 + Variable("x") will be treated as simply `2*x+2`. -The first argument to `minimize` and `maximize` is the objective function which must be either a single variable or a linear combination as in the example above. +The first argument to `minimize` and `maximize` is the objective function which must be either a single variable or a linear combination as in the examples above. ### Constraints @@ -97,7 +97,7 @@ There is no need to convert your constraints to the canonical form which uses on The return value of `minimize` and `maximize` is a 2-tuple containing the value of the objective function and a dictionary of variables and their values. -The functions may raise `pivotal.Infeasible` if the program is over-constrained (no solution exists) or `pivotal.Unbounded` if the program is under-constrained (the objective can be made arbitrarily small). +The functions may raise `pivotal.Infeasible` if the program is over-constrained (no solution exists) or `pivotal.Unbounded` if the program is under-constrained (the objective can be made arbitrarily small): ```python from pivotal import minimize, maximize, Variable, Infeasible @@ -123,13 +123,15 @@ except Infeasible: ## Limitations -- Currently, all variables are assumed to be positive +- Currently, all variables are assumed to be nonnegative i.e. x >= 0. ## TODO (Contributions welcome) - ✔️ Setting tolerance & max number of iterations - (WIP) Arbitrary variable bounds, e.g. `a <= x <= b` -- (WIP) Support for absolute values +- (WIP) Support for absolute values using Python's `abs()` + - in the objective function + - in constraints - MILP solver with branch & bound diff --git a/pivotal/__init__.py b/pivotal/__init__.py index e8a5621..87a4082 100644 --- a/pivotal/__init__.py +++ b/pivotal/__init__.py @@ -1,5 +1,6 @@ -from pivotal.api import Variable, maximize, minimize +from pivotal.api import maximize, minimize from pivotal.errors import Infeasible, Unbounded +from pivotal.expressions import Variable __version__ = "0.0.3" diff --git a/pivotal/api.py b/pivotal/api.py index dd4ee36..abc06bf 100644 --- a/pivotal/api.py +++ b/pivotal/api.py @@ -1,382 +1,31 @@ -import itertools import math +from typing import TypeAlias -import numpy as np +from pivotal.expressions import Constraint, Expression, get_variable_names +from pivotal.simplex import ProgramType, solve -from pivotal.simplex import LinearProgram as LP -from pivotal.simplex import canonicalize, solve +ResultType: TypeAlias = tuple[float, dict[str, float]] -class Expression: - pass +def optimize( + _type: ProgramType, objective: Expression, constraints: list[Constraint], max_iterations: int, tolerance: float +) -> ResultType: + value, variables = solve(_type, objective, constraints, max_iterations=max_iterations, tolerance=tolerance) + value = value if _type == ProgramType.MIN else -value -class Abs(Expression): - def __init__(self, arg, sign=1): - self.arg = arg - self.sign = sign + user_defined_vars = get_variable_names((objective, *constraints)) + variables = variables[: len(user_defined_vars)] + return value, ({user_defined_vars[i]: v for i, v in enumerate(variables)}) - def __abs__(self): - return Abs(self.arg) - def __neg__(self): - return Abs(self.arg, -self.sign) +def minimize( + objective: Expression, constraints: list[Constraint], *, max_iterations: int = math.inf, tolerance: float = 1e-6 +) -> ResultType: + return optimize(ProgramType.MIN, objective, constraints, max_iterations, tolerance) - def __add__(self, other): - match other: - case Abs(): - return Sum(self, other) - case Sum(elts=elts): - return Sum(self, *elts) - case int() | float(): - if other == 0: - return self - return Sum(self, other) - case Variable() | int() | float(): - return Sum(self, other) - case _: - return NotImplemented - def __sub__(self, other): - return self.__add__(-other) - - def __radd__(self, other): - match other: - case int() | float(): - if other == 0: - return self - return Sum(other, self) - case _: - return NotImplemented - - def __rsub__(self, other): - return -self.__radd__(other) - - def __mul__(self, other): - match other: - case int() | float(): - if other == 1: - return self - elif other == 0: - return 0 - elif other > 0: - return Abs(other*self.arg) - else: - return Abs(abs(other)*self.arg, sign=-1) - case _: - return NotImplemented - - __rmul__ = __mul__ - - def __str__(self): - sign = "-" if self.sign == -1 else "" - return f"{sign}|{self.arg}|" - - -class Comparable: - def __eq__(self, other): - return Eq(self, other) - - def __ge__(self, other): - return Gte(self, other) - - def __le__(self, other): - return Lte(self, other) - - -class Variable(Expression, Comparable): - id_iter = itertools.count() - - def __init__(self, name: str | None = None, coeff=1.0, lb=0, ub=0): - if name is None: - id = next(self.id_iter) - self.name = f"_{id}" - # TODO: underscored variable names - else: - self.name = name - self.coeff = coeff - - def __abs__(self): - return Abs(self) - - def __neg__(self): - return Variable(name=self.name, coeff=-self.coeff) - - def __add__(self, other): - match other: - case Variable(name=name, coeff=coeff) if name == self.name: - if self.coeff+coeff == 0: - return 0 - return Variable(name=name, coeff=self.coeff+coeff) - case int() | float() if other == 0: - return self - case Variable() | Abs() | int() | float(): - return Sum(self, other) - case Sum(elts=elts): - return Sum(self, *elts) - case _: - return NotImplemented - - def __radd__(self, other): - match other: - case int() | float(): - if other == 0: - return self - return Sum(other, self) - case _: - return NotImplemented - - def __sub__(self, other): - return self.__add__(-other) - - def __rsub__(self, other): - return -self.__radd__(other) - - def __mul__(self, other): - match other: - case int() | float(): - if other == 0: - return 0 - elif other == 1: - return self - return Variable(self.name, coeff=self.coeff * other) - case _: - return NotImplemented - - __rmul__ = __mul__ - - def __repr__(self): - if self.coeff == 1.0: - return self.name - elif self.coeff == -1.0: - return f"-{self.name}" - else: - return f"{self.coeff}*{self.name}" - - -def simplify_sum(elts): - c = 0 - variables = {} - for expr in elts: - match expr: - case int() | float(): - c += expr - case Variable(name=name): - if name in variables: - variables[name] += expr - else: - variables[name] = expr - case _: - raise Exception("cannot happen") - return variables, c - - -def simplify(expr): - match expr: - case int() | float(): - return {}, expr - case Variable(): - return {expr.name: expr}, 0 - case Sum(elts=elts): - return simplify_sum(elts) - case Constraint(): - return expr.normalize() - - -class Sum(Expression, Comparable): - def __init__(self, *elts): - self.elts = elts - - def __abs__(self): - return Abs(self) - - def __neg__(self): - return Sum(*(-expr for expr in self.elts)) - - def __add__(self, other): - match other: - case Sum(elts=elts): - return Sum(*self.elts, *elts) - case int() | float() if other == 0: - return self - case Variable() | Abs() | int() | float(): - return Sum(*self.elts, other) - case _: - return NotImplemented - - def __sub__(self, other): - return self.__add__(-other) - - def __radd__(self, other): - match other: - case int() | float() if other == 0: - if other == 0: - return self - return Sum(other, *self.elts) - case _: - return NotImplemented - - def __rsub__(self, other): - return -self.__radd__(other) - - def __mul__(self, other): - match other: - case int() | float(): - if other == 0: - return 0 - elif other == 1: - return self - return Sum(*(other*elt for elt in self.elts)) - case _: - return NotImplemented - - __rmul__ = __mul__ - - def __repr__(self): - # out = repr(self.elts[0]) - # for expr in self.elts[1:]: - # if - elts = " + ".join(repr(x) for x in self.elts) - return f"{elts}" - - -class Constraint: - def __init__(self, left, right): - self.left = left - self.right = right - - def normalize(self): - vars_left, c_left = simplify(self.left) - vars_right, c_right = simplify(self.right) - - c = c_right - c_left - variables = vars_left - for name in vars_right: - if name in variables: - variables[name] -= vars_right[name] - else: - variables[name] = -vars_right[name] - - return variables, c - - -class Eq(Constraint): - def __str__(self): - return f"{self.left} = {self.right}" - - -class Gte(Constraint): - def __str__(self): - return f"{self.left} >= {self.right}" - - -class Lte(Constraint): - def __str__(self): - return f"{self.left} <= {self.right}" - - -_type = { - Eq: 0, - Gte: 1, - Lte: -1, -} - - -class LinearProgram: - MIN = 'min' - MAX = 'max' - - def __init__(self, type, objective, constraints): - self.type = type - self.objective = objective - self.constraints = constraints - - def validate(self): - ... - # TODO: check for nesting - - # for elt in self.objective: - # if not isinstance(elt, Abs): - # continue - # if self.type == LinearProgram.MIN and elt.sign == -1: - # raise Exception("Cannot minimize a negative absolute value") - # if self.type == LinearProgram.MAX and elt.sign == 1: - # raise Exception("Cannot maximize a positive absolute value") - - # for constraint in self.constraints: - # for elt in self.objective: - # if not isinstance(elt, Abs): - # continue - # if self.type == LinearProgram.MIN and elt.sign == -1: - # raise Exception("Cannot minimize a negative absolute value") - # if self.type == LinearProgram.MAX and elt.sign == 1: - # raise Exception("Cannot maximize a positive absolute value") - - def as_matrix(self): - variables = sorted(list(self._get_variables())) - n_vars = len(variables) - n_constraints = len(self.constraints) - A = np.zeros((n_constraints, n_vars), dtype=float) - b = np.zeros(n_constraints, dtype=float) - c = np.zeros(n_vars, dtype=float) - constraint_types = np.zeros(n_constraints, dtype=int) - - _vars, _ = simplify(self.objective) - for name in _vars: - i = variables.index(name) - c[i] = _vars[name].coeff - - for row, constraint in enumerate(self.constraints): - _vars, _b = simplify(constraint) - for name in _vars: - column = variables.index(name) - A[row, column] = _vars[name].coeff - b[row] = _b - constraint_types[row] = _type[type(constraint)] - - return self.type, A, b, c, constraint_types - - def _get_variables(self): - variables = set() - _vars, _ = simplify(self.objective) - variables |= set(_vars.keys()) - for constraint in self.constraints: - _vars, _ = simplify(constraint) - variables |= set(_vars.keys()) - return variables - - @classmethod - def minimize(cls, objective): - return cls(LinearProgram.MIN, objective, []) - - @classmethod - def maximize(cls, objective): - return cls(LinearProgram.MAX, objective, []) - - def such_that(self, *constraints): - self.constraints = constraints - return self - - def optimize(self, max_iterations, tolerance): - all_vars = sorted(list(self._get_variables())) - A, b, c = canonicalize(*self.as_matrix()) - value, variables = solve(LP(A, b, c, {}), max_iterations=max_iterations, tolerance=tolerance) - value = value if self.type == LinearProgram.MIN else -value - variables = variables[:len(all_vars)] - return value, ({all_vars[i]: v for i, v in enumerate(variables)}) - - def __str__(self): - constraints = "\n".join((f" {constraint}" for constraint in self.constraints)) - return f"{self.type} {self.objective}\ns.t.\n{constraints}" - - -def minimize(objective: Expression, constraints: list[Constraint], - *, max_iterations=math.inf, tolerance=1e-6) -> tuple[float, dict[str, float]]: - return LinearProgram.minimize(objective).such_that(*constraints).optimize(max_iterations=max_iterations, - tolerance=tolerance) - - -def maximize(objective: Expression, constraints: list[Constraint], - *, max_iterations=math.inf, tolerance=1e-6) -> tuple[float, dict[str, float]]: - return LinearProgram.maximize(objective).such_that(*constraints).optimize(max_iterations=max_iterations, - tolerance=tolerance) +def maximize( + objective: Expression, constraints: list[Constraint], *, max_iterations: int = math.inf, tolerance: float = 1e-6 +) -> ResultType: + return optimize(ProgramType.MAX, objective, constraints, max_iterations, tolerance) diff --git a/pivotal/errors.py b/pivotal/errors.py index 620a498..170a87f 100644 --- a/pivotal/errors.py +++ b/pivotal/errors.py @@ -1,16 +1,6 @@ - - class Unbounded(Exception): pass class Infeasible(Exception): pass - - -class InvalidConstraint(Exception): - pass - - -class InvalidObjective(Exception): - pass diff --git a/pivotal/expressions.py b/pivotal/expressions.py new file mode 100644 index 0000000..5b89a06 --- /dev/null +++ b/pivotal/expressions.py @@ -0,0 +1,255 @@ +import itertools +from collections import defaultdict +from typing import Literal, Self + + +class Expression: + """Base class for all expressions.""" + + +class Constraint: + def __init__(self, left: Expression, right: Expression) -> None: + self.left = left + self.right = right + + +class Equal(Constraint): + def __repr__(self) -> str: + return f"{self.left} = {self.right}" + + +class GreaterOrEqual(Constraint): + def __repr__(self) -> str: + return f"{self.left} >= {self.right}" + + +class LessOrEqual(Constraint): + def __repr__(self) -> str: + return f"{self.left} <= {self.right}" + + +class ComparableMixin: + def __eq__(self, other: Expression) -> Equal: + return Equal(self, other) + + def __ge__(self, other: Expression) -> GreaterOrEqual: + return GreaterOrEqual(self, other) + + def __le__(self, other: Expression) -> LessOrEqual: + return LessOrEqual(self, other) + + +class Variable(Expression, ComparableMixin): + _id_iter = itertools.count() + + def __init__(self, name: str | None = None, coeff: float = 1.0) -> None: + if name is None: + _id = next(self._id_iter) + self.name = f"_{_id}" + else: + self.name = name + self.coeff = coeff + + def __abs__(self) -> "Abs": + return Abs(self) + + def __neg__(self) -> Self: + return Variable(name=self.name, coeff=-self.coeff) + + def __add__(self, other: "float | Self | Sum | Abs") -> "float | Self | Sum": + match other: + case Variable(name=name, coeff=coeff) if name == self.name: + if self.coeff + coeff == 0: + return 0 + return Variable(name=name, coeff=self.coeff + coeff) + case int() | float() if other == 0: + return self + case Variable() | Abs() | int() | float(): + return Sum(self, other) + case Sum(elts=elts): + return Sum(self, *elts) + case _: + return NotImplemented + + def __radd__(self, other: float) -> "Self | Sum": + match other: + case int() | float(): + if other == 0: + return self + return Sum(other, self) + case _: + return NotImplemented + + def __sub__(self, other: "float | Self | Sum | Abs") -> "float | Self | Sum": + return self.__add__(-other) + + def __rsub__(self, other: float) -> "Self | Sum": + return (-self).__radd__(other) + + def __mul__(self, other: float) -> "float | Self": + match other: + case int() | float(): + if other == 0: + return 0 + if other == 1: + return self + return Variable(self.name, coeff=self.coeff * other) + case _: + return NotImplemented + + __rmul__ = __mul__ + + def __repr__(self) -> str: + if self.coeff == 1.0: + return self.name + if self.coeff == -1.0: + return f"-{self.name}" + return f"{self.coeff}*{self.name}" + + +class Sum(Expression, ComparableMixin): + def __init__(self, *elts: list[Expression]) -> None: + self.elts = elts + + def __abs__(self) -> "Abs": + return Abs(self) + + def __neg__(self) -> Self: + return Sum(*(-expr for expr in self.elts)) + + def __add__(self, other: "float | Variable | Self | Abs") -> Self: + match other: + case Sum(elts=elts): + return Sum(*self.elts, *elts) + case int() | float() if other == 0: + return self + case Variable() | Abs() | int() | float(): + return Sum(*self.elts, other) + case _: + return NotImplemented + + def __sub__(self, other: "float | Variable | Self | Abs") -> Self: + return self.__add__(-other) + + def __radd__(self, other: float) -> Self: + match other: + case int() | float() if other == 0: + if other == 0: + return self + return Sum(other, *self.elts) + case _: + return NotImplemented + + def __rsub__(self, other: float) -> Self: + return (-self).__radd__(other) + + def __mul__(self, other: float) -> float | Self: + match other: + case int() | float(): + if other == 0: + return 0 + if other == 1: + return self + return Sum(*(other * elt for elt in self.elts)) + case _: + return NotImplemented + + __rmul__ = __mul__ + + def __repr__(self) -> str: + return " + ".join(repr(x) for x in self.elts) + + +class Abs(Expression): + def __init__(self, arg: Expression, sign: Literal[-1, 1] = 1) -> None: + self.arg = arg + self.sign = sign + + def __abs__(self) -> Self: + return Abs(self.arg) + + def __neg__(self) -> Self: + return Abs(self.arg, -self.sign) + + def __add__(self, other: float | Variable | Self | Sum) -> Self | Sum: + match other: + case Abs(): + return Sum(self, other) + case Sum(elts=elts): + return Sum(self, *elts) + case int() | float(): + if other == 0: + return self + return Sum(self, other) + case Variable() | int() | float(): + return Sum(self, other) + case _: + return NotImplemented + + def __sub__(self, other: float | Variable | Self | Sum) -> Self | Sum: + return self.__add__(-other) + + def __radd__(self, other: float) -> Self | Sum: + match other: + case int() | float(): + if other == 0: + return self + return Sum(other, self) + case _: + return NotImplemented + + def __rsub__(self, other: float) -> Self | Sum: + return (-self).__radd__(other) + + def __mul__(self, other: float) -> float | Self: + match other: + case int() | float(): + if other == 1: + return self + if other == 0: + return 0 + if other > 0: + return Abs(other * self.arg) + return Abs(abs(other) * self.arg, sign=-1) + case _: + return NotImplemented + + __rmul__ = __mul__ + + def __repr__(self) -> str: + sign = "-" if self.sign == -1 else "" + return f"{sign}|{self.arg}|" + + +def get_variable_coeffs(elem: Expression | Constraint) -> tuple[dict[str, float], float]: + c = 0 + variables = defaultdict(float) + match elem: + case int() | float() as const: + c += const + case Variable(name=name, coeff=coeff): + variables[name] += coeff + case Sum(elts=elts): + for expr in elts: + _variables, _c = get_variable_coeffs(expr) + c += _c + for v in _variables: + variables[v] += _variables[v] + case Constraint(left=left, right=right): + vars_left, c_left = get_variable_coeffs(left) + vars_right, c_right = get_variable_coeffs(right) + + c += c_right - c_left + for v in vars_left: + variables[v] += vars_left[v] + for v in vars_right: + variables[v] -= vars_right[v] + return variables, c + + +def get_variable_names(elems: list[Expression | Constraint]) -> list[str]: + variables = set() + for elem in elems: + coeffs, _ = get_variable_coeffs(elem) + variables |= set(coeffs.keys()) + return sorted(variables) diff --git a/pivotal/simplex.py b/pivotal/simplex.py index 3e2633f..56ce13b 100644 --- a/pivotal/simplex.py +++ b/pivotal/simplex.py @@ -1,32 +1,62 @@ import math +import warnings +from enum import Enum, auto +from typing import Literal import numpy as np from pivotal.errors import Infeasible, Unbounded +from pivotal.expressions import (Constraint, Equal, Expression, GreaterOrEqual, LessOrEqual, get_variable_coeffs, + get_variable_names) + + +class ProgramType(Enum): + MIN = auto() + MAX = auto() + + +def suppress_divide_by_zero_warning(fn) -> None: + def _fn_suppressed(*args, **kwargs): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + return fn(*args, **kwargs) + + return _fn_suppressed class Pivots: - def __init__(self, cr): + """ + Helper class to manage the positions of pivots in the tableau. + + The pivots can be accessed by row or column index using either `rc` or `cr` respectively. + The row and columns indices are based on the full M matrix. If you want to get the position + in the A matrix, you need to subtract 1 from the row index. + """ + + def __init__(self, cr: dict[int, int]) -> None: self.cr = cr self.rc = {cr[c]: c for c in cr} - def get(self, *, row=None, column=None): + def get(self, *, row: int | None = None, column: int | None = None) -> int: if row is not None: return self.rc[row] - else: - return self.cr[column] + return self.cr[column] + + def get_pivot(self, *, row: int | None = None, column: int | None = None) -> tuple[int, int]: + if row is not None: + return row, self.rc[row] + return self.cr[column], column - def set(self, *, row, column): + def set(self, *, row: int, column: int) -> None: self.rc[row] = column self.cr[column] = row - def has(self, *, row=None, column=None): + def has(self, *, row: int | None = None, column: int | None = None) -> bool: if row is not None: return row in self.rc - else: - return column in self.cr + return column in self.cr - def delete(self, *, row=None, column=None): + def delete(self, *, row: int | None = None, column: int | None = None) -> None: if row is not None: column = self.rc[row] del self.rc[row] @@ -36,210 +66,350 @@ def delete(self, *, row=None, column=None): del self.rc[row] del self.cr[column] - def __str__(self): - return str(self.rc) + def __repr__(self) -> str: + return repr(self.rc) -class LinearProgram: - def __init__(self, A, b, c, pivots: dict[int, int] | Pivots): - self.M = np.block([ - [np.atleast_2d(c), np.atleast_2d(0)], - [A, np.atleast_2d(b).T] - ]) - self.pivots = pivots if isinstance(pivots, Pivots) else Pivots(pivots) +class Tableau: + def __init__( + self, A: np.ndarray, b: np.ndarray, c: np.ndarray, pivots: Pivots | None = None, *, tolerance: float = 1e-6 + ) -> None: + self.M = np.block([[np.atleast_2d(c), np.atleast_2d(0)], [A, np.atleast_2d(b).T]]) + self.pivots = pivots if pivots else Pivots({}) + self.tolerance = tolerance @property - def n_vars(self): + def n_vars(self) -> int: return self.A.shape[1] @property - def n_constraints(self): + def n_constraints(self) -> int: return self.A.shape[0] @property - def A(self): + def A(self) -> np.ndarray: return self.M[1:, :-1] @property - def b(self): + def b(self) -> np.ndarray: return self.M[1:, -1] @property - def c(self): + def c(self) -> np.ndarray: return self.M[0, :-1] @property - def value(self): + def value(self) -> float: + """The value of the objective function.""" return -self.M[0, -1] - def __str__(self): - return f"{self.M}\n{self.pivots}" + @property + def solution(self) -> np.ndarray: + solution = np.zeros(self.n_vars, dtype=float) + for col in range(self.n_vars): + if self.pivots.has(column=col): + row = self.pivots.get(column=col) + solution[col] = self.b[row - 1] + return solution + def __repr__(self) -> str: + return f"Tableau:\n{self.M}\nPivots:\n{self.pivots}" -def zero_out_cj(program, column, *, tolerance): - row = program.pivots.get(column=column) - cj = program.c[column] - if is_zero(cj, tolerance): - return - program.M[0] -= cj*program.M[row] +def solve( + _type: ProgramType, objective: Expression, constraints: list[Constraint], *, max_iterations: int, tolerance: float +) -> tuple[float, np.ndarray]: + """ + The main entrypoint to the simplex algorithm. -def do_pivot(program, pivot): - row, column = pivot - program.M[row] /= program.M[row, column] + This function accepts a high-level representation of the LP and returns the optimal value and the solution. + """ + program = as_tableau(_type, objective, constraints, tolerance=tolerance) + return _solve(program, max_iterations=max_iterations, tolerance=tolerance) - for i in range(1, program.n_constraints+1): - if i == row: - continue - program.M[i] -= program.M[i, column]*program.M[row] +def _solve(program: Tableau, *, max_iterations: int, tolerance: float) -> tuple[float, np.ndarray]: + aux_program = create_phase_one_program(program) + run_simplex(aux_program) -def find_pivot(program, *, tolerance): - for j in range(program.n_vars): - if program.pivots.has(column=j): - continue - if program.c[j] < -tolerance: - min_value = math.inf - min_i = None - for i in range(1, program.n_constraints + 1): - if is_zero(program.M[i, j], tolerance): - continue - frac = program.M[i, -1]/program.M[i, j] - if program.M[i, j] > tolerance and (frac-min_value) < -tolerance: - min_value = frac - min_i = i - if min_i is not None: - return min_i, j - return None + # The auxiliary problem minimizes the sum of slack variables. + # If the sum is not zero, the original problem is infeasible. + if not is_zero(aux_program.value, tolerance): + msg = "The program is infeasible, try removing some constraints." + raise Infeasible(msg) + + if has_slack_variables_in_basis(aux_program, program.n_vars): + # Either the problem is degenerate (one or more original variables is zero) + # or there are redundant constraints (or both). + # Before we proceed with the second phase, we need to remove the redundant constraints and + # move the slack variables out of the base. + remove_redundant_constraints(aux_program, program.n_vars, tolerance) + move_pivots_from_slack_variables(aux_program, program.n_vars, tolerance) + + new_prog = create_phase_two_program(program, aux_program) + run_simplex(new_prog, max_iterations=max_iterations) + return new_prog.value, new_prog.solution + + +def as_tableau( + _type: ProgramType, objective: Expression, constraints: list[Constraint], *, tolerance: float +) -> Tableau: + """ + Convert the high-level representation of the LP to a Tableau backed by a numpy matrix. + + The LP is converted to its canonical form in the process. + """ + variables = get_variable_names((objective, *constraints)) + n_vars = len(variables) + n_constraints = len(constraints) + A = np.zeros((n_constraints, n_vars), dtype=float) + b = np.zeros(n_constraints, dtype=float) + c = np.zeros(n_vars, dtype=float) + constraint_types = np.zeros(n_constraints, dtype=int) + + # Create the objective function row + coeffs, _ = get_variable_coeffs(objective) + for name, coeff in coeffs.items(): + i = variables.index(name) + c[i] = coeff + + constraint_map = { + Equal: 0, + GreaterOrEqual: 1, + LessOrEqual: -1, + } + # Create the constraints + for row, constraint in enumerate(constraints): + coeffs, _b = get_variable_coeffs(constraint) + for name, coeff in coeffs.items(): + column = variables.index(name) + A[row, column] = coeff + b[row] = _b + constraint_types[row] = constraint_map[type(constraint)] + + A, b, c = canonicalize(_type, A, b, c, constraint_types) + return Tableau(A, b, c, tolerance=tolerance) + + +def canonicalize( + _type: ProgramType, A: np.ndarray, b: np.ndarray, c: np.ndarray, constraint_types: list[Literal[-1, 0, 1]] +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Convert the LP to the canonical form: + + min c^T*x + s.t. + Ax <= b + x >= 0 + b >= 0 + """ + if _type == ProgramType.MAX: + c *= -1 + # convert <= to >= + sel = constraint_types == -1 + b[sel] *= -1 + A[sel] *= -1 + constraint_types[sel] = 1 -def get_solution(program): - n = program.n_vars - solution = np.zeros(n, dtype=float) - for j in range(n): - if program.pivots.has(column=j): - row = program.pivots.get(column=j) - solution[j] = program.b[row-1] - return solution + # convert >= to == + sel = constraint_types == 1 + n = np.sum(sel) + if n > 0: + A = np.hstack((A, np.zeros((A.shape[0], n)))) + c = np.hstack((c, np.zeros(n))) + A[sel, -n:] = -np.eye(n) + constraint_types[sel] = 0 + # convert -b to b + sel = b < 0 + b[sel] *= -1 + A[sel] *= -1 -def run_simplex(program, *, max_iterations=math.inf, tolerance=1e-6): - iterations = 0 - while (pivot := find_pivot(program, tolerance=tolerance)) is not None: - old_pivot = program.pivots.get(row=pivot[0]) - do_pivot(program, pivot) - program.pivots.delete(column=old_pivot) - program.pivots.set(row=pivot[0], column=pivot[1]) - zero_out_cj(program, pivot[1], tolerance=tolerance) + return A, b, c - iterations += 1 - if iterations > max_iterations: - return - if np.any(program.c < -tolerance): - raise Unbounded("The program is unbounded, try adding more constraints.") +def has_slack_variables_in_basis(aux_program: Tableau, n_original_vars: int) -> bool: + return any(j >= n_original_vars for j in aux_program.pivots.cr) -def solve(program, *, max_iterations, tolerance): - aux_program = aux_problem(program, tolerance=tolerance) - run_simplex(aux_program, tolerance=tolerance) +def create_phase_one_program(program: Tableau) -> Tableau: + """ + Create the auxiliary problem for the first phase of the simplex algorithm. - if not is_zero(aux_program.value, tolerance): - raise Infeasible("The program is infeasible, try removing some constraints.") - - if any(j >= program.n_vars for j in aux_program.pivots.cr): - zeros = is_zero(aux_program.A[:, :program.n_vars], tolerance) - sel = np.all(zeros, axis=1) - sel = np.concatenate(([False], sel)) - - if np.any(sel): - aux_program.M = aux_program.M[~sel] - for i, s in enumerate(sel): - if s: - col = aux_program.pivots.get(row=i) - aux_program.pivots.delete(column=col) - - pivots = [] - for col in aux_program.pivots.cr: - if col >= program.n_vars: - pivots.append((aux_program.pivots.get(column=col), col)) - - for p in pivots: - row, col = p - for candidate_col in range(program.n_vars): - if candidate_col in aux_program.pivots.cr: - continue - if is_zero(aux_program.M[row, candidate_col], tolerance): - continue - do_pivot(aux_program, (row, candidate_col)) - aux_program.pivots.delete(column=col) - aux_program.pivots.set(row=row, column=candidate_col) - break - - A = aux_program.A[:, :program.n_vars] - b = np.copy(aux_program.b) - c = np.copy(program.c) - pivots = aux_program.pivots - new_prog = LinearProgram(A, b, c, pivots) + This adds additional slack variables to create a valid basis. + The sum of the slack variables is then minimized. - for column in pivots.cr: - if not is_zero(new_prog.c[column], tolerance): - zero_out_cj(new_prog, column, tolerance=tolerance) + For example, this problem: - run_simplex(new_prog, max_iterations=max_iterations, tolerance=tolerance) - return new_prog.value, get_solution(new_prog) + min 2*x + y + 3*z s.t. + x - y == 4 + y + 2*z == 2 + is converted to: -def aux_problem(program, *, tolerance): + min s_1 + s_2 s.t. + x - y + s_1 == 4 + y + 2*z + s_2 == 2 + """ m = program.n_constraints A = np.copy(np.c_[program.A, np.eye(m)]) b = np.copy(program.b) c = np.concatenate((np.zeros(program.n_vars), np.ones(m))) - pivots = {program.n_vars+i: i+1 for i in range(m)} + pivots = Pivots({program.n_vars + i: i + 1 for i in range(m)}) - aux_program = LinearProgram(A, b, c, pivots) - for j in range(program.n_vars, program.n_vars+m): - zero_out_cj(aux_program, j, tolerance=tolerance) + aux_program = Tableau(A, b, c, pivots) + for j in range(program.n_vars, program.n_vars + m): + zero_out_cj(aux_program, j) return aux_program -def is_zero(v, tol): - return np.isclose(v, 0, rtol=0, atol=tol) +def create_phase_two_program(original_program: Tableau, aux_program: Tableau) -> Tableau: + """ + Create the phase two program using the optimized phase one program. + To do so, we discard the slack variables used by the auxiliary problem and + restore the original objective function. + """ + # Discard slack variables + A = aux_program.A[:, : original_program.n_vars] + b = np.copy(aux_program.b) + # Use the orginal objective function + c = np.copy(original_program.c) + pivots = aux_program.pivots + new_prog = Tableau(A, b, c, pivots) -def canonicalize(type, A, b, c, constraint_types): - if type == 'max': - c *= -1 + # The simplex algorithm expects reduced costs to be zeroed out + for column in pivots.cr: + zero_out_cj(new_prog, column) + + return new_prog + + +def remove_redundant_constraints(aux_program: Tableau, n_original_vars: int, tolerance: float) -> None: + """ + Remove redundant constraints from the auxiliary problem. + + A constraint is redundant if all the coefficients of the original variables are zero. + For example in this tableau, the last constraint is redundant: + + x1 x2 x3 s1 s2 + | 1 0 0 | 0 0 | + | 0 1 0 | 0 0 | + | 0 0 0 | 0 1 | + + """ + zeros = is_zero(aux_program.A[:, :n_original_vars], tolerance) + sel = np.all(zeros, axis=1) + sel = np.concatenate(([False], sel)) + + # Delete pivots in the same rows as the redundant constraints + aux_program.M = aux_program.M[~sel] + for i in sel.nonzero()[0]: + col = aux_program.pivots.get(row=i) + aux_program.pivots.delete(column=col) + + +def move_pivots_from_slack_variables(aux_program: Tableau, n_original_vars: int, tolerance: float) -> None: + """ + Move pivots from the slack variables to the original variables. + + This happens when the original solution is degenrate (one or more original variables is zero). + For example, in this tableau, the pivot is in the slack variable s2: + + x1 x2 x3 s1 s2 + | 1 0 0 | 0 0 | + | 0 1 4 | 0 0 | + | 0 0 2 | 0 1 | + + After moving the pivot, the tableau becomes: + + x1 x2 x3' s1 s2' + | 1 0 0 | 0 0 | + | 0 1 0 | 0 -2 | + | 0 0 1 | 0 1/2 | + + """ + # Find the pivots in the slack variables + pivots = [aux_program.pivots.get_pivot(column=col) for col in aux_program.pivots.cr if (col >= n_original_vars)] + + # Move pivots out of the slack variables + for p in pivots: + row, col = p + # Search through the original variables + for candidate_col in range(n_original_vars): + # Already in the base + if candidate_col in aux_program.pivots.cr: + continue + # The new pivot location is zero, not possible + if is_zero(aux_program.M[row, candidate_col], tolerance): + continue + # Set as the new pivot + aux_program.pivots.delete(column=col) + aux_program.pivots.set(row=row, column=candidate_col) + # Update the tableau + normalize_pivot(aux_program, (row, candidate_col)) + zero_out_cj(aux_program, candidate_col) + break + + +def zero_out_cj(program: Tableau, column: int) -> None: + """Make the reduced cost in the given column zero by subtracting a multiple of the pivot row.""" + row = program.pivots.get(column=column) + cj = program.c[column] + if is_zero(cj, program.tolerance): + return + program.M[0] -= cj * program.M[row] - # convert <= to >= - sel = constraint_types == -1 - b[sel] *= -1 - A[sel] *= -1 - constraint_types[sel] = 1 - # convert >= to == - sel = constraint_types == 1 - n = np.sum(sel) +def normalize_pivot(program: Tableau, pivot: tuple[int, int]) -> None: + """Make the pivot equal to 1 by dividing the row by the pivot value.""" + row, column = pivot + program.M[row] /= program.M[row, column] - if n > 0: - print(A, np.zeros(A.shape[0], n)) - A = np.hstack((A, np.zeros((A.shape[0], n)))) - c = np.hstack((c, np.zeros(n))) - A[sel, -n:] = -np.eye(n) - constraint_types[sel] = 0 + for i in range(1, program.n_constraints + 1): + if i == row: + continue + program.M[i] -= program.M[i, column] * program.M[row] - # TODO: variable bounds - # for j in range(A.shape[1]): - # AA = np.r_[AA, np.zeros(A.shape[1])] - # b = np.c_[b, 0] - # AA[-1, j] = 1 - # AA[-1, -2] = -1 - # AA[-1, -1] = 1 - # c = np.c_[c, 0, 0] - # convert -b to b - sel = b < 0 - b[sel] *= -1 - A[sel] *= -1 +@suppress_divide_by_zero_warning +def find_pivot(program: Tableau) -> tuple[int, int] | None: + """Find a new pivot which will improve the objective function.""" + for j in range(program.n_vars): + if program.pivots.has(column=j): + continue + # cj must be negative + if program.c[j] < -program.tolerance: + # Pivot candidates must positive + positive = program.M[1:, j] > program.tolerance + frac = program.M[1:, -1] / program.M[1:, j] + frac[~positive] = np.inf + # Pick a pivot which minimizes b_i/p_ij + min_i = np.argmin(frac) + if frac[min_i] != np.inf: + return min_i + 1, j + return None - return A, b, c + +def run_simplex(program: Tableau, *, max_iterations: int | None = math.inf) -> None: + """Run the simplex algorithm starting from a feasible solution.""" + iterations = 0 + while (pivot := find_pivot(program)) is not None: + row, column = pivot + program.pivots.delete(column=program.pivots.get(row=row)) + program.pivots.set(row=row, column=column) + normalize_pivot(program, pivot) + zero_out_cj(program, column) + + iterations += 1 + if iterations > max_iterations: + return + + if np.any(program.c < -program.tolerance): + msg = "The program is unbounded, try adding more constraints." + raise Unbounded(msg) + + +def is_zero(v: float, tol: float) -> bool: + return np.isclose(v, 0, rtol=0, atol=tol) diff --git a/pyproject.toml b/pyproject.toml index b1ec9ae..218adf3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,8 +1,6 @@ [project] name = "pivotal-solver" -authors = [ - {name = "Tomas Roun", email = "tomas.roun8@gmail.com"}, -] +authors = [{ name = "Tomas Roun", email = "tomas.roun8@gmail.com" }] description = "High-level Linear Programming solver using the Simplex algorithm" readme = "README.md" requires-python = ">=3.10" @@ -18,7 +16,7 @@ keywords = [ "simplex", "simplex algorithm", ] -license = {file = "LICENSE"} +license = { file = "LICENSE" } classifiers = [ "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", @@ -27,32 +25,46 @@ classifiers = [ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Topic :: Scientific/Engineering", - "Topic :: Software Development :: Libraries" + "Topic :: Software Development :: Libraries", ] dependencies = ["numpy"] dynamic = ["version"] [project.optional-dependencies] -dev = [ - "pytest", - "flake8", - "isort", - "pre-commit" -] +dev = ["pytest", "ruff", "isort", "pre-commit"] [project.urls] Homepage = "https://github.com/tomasr8/pivotal" Github = "https://github.com/tomasr8/pivotal" [tool.isort] -line_length=120 -lines_after_imports=2 +line_length = 120 +lines_after_imports = 2 + +[tool.ruff.lint] +ignore = [ + "D", + "ANN002", + "ANN003", + "ANN101", + "RET503", + "COM812", + "ISC001", + "N802", + "N803", + "N806", + "N818", + "S101", +] + +[tool.ruff] +line-length = 120 [tool.setuptools] packages = ["pivotal"] [tool.setuptools.dynamic] -version = {attr = "pivotal.__version__"} +version = { attr = "pivotal.__version__" } [build-system] requires = ["setuptools", "setuptools-scm"] diff --git a/tests/test_api.py b/tests/test_api.py index ce8701e..8181083 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,9 +1,8 @@ -import math - import numpy as np from numpy.testing import assert_allclose -from pivotal.api import LinearProgram, Variable +from pivotal.api import maximize, minimize +from pivotal.expressions import Variable def assert_solution_almost_equal(expected, actual): @@ -31,47 +30,47 @@ def test_variables(): assert 0 + x is x assert x - 0 is x - assert 1*x is x - assert x*1 is x + assert 1 * x is x + assert x * 1 is x assert (-x).coeff == -1 - assert (2*x).coeff == 2 - assert (-2*x).coeff == -2 - assert (x*2).coeff == 2 - assert (x*-2).coeff == -2 + assert (2 * x).coeff == 2 + assert (-2 * x).coeff == -2 + assert (x * 2).coeff == 2 + assert (x * -2).coeff == -2 - assert (5*x).name == x.name + assert (5 * x).name == x.name def test_sum(): x = Variable() y = Variable() - assert_array_equal((x+1).elts, (x, 1)) - assert_array_equal((x-3).elts, (x, -3)) - assert_array_equal((1+x).elts, (1, x)) - assert_array_equal((3+x).elts, (3, x)) + assert_array_equal((x + 1).elts, (x, 1)) + assert_array_equal((x - 3).elts, (x, -3)) + assert_array_equal((1 + x).elts, (1, x)) + assert_array_equal((3 + x).elts, (3, x)) - s = 1-x + s = 1 - x assert_equal(s.elts[0], 1) assert s.elts[1].coeff == -1 - s = -x+1 + s = -x + 1 assert s.elts[0].coeff == -1 assert_equal(s.elts[1], 1) - s = 1+x+2 + s = 1 + x + 2 assert_equal(s.elts[0], 1) assert s.elts[1].coeff == 1 assert_equal(s.elts[2], 2) - s = x+y + s = x + y assert s.elts[0].coeff == 1 assert s.elts[1].coeff == 1 assert s.elts[0].name == x.name assert s.elts[1].name == y.name - s = x-y+3*x + s = x - y + 3 * x assert s.elts[0].coeff == 1 assert s.elts[1].coeff == -1 assert s.elts[2].coeff == 3 @@ -79,25 +78,25 @@ def test_sum(): assert s.elts[1].name == y.name assert s.elts[2].name == x.name - s = 5*(x+2*y) + s = 5 * (x + 2 * y) assert s.elts[0].coeff == 5 assert s.elts[1].coeff == 10 assert s.elts[0].name == x.name assert s.elts[1].name == y.name - s = (x+2*y)*5 + s = (x + 2 * y) * 5 assert s.elts[0].coeff == 5 assert s.elts[1].coeff == 10 assert s.elts[0].name == x.name assert s.elts[1].name == y.name - s = x+y - assert 1*s is s - assert s*1 is s + s = x + y + assert 1 * s is s + assert s * 1 is s x1 = Variable("x") x2 = Variable("x") - assert (x1+3*x2).coeff == 4 + assert (x1 + 3 * x2).coeff == 4 def test_abs(): @@ -118,69 +117,68 @@ def test_mixing_task(): m = Variable("m") z = Variable("z") - program = LinearProgram.minimize( - 25*b + 50*m + 80*z - ).such_that( - 2*b + m + z >= 8, - 2*b + 6*m + z >= 16, - b + 3*m + 6*z >= 8, + solution = minimize( + 25 * b + 50 * m + 80 * z, + ( + 2 * b + m + z >= 8, + 2 * b + 6 * m + z >= 16, + b + 3 * m + 6 * z >= 8, + ), ) - assert_solution_almost_equal((160, [3.2, 1.6, 0]), program.optimize(max_iterations=math.inf, tolerance=1e-6)) + assert_solution_almost_equal((160, [3.2, 1.6, 0]), solution) def test_assignment_problem(): # Assignment problem # technically ILP, but the LP relaxation is always optimal # opt.pdf page 149, section 11.4.1 - C = np.array([ - [1, 2, 3], - [2, 3, 1], - [10, 2, 1] - ], dtype=float) + C = np.array([[1, 2, 3], [2, 3, 1], [10, 2, 1]], dtype=float) X = [[Variable(f"x_{i}{j}") for j in range(3)] for i in range(3)] - program = LinearProgram.minimize( - sum(sum(C[i, j] * X[i][j] for j in range(3)) for i in range(3)) - ).such_that( - sum(X[0][j] for j in range(3)) == 1, - sum(X[1][j] for j in range(3)) == 1, - sum(X[2][j] for j in range(3)) == 1, - sum(X[i][0] for i in range(3)) == 1, - sum(X[i][1] for i in range(3)) == 1, - sum(X[i][2] for i in range(3)) == 1, + solution = minimize( + sum(sum(C[i, j] * X[i][j] for j in range(3)) for i in range(3)), + ( + sum(X[0][j] for j in range(3)) == 1, + sum(X[1][j] for j in range(3)) == 1, + sum(X[2][j] for j in range(3)) == 1, + sum(X[i][0] for i in range(3)) == 1, + sum(X[i][1] for i in range(3)) == 1, + sum(X[i][2] for i in range(3)) == 1, + ), ) - assert_solution_almost_equal((4, [1, 0, 0, 0, 0, 1, 0, 1, 0]), - program.optimize(max_iterations=math.inf, tolerance=1e-6)) + assert_solution_almost_equal((4, [1, 0, 0, 0, 0, 1, 0, 1, 0]), solution) def test_redundant_constraints(): # OCW, page 82 X = [Variable(f"x{i+1}") for i in range(3)] - program = LinearProgram.maximize( - X[0] + 2*X[1] - X[2] - ).such_that( - 2*X[0] - X[1] + X[2] == 12, - -X[0] + 2*X[1] + X[2] == 10, - X[0] + X[1] + 2*X[2] == 22, + solution = maximize( + X[0] + 2 * X[1] - X[2], + ( + 2 * X[0] - X[1] + X[2] == 12, + -X[0] + 2 * X[1] + X[2] == 10, + X[0] + X[1] + 2 * X[2] == 22, + ), ) - assert_solution_almost_equal((98/3, [34/3, 32/3, 0]), program.optimize(max_iterations=math.inf, tolerance=1e-6)) + assert_solution_almost_equal((98 / 3, [34 / 3, 32 / 3, 0]), solution) def test_degenerate_solution(): # OCW, page 84 X = [Variable(f"x{i+1}") for i in range(4)] - program = LinearProgram.minimize( - X[0] + X[1] + 3*X[2] - ).such_that( - X[0] + 5*X[1] + X[2] + X[3] == 7, - X[0] - X[1] + X[2] == 5, - 0.5*X[0] - 2*X[1] + X[2] == 5, + solution = minimize( + X[0] + X[1] + 3 * X[2], + ( + X[0] + 5 * X[1] + X[2] + X[3] == 7, + X[0] - X[1] + X[2] == 5, + 0.5 * X[0] - 2 * X[1] + X[2] == 5, + ), ) - assert_solution_almost_equal((15, [0, 0, 5, 2]), program.optimize(max_iterations=math.inf, tolerance=1e-6)) + assert_solution_almost_equal((15, [0, 0, 5, 2]), solution) diff --git a/tests/test_simplex.py b/tests/test_simplex.py index 04c221f..d5a39d8 100644 --- a/tests/test_simplex.py +++ b/tests/test_simplex.py @@ -4,123 +4,125 @@ import pytest from numpy.testing import assert_array_equal -from pivotal.simplex import LinearProgram, canonicalize, do_pivot, find_pivot, get_solution, solve, zero_out_cj +from pivotal.simplex import Pivots, Tableau, _solve, canonicalize, find_pivot, normalize_pivot, zero_out_cj def test_zero_out_cj_1(): A = np.eye(2) b = np.ones(2) c = np.zeros(2) - program = LinearProgram(A, b, c, {0: 1, 1: 2}) + program = Tableau(A, b, c, Pivots({0: 1, 1: 2})) M = np.copy(program.M) - zero_out_cj(program, 0, tolerance=1e-6) + zero_out_cj(program, 0) assert_array_equal(M, program.M) - zero_out_cj(program, 1, tolerance=1e-6) + zero_out_cj(program, 1) assert_array_equal(M, program.M) def test_zero_out_cj_2(): A = np.eye(2) b = np.ones(2) - c = 2*np.ones(2) - program = LinearProgram(A, b, c, {0: 1, 1: 2}) - zero_out_cj(program, 0, tolerance=1e-6) - assert_array_equal(program.M, np.array([ - [0, 2, -2], - [1, 0, 1], - [0, 1, 1], - ])) - zero_out_cj(program, 1, tolerance=1e-6) - assert_array_equal(program.M, np.array([ - [0, 0, -4], - [1, 0, 1], - [0, 1, 1], - ])) - - -def test_do_pivot_1(): - A = np.array([ - [1, 0], - [0, 1] - ], dtype=float) + c = 2 * np.ones(2) + program = Tableau(A, b, c, Pivots({0: 1, 1: 2})) + zero_out_cj(program, 0) + assert_array_equal( + program.M, + np.array( + [ + [0, 2, -2], + [1, 0, 1], + [0, 1, 1], + ] + ), + ) + zero_out_cj(program, 1) + assert_array_equal( + program.M, + np.array( + [ + [0, 0, -4], + [1, 0, 1], + [0, 1, 1], + ] + ), + ) + + +def test_normalize_pivot_1(): + A = np.array([[1, 0], [0, 1]], dtype=float) b = [3, 5] c = [1, 1] - program = LinearProgram(A, b, c, {}) - do_pivot(program, (1, 0)) - - assert_array_equal(program.M, np.array([ - [1, 1, 0], - [1, 0, 3], - [0, 1, 5], - ])) - - -def test_do_pivot_2(): - A = np.array([ - [12, 3], - [0, 2] - ], dtype=float) + program = Tableau(A, b, c, Pivots({})) + normalize_pivot(program, (1, 0)) + + assert_array_equal( + program.M, + np.array( + [ + [1, 1, 0], + [1, 0, 3], + [0, 1, 5], + ] + ), + ) + + +def test_normalize_pivot_2(): + A = np.array([[12, 3], [0, 2]], dtype=float) b = [3, 5] c = [1, 1] - program = LinearProgram(A, b, c, {}) - do_pivot(program, (1, 1)) + program = Tableau(A, b, c, Pivots({})) + normalize_pivot(program, (1, 1)) - assert_array_equal(program.M, np.array([ - [1, 1, 0], - [4, 1, 1], - [-8, 0, 3], - ])) + assert_array_equal( + program.M, + np.array( + [ + [1, 1, 0], + [4, 1, 1], + [-8, 0, 3], + ] + ), + ) def test_find_pivot_no_pivot_cj_positive(): - A = np.array([ - [1, 0, 2], - [0, 1, 7] - ], dtype=float) + A = np.array([[1, 0, 2], [0, 1, 7]], dtype=float) b = [3, 5] c = [1, 1, 1] - program = LinearProgram(A, b, c, {}) - pivot = find_pivot(program, tolerance=1e-6) + program = Tableau(A, b, c, Pivots({})) + pivot = find_pivot(program) assert pivot is None def test_find_pivot_no_pivot_negative_column(): - A = np.array([ - [1, 0, -2], - [0, 1, -7] - ], dtype=float) + A = np.array([[1, 0, -2], [0, 1, -7]], dtype=float) b = [3, 5] c = [1, 1, -3] - program = LinearProgram(A, b, c, {}) - pivot = find_pivot(program, tolerance=1e-6) + program = Tableau(A, b, c, Pivots({})) + pivot = find_pivot(program) assert pivot is None def test_find_pivot(): - A = np.array([ - [1, 0, 2], - [0, 1, 7] - ], dtype=float) + A = np.array([[1, 0, 2], [0, 1, 7]], dtype=float) b = [3, 5] c = [1, 1, -2] - program = LinearProgram(A, b, c, {0: 1, 1: 2}) - pivot = find_pivot(program, tolerance=1e-6) + program = Tableau(A, b, c, Pivots({0: 1, 1: 2})) + pivot = find_pivot(program) assert_array_equal(pivot, [2, 2]) def test_get_solution_1(): - A = np.array([ - [1, 0], - [0, 1] - ], dtype=float) + A = np.array([[1, 0], [0, 1]], dtype=float) b = [3, 5] c = [1, 1] - program = LinearProgram(A, b, c, {0: 1, 1: 2}) - solution = get_solution(program) + program = Tableau(A, b, c, Pivots({0: 1, 1: 2})) + solution = program.solution assert_array_equal(solution, [3, 5]) - program = LinearProgram(A, b, c, {1: 2, 0: 1}) - solution = get_solution(program) + program = Tableau(A, b, c, Pivots({1: 2, 0: 1})) + solution = program.solution assert_array_equal(solution, [3, 5]) @@ -130,7 +132,7 @@ def test_canonicalize_1(): c = np.array([3, 4]) types = np.array([0, 0]) - Ac, bc, cc = canonicalize('min', A, b, c, types) + Ac, bc, cc = canonicalize("min", A, b, c, types) assert_array_equal(A, Ac) assert_array_equal(b, bc) assert_array_equal(c, cc) @@ -142,10 +144,8 @@ def test_canonicalize_2(): c = np.array([3, 4]) types = np.array([0]) - Ac, bc, cc = canonicalize('min', A, b, c, types) - assert_array_equal(Ac, [ - [-1, -1] - ]) + Ac, bc, cc = canonicalize("min", A, b, c, types) + assert_array_equal(Ac, [[-1, -1]]) assert_array_equal(bc, [2]) assert_array_equal(cc, [3, 4]) @@ -156,10 +156,8 @@ def test_canonicalize_3(): c = np.array([3, 4]) types = np.array([1]) - Ac, bc, cc = canonicalize('min', A, b, c, types) - assert_array_equal(Ac, [ - [1, 1, -1] - ]) + Ac, bc, cc = canonicalize("min", A, b, c, types) + assert_array_equal(Ac, [[1, 1, -1]]) assert_array_equal(bc, [2]) assert_array_equal(cc, [3, 4, 0]) @@ -170,24 +168,25 @@ def test_canonicalize_4(): c = np.array([3, 4]) types = np.array([-1]) - Ac, bc, cc = canonicalize('min', A, b, c, types) - assert_array_equal(Ac, [ - [1, 1, 1] - ]) + Ac, bc, cc = canonicalize("min", A, b, c, types) + assert_array_equal(Ac, [[1, 1, 1]]) assert_array_equal(bc, [2]) assert_array_equal(cc, [3, 4, 0]) def test_program_1(): - A = np.array([ - [1, 3], - [4, 1], - ], dtype=float) + A = np.array( + [ + [1, 3], + [4, 1], + ], + dtype=float, + ) b = np.array([6, 5], dtype=float) c = np.array([1, 2], dtype=float) - value, X = solve(LinearProgram(A, b, c, {}), max_iterations=math.inf, tolerance=1e-6) + value, X = _solve(Tableau(A, b, c, Pivots({})), max_iterations=math.inf, tolerance=1e-6) X_true = [0.81818182, 1.72727273] assert value == pytest.approx(4.27272727) for x, xt in zip(X, X_true): @@ -195,15 +194,18 @@ def test_program_1(): def test_program_2(): - A = np.array([ - [3, 2, 1], - [1, 2, 2], - ], dtype=float) + A = np.array( + [ + [3, 2, 1], + [1, 2, 2], + ], + dtype=float, + ) b = np.array([10, 15], dtype=float) c = np.array([-20, -30, -40], dtype=float) - value, X = solve(LinearProgram(A, b, c, {}), max_iterations=math.inf, tolerance=1e-6) + value, X = _solve(Tableau(A, b, c, Pivots({})), max_iterations=math.inf, tolerance=1e-6) X_true = [1, 0, 7] assert value == pytest.approx(-300) for x, xt in zip(X, X_true):