From 63b80dbbbb8db651b70f010e68211d2ef3b3101a Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 12 Mar 2024 14:36:11 +0900 Subject: [PATCH 1/2] RFC: linprog_simplex: Separate out Phase 1 --- quantecon/optimize/linprog_simplex.py | 35 +++++++++++++++++++++------ 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/quantecon/optimize/linprog_simplex.py b/quantecon/optimize/linprog_simplex.py index 098854437..b10f44170 100644 --- a/quantecon/optimize/linprog_simplex.py +++ b/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,32 @@ 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`. + + """ + success, status, num_iter_1 = \ + solve_tableau(tableau, basis, max_iter, skip_aux=False, + piv_options=piv_options) + + if success and tableau[-1, -1] > piv_options.fea_tol: # Infeasible + success = False + status = 2 + + return success, status, num_iter_1 + + @jit(nopython=True, cache=True) def _pivot_col(tableau, skip_aux, piv_options): """ From 85ca2a1ad20847c968384fbfcc8b93e2a771a7f6 Mon Sep 17 00:00:00 2001 From: Daisuke Oyama Date: Tue, 12 Mar 2024 16:25:39 +0900 Subject: [PATCH 2/2] FIX: linprog_simplex: Fix bug in Phase 1 Fix #725 --- quantecon/optimize/linprog_simplex.py | 20 ++++++++++++++++++- .../optimize/tests/test_linprog_simplex.py | 12 +++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/quantecon/optimize/linprog_simplex.py b/quantecon/optimize/linprog_simplex.py index b10f44170..73b8f7371 100644 --- a/quantecon/optimize/linprog_simplex.py +++ b/quantecon/optimize/linprog_simplex.py @@ -452,13 +452,31 @@ def solve_phase_1(tableau, basis, max_iter=10**6, piv_options=PivOptions()): 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 success and tableau[-1, -1] > piv_options.fea_tol: # Infeasible + 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 diff --git a/quantecon/optimize/tests/test_linprog_simplex.py b/quantecon/optimize/tests/test_linprog_simplex.py index 6f17ed645..50bc92576 100644 --- a/quantecon/optimize/tests/test_linprog_simplex.py +++ b/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)