Python/linear_programming/simplex.py at 6d20b1fc2129c01d0a03142fdcc7cd24c9d23a83 · MeriemBchini/Python · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
from typing import Dict
import numpy as np


class Tableau:
    """
    Simplex algorithm implementation using tableau form.

    Solves linear programming problems of the form:
    maximize c^T x subject to Ax <= b and x >= 0
    """

    maxiter = 100

    def __init__(
        self, tableau: np.ndarray, n_vars: int, n_artificial_vars: int
    ) -> None:
        tableau = tableau.astype("float64")

        if tableau.dtype != "float64":
            raise TypeError("Tableau must have type float64")

        if not (tableau[:, -1] >= 0).all():
            raise ValueError("RHS must be > 0")

        if n_vars < 2 or n_artificial_vars < 0:
            raise ValueError("Invalid number of variables")

        self.tableau = tableau
        self.n_rows, n_cols = tableau.shape

        self.n_vars = n_vars
        self.n_artificial_vars = n_artificial_vars

        self.n_stages = (self.n_artificial_vars > 0) + 1
        self.n_slack = n_cols - self.n_vars - self.n_artificial_vars - 1

        self.objectives = ["max"]
        if self.n_artificial_vars:
            self.objectives.append("min")

        self.col_titles = self.generate_col_titles()

        self.stop_iter = False

    def generate_col_titles(self) -> list[str]:
        string_starts = ["x", "s"]
        sizes = [self.n_vars, self.n_slack]

        titles: list[str] = []

        for i in range(2):
            for j in range(sizes[i]):
                titles.append(string_starts[i] + str(j + 1))

        titles.append("RHS")
        return titles

    def find_pivot(self) -> tuple[int, int]:
        objective = self.objectives[-1]

        sign = (objective == "min") - (objective == "max")
        col_idx = int(np.argmax(sign * self.tableau[0, :-1]))

        if sign * self.tableau[0, col_idx] <= 0:
            self.stop_iter = True
            return 0, 0

        s = slice(self.n_stages, self.n_rows)

        dividend = self.tableau[s, -1]
        divisor = self.tableau[s, col_idx]

        nans = np.full(self.n_rows - self.n_stages, np.nan)

        quotients = np.divide(dividend, divisor, out=nans, where=divisor > 0)

        row_idx = int(np.nanargmin(quotients) + self.n_stages)
        return row_idx, col_idx

    def pivot(self, row_idx: int, col_idx: int) -> np.ndarray:
        tableau = self.tableau

        piv_row = tableau[row_idx].copy()
        piv_row /= piv_row[col_idx]

        tableau -= tableau[:, col_idx][:, None] * piv_row
        tableau[row_idx] = piv_row

        self.tableau = tableau
        return tableau

    def change_stage(self) -> np.ndarray:
        self.objectives.pop()

        if not self.objectives:
            return self.tableau

        s = slice(-self.n_artificial_vars - 1, -1)

        self.tableau = np.delete(self.tableau, s, axis=1)
        self.tableau = np.delete(self.tableau, 0, axis=0)

        self.n_stages = 1
        self.n_rows -= 1
        self.n_artificial_vars = 0
        self.stop_iter = False

        return self.tableau

    def run_simplex(self) -> Dict[str, float]:
        """
        Run simplex algorithm until optimal solution is found.

        >>> t = Tableau(np.array([
        ... [-1, -1, 0, 0, 0],
        ... [1, 3, 1, 0, 4],
        ... [3, 1, 0, 1, 4]
        ... ], dtype="float64"), 2, 0)
        >>> result = t.run_simplex()
        >>> result["P"]
        2.0
        """

        for iteration in range(Tableau.maxiter):
            if not self.objectives:
                return self.interpret_tableau()

            row_idx, col_idx = self.find_pivot()

            if self.stop_iter:
                self.tableau = self.change_stage()
            else:
                self.tableau = self.pivot(row_idx, col_idx)

        raise ValueError(
            "Simplex did not converge.\n"
            f"- Iterations: {Tableau.maxiter}\n"
            f"- Remaining objectives: {self.objectives}"
        )

    def interpret_tableau(self) -> Dict[str, float]:
        output: Dict[str, float] = {"P": float(abs(self.tableau[0, -1]))}

        for i in range(self.n_vars):
            nonzero = np.nonzero(self.tableau[:, i])[0]

            if len(nonzero) == 1:
                row = nonzero[0]
                if self.tableau[row, i] == 1:
                    output[self.col_titles[i]] = float(self.tableau[row, -1])

        return output


if __name__ == "__main__":
    import doctest

    doctest.testmod()