operations_research_60016/simplex.py at main · ohmdelta/operations_research_60016 · 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
import numpy as np


def swapIn(tableu: np.ndarray, indIn: int ,indOut: int):

    basic = tableu[1:,0]
    basic[indOut] = indIn

    extMatrix = tableu[:,1:]
    extMatrix[indOut + 1] /= extMatrix[indOut + 1, indIn]

    for i in range(len(tableu)):
        if i == indOut + 1:
            continue
        else:
            extMatrix[i] -= extMatrix[indOut + 1] * extMatrix[i, indIn]


def simplex(tableu: np.ndarray, zeroIndexed=True, printIntermediateTableu=False, blandsRule=False, err: float = 1e-9) -> np.ndarray:

    z = tableu[0, 1:-1]
    basic = tableu[1:, 0]

    cost = tableu[:, -1]

    matrix = tableu[1:, 1:-1]

    if printIntermediateTableu: print(np.c_[tableu[:,0] + int(zeroIndexed),tableu[:,1:]])

    if np.all(z <= 0):
        # overallCost = tableu[0,-1]
        return basic + int(zeroIndexed), cost[1:] , cost[0]
    else:
        if not blandsRule:
            indIn = np.argsort(z[:-1])[::-1]
            for ind in indIn:
                arr = [(cost[k+1] / i if i > 0 else np.inf) for k, i in enumerate(matrix[:, ind])]
                if z[:-1][ind] <= err:
                    return basic + int(zeroIndexed), cost[1:] , cost[0]
                if np.all(np.array(arr) == np.inf):
                    continue
                    # return basic + int(zeroIndexed), cost[1:] , cost[0]
                else:
                    break
            if z[:-1][ind] <= err:
                return basic + int(zeroIndexed), cost[1:] , cost[0]

            indOut = np.argmin(arr)
            indIn = ind

            swapIn(tableu, indIn, indOut)

            return simplex(tableu, zeroIndexed, printIntermediateTableu)

        else:
            indIn = range(len(z[:-1]))
            for ind in indIn:
                arr = [(cost[k+1] / i if i > 0 else np.inf) for k, i in enumerate(matrix[:, ind])]
                if z[:-1][ind] <= err or ind in tableu[1:,0]:
                    continue
                if np.all(np.array(arr) == np.inf):
                    continue
                    # return basic + int(zeroIndexed), cost[1:] , cost[0]
                else:
                    break

            if z[:-1][ind] <= err:
                return basic + int(zeroIndexed), cost[1:] , cost[0]

            indOut = np.argmin(arr)
            indIn = ind
            swapIn(tableu, indIn, indOut)

            return simplex(tableu, zeroIndexed, printIntermediateTableu, blandsRule)


np.set_printoptions(precision=4, suppress=True, linewidth=10000)

if __name__ == "__main__":
    pass