#! /usr/bin/env python3
import math
import sys
# ! /usr/bin/env python3
import numpy as np
class Simplex:
def __init__(self, A, b, c, B_idx, N_idx=None):
self.A = np.matrix(A)
self.b = np.matrix(b)
self.c = np.matrix(c)
self.B_idx = B_idx
if N_idx is None:
self.N_idx = [i for i in range(self.A.shape[1]) if i not in self.B_idx]
else:
self.N_idx = N_idx
def getB(self):
return self.A[:, self.B_idx]
def getCb(self):
return self.c[:, self.B_idx]
def getCn(self):
return self.c[:, self.N_idx]
def getN(self):
return self.A[:, self.N_idx]
def calcZvar(self):
InvB = np.linalg.inv(self.getB())
N = self.getN()
Cb = self.getCb()
Cb_InvB = Cb @ InvB
Cb_InvB_b = Cb_InvB @ self.b.T
Cb_InvB_N = Cb_InvB @ N
bracadddinget_term = Cb_InvB_N - self.getCn()
return Cb_InvB_b, -1 * bracadddinget_term
def getVector(self):
X = self.computeX()
res = [0.0] * self.c.shape[1]
index = 0
for j in self.B_idx:
res[j] = X[index, 0]
index += 1
return res
def isOptimal(self, z):
return all([x <= 0 for x in np.nditer(z)])
def enteringVariable(self, z):
# pri('CHOOSE ENTERING')nt
max_ascent = -1 * np.inf
index = -1
it = np.nditer(z.T, flags=['f_index'])
while not it.finished:
# print(it[0])
if it[0] > max_ascent:
max_ascent = it[0]
index = it.index
it.iternext()
if index != -1:
InvB = np.linalg.inv(self.getB())
val = InvB @ self.A[:, self.N_idx[index]]
return index, val
else:
return None, None
def computeX(self):
InvB = np.linalg.inv(self.getB())
InvB_b = InvB @ self.b.T
return InvB_b
def exitingVariable(self, entering_pos, ev):
X = self.computeX()
# print('CHOOSE EXITING')
min_value = np.inf
index = -1
for idx, cr in enumerate(np.nditer([X, ev])):
x, y = cr[0], cr[1]
# print(x, y)
if y <= 0:
continue
ratio = x / y
if ratio < min_value:
min_value = ratio
index = idx
# print('ratio={}'.format(min_value))
if index != -1:
return index
else:
return None, None
def do(self, verbose=False):
while True:
z, z_var = self.calcZvar()
# print(z, z_var)
if self.isOptimal(z_var):
if verbose:
print("solution")
return z, self.getVector(), self.B_idx
entering_pos, val = self.enteringVariable(z_var)
if entering_pos is None:
# print('Entering None')
# print(self.getVector())
return None, None, None
exiting_pos = self.exitingVariable(entering_pos, val)
if exiting_pos is None:
# print('Exiting None')
# print(self.getVector())
return None, None, None
self.B_idx[exiting_pos] = self.N_idx[entering_pos]
self.N_idx = [i for i in range(self.A.shape[1]) if i not in self.B_idx]
def calc_simplex_max(orig_A, orig_b, orig_c):
A = []
for l in orig_A:
A.append(l.copy())
b = orig_b.copy()
c = [0.0] * len(A[0])
base = [i + len(A[0]) for i in range(len(b))]
for i, row in enumerate(A):
for j in range(len(b)):
if i == j:
row.append(1.0)
else:
row.append(0.0)
for _ in range(len(b)):
c.append(-1.0)
# print(A)
# print(b)
# print(c)
# print(base)
simplex = Simplex(A, b, c, base)
z, res, base = simplex.do()
if z is None:
return None, None
if any([x >= len(orig_c) for x in base]):
return None, None
A = []
for l in orig_A:
A.append(l.copy())
b = orig_b.copy()
c = orig_c.copy()
# c = simplex.c
# notbase = [x for x in notbase if x < len(orig_c)]
# print(A)
# print(b)
# print(c)
# print(base)
# print(notbase)
simplex = Simplex(A, b, c, base)
z, res, _ = simplex.do()
return z, res
def isAlmostEqual(x, y, epsilon=10 ** (-8)):
return abs(x - y) <= epsilon
def isAlmostEqualLists(x, y, epsilon=10 ** (-8)):
return all(isAlmostEqual(x_i, y_i, epsilon) for (x_i, y_i) in zip(x, y))
def isIntList(list):
rounded_list = [round(x) for x in list]
return isAlmostEqualLists(list, rounded_list)
def createBranch(x, index, A, b):
A1 = [condition + [0.] for condition in A]
new_A1 = [0.] * len(A1[0])
new_A1[index] = 1.
new_A1[-1] = 1.
A1.append(new_A1)
b1 = b.copy()
b1.append(math.floor(x[index]))
# print('Created: {} {}'.format(A1, b1))
A2 = [condition + [0.] for condition in A]
new_A2 = [0.] * len(A1[0])
# new_A2[index] = -1.
new_A2[index] = 1.
new_A2[-1] = -1.
A2.append(new_A2)
b2 = b.copy()
b2.append(math.floor(x[index]) + 1.)
# print('Created: {} {}'.format(A2, b2))
return A1, b1, A2, b2
def process(A, b, c, q, xs):
# print('')
# print(A)
# print(b)
# print(c)
try:
f_val, x_val = calc_simplex_max(A, b, c)
except IndexError:
print(len(A[0]), len(b), len(c))
raise IndexError
# print(f_val, x_val)
if f_val is None:
return
if isIntList(x_val):
q.append((f_val, x_val, A, b))
xs.append((f_val, x_val))
else:
q.append((f_val, x_val, A, b))
def branchsAndBounds(A, b, c, iterations_max=30):
k = 0
f_val, x_val = calc_simplex_max(A, b, c)
if f_val is None:
return None, None, None
q = [(f_val, x_val.copy(), A, b)]
xs = []
if isIntList(x_val):
return f_val, x_val, k
k = 0
while q:
# print("Q len = {}".format(len(q)))
x_values = [x[1] for x in q]
f_values = [x[0] for x in q]
max_f_value_index = f_values.index(max(f_values))
_, new_x, c_A, c_b = q.pop(max_f_value_index)
frac_parts = [i for i, x in enumerate(new_x) if divmod(x, 1)[1] > 0]
if not frac_parts:
continue
index = frac_parts[0]
A1, b1, A2, b2 = createBranch(new_x, index, c_A, c_b)
process(A1, b1, c + [0.0] * (len(A1[0]) - len(A[0])), q, xs)
process(A2, b2, c + [0.0] * (len(A2[0]) - len(A[0])), q, xs)
if k > iterations_max:
break
k += 1
# sys.stdin.readline()
f_values = [x[0] for x in xs]
print(f_values)
max_val = max(f_values)
max_x = xs[f_values.index(max(f_values))][1]
return max_val, max_x, k
def main():
# A = [[1, 3, 4, -1], [2, 5, 1, 1]]
# b = [23, 14]
# c = [4, 2, -1, 0]
A = [[1.0, 3.0, 4.0, 0.0], [2.0, 5.0, 7.0, 1.0]]
b = [28.0, 64.0]
c = [3.0, 4.0, -10.0, 0.0]
print(branchsAndBounds(A, b, c))
if __name__ == '__main__':
main()