from functools import reduce
def show(arr):
return '\n'.join(' | '.join(str(x) for x in row) for row in arr)
class Triplet:
def __init__(self, x, y, z=0):
self.x = x
self.y = y
self.z = z
def __copy__(self):
new_object = type(self)(self.x, self.y, self.z)
return new_object
def __str__(self):
return f"{self.x} {'+' if self.y >= 0 else '-'} {abs(self.y)}·ε"
def __neg__(self, other):
return Triplet(-self.x, -self.y)
def is_pos(self):
return self > Triplet(0, 0)
def __add__(self, other):
return Triplet(self.x + other.x, self.y + other.y, self.z + other.z)
def __sub__(self, other):
return Triplet(self.x - other.x, self.y - other.y, self.z + other.z)
def __lt__(self, other):
return (self.x, self.y, self.z) < (other.x, other.y, other.z)
def __le__(self, other):
return (self.x, self.y, self.z) <= (other.x, other.y, other.z)
def __gt__(self, other):
return (self.x, self.y, self.z) > (other.x, other.y, other.z)
def __ge__(self, other):
return (self.x, self.y, self.z) > (other.x, other.y, other.z)
def __mul__(self, k):
assert isinstance(k, int), "Multiplication undefined"
return Triplet(k * self.x, k * self.y, k * self.z)
def __rmul__(self, k):
assert isinstance(k, int), "Multiplication undefined"
return Triplet(k * self.x, k * self.y, k * self.z)
@classmethod
def Triplet_sum(cls, Triplets):
return reduce(lambda acc, x: acc + x, Triplets, Triplet(0, 0))
class Simplex:
def __init__(self, a, b, c):
self.m, self.n = len(a), len(b)
self.costs = c
self.a = [Triplet(x, 1) for i, x in enumerate(a)]
self.b = [Triplet(x, self.m) if i == self.n - 1 else Triplet(x, 0) for i, x in enumerate(b)]
def make_initial_solution(self):
arr = [[Triplet(0, 0) for _ in range(self.n + 1)] for _ in range(self.m + 1)]
for i in range(self.m):
arr[i][self.n] = self.a[i]
for j in range(self.n):
arr[self.m][j] = self.b[j]
i, j = 0, 0
num_pos = 0
inv_pos = {}
pos, row, col = [], {i: set() for i in range(self.m)}, {j: set() for j in range(self.n)},
while i < self.m or j < self.n:
sum_row = arr[i][self.n] - Triplet.Triplet_sum([arr[i][k] for k in range(j)]) if i < self.m else Triplet(0, 0)
sum_col = arr[self.m][j] - Triplet.Triplet_sum([arr[k][j] for k in range(i)]) if j < self.n else Triplet(0, 0)
if not sum_row.is_pos():
i += 1
elif not sum_col.is_pos():
j += 1
else:
if sum_row < sum_col:
arr[i][j] = sum_row
pos.append((i, j))
inv_pos[(i, j)] = num_pos
num_pos += 1
row[i].add((i, j))
col[j].add((i, j))
i += 1
else:
arr[i][j] = sum_col
pos.append((i, j))
inv_pos[(i, j)] = num_pos
num_pos += 1
row[i].add((i, j))
col[j].add((i, j))
j += 1
arr[self.m][self.n] = Triplet.Triplet_sum([arr[self.m][k] for k in range(self.n)])
print("Init :")
print(pos)
print('\n\n')
return arr, pos, inv_pos, row, col, float("inf")
def step(self, arr, pos, inv_pos, row, col, theta):
x, y = pos[0]
u, v, visited = [None for _ in range(self.m)], [None for _ in range(self.n)], set([(x, y)])
lam, mu = {i: {} for i in range(self.m)}, {i: {} for i in range(self.n)}
lam[x][0] = 1
u[x] = self.costs[x][y]
v[y] = 0
queue = [(x, y)]
# print("start : " ,x,y)
while queue:
i, j = queue.pop()
for a, b in row[i]:
if not (a, b) in visited:
# print("v before:", v, b, j)
if v[j] is not None:
v[b] = v[j] + self.costs[a][b] - self.costs[i][j]
i1, i2 = inv_pos[(i, j)], inv_pos[(a, b)]
d = dict(mu[j])
d[i2] = d.get(i2, 0) + 1
d[i1] = d.get(i1, 0) - 1
mu[b] = d
# elif v[j] is None:
# v[j] = v[b] - (self.costs[a][b] - self.costs[i][j])
# i1, i2 = inv_pos[(i, j)], inv_pos[(a, b)]
# d = dict(mu[j])
# d[i2] = d.get(i2, 0) - 1
# d[i1] = d.get(i1, 0) + 1
# mu[b] = d
# print("v after:", v)
visited.add((a, b))
queue.append((a, b))
for a, b in col[j]:
if not (a, b) in visited:
# print("u before :", u, i, a)
if u[i] is not None:
# print("new visited (col):", a, b)
u[a] = u[i] + self.costs[a][b] - self.costs[i][j]
i1, i2 = inv_pos[(i, j)], inv_pos[(a, b)]
# print("i1, i2 :", i1, i2)
d = dict(lam[i])
d[i2] = d.get(i2, 0) + 1
d[i1] = d.get(i1, 0) - 1
lam[a] = d
# elif u[i] is None:
# # print("new visited (col):", a, b)
# u[i] = u[a] - (self.costs[a][b] - self.costs[i][j])
# i1, i2 = inv_pos[(i, j)], inv_pos[(a, b)]
# # print("i1, i2 :", i1, i2)
# d = dict(lam[a])
# d[i2] = d.get(i2, 0) - 1
# d[i1] = d.get(i1, 0) + 1
# lam[i] = d
# print("u after:", u)
visited.add((a, b))
queue.append((a, b))
indirect_cost = [[u[i] + v[j] for j in range(self.n)] for i in range(self.m)]
delta = [[u[i] + v[j] - self.costs[i][j] for j in range(self.n)] for i in range(self.m)]
print("u :", u)
print("v :", v)
print("delta :", delta)
print("indirect cost", indirect_cost)
max_delta = 0
for i in range(self.m):
for j in range(self.n):
delta = u[i] + v[j] - self.costs[i][j]
if delta > max_delta:
max_delta = delta
i0, j0 = i, j
print("lambda : ", lam)
print("mu : ", mu)
if max_delta == 0:
return (arr, pos, inv_pos, row, col, max_delta)
print("max delta :", max_delta)
print("Position of max delta :", i0, j0)
##############
theta = Triplet(float('inf'), float('inf'), float('inf'))
theta_idx = set()
print("nu :", [lam[i0].get(k, 0) + mu[j0].get(k, 0) for k in range(self.n + self.m - 1)])
x_rem, y_rem = 0, 0
for k in range(self.n + self.m - 1):
nu = lam[i0].get(k, 0) + mu[j0].get(k, 0)
x, y = pos[k]
if nu == 1:
if arr[x][y] < theta:
theta = arr[x][y].__copy__()
x_rem, y_rem = x, y
theta_idx.add((1, (x, y)))
elif nu == -1:
theta_idx.add((-1, (x, y)))
print("theta idx:", theta_idx)
print("Before remove", show(arr), sep='\n')
print("Position to remove :", x_rem, y_rem)
print("theta :", theta)
for v, (x, y) in theta_idx:
arr[x][y] -= v * theta
arr[i0][j0] += theta
idx = inv_pos[(x_rem, y_rem)]
pos[idx] = (i0, j0)
del inv_pos[(x_rem, y_rem)]
inv_pos[(i0, j0)] = idx
row[x_rem].remove((x_rem, y_rem))
col[y_rem].remove((x_rem, y_rem))
row[i0].add((i0, j0))
col[j0].add((i0, j0))
# print("theta :", theta)
print("pos :", pos)
print("inv pos :", inv_pos)
print("row :", row)
print("col :", col)
print("After remove", show(arr), sep='\n')
return arr, pos, inv_pos, row, col, max_delta
def solve(self):
arr, pos, inv_pos, row, col, max_delta = self.make_initial_solution()
while max_delta > 0:
arr, pos, inv_pos, row, col, max_delta = self.step(arr, pos, inv_pos, row, col, max_delta)
print('\n\n')
ans = 0
for i in range(self.m):
for j in range(self.n):
ans += self.costs[i][j] * arr[i][j].x
return ans
a = [1, 5, 7]
b = [3, 3, 3, 2, 2]
c = [[3, 2, 1, 2, 3], [5, 4, 3, -1, 1], [0, 2, 3, 4, 5]]
a = [31, 16]
b = [14, 17, 16]
c = [
[41, 18, 0],
[4, 16, 37]
]
simp = Simplex(a, b, c)
r = simp.solve()
print(r)
To embed this program on your website, copy the following code and paste it into your website's HTML: