import numpy as np
np.set_printoptions(suppress=True)
alpha = 0.99995
eps = 0.00001
beta = 0.1
def make_mat(x, sigma, a):
n, m = len(a), len(a[0])
mat = np.array(a)
at = np.transpose(mat)
a1 = np.diagflat(-sigma / x)
return np.block([[a1, at], [mat, np.zeros((n, n))]])
def interior_point_solve(a, b, c):
n, m = len(a), len(a[0])
a_t = np.transpose(a)
x = np.ones(m, dtype=float)
pi = np.ones(n, dtype=float)
sigma = np.ones(m, dtype=float)
mu = beta * sum(x * sigma) / m
while max(abs(x[j] * sigma[j] - mu) for j in range(m)) > eps:
mu = beta * sum(x * sigma) / m
fst = c - np.transpose(a) @ pi - mu / x
snd = b - a @ x
B = np.concatenate((fst, snd))
A = make_mat(x, sigma, a)
x0 = np.linalg.solve(A, B)
delta_x, delta_pi = x0[:m], x0[m:]
delta_sigma = - np.array(sigma) - (sigma * delta_x - mu) / x #[(sigma[i] * delta_x[i] - mu) / x[i] for i in range(m)]
theta_x = min((x[j] / -delta_x[j] for j in range(m) if delta_x[j] < 0), default=1)
phi_s = min((sigma[i] / -delta_sigma[i] for i in range(m) if delta_sigma[i] < 0), default=1)
theta = min(1, alpha * theta_x)
phi = min(1, alpha * phi_s)
# print("total cost :", sum(c * x))
x += theta * delta_x
sigma += phi * delta_sigma
pi += phi * delta_pi
return round(sum(c * x))
def make_array(suppliers, consumers):
s, c = len(suppliers), len(consumers)
n, m = s + c, s * c
dct = {}
for i in range(s):
d = {}
for k in range(i * c, (i + 1) * c):
d[k] = 1
d[m] = suppliers[i]
dct[i] = d
for i in range(c):
d = {}
for k in range(i, s * c, c):
d[k] = 1
d[m] = consumers[i]
dct[s + i] = d
return (dct, n, m + 1)
def sparse_gauss(arr, n, m):
h, k = 0, 0
while h < n and k < m:
idx = next((i for i in range(h, n) if arr[i].get(k, 0) != 0), None)
if idx is None:
k += 1
else:
cpy = dict(arr[idx])
arr[idx] = arr[h]
arr[h] = dict(cpy)
for i in range(h + 1, n):
f = arr[i].get(k, 0) // arr[h][k]
arr[i][k] = 0
for j in range(k + 1, m):
if j in arr[h]:
arr[i][j] = arr[i].get(j, 0) - arr[h][j] * f
h += 1
k += 1
ans = {}
for i in range(n - 1):
ans[i] = {k: v for k, v in arr[i].items() if v != 0 or k == m - 1}
return ans
def solve_eqs(suppliers, consumers):
arr, n, m = make_array(suppliers, consumers)
red_arr = sparse_gauss(arr, n, m)
frees, constrained = set(), set()
keys = list(reversed(red_arr.keys()))
sol = {}
for key in keys:
d = red_arr[key]
ks = sorted([k for k in d.keys() if k != m - 1])
if len(ks) > 1:
nk = ks.pop(0)
dct = {m - 1: d[m - 1] // d[nk]}
constrained.add(nk)
for k in ks:
if k != m - 1 and not k in constrained:
frees.add(k)
if not k in sol:
sol[k] = {k: 1, m - 1: 0}
for other in ks:
for k, v in sol[other].items():
dct[k] = dct.get(k, 0) - d.get(other, 0) * (v // d[nk])
sol[nk] = dct
return sol
def minimum_transportation_price(suppliers, consumers, _costs):
costs = sum(_costs,[])
sol = solve_eqs(suppliers, consumers)
print("Eqns before adding slack variables :", sol)
cost = {}
for i, c in enumerate(costs):
for k, v in sol[i].items():
cost[k] = cost.get(k, 0) + c * v
INIT_VAR = len(sol)
VMAX = cost[INIT_VAR]
NUM_VAR = INIT_VAR
for k, d in sol.items():
if not k in d:
NUM_VAR += 1
# print(sol)
# print(INIT_VAR, NUM_VAR)
_cost_ = np.zeros((NUM_VAR,))
for k, v in cost.items():
if k != INIT_VAR:
_cost_[k] = v
a, b = [], []
slack = 0
for k in sorted(sol.keys()):
d = sol[k]
row = [0 for _ in range(NUM_VAR)]
if not k in d:
sgn = 1 if d[INIT_VAR] > 0 else -1
for key, value in d.items():
if key == INIT_VAR:
b.append(sgn * value)
else:
row[key] = -sgn * value
row[INIT_VAR + slack] = sgn
slack += 1
a.append(row)
def to_dense(dct, n, m):
return [[dct[i].get(j, 0) if i in dct else 0 for j in range(m)] for i in range(n)]
def make_sparse_params(sol, init_var, num_var):
# getting a_dct
a_dct = {}
b_dct = {}
slack_dct = 0
for k in sorted(sol.keys()):
d = sol[k]
row_dct = {}
if not k in d:
# skip condition xi >= 0
sgn = 1 if d[init_var] > 0 else -1
for key, value in d.items():
# for this value of key, this is the second member of the condition
if key == init_var:
b_dct[slack_dct] = sgn * value
else:
row_dct[key] = row_dct.get(key, 0) - sgn * value
row_dct[init_var + slack_dct] = sgn
a_dct[slack_dct] = {k: v for k, v in row_dct.items() if v != 0}
slack_dct += 1
print("a_dct", a_dct)
print("b_dct", b_dct)
print("height", slack_dct)
print("width", num_var)
# transpose of a
a_dct_tr = {}
for i in a_dct:
for j, v in a_dct[i].items():
if j in a_dct_tr:
a_dct_tr[j][i] = a_dct[i][j]
else:
a_dct_tr[j] = {i: a_dct[i][j]}
print("transpose ", a_dct_tr)
A_dct = {}
for i in a_dct:
A_dct[num_var + i] = dict(a_dct[i])
for i in a_dct_tr:
for j, v in a_dct_tr[i].items():
if i in A_dct:
A_dct[i][j + num_var] = a_dct_tr[i][j]
else:
A_dct[i] = {j + num_var: a_dct_tr[i][j]}
size = slack_dct + num_var
# print(A_dct)
print(np.array(to_dense(A_dct, size, size)))
return A_dct, b_dct
make_sparse_params(sol, INIT_VAR, NUM_VAR)
A_eq = np.array(a)
print("A_eq :\n", A_eq, sep='')
# print(len(A_eq), len(A_eq[0]))
B_eq = np.array(b)
print("B_eq :\n", B_eq, sep='')
# print(_cost_)
r = interior_point_solve(A_eq, B_eq, _cost_)
return r + VMAX
tests = (
([10, 7, 13], [6, 20, 4], [[4, 12, 3], [20, 1, 6], [7, 0, 5]], 43),
([8, 15, 21], [8, 36], [[9, 16], [7, 13], [25, 1]], 288),
([31, 16], [14, 17, 16], [[41, 18, 0], [4, 16, 37]], 358),
(
[10, 20, 20],
[5, 25, 10, 10],
[[2, 5, 3, 0], [3, 4, 1, 4], [2, 6, 5, 2]],
150,
),
(
[13, 44, 27, 39, 17],
[28, 12, 30, 17, 19, 34],
[
[6, 6, 12, 8, 13, 13],
[7, 20, 5, 16, 11, 16],
[4, 6, 19, 0, 2, 18],
[1, 16, 6, 11, 8, 11],
[5, 6, 11, 1, 6, 14],
],
759,
),
)
tests = (
([10, 7, 13], [6, 20, 4], [[4, 12, 3], [20, 1, 6], [7, 0, 5]], 43),
([8, 15, 21], [8, 36], [[9, 16], [7, 13], [25, 1]], 288),
([31, 16], [14, 17, 16], [[41, 18, 0], [4, 16, 37]], 358),
(
[10, 20, 20],
[5, 25, 10, 10],
[[2, 5, 3, 0], [3, 4, 1, 4], [2, 6, 5, 2]],
150,
),
(
[13, 44, 27, 39, 17],
[28, 12, 30, 17, 19, 34],
[
[6, 6, 12, 8, 13, 13],
[7, 20, 5, 16, 11, 16],
[4, 6, 19, 0, 2, 18],
[1, 16, 6, 11, 8, 11],
[5, 6, 11, 1, 6, 14],
],
759,
),
)
for suppliers, consumers, costs, res in tests:
sol = minimum_transportation_price(suppliers, consumers, costs)
print(sol, res)
To embed this program on your website, copy the following code and paste it into your website's HTML: