Initial commit (uploading existing code)

Author: SolidityD

.gitignore

@@ -0,0 +1 @@
+__pychache__/

cuttingplanes.py

@@ -0,0 +1,90 @@
+import simplex
+import numpy as np
+import math
+
+EPSILON = 0.0001
+
+def solve(c, b, A, debug_output=False, blands_rule=False):
+    n = len(c)
+    m = len(b)
+
+    simplex_tableau = simplex.make_tableau(c, b, A)
+
+    return solve_tableau(n, m, simplex_tableau, debug_output, blands_rule)
+
+def solve_tableau(n, m, simplex_tableau, debug_output=False, blands_rule=False):
+    while True:    
+        if debug_output:
+            print("Current simplex + cutting planes tableau :")
+            print(simplex_tableau)
+
+        # Solving the relaxed linear program
+        solved_simplex_tableau = simplex.solve_tableau(n, m, simplex_tableau.copy(), False, blands_rule)
+
+        if debug_output:
+            print("Solved simplex + cutting planes tableau :")
+            print(solved_simplex_tableau)
+
+        # Finding a non integer basic variable
+        basic_variables = simplex.get_basic_variables(solved_simplex_tableau)
+        basic_variable_row = -1
+        basic_variable_column = -1
+        for row, column, value in basic_variables:
+            # This is how we test if a value is close to an integer (we have to because of floating point precision errors)
+            if not -EPSILON < value - round(value) < EPSILON:
+                basic_variable_row = row
+                basic_variable_column = column
+                break
+        else:
+            # We found a solution
+            if debug_output:
+                print("Optimal cutting planes + simplex tableau found")
+            simplex_tableau = solved_simplex_tableau
+            break
+
+        # Creating a linear constraint (used inequality : x_i + sum_j(floor(a_i,j)x_j) <= floor(b_i))
+        constraint_row = np.zeros(n + m)
+        for c in range(n + m):
+            if c != basic_variable_column:
+                constraint_row[c] = math.floor(solved_simplex_tableau[basic_variable_row, c])
+        constraint_row[basic_variable_column] = 1.0
+        constraint_b_value = math.floor(solved_simplex_tableau[basic_variable_row, -1])
+
+        simplex_tableau = simplex.add_constraint(simplex_tableau, constraint_row, constraint_b_value)
+        m += 1
+    
+    return simplex_tableau
+
+# Test case : max x, -x <= 0
+# Output should be : Unbounded tableau exception
+
+# Test case : max x, x <= 1.5
+# Output should be x = 1
+
+# Test case : max x + y, x + y/3 <= 3, x/3 + y <= 3
+# Output should be x = 2, y = 2
+
+# Test case : max y, -x + y <= 1, 3x + 2y <= 12, 2x + 3y <= 12
+# Output should be x = 1, y = 2 or x = 2, y = 2
+    
+if __name__ == '__main__':
+    print("Cutting planes + simplex algorithm")
+    print("Maximize <c, x> subject to constraints Ax <= b, x_i >= 0, x_i integer, where A is an m*n matrix")
+    print("==================")
+
+    n = int(input("Size n (size of vector x) : "))
+    m = int(input("Size m (amount of constraints) : "))
+
+    print("c vector :")
+    c = np.array([float(input("c_" + str(i) + " : ")) for i in range(n)])
+    print("b vector :")
+    b = np.array([float(input("b_" + str(i) + " : ")) for i in range(m)])
+    print("A matrix :")
+    A = []
+    for i in range(m):
+        A.append([float(input("A_" + str(i) + "," + str(j) + " : ")) for j in range(n)])
+    A = np.array(A)
+
+    simplex_tableau = solve(c, b, A, True)
+
+    exit(0)

simplex.py

