OpenSoure_Repos/TheAlgorithms-Python
1
0
Fork 0
mirror of https://hub.njuu.cf/TheAlgorithms/Python.git synced 2023-10-11 13:06:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix simplex.py (#8843)

Browse Source 
* changes to accommodate special case

* changed n_slack calculation method

* fix precommit typehints

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* n_art_vars inputs

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix: docstrings and typehints

* fix: doctest issues when running code

* additional check and doctests

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix ruff

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix whitespace

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
This commit is contained in:
Ilkin Mengusoglu 2023-08-17 22:34:53 +01:00 committed by GitHub
parent f6b12420ce
commit a207187ddb
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 127 additions and 100 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

227
linear_programming/simplex.py
View File

@ -20,40 +20,60 @@ import numpy as np
class Tableau: class Tableau:
"""Operate on simplex tableaus """Operate on simplex tableaus
>>> t = Tableau(np.array([[-1,-1,0,0,-1],[1,3,1,0,4],[3,1,0,1,4.]]), 2) >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4]]), 2, 2)
Traceback (most recent call last):
...
TypeError: Tableau must have type float64
>>> Tableau(np.array([[-1,-1,0,0,-1],[1,3,1,0,4],[3,1,0,1,4.]]), 2, 2)
Traceback (most recent call last): Traceback (most recent call last):
... ...
ValueError: RHS must be > 0 ValueError: RHS must be > 0
>>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]), -2, 2)
Traceback (most recent call last):
...
ValueError: number of (artificial) variables must be a natural number
""" """
def __init__(self, tableau: np.ndarray, n_vars: int) -> None: # Max iteration number to prevent cycling
maxiter = 100
def __init__(
self, tableau: np.ndarray, n_vars: int, n_artificial_vars: int
) -> None:
if tableau.dtype != "float64":
raise TypeError("Tableau must have type float64")
# Check if RHS is negative # Check if RHS is negative
if np.any(tableau[:, -1], where=tableau[:, -1] < 0): if not (tableau[:, -1] >= 0).all():
raise ValueError("RHS must be > 0") raise ValueError("RHS must be > 0")
if n_vars < 2 or n_artificial_vars < 0:
raise ValueError(
"number of (artificial) variables must be a natural number"
)
self.tableau = tableau self.tableau = tableau
self.n_rows, _ = tableau.shape self.n_rows, n_cols = tableau.shape
# Number of decision variables x1, x2, x3... # Number of decision variables x1, x2, x3...
self.n_vars = n_vars self.n_vars, self.n_artificial_vars = n_vars, n_artificial_vars
# Number of artificial variables to be minimised
self.n_art_vars = len(np.where(tableau[self.n_vars : -1] == -1)[0])
# 2 if there are >= or == constraints (nonstandard), 1 otherwise (std) # 2 if there are >= or == constraints (nonstandard), 1 otherwise (std)
self.n_stages = (self.n_art_vars > 0) + 1 self.n_stages = (self.n_artificial_vars > 0) + 1
# Number of slack variables added to make inequalities into equalities # Number of slack variables added to make inequalities into equalities
self.n_slack = self.n_rows - self.n_stages self.n_slack = n_cols - self.n_vars - self.n_artificial_vars - 1
# Objectives for each stage # Objectives for each stage
self.objectives = ["max"] self.objectives = ["max"]
# In two stage simplex, first minimise then maximise # In two stage simplex, first minimise then maximise
if self.n_art_vars: if self.n_artificial_vars:
self.objectives.append("min") self.objectives.append("min")
self.col_titles = [""] self.col_titles = self.generate_col_titles()
# Index of current pivot row and column # Index of current pivot row and column
self.row_idx = None self.row_idx = None
@ -62,48 +82,39 @@ class Tableau:
# Does objective row only contain (non)-negative values? # Does objective row only contain (non)-negative values?
self.stop_iter = False self.stop_iter = False
@staticmethod def generate_col_titles(self) -> list[str]:
def generate_col_titles(*args: int) -> list[str]:
"""Generate column titles for tableau of specific dimensions """Generate column titles for tableau of specific dimensions
>>> Tableau.generate_col_titles(2, 3, 1) >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]),
['x1', 'x2', 's1', 's2', 's3', 'a1', 'RHS'] ... 2, 0).generate_col_titles()
['x1', 'x2', 's1', 's2', 'RHS']
>>> Tableau.generate_col_titles() >>> Tableau(np.array([[-1,-1,0,0,1],[1,3,1,0,4],[3,1,0,1,4.]]),
Traceback (most recent call last): ... 2, 2).generate_col_titles()
... ['x1', 'x2', 'RHS']
ValueError: Must provide n_vars, n_slack, and n_art_vars
>>> Tableau.generate_col_titles(-2, 3, 1)
Traceback (most recent call last):
...
ValueError: All arguments must be non-negative integers
""" """
if len(args) != 3: args = (self.n_vars, self.n_slack)
raise ValueError("Must provide n_vars, n_slack, and n_art_vars")
if not all(x >= 0 and isinstance(x, int) for x in args): # decision | slack
raise ValueError("All arguments must be non-negative integers") string_starts = ["x", "s"]
# decision | slack | artificial
string_starts = ["x", "s", "a"]
titles = [] titles = []
for i in range(3): for i in range(2):
for j in range(args[i]): for j in range(args[i]):
titles.append(string_starts[i] + str(j + 1)) titles.append(string_starts[i] + str(j + 1))
titles.append("RHS") titles.append("RHS")
return titles return titles
def find_pivot(self, tableau: np.ndarray) -> tuple[Any, Any]: def find_pivot(self) -> tuple[Any, Any]:
"""Finds the pivot row and column. """Finds the pivot row and column.
>>> t = Tableau(np.array([[-2,1,0,0,0], [3,1,1,0,6], [1,2,0,1,7.]]), 2) >>> Tableau(np.array([[-2,1,0,0,0], [3,1,1,0,6], [1,2,0,1,7.]]),
>>> t.find_pivot(t.tableau) ... 2, 0).find_pivot()
(1, 0) (1, 0)
""" """
objective = self.objectives[-1] objective = self.objectives[-1]
# Find entries of highest magnitude in objective rows # Find entries of highest magnitude in objective rows
sign = (objective == "min") - (objective == "max") sign = (objective == "min") - (objective == "max")
col_idx = np.argmax(sign * tableau[0, : self.n_vars]) col_idx = np.argmax(sign * self.tableau[0, :-1])
# Choice is only valid if below 0 for maximise, and above for minimise # Choice is only valid if below 0 for maximise, and above for minimise
if sign * self.tableau[0, col_idx] <= 0: if sign * self.tableau[0, col_idx] <= 0:
@ -117,15 +128,15 @@ class Tableau:
s = slice(self.n_stages, self.n_rows) s = slice(self.n_stages, self.n_rows)
# RHS # RHS
dividend = tableau[s, -1] dividend = self.tableau[s, -1]
# Elements of pivot column within slice # Elements of pivot column within slice
divisor = tableau[s, col_idx] divisor = self.tableau[s, col_idx]
# Array filled with nans # Array filled with nans
nans = np.full(self.n_rows - self.n_stages, np.nan) nans = np.full(self.n_rows - self.n_stages, np.nan)
# If element in pivot column is greater than zeron_stages, return # If element in pivot column is greater than zero, return
# quotient or nan otherwise # quotient or nan otherwise
quotients = np.divide(dividend, divisor, out=nans, where=divisor > 0) quotients = np.divide(dividend, divisor, out=nans, where=divisor > 0)
@ -134,18 +145,18 @@ class Tableau:
row_idx = np.nanargmin(quotients) + self.n_stages row_idx = np.nanargmin(quotients) + self.n_stages
return row_idx, col_idx return row_idx, col_idx
def pivot(self, tableau: np.ndarray, row_idx: int, col_idx: int) -> np.ndarray: def pivot(self, row_idx: int, col_idx: int) -> np.ndarray:
"""Pivots on value on the intersection of pivot row and column. """Pivots on value on the intersection of pivot row and column.
>>> t = Tableau(np.array([[-2,-3,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]), 2) >>> Tableau(np.array([[-2,-3,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]),
>>> t.pivot(t.tableau, 1, 0).tolist() ... 2, 2).pivot(1, 0).tolist()
... # doctest: +NORMALIZE_WHITESPACE ... # doctest: +NORMALIZE_WHITESPACE
[[0.0, 3.0, 2.0, 0.0, 8.0], [[0.0, 3.0, 2.0, 0.0, 8.0],
[1.0, 3.0, 1.0, 0.0, 4.0], [1.0, 3.0, 1.0, 0.0, 4.0],
[0.0, -8.0, -3.0, 1.0, -8.0]] [0.0, -8.0, -3.0, 1.0, -8.0]]
""" """
# Avoid changes to original tableau # Avoid changes to original tableau
piv_row = tableau[row_idx].copy() piv_row = self.tableau[row_idx].copy()
piv_val = piv_row[col_idx] piv_val = piv_row[col_idx]
@ -153,48 +164,47 @@ class Tableau:
piv_row *= 1 / piv_val piv_row *= 1 / piv_val
# Variable in pivot column becomes basic, ie the only non-zero entry # Variable in pivot column becomes basic, ie the only non-zero entry
for idx, coeff in enumerate(tableau[:, col_idx]): for idx, coeff in enumerate(self.tableau[:, col_idx]):
tableau[idx] += -coeff * piv_row self.tableau[idx] += -coeff * piv_row
tableau[row_idx] = piv_row self.tableau[row_idx] = piv_row
return tableau return self.tableau
def change_stage(self, tableau: np.ndarray) -> np.ndarray: def change_stage(self) -> np.ndarray:
"""Exits first phase of the two-stage method by deleting artificial """Exits first phase of the two-stage method by deleting artificial
rows and columns, or completes the algorithm if exiting the standard rows and columns, or completes the algorithm if exiting the standard
case. case.
>>> t = Tableau(np.array([ >>> Tableau(np.array([
... [3, 3, -1, -1, 0, 0, 4], ... [3, 3, -1, -1, 0, 0, 4],
... [2, 1, 0, 0, 0, 0, 0.], ... [2, 1, 0, 0, 0, 0, 0.],
... [1, 2, -1, 0, 1, 0, 2], ... [1, 2, -1, 0, 1, 0, 2],
... [2, 1, 0, -1, 0, 1, 2] ... [2, 1, 0, -1, 0, 1, 2]
... ]), 2) ... ]), 2, 2).change_stage().tolist()
>>> t.change_stage(t.tableau).tolist()
... # doctest: +NORMALIZE_WHITESPACE ... # doctest: +NORMALIZE_WHITESPACE
[[2.0, 1.0, 0.0, 0.0, 0.0, 0.0], [[2.0, 1.0, 0.0, 0.0, 0.0],
[1.0, 2.0, -1.0, 0.0, 1.0, 2.0], [1.0, 2.0, -1.0, 0.0, 2.0],
[2.0, 1.0, 0.0, -1.0, 0.0, 2.0]] [2.0, 1.0, 0.0, -1.0, 2.0]]
""" """
# Objective of original objective row remains # Objective of original objective row remains
self.objectives.pop() self.objectives.pop()
if not self.objectives: if not self.objectives:
return tableau return self.tableau
# Slice containing ids for artificial columns # Slice containing ids for artificial columns
s = slice(-self.n_art_vars - 1, -1) s = slice(-self.n_artificial_vars - 1, -1)
# Delete the artificial variable columns # Delete the artificial variable columns
tableau = np.delete(tableau, s, axis=1) self.tableau = np.delete(self.tableau, s, axis=1)
# Delete the objective row of the first stage # Delete the objective row of the first stage
tableau = np.delete(tableau, 0, axis=0) self.tableau = np.delete(self.tableau, 0, axis=0)
self.n_stages = 1 self.n_stages = 1
self.n_rows -= 1 self.n_rows -= 1
self.n_art_vars = 0 self.n_artificial_vars = 0
self.stop_iter = False self.stop_iter = False
return tableau return self.tableau
def run_simplex(self) -> dict[Any, Any]: def run_simplex(self) -> dict[Any, Any]:
"""Operate on tableau until objective function cannot be """Operate on tableau until objective function cannot be
@ -205,15 +215,29 @@ class Tableau:
ST: x1 + 3x2 <= 4 ST: x1 + 3x2 <= 4
3x1 + x2 <= 4 3x1 + x2 <= 4
>>> Tableau(np.array([[-1,-1,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]), >>> Tableau(np.array([[-1,-1,0,0,0],[1,3,1,0,4],[3,1,0,1,4.]]),
... 2).run_simplex() ... 2, 0).run_simplex()
{'P': 2.0, 'x1': 1.0, 'x2': 1.0} {'P': 2.0, 'x1': 1.0, 'x2': 1.0}
# Standard linear program with 3 variables:
Max: 3x1 + x2 + 3x3
ST: 2x1 + x2 + x3 ≤ 2
x1 + 2x2 + 3x3 ≤ 5
2x1 + 2x2 + x3 ≤ 6
>>> Tableau(np.array([
... [-3,-1,-3,0,0,0,0],
... [2,1,1,1,0,0,2],
... [1,2,3,0,1,0,5],
... [2,2,1,0,0,1,6.]
... ]),3,0).run_simplex() # doctest: +ELLIPSIS
{'P': 5.4, 'x1': 0.199..., 'x3': 1.6}
# Optimal tableau input: # Optimal tableau input:
>>> Tableau(np.array([ >>> Tableau(np.array([
... [0, 0, 0.25, 0.25, 2], ... [0, 0, 0.25, 0.25, 2],
... [0, 1, 0.375, -0.125, 1], ... [0, 1, 0.375, -0.125, 1],
... [1, 0, -0.125, 0.375, 1] ... [1, 0, -0.125, 0.375, 1]
... ]), 2).run_simplex() ... ]), 2, 0).run_simplex()
{'P': 2.0, 'x1': 1.0, 'x2': 1.0} {'P': 2.0, 'x1': 1.0, 'x2': 1.0}
# Non-standard: >= constraints # Non-standard: >= constraints
@ -227,7 +251,7 @@ class Tableau:
... [1, 1, 1, 1, 0, 0, 0, 0, 40], ... [1, 1, 1, 1, 0, 0, 0, 0, 40],
... [2, 1, -1, 0, -1, 0, 1, 0, 10], ... [2, 1, -1, 0, -1, 0, 1, 0, 10],
... [0, -1, 1, 0, 0, -1, 0, 1, 10.] ... [0, -1, 1, 0, 0, -1, 0, 1, 10.]
... ]), 3).run_simplex() ... ]), 3, 2).run_simplex()
{'P': 70.0, 'x1': 10.0, 'x2': 10.0, 'x3': 20.0} {'P': 70.0, 'x1': 10.0, 'x2': 10.0, 'x3': 20.0}
# Non standard: minimisation and equalities # Non standard: minimisation and equalities
@ -235,73 +259,76 @@ class Tableau:
ST: 2x1 + x2 = 12 ST: 2x1 + x2 = 12
6x1 + 5x2 = 40 6x1 + 5x2 = 40
>>> Tableau(np.array([ >>> Tableau(np.array([
... [8, 6, 0, -1, 0, -1, 0, 0, 52], ... [8, 6, 0, 0, 52],
... [1, 1, 0, 0, 0, 0, 0, 0, 0], ... [1, 1, 0, 0, 0],
... [2, 1, 1, 0, 0, 0, 0, 0, 12], ... [2, 1, 1, 0, 12],
... [2, 1, 0, -1, 0, 0, 1, 0, 12], ... [6, 5, 0, 1, 40.],
... [6, 5, 0, 0, 1, 0, 0, 0, 40], ... ]), 2, 2).run_simplex()
... [6, 5, 0, 0, 0, -1, 0, 1, 40.]
... ]), 2).run_simplex()
{'P': 7.0, 'x1': 5.0, 'x2': 2.0} {'P': 7.0, 'x1': 5.0, 'x2': 2.0}
# Pivot on slack variables
Max: 8x1 + 6x2
ST: x1 + 3x2 <= 33
4x1 + 2x2 <= 48
2x1 + 4x2 <= 48
x1 + x2 >= 10
x1 >= 2
>>> Tableau(np.array([
... [2, 1, 0, 0, 0, -1, -1, 0, 0, 12.0],
... [-8, -6, 0, 0, 0, 0, 0, 0, 0, 0.0],
... [1, 3, 1, 0, 0, 0, 0, 0, 0, 33.0],
... [4, 2, 0, 1, 0, 0, 0, 0, 0, 60.0],
... [2, 4, 0, 0, 1, 0, 0, 0, 0, 48.0],
... [1, 1, 0, 0, 0, -1, 0, 1, 0, 10.0],
... [1, 0, 0, 0, 0, 0, -1, 0, 1, 2.0]
... ]), 2, 2).run_simplex() # doctest: +ELLIPSIS
{'P': 132.0, 'x1': 12.000... 'x2': 5.999...}
""" """
# Stop simplex algorithm from cycling. # Stop simplex algorithm from cycling.
for _ in range(100): for _ in range(Tableau.maxiter):
# Completion of each stage removes an objective. If both stages # Completion of each stage removes an objective. If both stages
# are complete, then no objectives are left # are complete, then no objectives are left
if not self.objectives: if not self.objectives:
self.col_titles = self.generate_col_titles(
self.n_vars, self.n_slack, self.n_art_vars
)
# Find the values of each variable at optimal solution # Find the values of each variable at optimal solution
return self.interpret_tableau(self.tableau, self.col_titles) return self.interpret_tableau()
row_idx, col_idx = self.find_pivot(self.tableau) row_idx, col_idx = self.find_pivot()
# If there are no more negative values in objective row # If there are no more negative values in objective row
if self.stop_iter: if self.stop_iter:
# Delete artificial variable columns and rows. Update attributes # Delete artificial variable columns and rows. Update attributes
self.tableau = self.change_stage(self.tableau) self.tableau = self.change_stage()
else: else:
self.tableau = self.pivot(self.tableau, row_idx, col_idx) self.tableau = self.pivot(row_idx, col_idx)
return {} return {}
def interpret_tableau( def interpret_tableau(self) -> dict[str, float]:
self, tableau: np.ndarray, col_titles: list[str]
) -> dict[str, float]:
"""Given the final tableau, add the corresponding values of the basic """Given the final tableau, add the corresponding values of the basic
decision variables to the `output_dict` decision variables to the `output_dict`
>>> tableau = np.array([ >>> Tableau(np.array([
... [0,0,0.875,0.375,5], ... [0,0,0.875,0.375,5],
... [0,1,0.375,-0.125,1], ... [0,1,0.375,-0.125,1],
... [1,0,-0.125,0.375,1] ... [1,0,-0.125,0.375,1]
... ]) ... ]),2, 0).interpret_tableau()
>>> t = Tableau(tableau, 2)
>>> t.interpret_tableau(tableau, ["x1", "x2", "s1", "s2", "RHS"])
{'P': 5.0, 'x1': 1.0, 'x2': 1.0} {'P': 5.0, 'x1': 1.0, 'x2': 1.0}
""" """
# P = RHS of final tableau # P = RHS of final tableau
output_dict = {"P": abs(tableau[0, -1])} output_dict = {"P": abs(self.tableau[0, -1])}
for i in range(self.n_vars): for i in range(self.n_vars):
# Gives ids of nonzero entries in the ith column # Gives indices of nonzero entries in the ith column
nonzero = np.nonzero(tableau[:, i]) nonzero = np.nonzero(self.tableau[:, i])
n_nonzero = len(nonzero[0]) n_nonzero = len(nonzero[0])
# First entry in the nonzero ids # First entry in the nonzero indices
nonzero_rowidx = nonzero[0][0] nonzero_rowidx = nonzero[0][0]
nonzero_val = tableau[nonzero_rowidx, i] nonzero_val = self.tableau[nonzero_rowidx, i]
# If there is only one nonzero value in column, which is one # If there is only one nonzero value in column, which is one
if n_nonzero == nonzero_val == 1: if n_nonzero == 1 and nonzero_val == 1:
rhs_val = tableau[nonzero_rowidx, -1] rhs_val = self.tableau[nonzero_rowidx, -1]
output_dict[col_titles[i]] = rhs_val output_dict[self.col_titles[i]] = rhs_val
# Check for basic variables
for title in col_titles:
# Don't add RHS or slack variables to output dict
if title[0] not in "R-s-a":
output_dict.setdefault(title, 0)
return output_dict return output_dict