added doctests

Author: shoresj

linear_programming/interior_point_method.py

@@ -1,11 +1,3 @@
-"""
-Python implementation of the Primal-Dual Interior-Point Method for solving linear
-programs with - `>=`, `<=`, and `=` constraints and - each variable `x1, x2, ... >= 0`.
-
-Resources:
-https://en.wikipedia.org/wiki/Interior-point_method
-"""
-
 import numpy as np
 
 
@@ -15,11 +7,14 @@ class InteriorPointMethod:
 
     Attributes:
         objective_coefficients (np.ndarray): Coefficient matrix for the objective
-function.
+            function.
         constraint_matrix (np.ndarray): Constraint matrix.
         constraint_bounds (np.ndarray): Constraint bounds.
         tol (float): Tolerance for stopping criterion.
         max_iter (int): Maximum number of iterations.
+
+    Methods:
+        solve: Solve the linear programming problem.
     """
 
     def __init__(
@@ -40,18 +35,53 @@ def __init__(
             raise ValueError("Invalid input for the linear programming problem.")
 
     def _is_valid_input(self) -> bool:
-        """Validate the input for the linear programming problem."""
+        """
+        Validate the input for the linear programming problem.
+
+        Returns:
+            bool: True if input is valid, False otherwise.
+
+        >>> objective_coefficients = np.array([1, 2])
+        >>> constraint_matrix = np.array([[1, 1], [1, -1]])
+        >>> constraint_bounds = np.array([2, 0])
+        >>> ipm = InteriorPointMethod(objective_coefficients, constraint_matrix,
+constraint_bounds)
+        >>> ipm._is_valid_input()
+        True
+        >>> constraint_bounds = np.array([2, 0, 1])
+        >>> ipm = InteriorPointMethod(objective_coefficients, constraint_matrix,
+constraint_bounds)
+        >>> ipm._is_valid_input()
+        False
+        """
         return (
             self.constraint_matrix.shape[0] == self.constraint_bounds.shape[0]
             and self.constraint_matrix.shape[1] == self.objective_coefficients.shape[0]
         )
 
     def _convert_to_standard_form(self) -> tuple[np.ndarray, np.ndarray]:
-        """Convert constraints to standard form by adding slack and surplus
-variables."""
+        """
+        Convert constraints to standard form by adding slack variables.
+
+        Returns:
+            tuple: A tuple of the standard form constraint matrix and objective
+coefficients.
+
+        >>> objective_coefficients = np.array([1, 2])
+        >>> constraint_matrix = np.array([[1, 1], [1, -1]])
+        >>> constraint_bounds = np.array([2, 0])
+        >>> ipm = InteriorPointMethod(objective_coefficients, constraint_matrix,
+constraint_bounds)
+        >>> a_standard, c_standard = ipm._convert_to_standard_form()
+        >>> a_standard
+        array([[ 1.,  1.,  1.,  0.],
+               [ 1., -1.,  0.,  1.]])
+        >>> c_standard
+        array([1., 2., 0., 0.])
+        """
         (m, n) = self.constraint_matrix.shape
-        slack_surplus = np.eye(m)
-        a_standard = np.hstack([self.constraint_matrix, slack_surplus])
+        slack = np.eye(m)
+        a_standard = np.hstack([self.constraint_matrix, slack])
         c_standard = np.hstack([self.objective_coefficients, np.zeros(m)])
         return a_standard, c_standard
 
@@ -60,15 +90,24 @@ def solve(self) -> tuple[np.ndarray, float]:
         Solve problem with Primal-Dual Interior-Point Method.
 
         Returns:
-        tuple: A tuple with optimal soln and the optimal value.
+            tuple: A tuple with optimal solution and the optimal value.
+
+        >>> objective_coefficients = np.array([1, 2])
+        >>> constraint_matrix = np.array([[1, 1], [1, -1]])
+        >>> constraint_bounds = np.array([2, 0])
+        >>> ipm = InteriorPointMethod(objective_coefficients, constraint_matrix,
+constraint_bounds)
+        >>> solution, value = ipm.solve()
+        >>> np.isclose(value, np.dot(objective_coefficients, solution))
+        True
         """
         a, c = self._convert_to_standard_form()
         m, n = a.shape
         x = np.ones(n)
         s = np.ones(n)
         y = np.ones(m)
 
-        for _ in range(self.max_iter):
+        for iteration in range(self.max_iter):
             x_diag = np.diag(x)
             s_diag = np.diag(s)
 
@@ -99,23 +138,35 @@ def solve(self) -> tuple[np.ndarray, float]:
             r = np.hstack([-r2, -r1, -r3 + mu * np.ones(n)])
 
             # Solve the KKT system
-            delta = np.linalg.solve(kkt, r)
+            try:
+                delta = np.linalg.solve(kkt, r)
+            except np.linalg.LinAlgError:
+                print("KKT matrix is singular, switching to least squares solution")
+                delta = np.linalg.lstsq(kkt, r, rcond=None)[0]
 
             dx = delta[:n]
             dy = delta[n : n + m]
             ds = delta[n + m :]
 
             # Step size
-            alpha_primal = min(1, 0.99 * min(-x[dx < 0] / dx[dx < 0], default=1))
-            alpha_dual = min(1, 0.99 * min(-s[ds < 0] / ds[ds < 0], default=1))
+            alpha_primal = min(
+                1, 0.99 * min([-x[i] / dx[i] for i in range(n) if dx[i] < 0], default=1)
+            )
+            alpha_dual = min(
+                1, 0.99 * min([-s[i] / ds[i] for i in range(n) if ds[i] < 0], default=1)
+            )
 
             # Update variables
             x += alpha_primal * dx
             y += alpha_dual * dy
             s += alpha_dual * ds
 
-        optimal_value = np.dot(c, x)
-        return x, optimal_value
+            print(f"Iteration {iteration}: x = {x}, s = {s}, y = {y}")
+
+        # Extract the solution (remove slack variables)
+        original_vars = x[: self.objective_coefficients.shape[0]]
+        optimal_value = np.dot(self.objective_coefficients, original_vars)
+        return original_vars, optimal_value
 
 
 if __name__ == "__main__":
@@ -127,5 +178,10 @@ def solve(self) -> tuple[np.ndarray, float]:
         objective_coefficients, constraint_matrix, constraint_bounds
     )
     solution, value = ipm.solve()
-    print("Optimal solution:", solution)
-    print("Optimal value:", value)
+    print("Final optimal solution:", solution)
+    print("Final optimal value:", value)
+
+    # Verify the solution
+    print("\nVerification:")
+    print("Objective value calculation matches final optimal value:")
+    print(np.isclose(value, np.dot(objective_coefficients, solution)))