Skip to content

Commit

Permalink
Merge pull request #727 from QuantEcon/fix-linprog
Browse files Browse the repository at this point in the history
FIX: linprog_simplex: Fix bug in Phase 1
  • Loading branch information
oyamad authored Mar 14, 2024
2 parents d7c2f52 + 85ca2a1 commit 0adb3ec
Show file tree
Hide file tree
Showing 2 changed files with 58 additions and 7 deletions.
53 changes: 46 additions & 7 deletions quantecon/optimize/linprog_simplex.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
"""
Expand Down
12 changes: 12 additions & 0 deletions quantecon/optimize/tests/test_linprog_simplex.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

0 comments on commit 0adb3ec

Please sign in to comment.