FIX: linprog_simplex: Fix bug in Phase 1

quantecon/optimize/linprog_simplex.py

@@ -163,14 +163,9 @@ def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
 
     # Phase 1
     success, status, num_iter_1 = \
-        solve_tableau(tableau, basis, max_iter, skip_aux=False,
-                      piv_options=piv_options)
+        solve_phase_1(tableau, basis, max_iter, piv_options=piv_options)
     num_iter += num_iter_1
-    if not success:  # max_iter exceeded
-        return SimplexResult(x, lambd, fun, success, status, num_iter)
-    if tableau[-1, -1] > piv_options.fea_tol:  # Infeasible
-        success = False
-        status = 2
+    if not success:
         return SimplexResult(x, lambd, fun, success, status, num_iter)
 
     # Modify the criterion row for Phase 2
@@ -442,6 +437,50 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True,
     return success, status, num_iter
 
 
+@jit(nopython=True, cache=True)
+def solve_phase_1(tableau, basis, max_iter=10**6, piv_options=PivOptions()):
+    """
+    Perform the simplex algorithm for Phase 1 on a given tableau in
+    canonical form, by calling `solve_tableau` with `skip_aux=False`.
+
+    Parameters
+    ----------
+    See `solve_tableau`.
+
+    Returns
+    -------
+    See `solve_tableau`.
+
+    """
+    L = tableau.shape[0] - 1
+    nm = tableau.shape[1] - (L+1)  # n + m
+
+    success, status, num_iter_1 = \
+        solve_tableau(tableau, basis, max_iter, skip_aux=False,
+                      piv_options=piv_options)
+
+    if not success:  # max_iter exceeded
+        return success, status, num_iter_1
+    if tableau[-1, -1] > piv_options.fea_tol:  # Infeasible
+        success = False
+        status = 2
+        return success, status, num_iter_1
+
+    # Check artifial variables have been eliminated
+    tol_piv = piv_options.tol_piv
+    for i in range(L):
+        if basis[i] >= nm:  # Artifial variable not eliminated
+            for j in range(nm):
+                if tableau[i, j] < -tol_piv or \
+                   tableau[i, j] > tol_piv:  # Treated nonzero
+                    _pivoting(tableau, j, i)
+                    basis[i] = j
+                    num_iter_1 += 1
+                    break
+
+    return success, status, num_iter_1
+
+
 @jit(nopython=True, cache=True)
 def _pivot_col(tableau, skip_aux, piv_options):
     """

quantecon/optimize/tests/test_linprog_simplex.py

@@ -239,3 +239,15 @@ def test_bug_8662(self):
         desired_fun = -36.0000000000
         res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
         _assert_success(res, c, b_ub=b_ub, desired_fun=desired_fun)
+
+
+class TestLinprogSimplex:
+    def test_phase_1_bug_725(self):
+        # Identified a bug in Phase 1
+        # https://github.com/QuantEcon/QuantEcon.py/issues/725
+        c = np.array([-4.09555556, 4.59044444])
+        A_ub = np.array([[1, 0.1], [-1, -0.1], [1, 1]])
+        b_ub = np.array([9.1, -0.1, 0.1])
+        desired_x = [0.1, 0]
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_success(res, c, b_ub=b_ub, desired_x=desired_x)