@@ -0,0 +1,128 @@
+import numpy as np
+
+def make_tableau(c, b, A):
+    n = len(c)
+    m = len(b)
+
+    simplex_tableau = np.hstack((A, np.identity(m), np.zeros((m, 1)), np.array([b]).T))
+    simplex_tableau = np.vstack((simplex_tableau, np.array([np.hstack((-c, np.zeros(m), np.ones(1), np.zeros(1)))])))
+
+    return simplex_tableau
+
+def solve(c, b, A, debug_output=False, blands_rule=False):
+    n = len(c)
+    m = len(b)
+
+    simplex_tableau = make_tableau(c, b, A)
+
+    return solve_tableau(n, m, simplex_tableau, debug_output, blands_rule)
+
+def solve_tableau(n, m, simplex_tableau, debug_output=False, blands_rule=False):
+    if debug_output:
+        print("Simplex tableau :")
+        print(simplex_tableau)
+
+    while True:
+        # Column selection
+        # Equivalent to searching for the variable in the base that would increase the output function the most
+        optimal = True
+        column = 0
+        if blands_rule:
+            for i, e in enumerate(simplex_tableau[-1]):
+                if e < 0:
+                    optimal = False
+                    column = i
+                    break
+        else:
+            for i, e in enumerate(simplex_tableau[-1]):
+                if e < 0:
+                    optimal = False
+                    if e < simplex_tableau[-1, column]:
+                        column = i
+        
+        if optimal:
+            if debug_output:
+                print("Optimal simplex tableau found")
+            break
+
+        # Row selection
+        # Equivalent to searching for the variable outside the base that constraints the augmentation of the output function the most (smallest b_i / leaving variable coefficient)
+        row = -1
+        for i in range(m):
+            # Only positive entries are considered
+            if simplex_tableau[i, column] <= 0:
+                continue
+            if row == -1 or simplex_tableau[row, -1] / simplex_tableau[row, column] > simplex_tableau[i, -1] / simplex_tableau[i, column]:
+                row = i
+        if row == -1:
+            raise Exception(f"Unbounded solutions in tableau (selected column : {column})", simplex_tableau)
+    
+        pivot = simplex_tableau[row, column]
+        simplex_tableau[row] /= pivot
+    
+        for r in range(len(simplex_tableau)):
+            if r != row:
+                simplex_tableau[r] -= simplex_tableau[row] * simplex_tableau[r, column] 
+    
+        if debug_output:
+            print("Current simplex tableau :")
+            print(simplex_tableau)
+
+    return simplex_tableau
+
+# Returns an array of basic variables (row, column, value)
+def get_basic_variables(simplex_tableau):
+    x = []
+
+    # We are excluding the last two columns as they do not contain variables
+    for i in range(len(simplex_tableau.T) - 2):
+        nonzeros = 0
+        nonones = 0
+        oneidx = -1
+        for cidx, c in enumerate(simplex_tableau.T[i]):
+            if c != 0:
+                nonzeros += 1
+            if c != 1:
+                nonones += 1
+            if c == 1:
+                oneidx = cidx
+    
+        if nonzeros == 1 and nonones == len(simplex_tableau.T[i]) - 1:
+            # This is a basic variable
+            x.append((oneidx, i, simplex_tableau[oneidx, -1]))
+    
+    return x
+
+def add_constraint(simplex_tableau, constraint_row, constraint_b_value):
+    rows, columns = simplex_tableau.shape
+
+    constraint_row = np.concatenate((constraint_row, [1.0, 0.0, constraint_b_value]))
+
+    simplex_tableau = np.insert(simplex_tableau, -2, np.zeros(rows), axis=1)
+    simplex_tableau = np.insert(simplex_tableau, -1, constraint_row, axis=0)
+
+    slack_variable = np.zeros(rows)
+    slack_variable[-2] = 1
+    slack_variable = np.array([slack_variable])
+
+    return simplex_tableau
+
+if __name__ == '__main__':
+    print("Simplex algorithm")
+    print("Maximize <c, x> subject to constraints Ax <= b, x_i >= 0, where A is an m*n matrix")
+    print("==================")
+
+    n = int(input("Size n (size of vector x) : "))
+    m = int(input("Size m (amount of constraints) : "))
+
+    print("c vector :")
+    c = np.array([float(input("c_" + str(i) + " : ")) for i in range(n)])
+    print("b vector :")
+    b = np.array([float(input("b_" + str(i) + " : ")) for i in range(m)])
+    print("A matrix :")
+    A = []
+    for i in range(m):
+        A.append([float(input("A_" + str(i) + "," + str(j) + " : ")) for j in range(n)])
+    A = np.array(A)
+
+    solve(c, b, A, True)