simplex and IPA choice

Author: Arthod

__pycache__/linear_programming.cpython-39.pyc

@@ -24,8 +24,8 @@ def interior_point_iteration(A, c, x_initial, alpha):
     x = D @ x_tilde
     return x
 
-def interior_point(A, c, x_initial, alpha=0.5, optimal_tol=10e-5):
-    iterations = 0
+def interior_point(A, c, x_initial, alpha=0.5, optimal_tol=10e-6):
+    iterations_count = 0
     path = [x_initial]
 
     x_prev = -10 * np.ones(c.size)
@@ -35,8 +35,8 @@ def interior_point(A, c, x_initial, alpha=0.5, optimal_tol=10e-5):
         x = interior_point_iteration(A, c, x, alpha)
 
         path.append(x)
-        iterations += 1
+        iterations_count += 1
 
-    return x, path, iterations
+    return x, path, iterations_count
 
 

interior_point.py

@@ -6,7 +6,7 @@
 import matplotlib.pyplot as plt
 
 import interior_point
-from simplex import simplex
+import simplex
 
 
 GREATER_EQUAL = -1
@@ -105,6 +105,10 @@ def copy(self) -> "LinearProgrammingProblem":
         lp.is_maximizing = self.is_maximizing
         lp.objective = self.objective.copy()
 
+        lp._A_updated = False
+        lp._b_updated = False
+        lp._c_updated = False
+
         return lp
 
 
@@ -220,22 +224,54 @@ def compute_feasible_basis(self):
     def is_feasible(lp: "LinearProgrammingProblem", x: np.array, tol=10e-5) -> bool:
         return all(np.isclose(lp.A @ x, lp.b, atol=tol))
 
-    def solve(self, method=SOLVER_SIMPLEX):
-        if (method == SOLVER_SIMPLEX):
-            lp_slacked = self.slacken_problem()
+    def solve(self, method=None):
+        lp_slacked = self.slacken_problem()
+        # Decide which solver to use
+        # Simplex?
+        if (method == None):
+            x_zero = simplex.zero_point_solution(lp_slacked.A, lp_slacked.b, lp_slacked.c)
+            if (simplex.is_feasible(lp_slacked.A, x_zero, lp_slacked.b)):
+                # If zero point solution is feasible, use simplex
+                method = SOLVER_SIMPLEX
+
+        if (method == None):
+            # Solve using interior point algorithm
+            x_ones = np.ones(len(self.variables))
+            x_initial = np.hstack((x_ones, self.b - (self.A @ x_ones)))
+            if (self.is_feasible(self, x_initial)):
+                # IPA is feasible
+                method = SOLVER_INTERIOR_POINT_METHOD
 
-            x_sol = simplex(lp_slacked.A, lp_slacked.b, lp_slacked.c, lp_slacked.vars_slack_amount)
+
+        x_sol = None
+        iterations_count = None
+
+        if (method == SOLVER_SIMPLEX):
+            # Solve using simplex
+            x_sol, iterations_count = simplex.simplex(lp_slacked.A, lp_slacked.b, lp_slacked.c, lp_slacked.vars_slack_amount)
 
             # Cut slack variables
             x_sol = x_sol[:lp_slacked.c.size - lp_slacked.vars_slack_amount]
-            print(x_sol)
 
-            return x_sol
+        elif (method == SOLVER_INTERIOR_POINT_METHOD):
+            # Solve using interior point algorithm
+            x_ones = np.ones(len(self.variables))
+            x_initial = np.hstack((x_ones, self.b - (self.A @ x_ones)))
+
+            #x_initial = np.array([1 for _ in range(len(lp_slacked.variables))])
+            x_sol, path, iterations_count = interior_point.interior_point(lp_slacked.A, lp_slacked.c, x_initial)
+
+            # Cut slack variables
+            x_sol = x_sol[:lp_slacked.c.size - lp_slacked.vars_slack_amount]
+
+        else:
+            print("No solution found")
+            return
 
-        if (method == SOLVER_INTERIOR_POINT_METHOD):
-            pass
+        print(f"Solving using {method}")
+        print(f"Solved in {iterations_count} iterations")
 
-        return "x"
+        return x_sol
 
 
     def plot_solution_path(self):

linear_programming.py

@@ -12,14 +12,17 @@
 """
 
 
-def simplex(A: np.array, b: np.array, c: np.array, slacks_amt: int) -> np.array:
+def simplex(A: np.array, b: np.array, c: np.array, slacks_amt: int) -> "tuple[np.array, int]":
     # https://sdu.itslearning.com/ContentArea/ContentArea.aspx?LocationID=17727&LocationType=1&ElementID=589350
     #x_current = basic_feasible_solution(A, b, c)
     #N = [i for i in range(c.size - slacks_amt)]   # Not In basis
     B = [i + (c.size - slacks_amt) for i in range(slacks_amt)]    # In basis
+    # Zero basic solution
+
+    iterations_count = 0
 
     while True:
-        tableau(A, b, c)
+        iterations_count += 1
         ## Choosing a pivot
         # Pick a non-negative column    TODO
         positives = np.where(c > 0)
@@ -73,19 +76,19 @@ def simplex(A: np.array, b: np.array, c: np.array, slacks_amt: int) -> np.array:
     x_sol = np.zeros(c.size)
     for i in range(len(B)):
         x_sol[B[i]] = b[i]
-    return np.array(x_sol)
+    return np.array(x_sol), iterations_count
 
 def row_mult(A: np.array, row_index: int, scalar: int) -> np.array:
-    # Row operation: R1 = scalar * R1
-    # Inspiration: (TODO - optimize)
+    # Row operation: R1 = scalar * R1 (TODO - optimize)
+    # Inspiration:
     # https://personal.math.ubc.ca/~pwalls/math-python/linear-algebra/solving-linear-systems/
     I = np.eye(A.shape[0])
     I[row_index, row_index] = scalar
     return I @ A
 
 def row_add(A: np.array, row_main_index: int, row_other_index: int, scalar: int) -> np.array:
-    # Row operation: R1 = R1 + scalar * R2
-    # Inspiration: (TODO - optimize)
+    # Row operation: R1 = R1 + scalar * R2 (TODO - optimize)
+    # Inspiration:
     # https://personal.math.ubc.ca/~pwalls/math-python/linear-algebra/solving-linear-systems/
 
     I = np.eye(A.shape[0])
@@ -109,8 +112,7 @@ def zero_point_solution(A: np.array, b: np.array, c: np.array) -> np.array:
     slack_sign = [A[i][num_decision+i] for i in range(len(b))]
     slack_values = [b[i]/slack_sign[i] for i in range(b.size)]
     zero_point = np.hstack((np.zeros(num_decision), slack_values))
-    print("zero_point:", zero_point)
-    print("Zero point solution is feasible:", is_feasible(A, zero_point, b))
+    
     return zero_point