Merge pull request scipy#4762 from Andreea-G/linprog_artificial

BUGS: Fixes scipy#4746 and scipy#4594: linprog returns solution violating equality constraint and IndexError when a callback is provided

Author: argriffing

scipy/optimize/_linprog.py

@@ -305,7 +305,6 @@ def _solve_simplex(T, n, basis, maxiter=1000, phase=2, callback=None,
     """
     nit = nit0
     complete = False
-    solution = np.zeros(T.shape[1]-1, dtype=np.float64)
 
     if phase == 1:
         m = T.shape[0]-2
@@ -314,6 +313,37 @@ def _solve_simplex(T, n, basis, maxiter=1000, phase=2, callback=None,
     else:
         raise ValueError("Argument 'phase' to _solve_simplex must be 1 or 2")
 
+    if phase == 2:
+        # Check if any artificial variables are still in the basis.
+        # If yes, check if any coefficients from this row and a column
+        # corresponding to one of the non-artificial variable is non-zero.
+        # If found, pivot at this term. If not, start phase 2.
+        # Do this for all artificial variables in the basis.
+        # Ref: "An Introduction to Linear Programming and Game Theory"
+        # by Paul R. Thie, Gerard E. Keough, 3rd Ed,
+        # Chapter 3.7 Redundant Systems (pag 102)
+        for pivrow in [row for row in range(basis.size)
+                       if basis[row] > T.shape[1] - 2]:
+            non_zero_row = [col for col in range(T.shape[1] - 1)
+                            if T[pivrow, col] != 0]
+            if len(non_zero_row) > 0:
+                pivcol = non_zero_row[0]
+                # variable represented by pivcol enters
+                # variable in basis[pivrow] leaves
+                basis[pivrow] = pivcol
+                pivval = T[pivrow][pivcol]
+                T[pivrow, :] = T[pivrow, :] / pivval
+                for irow in range(T.shape[0]):
+                    if irow != pivrow:
+                        T[irow, :] = T[irow, :] - T[pivrow, :]*T[irow, pivcol]
+                nit += 1
+
+    if len(basis[:m]) == 0:
+        solution = np.zeros(T.shape[1] - 1, dtype=np.float64)
+    else:
+        solution = np.zeros(max(T.shape[1] - 1, max(basis[:m]) + 1),
+                            dtype=np.float64)
+
     while not complete:
         # Find the pivot column
         pivcol_found, pivcol = _pivot_col(T, tol, bland)

scipy/optimize/tests/test_linprog.py

@@ -264,7 +264,13 @@ def test_network_flow_limited_capacity():
             [0, p, p, 0, n],
             [0, 0, 0, p, p]]
     b_eq = [-4, 0, 0, 4]
-    res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
+    # Including the callback here ensures the solution can be
+    # calculated correctly, even when phase 1 terminated
+    # with some of the artificial variables as pivots
+    # (i.e. basis[:m] contains elements corresponding to
+    # the artificial variables)
+    res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
+                  callback=lambda x, **kwargs: None)
     _assert_success(res, desired_fun=14)
 
 
@@ -301,14 +307,16 @@ def test_enzo_example():
 def test_enzo_example_b():
     # rescued from https://github.com/scipy/scipy/pull/218
     c = [2.8, 6.3, 10.8, -2.8, -6.3, -10.8]
-    A_eq = [
-            [-1, -1, -1, 0, 0, 0],
+    A_eq = [[-1, -1, -1, 0, 0, 0],
             [0, 0, 0, 1, 1, 1],
             [1, 0, 0, 1, 0, 0],
             [0, 1, 0, 0, 1, 0],
             [0, 0, 1, 0, 0, 1]]
     b_eq = [-0.5, 0.4, 0.3, 0.3, 0.3]
-    res = linprog(c=c, A_eq=A_eq, b_eq=b_eq)
+    # Including the callback here ensures the solution can be
+    # calculated correctly.
+    res = linprog(c=c, A_eq=A_eq, b_eq=b_eq,
+                  callback=lambda x, **kwargs: None)
     _assert_success(res, desired_fun=-1.77,
                     desired_x=[0.3, 0.2, 0.0, 0.0, 0.1, 0.3])
 
@@ -426,5 +434,21 @@ def test_invalid_inputs():
     assert_raises(ValueError, linprog, [1,2], A_ub=np.zeros((1,1,3)), b_eq=1)
 
 
+def test_basic_artificial_vars():
+    # Test if linprog succeeds when at the end of Phase 1 some artificial
+    # variables remain basic, and the row in T corresponding to the
+    # artificial variables is not all zero.
+    c = np.array([-0.1, -0.07, 0.004, 0.004, 0.004, 0.004])
+    A_ub = np.array([[1.0, 0, 0, 0, 0, 0], [-1.0, 0, 0, 0, 0, 0],
+                     [0, -1.0, 0, 0, 0, 0], [0, 1.0, 0, 0, 0, 0],
+                     [1.0, 1.0, 0, 0, 0, 0]])
+    b_ub = np.array([3.0, 3.0, 3.0, 3.0, 20.0])
+    A_eq = np.array([[1.0, 0, -1, 1, -1, 1], [0, -1.0, -1, 1, -1, 1]])
+    b_eq = np.array([0, 0])
+    res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
+                  callback=lambda x, **kwargs: None)
+    _assert_success(res, desired_fun=0, desired_x=np.zeros_like(c))
+
+
 if __name__ == '__main__':
     run_module_suite()