From 1f1ad4be5ff69ed9aa0b46174ed9631abf5541cf Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sun, 1 Jul 2018 15:38:56 +0900
Subject: [PATCH 1/6] MAINT: Move pivoting routines to `optimize.pivoting`

---
 docs/source/optimize.rst              |   1 +
 docs/source/optimize/pivoting.rst     |   7 ++
 quantecon/game_theory/lemke_howson.py | 150 +-----------------------
 quantecon/optimize/pivoting.py        | 158 ++++++++++++++++++++++++++
 4 files changed, 169 insertions(+), 147 deletions(-)
 create mode 100644 docs/source/optimize/pivoting.rst
 create mode 100644 quantecon/optimize/pivoting.py

diff --git a/docs/source/optimize.rst b/docs/source/optimize.rst
index 564fc8ab6..049b722fe 100644
--- a/docs/source/optimize.rst
+++ b/docs/source/optimize.rst
@@ -6,4 +6,5 @@ Optimize
 
    optimize/nelder_mead
    optimize/root_finding
+   optimize/pivoting
    optimize/scalar_maximization
diff --git a/docs/source/optimize/pivoting.rst b/docs/source/optimize/pivoting.rst
new file mode 100644
index 000000000..65e5a5e15
--- /dev/null
+++ b/docs/source/optimize/pivoting.rst
@@ -0,0 +1,7 @@
+pivoting
+========
+
+.. automodule:: quantecon.optimize.pivoting
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/quantecon/game_theory/lemke_howson.py b/quantecon/game_theory/lemke_howson.py
index f8a05a33d..1f01d23e6 100644
--- a/quantecon/game_theory/lemke_howson.py
+++ b/quantecon/game_theory/lemke_howson.py
@@ -7,10 +7,7 @@
 import numpy as np
 from numba import jit
 from .utilities import NashResult
-
-
-TOL_PIV = 1e-10
-TOL_RATIO_DIFF = 1e-15
+from ..optimize.pivoting import _pivoting, _lex_min_ratio_test
 
 
 def lemke_howson(g, init_pivot=0, max_iter=10**6, capping=None,
@@ -398,8 +395,8 @@ def _lemke_howson_tbl(tableaux, bases, init_pivot, max_iter):
     while True:
         for pl in pls:
             # Determine the leaving variable
-            row_min = _lex_min_ratio_test(tableaux[pl], pivot,
-                                          slack_starts[pl], argmins)
+            _, row_min = _lex_min_ratio_test(tableaux[pl], pivot,
+                                             slack_starts[pl], argmins)
 
             # Pivoting step: modify tableau in place
             _pivoting(tableaux[pl], pivot, row_min)
@@ -421,147 +418,6 @@ def _lemke_howson_tbl(tableaux, bases, init_pivot, max_iter):
     return converged, num_iter
 
 
-@jit(nopython=True, cache=True)
-def _pivoting(tableau, pivot, pivot_row):
-    """
-    Perform a pivoting step. Modify `tableau` in place.
-
-    Parameters
-    ----------
-    tableau : ndarray(float, ndim=2)
-        Array containing the tableau.
-
-    pivot : scalar(int)
-        Pivot.
-
-    pivot_row : scalar(int)
-        Pivot row index.
-
-    Returns
-    -------
-    tableau : ndarray(float, ndim=2)
-        View to `tableau`.
-
-    """
-    nrows, ncols = tableau.shape
-
-    pivot_elt = tableau[pivot_row, pivot]
-    for j in range(ncols):
-        tableau[pivot_row, j] /= pivot_elt
-
-    for i in range(nrows):
-        if i == pivot_row:
-            continue
-        multiplier = tableau[i, pivot]
-        if multiplier == 0:
-            continue
-        for j in range(ncols):
-            tableau[i, j] -= tableau[pivot_row, j] * multiplier
-
-    return tableau
-
-
-@jit(nopython=True, cache=True)
-def _min_ratio_test_no_tie_breaking(tableau, pivot, test_col,
-                                    argmins, num_candidates):
-    """
-    Perform the minimum ratio test, without tie breaking, for the
-    candidate rows in `argmins[:num_candidates]`. Return the number
-    `num_argmins` of the rows minimizing the ratio and store thier
-    indices in `argmins[:num_argmins]`.
-
-    Parameters
-    ----------
-    tableau : ndarray(float, ndim=2)
-        Array containing the tableau.
-
-    pivot : scalar(int)
-        Pivot.
-
-    test_col : scalar(int)
-        Index of the column used in the test.
-
-    argmins : ndarray(int, ndim=1)
-        Array containing the indices of the candidate rows. Modified in
-        place to store the indices of minimizing rows.
-
-    num_candidates : scalar(int)
-        Number of candidate rows in `argmins`.
-
-    Returns
-    -------
-    num_argmins : scalar(int)
-        Number of minimizing rows.
-
-    """
-    ratio_min = np.inf
-    num_argmins = 0
-
-    for k in range(num_candidates):
-        i = argmins[k]
-        if tableau[i, pivot] <= TOL_PIV:  # Treated as nonpositive
-            continue
-        ratio = tableau[i, test_col] / tableau[i, pivot]
-        if ratio > ratio_min + TOL_RATIO_DIFF:  # Ratio large for i
-            continue
-        elif ratio < ratio_min - TOL_RATIO_DIFF:  # Ratio smaller for i
-            ratio_min = ratio
-            num_argmins = 1
-        else:  # Ratio equal
-            num_argmins += 1
-        argmins[num_argmins-1] = i
-
-    return num_argmins
-
-
-@jit(nopython=True, cache=True)
-def _lex_min_ratio_test(tableau, pivot, slack_start, argmins):
-    """
-    Perform the lexico-minimum ratio test.
-
-    Parameters
-    ----------
-    tableau : ndarray(float, ndim=2)
-        Array containing the tableau.
-
-    pivot : scalar(int)
-        Pivot.
-
-    slack_start : scalar(int)
-        First index for the slack variables.
-
-    argmins : ndarray(int, ndim=1)
-        Empty array used to store the row indices. Its length must be no
-        smaller than the number of the rows of `tableau`.
-
-    Returns
-    -------
-    row_min : scalar(int)
-        Index of the row with the lexico-minimum ratio.
-
-    """
-    nrows = tableau.shape[0]
-    num_candidates = nrows
-
-    # Initialize `argmins`
-    for i in range(nrows):
-        argmins[i] = i
-
-    num_argmins = _min_ratio_test_no_tie_breaking(tableau, pivot, -1,
-                                                  argmins, num_candidates)
-    if num_argmins == 1:
-        return argmins[0]
-
-    for j in range(slack_start, slack_start+nrows):
-        if j == pivot:
-            continue
-        num_argmins = _min_ratio_test_no_tie_breaking(tableau, pivot, j,
-                                                      argmins, num_argmins)
-        if num_argmins == 1:
-            break
-    return argmins[0]
-
-
 @jit(nopython=True, cache=True)
 def _get_mixed_actions(tableaux, bases):
     """
diff --git a/quantecon/optimize/pivoting.py b/quantecon/optimize/pivoting.py
new file mode 100644
index 000000000..18b312710
--- /dev/null
+++ b/quantecon/optimize/pivoting.py
@@ -0,0 +1,158 @@
+"""
+Contain pivoting routines commonly used in the Simplex Algorithm and
+Lemke-Howson Algorithm routines.
+
+"""
+import numpy as np
+from numba import jit
+
+
+TOL_PIV = 1e-10
+TOL_RATIO_DIFF = 1e-15
+
+
+@jit(nopython=True, cache=True)
+def _pivoting(tableau, pivot, pivot_row):
+    """
+    Perform a pivoting step. Modify `tableau` in place.
+
+    Parameters
+    ----------
+    tableau : ndarray(float, ndim=2)
+        Array containing the tableau.
+
+    pivot : scalar(int)
+        Pivot.
+
+    pivot_row : scalar(int)
+        Pivot row index.
+
+    Returns
+    -------
+    tableau : ndarray(float, ndim=2)
+        View to `tableau`.
+
+    """
+    nrows, ncols = tableau.shape
+
+    pivot_elt = tableau[pivot_row, pivot]
+    for j in range(ncols):
+        tableau[pivot_row, j] /= pivot_elt
+
+    for i in range(nrows):
+        if i == pivot_row:
+            continue
+        multiplier = tableau[i, pivot]
+        if multiplier == 0:
+            continue
+        for j in range(ncols):
+            tableau[i, j] -= tableau[pivot_row, j] * multiplier
+
+    return tableau
+
+
+@jit(nopython=True, cache=True)
+def _min_ratio_test_no_tie_breaking(tableau, pivot, test_col,
+                                    argmins, num_candidates):
+    """
+    Perform the minimum ratio test, without tie breaking, for the
+    candidate rows in `argmins[:num_candidates]`. Return the number
+    `num_argmins` of the rows minimizing the ratio and store thier
+    indices in `argmins[:num_argmins]`.
+
+    Parameters
+    ----------
+    tableau : ndarray(float, ndim=2)
+        Array containing the tableau.
+
+    pivot : scalar(int)
+        Pivot.
+
+    test_col : scalar(int)
+        Index of the column used in the test.
+
+    argmins : ndarray(int, ndim=1)
+        Array containing the indices of the candidate rows. Modified in
+        place to store the indices of minimizing rows.
+
+    num_candidates : scalar(int)
+        Number of candidate rows in `argmins`.
+
+    Returns
+    -------
+    num_argmins : scalar(int)
+        Number of minimizing rows.
+
+    """
+    ratio_min = np.inf
+    num_argmins = 0
+
+    for k in range(num_candidates):
+        i = argmins[k]
+        if tableau[i, pivot] <= TOL_PIV:  # Treated as nonpositive
+            continue
+        ratio = tableau[i, test_col] / tableau[i, pivot]
+        if ratio > ratio_min + TOL_RATIO_DIFF:  # Ratio large for i
+            continue
+        elif ratio < ratio_min - TOL_RATIO_DIFF:  # Ratio smaller for i
+            ratio_min = ratio
+            num_argmins = 1
+        else:  # Ratio equal
+            num_argmins += 1
+        argmins[num_argmins-1] = i
+
+    return num_argmins
+
+
+@jit(nopython=True, cache=True)
+def _lex_min_ratio_test(tableau, pivot, slack_start, argmins):
+    """
+    Perform the lexico-minimum ratio test.
+
+    Parameters
+    ----------
+    tableau : ndarray(float, ndim=2)
+        Array containing the tableau.
+
+    pivot : scalar(int)
+        Pivot.
+
+    slack_start : scalar(int)
+        First index for the slack variables.
+
+    argmins : ndarray(int, ndim=1)
+        Empty array used to store the row indices. Its length must be no
+        smaller than the number of the rows of `tableau`.
+
+    Returns
+    -------
+    found : bool
+        False if there is no positive entry in the pivot column.
+
+    row_min : scalar(int)
+        Index of the row with the lexico-minimum ratio.
+
+    """
+    nrows = tableau.shape[0]
+    num_candidates = nrows
+
+    found = False
+
+    # Initialize `argmins`
+    for i in range(nrows):
+        argmins[i] = i
+
+    num_argmins = _min_ratio_test_no_tie_breaking(tableau, pivot, -1,
+                                                  argmins, num_candidates)
+    if num_argmins == 1:
+        found = True
+    elif num_argmins >= 2:
+        for j in range(slack_start, slack_start+nrows):
+            if j == pivot:
+                continue
+            num_argmins = _min_ratio_test_no_tie_breaking(tableau, pivot, j,
+                                                          argmins, num_argmins)
+            if num_argmins == 1:
+                found = True
+                break
+    return found, argmins[0]

From 85e4c2724fe15010419b844ef5e69abdfdf44156 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Mon, 9 Jul 2018 16:55:57 +0900
Subject: [PATCH 2/6] EHN: Add linprog_simplex

---
 quantecon/optimize/__init__.py                |   2 +-
 quantecon/optimize/linprog_simplex.py         | 238 +++++++++++++++++
 .../optimize/tests/test_linprog_simplex.py    | 251 ++++++++++++++++++
 3 files changed, 490 insertions(+), 1 deletion(-)
 create mode 100644 quantecon/optimize/linprog_simplex.py
 create mode 100644 quantecon/optimize/tests/test_linprog_simplex.py

diff --git a/quantecon/optimize/__init__.py b/quantecon/optimize/__init__.py
index 9b5a2ef53..372aef025 100644
--- a/quantecon/optimize/__init__.py
+++ b/quantecon/optimize/__init__.py
@@ -2,7 +2,7 @@
 """
 Initialization of the optimize subpackage
 """
-
+from .linprog_simplex import linprog_simplex, solve_tableau, get_solution
 from .scalar_maximization import brent_max
 from .nelder_mead import nelder_mead
 from .root_finding import newton, newton_halley, newton_secant, bisect, brentq
diff --git a/quantecon/optimize/linprog_simplex.py b/quantecon/optimize/linprog_simplex.py
new file mode 100644
index 000000000..2c37083ce
--- /dev/null
+++ b/quantecon/optimize/linprog_simplex.py
@@ -0,0 +1,238 @@
+"""
+Contain a linear programming solver routine based on the Simplex Method.
+
+"""
+from collections import namedtuple
+import numpy as np
+from numba import jit
+from .pivoting import _pivoting, _lex_min_ratio_test
+
+
+FEA_TOL = 1e-6
+
+
+SimplexResult = namedtuple(
+    'SimplexResult', ['x', 'lambd', 'fun', 'success', 'status', 'num_iter']
+)
+
+
+@jit(nopython=True, cache=True)
+def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
+                    A_eq=np.empty((0, 0)), b_eq=np.empty((0,)), max_iter=10**6,
+                    tableau=None, basis=None, x=None, lambd=None):
+    n, m, k = c.shape[0], A_ub.shape[0], A_eq.shape[0]
+    L = m + k
+
+    if tableau is None:
+        tableau = np.empty((L+1, n+m+L+1))
+    if basis is None:
+        basis = np.empty(L, dtype=np.int_)
+    if x is None:
+        x = np.empty(n)
+    if lambd is None:
+        lambd = np.empty(L)
+
+    num_iter = 0
+    fun = -np.inf
+
+    b_signs = np.empty(L, dtype=np.bool_)
+    for i in range(m):
+        b_signs[i] = True if b_ub[i] >= 0 else False
+    for i in range(k):
+        b_signs[m+i] = True if b_eq[i] >= 0 else False
+
+    # Construct initial tableau for Phase 1
+    _initialize_tableau(A_ub, b_ub, A_eq, b_eq, tableau, basis)
+
+    # Phase 1
+    success, status, num_iter_1 = \
+        solve_tableau(tableau, basis, max_iter, skip_aux=False)
+    num_iter += num_iter_1
+    if not success:  # max_iter exceeded
+        return SimplexResult(x, lambd, fun, success, status, num_iter)
+    if tableau[-1, -1] > FEA_TOL:  # Infeasible
+        success = False
+        status = 2
+        return SimplexResult(x, lambd, fun, success, status, num_iter)
+
+    # Modify the criterion row for Phase 2
+    _set_criterion_row(c, basis, tableau)
+
+    # Phase 2
+    success, status, num_iter_2 = \
+        solve_tableau(tableau, basis, max_iter-num_iter, skip_aux=True)
+    num_iter += num_iter_2
+    fun = get_solution(tableau, basis, x, lambd, b_signs)
+
+    return SimplexResult(x, lambd, fun, success, status, num_iter)
+
+
+@jit(nopython=True, cache=True)
+def _initialize_tableau(A_ub, b_ub, A_eq, b_eq, tableau, basis):
+    m, k = A_ub.shape[0], A_eq.shape[0]
+    L = m + k
+    n = tableau.shape[1] - (m+L+1)
+
+    for i in range(m):
+        for j in range(n):
+            tableau[i, j] = A_ub[i, j]
+    for i in range(k):
+        for j in range(n):
+            tableau[m+i, j] = A_eq[i, j]
+
+    tableau[:L, n:-1] = 0
+
+    for i in range(m):
+        tableau[i, -1] = b_ub[i]
+        if tableau[i, -1] < 0:
+            for j in range(n):
+                tableau[i, j] *= -1
+            tableau[i, n+i] = -1
+            tableau[i, -1] *= -1
+        else:
+            tableau[i, n+i] = 1
+        tableau[i, n+m+i] = 1
+    for i in range(k):
+        tableau[m+i, -1] = b_eq[i]
+        if tableau[m+i, -1] < 0:
+            for j in range(n):
+                tableau[m+i, j] *= -1
+            tableau[m+i, -1] *= -1
+        tableau[m+i, n+m+m+i] = 1
+
+    tableau[-1, :] = 0
+    for i in range(L):
+        for j in range(n+m):
+            tableau[-1, j] += tableau[i, j]
+        tableau[-1, -1] += tableau[i, -1]
+
+    for i in range(L):
+        basis[i] = n+m+i
+
+    return tableau, basis
+
+
+@jit(nopython=True, cache=True)
+def _set_criterion_row(c, basis, tableau):
+    n = c.shape[0]
+    L = basis.shape[0]
+
+    for j in range(n):
+        tableau[-1, j] = c[j]
+    tableau[-1, n:] = 0
+
+    for i in range(L):
+        multiplier = tableau[-1, basis[i]]
+        for j in range(tableau.shape[1]):
+            tableau[-1, j] -= tableau[i, j] * multiplier
+
+    return tableau
+
+
+@jit(nopython=True, cache=True)
+def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
+    """
+    Perform the simplex algorithm on a given tableau in canonical form.
+
+    Used to solve a linear program in the following form:
+
+        maximize:     c @ x
+
+        subject to:   A_ub @ x <= b_ub
+                      A_eq @ x == b_eq
+                      x >= 0
+
+    where A_ub is of shape (m, n) and A_eq is of shape (k, n). Thus,
+    `tableau` is of shape (L+1, n+m+L+1), where L=m+k, and
+
+    * `tableau[np.arange(L), :][:, basis]` must be an identity matrix,
+      and
+    * the elements of `tableau[:-1, -1]` must be nonnegative.
+
+    Parameters
+    ----------
+    tableau : ndarray(float, ndim=2)
+        ndarray of shape (L+1, n+m+L+1) containing the tableau. Modified
+        in place.
+
+    basis : ndarray(int, ndim=1)
+        ndarray of shape (L,) containing the basic variables. Modified
+        in place.
+
+    max_iter : scalar(int), optional(default=10**6)
+        Maximum number of pivoting steps.
+
+    skip_aux : bool, optional(default=True)
+        Whether to skip the coefficients of the auxiliary (or
+        artificial) variables in pivot column selection.
+
+    """
+    L = tableau.shape[0] - 1
+
+    # Array to store row indices in lex_min_ratio_test
+    argmins = np.empty(L, dtype=np.int_)
+
+    success = False
+    status = 1
+    num_iter = 0
+
+    while num_iter < max_iter:
+        num_iter += 1
+
+        pivcol_found, pivcol = _pivot_col(tableau, skip_aux)
+
+        if not pivcol_found:  # Optimal
+            success = True
+            status = 0
+            break
+
+        aux_start = tableau.shape[1] - L - 1
+        pivrow_found, pivrow = _lex_min_ratio_test(tableau[:-1, :], pivcol,
+                                                   aux_start, argmins)
+
+        if not pivrow_found:  # Unbounded
+            success = False
+            status = 3
+            break
+
+        _pivoting(tableau, pivcol, pivrow)
+        basis[pivrow] = pivcol
+
+    return success, status, num_iter
+
+
+@jit(nopython=True, cache=True)
+def _pivot_col(tableau, skip_aux):
+    L = tableau.shape[0] - 1
+    criterion_row_stop = tableau.shape[1] - 1
+    if skip_aux:
+        criterion_row_stop -= L
+
+    found = False
+    pivcol = -1
+    coeff = FEA_TOL
+    for j in range(criterion_row_stop):
+        if tableau[-1, j] > coeff:
+            coeff = tableau[-1, j]
+            pivcol = j
+            found = True
+
+    return found, pivcol
+
+
+@jit(nopython=True, cache=True)
+def get_solution(tableau, basis, x, lambd, b_signs):
+    n, L = x.size, lambd.size
+    aux_start = tableau.shape[1] - L - 1
+
+    x[:] = 0
+    for i in range(L):
+        if basis[i] < n:
+            x[basis[i]] = tableau[i, -1]
+    for j in range(L):
+        lambd[j] = tableau[-1, aux_start+j]
+        if lambd[j] != 0 and b_signs[j]:
+            lambd[j] *= -1
+    fun = tableau[-1, -1] * (-1)
+
+    return fun
diff --git a/quantecon/optimize/tests/test_linprog_simplex.py b/quantecon/optimize/tests/test_linprog_simplex.py
new file mode 100644
index 000000000..d5006b1f0
--- /dev/null
+++ b/quantecon/optimize/tests/test_linprog_simplex.py
@@ -0,0 +1,251 @@
+"""
+Tests for linprog_simplex
+
+"""
+import numpy as np
+from numpy.testing import assert_, assert_allclose, assert_equal
+
+from quantecon.optimize import linprog_simplex
+
+
+# Assert functions from scipy/optimize/tests/test_linprog.py
+
+def _assert_infeasible(res):
+    # res: linprog result object
+    assert_(not res.success, "incorrectly reported success")
+    assert_equal(res.status, 2, "failed to report infeasible status")
+
+
+def _assert_unbounded(res):
+    # res: linprog result object
+    assert_(not res.success, "incorrectly reported success")
+    assert_equal(res.status, 3, "failed to report unbounded status")
+
+
+def _assert_success(res, c, b_ub=np.array([]), b_eq=np.array([]),
+                    desired_fun=None, desired_x=None, rtol=1e-15, atol=1e-15):
+    if not res.success:
+        msg = "linprog status {0}".format(res.status)
+        raise AssertionError(msg)
+
+    assert_equal(res.status, 0)
+    if desired_fun is not None:
+        assert_allclose(res.fun, desired_fun,
+                        err_msg="converged to an unexpected objective value",
+                        rtol=rtol, atol=atol)
+    if desired_x is not None:
+        assert_allclose(res.x, desired_x,
+                        err_msg="converged to an unexpected solution",
+                        rtol=rtol, atol=atol)
+
+    fun_p = c @ res.x
+    fun_d = np.concatenate([b_ub, b_eq]) @ res.lambd
+    assert_allclose(fun_p, fun_d, rtol=rtol, atol=atol)
+    assert_allclose(res.fun, fun_d, rtol=rtol, atol=atol)
+
+
+class TestLinprogSimplexScipy:
+    # Test cases from scipy/optimize/tests/test_linprog.py
+
+    def test_infeasible(self):
+        # Test linprog response to an infeasible problem
+        c = np.array([-1, -1]) * (-1)
+        A_ub = [[1, 0],
+                [0, 1],
+                [-1, -1]]
+        b_ub = [2, 2, -5]
+        c, A_ub, b_ub = map(np.asarray, [c, A_ub, b_ub])
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_infeasible(res)
+
+    def test_unbounded(self):
+        # Test linprog response to an unbounded problem
+        c = np.array([1, 1])
+        A_ub = [[-1, 1],
+                [-1, -1]]
+        b_ub = [-1, -2]
+        c, A_ub, b_ub = map(np.asarray, [c, A_ub, b_ub])
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_unbounded(res)
+
+    def test_nontrivial_problem(self):
+        # Test linprog for a problem involving all constraint types,
+        # negative resource limits, and rounding issues.
+        c = np.array([-1, 8, 4, -6]) * (-1)
+        A_ub = [[-7, -7, 6, 9],
+                [1, -1, -3, 0],
+                [10, -10, -7, 7],
+                [6, -1, 3, 4]]
+        b_ub = [-3, 6, -6, 6]
+        A_eq = [[-10, 1, 1, -8]]
+        b_eq = [-4]
+        c, A_ub, b_ub, A_eq, b_eq = \
+            map(np.asarray, [c, A_ub, b_ub, A_eq, b_eq])
+        desired_fun = 7083 / 1391 * (-1)
+        desired_x = [101 / 1391, 1462 / 1391, 0, 752 / 1391]
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)
+        _assert_success(res, c, b_ub=b_ub, b_eq=b_eq, desired_fun=desired_fun,
+                        desired_x=desired_x)
+
+    def test_network_flow(self):
+        # A network flow problem with supply and demand at nodes
+        # and with costs along directed edges.
+        # https://www.princeton.edu/~rvdb/542/lectures/lec10.pdf
+        c = np.array([2, 4, 9, 11, 4, 3, 8, 7, 0, 15, 16, 18]) * (-1)
+        n, p = -1, 1
+        A_eq = [
+            [n, n, p, 0, p, 0, 0, 0, 0, p, 0, 0],
+            [p, 0, 0, p, 0, p, 0, 0, 0, 0, 0, 0],
+            [0, 0, n, n, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, p, p, 0, 0, p, 0],
+            [0, 0, 0, 0, n, n, n, 0, p, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, n, n, 0, 0, p],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, n, n, n]]
+        b_eq = [0, 19, -16, 33, 0, 0, -36]
+        c, A_eq, b_eq = map(np.asarray, [c, A_eq, b_eq])
+        desired_fun = -755
+        res = linprog_simplex(c, A_eq=A_eq, b_eq=b_eq)
+        _assert_success(res, c, b_eq=b_eq, desired_fun=desired_fun)
+
+    def test_basic_artificial_vars(self):
+        # 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]) * (-1)
+        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_simplex(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)
+        _assert_success(res, c, b_ub=b_ub, b_eq=b_eq, desired_fun=0,
+                        desired_x=np.zeros_like(c))
+
+    def test_bug_5400(self):
+        # https://github.com/scipy/scipy/issues/5400
+        f = 1 / 9
+        g = -1e4
+        h = -3.1
+        A_ub = np.array([
+            [1, -2.99, 0, 0, -3, 0, 0, 0, -1, -1, 0, -1, -1, 1, 1, 0, 0, 0, 0],
+            [1, 0, -2.9, h, 0, -3, 0, -1, 0, 0, -1, 0, -1, 0, 0, 1, 1, 0, 0],
+            [1, 0, 0, h, 0, 0, -3, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, 1, 1],
+            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
+            [0, 1.99, -1, -1, 0, 0, 0, -1, f, f, 0, 0, 0, g, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, -1, -1, 0, 0, 0, -1, f, f, 0, g, 0, 0, 0, 0],
+            [0, -1, 1.9, 2.1, 0, 0, 0, f, -1, -1, 0, 0, 0, 0, 0, g, 0, 0, 0],
+            [0, 0, 0, 0, -1, 2, -1, 0, 0, 0, f, -1, f, 0, 0, 0, g, 0, 0],
+            [0, -1, -1, 2.1, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, 0, 0, g, 0],
+            [0, 0, 0, 0, -1, -1, 2, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, g]])
+
+        b_ub = np.array([
+            0.0, 0, 0, 100, 100, 100, 100, 100, 100, 900, 900, 900, 900, 900,
+            900, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
+
+        c = np.array([-1.0, 1, 1, 1, 1, 1, 1, 1, 1,
+                      1, 1, 1, 1, 0, 0, 0, 0, 0, 0]) * (-1)
+
+        desired_fun = 106.63507541835018
+
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_success(res, c, b_ub=b_ub, desired_fun=desired_fun,
+                        rtol=1e-8, atol=1e-8)
+
+    def test_issue_8174_stackoverflow(self):
+        # Test supplementary example from issue 8174.
+        # https://github.com/scipy/scipy/issues/8174
+        # https://stackoverflow.com/questions/47717012/
+        c = np.array([1, 0, 0, 0, 0, 0, 0]) * (-1)
+        A_ub = -np.identity(7)
+        b_ub = np.full(7, -2)
+        A_eq = np.array([
+            [1, 1, 1, 1, 1, 1, 0],
+            [0.3, 1.3, 0.9, 0, 0, 0, -1],
+            [0.3, 0, 0, 0, 0, 0, -2/3],
+            [0, 0.65, 0, 0, 0, 0, -1/15],
+            [0, 0, 0.3, 0, 0, 0, -1/15]
+        ])
+        b_eq = np.array([100, 0, 0, 0, 0])
+
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq)
+        _assert_success(res, c, b_ub=b_ub, b_eq=b_eq)
+
+    def test_linprog_cyclic_recovery(self):
+        # Test linprogs recovery from cycling using the Klee-Minty problem
+        # Klee-Minty  https://www.math.ubc.ca/~israel/m340/kleemin3.pdf
+        c = np.array([100, 10, 1])
+        A_ub = [[1, 0, 0],
+                [20, 1, 0],
+                [200, 20, 1]]
+        b_ub = [1, 100, 10000]
+        c, A_ub, b_ub = map(np.asarray, [c, A_ub, b_ub])
+        desired_x = [0, 0, 10000]
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_success(res, c, b_ub=b_ub, desired_x=desired_x)
+
+    def test_linprog_cyclic_bland(self):
+        # Test the effect of Bland's rule on a cycling problem
+        c = np.array([-10, 57, 9, 24.]) * (-1)
+        A_ub = np.array([[0.5, -5.5, -2.5, 9],
+                         [0.5, -1.5, -0.5, 1],
+                         [1, 0, 0, 0]])
+        b_ub = np.array([0, 0, 1])
+        desired_x = [1, 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)
+
+    def test_linprog_cyclic_bland_bug_8561(self):
+        # Test that pivot row is chosen correctly when using Bland's rule
+        c = np.array([7, 0, -4, 1.5, 1.5]) * (-1)
+        A_ub = np.array([
+            [4, 5.5, 1.5, 1.0, -3.5],
+            [1, -2.5, -2, 2.5, 0.5],
+            [3, -0.5, 4, -12.5, -7],
+            [-1, 4.5, 2, -3.5, -2],
+            [5.5, 2, -4.5, -1, 9.5]])
+        b_ub = np.array([0, 0, 0, 0, 1])
+        # desired_x = [0, 0, 19, 16/3, 29/3]
+        res = linprog_simplex(c, A_ub=A_ub, b_ub=b_ub)
+        _assert_success(res, c, b_ub=b_ub)
+
+    def test_bug_8662(self):
+        # linprog simplex used to report inncorrect optimal results
+        # https://github.com/scipy/scipy/issues/8662
+        c = np.array([-10, 10, 6, 3]) * (-1)
+        A_ub = [[8, -8, -4, 6],
+                [-8, 8, 4, -6],
+                [-4, 4, 8, -4],
+                [3, -3, -3, -10]]
+        b_ub = [9, -9, -9, -4]
+        c, A_ub, b_ub = map(np.asarray, [c, A_ub, b_ub])
+        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)
+
+
+if __name__ == '__main__':
+    import sys
+    import nose
+
+    argv = sys.argv[:]
+    argv.append('--verbose')
+    argv.append('--nocapture')
+    nose.main(argv=argv, defaultTest=__file__)

From 9cea0e06350c3882c436744c4cd69d13132940ce Mon Sep 17 00:00:00 2001
From: Kyohei Okumura <major0422@gmail.com>
Date: Thu, 28 Mar 2019 23:06:40 +0900
Subject: [PATCH 3/6] DOC: Add docstrings for linprog_simplex

---
 docs/source/optimize.rst                 |   1 +
 docs/source/optimize/linprog_simplex.rst |   7 +
 quantecon/optimize/linprog_simplex.py    | 262 ++++++++++++++++++++++-
 3 files changed, 265 insertions(+), 5 deletions(-)
 create mode 100644 docs/source/optimize/linprog_simplex.rst

diff --git a/docs/source/optimize.rst b/docs/source/optimize.rst
index 049b722fe..d971c19cb 100644
--- a/docs/source/optimize.rst
+++ b/docs/source/optimize.rst
@@ -6,5 +6,6 @@ Optimize
 
    optimize/nelder_mead
    optimize/root_finding
+   optimize/linprog_simplex
    optimize/pivoting
    optimize/scalar_maximization
diff --git a/docs/source/optimize/linprog_simplex.rst b/docs/source/optimize/linprog_simplex.rst
new file mode 100644
index 000000000..f7ca4a388
--- /dev/null
+++ b/docs/source/optimize/linprog_simplex.rst
@@ -0,0 +1,7 @@
+linprog_simplex
+===============
+
+.. automodule:: quantecon.optimize.linprog_simplex
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/quantecon/optimize/linprog_simplex.py b/quantecon/optimize/linprog_simplex.py
index 2c37083ce..4787953a3 100644
--- a/quantecon/optimize/linprog_simplex.py
+++ b/quantecon/optimize/linprog_simplex.py
@@ -20,6 +20,89 @@
 def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
                     A_eq=np.empty((0, 0)), b_eq=np.empty((0,)), max_iter=10**6,
                     tableau=None, basis=None, x=None, lambd=None):
+    """
+    Solve a linear program in the following form by the simplex
+    algorithm (with the lexicographic pivoting rule):
+
+    maximize::
+
+        c @ x
+
+    subject to::
+
+        A_ub @ x <= b_ub
+        A_eq @ x == b_eq
+               x >= 0
+
+    Parameters
+    ----------
+    c : ndarray(float, ndim=1)
+        ndarray of shape (n,).
+
+    A_ub : ndarray(float, ndim=2), optional
+        ndarray of shape (m, n).
+
+    b_ub : ndarray(float, ndim=1), optional
+        ndarray of shape (m,).
+
+    A_eq : ndarray(float, ndim=2), optional
+        ndarray of shape (k, n).
+
+    b_eq : ndarray(float, ndim=1), optional
+        ndarray of shape (k,).
+
+    max_iter : int, optional(default=10**6)
+        Maximum number of iteration to perform.
+
+    tableau : ndarray(float, ndim=2), optional
+        Temporary ndarray of shape (L+1, n+m+L+1) to store the tableau,
+        where L=m+k. Modified in place.
+
+    basis : ndarray(int, ndim=1), optional
+        Temporary ndarray of shape (L,) to store the basic variables.
+        Modified in place.
+
+    x : ndarray(float, ndim=1), optional
+        Output ndarray of shape (n,) to store the primal solution.
+
+    lambd : ndarray(float, ndim=1), optional
+        Output ndarray of shape (L,) to store the dual solution.
+
+    Returns
+    -------
+    res : SimplexResult
+        namedtuple consisting of the fields:
+
+            x : ndarray(float, ndim=1)
+                ndarray of shape (n,) containing the primal solution.
+
+            lambd : ndarray(float, ndim=1)
+                ndarray of shape (L,) containing the dual solution.
+
+            fun : float
+                Value of the objective function.
+
+            success : bool
+                True if the algorithm succeeded in finding an optimal
+                solution.
+
+            status : int
+                An integer representing the exit status of the
+                optimization::
+
+                    0 : Optimization terminated successfully
+                    1 : Iteration limit reached
+                    2 : Problem appears to be infeasible
+                    3 : Problem apperas to be unbounded
+
+            num_iter : int
+                The number of iterations performed.
+
+    References
+    ----------
+    * K. C. Border, "The Gauss–Jordan and Simplex Algorithms," 2004.
+
+    """
     n, m, k = c.shape[0], A_ub.shape[0], A_eq.shape[0]
     L = m + k
 
@@ -69,6 +152,81 @@ def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
 
 @jit(nopython=True, cache=True)
 def _initialize_tableau(A_ub, b_ub, A_eq, b_eq, tableau, basis):
+    """
+    Initialize the `tableau` and `basis` arrays in place for Phase 1.
+
+    Suppose that the original linear program has the following form:
+
+    maximize::
+
+        c @ x
+
+    subject to::
+
+        A_ub @ x <= b_ub
+        A_eq @ x == b_eq
+               x >= 0
+
+    Let s be a vector of slack variables converting the inequality
+    constraint to an equality constraint so that the problem turns to be
+    the standard form:
+
+    maximize::
+
+        c @ x
+
+    subject to::
+
+        A_ub @ x + s == b_ub
+        A_eq @ x     == b_eq
+        x, s         >= 0
+
+    Then, let (z1, z2) be a vector of artificial variables for Phase 1:
+    we solve the following LP:
+
+    maximize::
+
+        -(1 @ z1 + 1 @ z2)
+
+    subject to::
+
+        A_ub @ x + s + z1 == b_ub
+        A_eq @ x + z2     == b_eq
+        x, s, z1, z2      >= 0
+
+    The tableau needs to be of shape (L+1, n+m+L+1), where L=m+k.
+
+    Parameters
+    ----------
+    A_ub : ndarray(float, ndim=2)
+        ndarray of shape (m, n).
+
+    b_ub : ndarray(float, ndim=1)
+        ndarray of shape (m,).
+
+    A_eq : ndarray(float, ndim=2)
+        ndarray of shape (k, n).
+
+    b_eq : ndarray(float, ndim=1)
+        ndarray of shape (k,).
+
+    tableau : ndarray(float, ndim=2)
+        Empty ndarray of shape (L+1, n+m+L+1) to store the tableau.
+        Modified in place.
+
+    basis : ndarray(int, ndim=1)
+        Empty ndarray of shape (L,) to store the basic variables.
+        Modified in place.
+
+    Returns
+    -------
+    tableau : ndarray(float, ndim=2)
+        View to `tableau`.
+
+    basis : ndarray(int, ndim=1)
+        View to `basis`.
+
+    """
     m, k = A_ub.shape[0], A_eq.shape[0]
     L = m + k
     n = tableau.shape[1] - (m+L+1)
@@ -114,6 +272,27 @@ def _initialize_tableau(A_ub, b_ub, A_eq, b_eq, tableau, basis):
 
 @jit(nopython=True, cache=True)
 def _set_criterion_row(c, basis, tableau):
+    """
+    Modify the criterion row of the tableau for Phase 2.
+
+    Parameters
+    ----------
+    c : ndarray(float, ndim=1)
+        ndarray of shape (n,).
+
+    basis : ndarray(int, ndim=1)
+        ndarray of shape (L,) containing the basis obtained by Phase 1.
+
+    tableau : ndarray(float, ndim=2)
+        ndarray of shape (L+1, n+m+L+1) containing the tableau obtained
+        by Phase 1. Modified in place.
+
+    Returns
+    -------
+    tableau : ndarray(float, ndim=2)
+        View to `tableau`.
+
+    """
     n = c.shape[0]
     L = basis.shape[0]
 
@@ -136,11 +315,15 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
 
     Used to solve a linear program in the following form:
 
-        maximize:     c @ x
+    maximize::
 
-        subject to:   A_ub @ x <= b_ub
-                      A_eq @ x == b_eq
-                      x >= 0
+        c @ x
+
+    subject to::
+
+        A_ub @ x <= b_ub
+        A_eq @ x == b_eq
+               x >= 0
 
     where A_ub is of shape (m, n) and A_eq is of shape (k, n). Thus,
     `tableau` is of shape (L+1, n+m+L+1), where L=m+k, and
@@ -159,13 +342,24 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
         ndarray of shape (L,) containing the basic variables. Modified
         in place.
 
-    max_iter : scalar(int), optional(default=10**6)
+    max_iter : int, optional(default=10**6)
         Maximum number of pivoting steps.
 
     skip_aux : bool, optional(default=True)
         Whether to skip the coefficients of the auxiliary (or
         artificial) variables in pivot column selection.
 
+    Returns
+    -------
+    success : bool
+        True if the algorithm succeeded in finding an optimal solution.
+
+    status : int
+        An integer representing the exit status of the optimization.
+
+    num_iter : int
+        The number of iterations performed.
+
     """
     L = tableau.shape[0] - 1
 
@@ -203,6 +397,33 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
 
 @jit(nopython=True, cache=True)
 def _pivot_col(tableau, skip_aux):
+    """
+    Choose the column containing the pivot element: the column
+    containing the maximum positive element in the last row of the
+    tableau.
+
+    `skip_aux` should be True in phase 1, and False in phase 2.
+
+    Parameters
+    ----------
+    tableau : ndarray(float, ndim=2)
+        ndarray of shape (L+1, n+m+L+1) containing the tableau.
+
+    skip_aux : bool
+        Whether to skip the coefficients of the auxiliary (or
+        artificial) variables in pivot column selection.
+
+    Returns
+    -------
+    found : bool
+        True iff there is a positive element in the last row of the
+        tableau (and then pivotting should be conducted).
+
+    pivcol : int
+        The index of column containing the pivot element. (-1 if `found
+        == False`.)
+
+    """
     L = tableau.shape[0] - 1
     criterion_row_stop = tableau.shape[1] - 1
     if skip_aux:
@@ -222,6 +443,37 @@ def _pivot_col(tableau, skip_aux):
 
 @jit(nopython=True, cache=True)
 def get_solution(tableau, basis, x, lambd, b_signs):
+    """
+    Fetch the optimal solution and value from an optimal tableau.
+
+    Parameters
+    ----------
+    tableau : ndarray(float, ndim=2)
+        ndarray of shape (L+1, n+m+L+1) containing the optimal tableau,
+        where L=m+k.
+
+    basis : ndarray(int, ndim=1)
+        Empty ndarray of shape (L,) to store the basic variables.
+        Modified in place.
+
+    x : ndarray(float, ndim=1)
+        Empty ndarray of shape (n,) to store the primal solution.
+        Modified in place.
+
+    lambd : ndarray(float, ndim=1)
+        Empty ndarray of shape (L,) to store the dual solution. Modified
+        in place.
+
+    b_signs : ndarray(bool, ndim=1)
+        ndarray of shape (L,) whose i-th element is True iff the i-th
+        element of the vector (b_ub, b_eq) is positive.
+
+    Returns
+    -------
+    fun : float
+        The optimal value.
+
+    """
     n, L = x.size, lambd.size
     aux_start = tableau.shape[1] - L - 1
 

From 5da62e7239b52821e637594b51fde359b93fb639 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sat, 30 Mar 2019 22:14:07 +0900
Subject: [PATCH 4/6] pivoting: Add tol options to `_lex_min_ratio_test`

---
 quantecon/optimize/pivoting.py | 53 +++++++++++++++++++++++++---------
 1 file changed, 39 insertions(+), 14 deletions(-)

diff --git a/quantecon/optimize/pivoting.py b/quantecon/optimize/pivoting.py
index 18b312710..21c8ca19f 100644
--- a/quantecon/optimize/pivoting.py
+++ b/quantecon/optimize/pivoting.py
@@ -12,7 +12,7 @@
 
 
 @jit(nopython=True, cache=True)
-def _pivoting(tableau, pivot, pivot_row):
+def _pivoting(tableau, pivot_col, pivot_row):
     """
     Perform a pivoting step. Modify `tableau` in place.
 
@@ -21,8 +21,8 @@ def _pivoting(tableau, pivot, pivot_row):
     tableau : ndarray(float, ndim=2)
         Array containing the tableau.
 
-    pivot : scalar(int)
-        Pivot.
+    pivot_col : scalar(int)
+        Pivot column index.
 
     pivot_row : scalar(int)
         Pivot row index.
@@ -35,14 +35,14 @@ def _pivoting(tableau, pivot, pivot_row):
     """
     nrows, ncols = tableau.shape
 
-    pivot_elt = tableau[pivot_row, pivot]
+    pivot_elt = tableau[pivot_row, pivot_col]
     for j in range(ncols):
         tableau[pivot_row, j] /= pivot_elt
 
     for i in range(nrows):
         if i == pivot_row:
             continue
-        multiplier = tableau[i, pivot]
+        multiplier = tableau[i, pivot_col]
         if multiplier == 0:
             continue
         for j in range(ncols):
@@ -53,7 +53,8 @@ def _pivoting(tableau, pivot, pivot_row):
 
 @jit(nopython=True, cache=True)
 def _min_ratio_test_no_tie_breaking(tableau, pivot, test_col,
-                                    argmins, num_candidates):
+                                    argmins, num_candidates,
+                                    tol_piv, tol_ratio_diff):
     """
     Perform the minimum ratio test, without tie breaking, for the
     candidate rows in `argmins[:num_candidates]`. Return the number
@@ -78,6 +79,13 @@ def _min_ratio_test_no_tie_breaking(tableau, pivot, test_col,
     num_candidates : scalar(int)
         Number of candidate rows in `argmins`.
 
+    tol_piv : scalar(float)
+        Pivot tolerance below which a number is considered to be
+        nonpositive.
+
+    tol_ratio_diff : scalar(float)
+        Tolerance to determine a tie between ratio values.
+
     Returns
     -------
     num_argmins : scalar(int)
@@ -89,12 +97,12 @@ def _min_ratio_test_no_tie_breaking(tableau, pivot, test_col,
 
     for k in range(num_candidates):
         i = argmins[k]
-        if tableau[i, pivot] <= TOL_PIV:  # Treated as nonpositive
+        if tableau[i, pivot] <= tol_piv:  # Treated as nonpositive
             continue
         ratio = tableau[i, test_col] / tableau[i, pivot]
-        if ratio > ratio_min + TOL_RATIO_DIFF:  # Ratio large for i
+        if ratio > ratio_min + tol_ratio_diff:  # Ratio large for i
             continue
-        elif ratio < ratio_min - TOL_RATIO_DIFF:  # Ratio smaller for i
+        elif ratio < ratio_min - tol_ratio_diff:  # Ratio smaller for i
             ratio_min = ratio
             num_argmins = 1
         else:  # Ratio equal
@@ -105,7 +113,8 @@ def _min_ratio_test_no_tie_breaking(tableau, pivot, test_col,
 
 
 @jit(nopython=True, cache=True)
-def _lex_min_ratio_test(tableau, pivot, slack_start, argmins):
+def _lex_min_ratio_test(tableau, pivot, slack_start, argmins,
+                        tol_piv=TOL_PIV, tol_ratio_diff=TOL_RATIO_DIFF):
     """
     Perform the lexico-minimum ratio test.
 
@@ -124,6 +133,14 @@ def _lex_min_ratio_test(tableau, pivot, slack_start, argmins):
         Empty array used to store the row indices. Its length must be no
         smaller than the number of the rows of `tableau`.
 
+    tol_piv : scalar(float), optional
+        Pivot tolerance below which a number is considered to be
+        nonpositive. Default value is {TOL_PIV}.
+
+    tol_ratio_diff : scalar(float), optional
+        Tolerance to determine a tie between ratio values. Default value
+        is {TOL_RATIO_DIFF}.
+
     Returns
     -------
     found : bool
@@ -142,17 +159,25 @@ def _lex_min_ratio_test(tableau, pivot, slack_start, argmins):
     for i in range(nrows):
         argmins[i] = i
 
-    num_argmins = _min_ratio_test_no_tie_breaking(tableau, pivot, -1,
-                                                  argmins, num_candidates)
+    num_argmins = _min_ratio_test_no_tie_breaking(
+        tableau, pivot, -1, argmins, num_candidates, tol_piv, tol_ratio_diff
+    )
     if num_argmins == 1:
         found = True
     elif num_argmins >= 2:
         for j in range(slack_start, slack_start+nrows):
             if j == pivot:
                 continue
-            num_argmins = _min_ratio_test_no_tie_breaking(tableau, pivot, j,
-                                                          argmins, num_argmins)
+            num_argmins = _min_ratio_test_no_tie_breaking(
+                tableau, pivot, j, argmins, num_argmins,
+                tol_piv, tol_ratio_diff
+            )
             if num_argmins == 1:
                 found = True
                 break
     return found, argmins[0]
+
+
+_lex_min_ratio_test.__doc__ = _lex_min_ratio_test.__doc__.format(
+    TOL_PIV=TOL_PIV, TOL_RATIO_DIFF=TOL_RATIO_DIFF
+)

From 49991e5471d29eeedb3c5287a0ba8f90f5f6f69d Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Sun, 23 Jun 2019 15:07:18 +0900
Subject: [PATCH 5/6] Add `PivOptions`

---
 quantecon/optimize/__init__.py        |  4 +-
 quantecon/optimize/linprog_simplex.py | 62 +++++++++++++++++++++------
 2 files changed, 52 insertions(+), 14 deletions(-)

diff --git a/quantecon/optimize/__init__.py b/quantecon/optimize/__init__.py
index 372aef025..63cb03ce9 100644
--- a/quantecon/optimize/__init__.py
+++ b/quantecon/optimize/__init__.py
@@ -2,7 +2,9 @@
 """
 Initialization of the optimize subpackage
 """
-from .linprog_simplex import linprog_simplex, solve_tableau, get_solution
+from .linprog_simplex import (
+    linprog_simplex, solve_tableau, get_solution, PivOptions
+)
 from .scalar_maximization import brent_max
 from .nelder_mead import nelder_mead
 from .root_finding import newton, newton_halley, newton_secant, bisect, brentq
diff --git a/quantecon/optimize/linprog_simplex.py b/quantecon/optimize/linprog_simplex.py
index 4787953a3..6c02892e2 100644
--- a/quantecon/optimize/linprog_simplex.py
+++ b/quantecon/optimize/linprog_simplex.py
@@ -8,17 +8,24 @@
 from .pivoting import _pivoting, _lex_min_ratio_test
 
 
-FEA_TOL = 1e-6
-
-
 SimplexResult = namedtuple(
     'SimplexResult', ['x', 'lambd', 'fun', 'success', 'status', 'num_iter']
 )
 
+FEA_TOL = 1e-6
+TOL_PIV = 1e-7
+TOL_RATIO_DIFF = 1e-13
+
+PivOptions = namedtuple(
+    'PivOptions', ['fea_tol', 'tol_piv', 'tol_ratio_diff']
+)
+PivOptions.__new__.__defaults__ = (FEA_TOL, TOL_PIV, TOL_RATIO_DIFF)
+
 
 @jit(nopython=True, cache=True)
 def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
-                    A_eq=np.empty((0, 0)), b_eq=np.empty((0,)), max_iter=10**6,
+                    A_eq=np.empty((0, 0)), b_eq=np.empty((0,)),
+                    max_iter=10**6, piv_options=PivOptions(),
                     tableau=None, basis=None, x=None, lambd=None):
     """
     Solve a linear program in the following form by the simplex
@@ -54,6 +61,19 @@ def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
     max_iter : int, optional(default=10**6)
         Maximum number of iteration to perform.
 
+    piv_options : PivOptions, optional
+        PivOptions namedtuple to set the following tolerance values:
+
+            fea_tol : float
+                Feasibility tolerance (default={FEA_TOL}).
+
+            tol_piv : float
+                Pivot tolerance (default={TOL_PIV}).
+
+            tol_ratio_diff : float
+                Tolerance used in the the lexicographic pivoting
+                (default={TOL_RATIO_DIFF}).
+
     tableau : ndarray(float, ndim=2), optional
         Temporary ndarray of shape (L+1, n+m+L+1) to store the tableau,
         where L=m+k. Modified in place.
@@ -129,11 +149,12 @@ 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)
+        solve_tableau(tableau, basis, max_iter, skip_aux=False,
+                      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] > FEA_TOL:  # Infeasible
+    if tableau[-1, -1] > piv_options.fea_tol:  # Infeasible
         success = False
         status = 2
         return SimplexResult(x, lambd, fun, success, status, num_iter)
@@ -143,13 +164,19 @@ def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
 
     # Phase 2
     success, status, num_iter_2 = \
-        solve_tableau(tableau, basis, max_iter-num_iter, skip_aux=True)
+        solve_tableau(tableau, basis, max_iter-num_iter, skip_aux=True,
+                      piv_options=piv_options)
     num_iter += num_iter_2
     fun = get_solution(tableau, basis, x, lambd, b_signs)
 
     return SimplexResult(x, lambd, fun, success, status, num_iter)
 
 
+linprog_simplex.__doc__ = linprog_simplex.__doc__.format(
+    FEA_TOL=FEA_TOL, TOL_PIV=TOL_PIV, TOL_RATIO_DIFF=TOL_RATIO_DIFF
+)
+
+
 @jit(nopython=True, cache=True)
 def _initialize_tableau(A_ub, b_ub, A_eq, b_eq, tableau, basis):
     """
@@ -309,7 +336,8 @@ def _set_criterion_row(c, basis, tableau):
 
 
 @jit(nopython=True, cache=True)
-def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
+def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True,
+                  piv_options=PivOptions()):
     """
     Perform the simplex algorithm on a given tableau in canonical form.
 
@@ -349,6 +377,9 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
         Whether to skip the coefficients of the auxiliary (or
         artificial) variables in pivot column selection.
 
+    piv_options : PivOptions, optional
+        PivOptions namedtuple to set the tolerance values.
+
     Returns
     -------
     success : bool
@@ -373,7 +404,7 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
     while num_iter < max_iter:
         num_iter += 1
 
-        pivcol_found, pivcol = _pivot_col(tableau, skip_aux)
+        pivcol_found, pivcol = _pivot_col(tableau, skip_aux, piv_options)
 
         if not pivcol_found:  # Optimal
             success = True
@@ -381,8 +412,10 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
             break
 
         aux_start = tableau.shape[1] - L - 1
-        pivrow_found, pivrow = _lex_min_ratio_test(tableau[:-1, :], pivcol,
-                                                   aux_start, argmins)
+        pivrow_found, pivrow = _lex_min_ratio_test(
+            tableau[:-1, :], pivcol, aux_start, argmins,
+            piv_options.tol_piv, piv_options.tol_ratio_diff
+        )
 
         if not pivrow_found:  # Unbounded
             success = False
@@ -396,7 +429,7 @@ def solve_tableau(tableau, basis, max_iter=10**6, skip_aux=True):
 
 
 @jit(nopython=True, cache=True)
-def _pivot_col(tableau, skip_aux):
+def _pivot_col(tableau, skip_aux, piv_options):
     """
     Choose the column containing the pivot element: the column
     containing the maximum positive element in the last row of the
@@ -413,6 +446,9 @@ def _pivot_col(tableau, skip_aux):
         Whether to skip the coefficients of the auxiliary (or
         artificial) variables in pivot column selection.
 
+    piv_options : PivOptions, optional
+        PivOptions namedtuple to set the tolerance values.
+
     Returns
     -------
     found : bool
@@ -431,7 +467,7 @@ def _pivot_col(tableau, skip_aux):
 
     found = False
     pivcol = -1
-    coeff = FEA_TOL
+    coeff = piv_options.fea_tol
     for j in range(criterion_row_stop):
         if tableau[-1, j] > coeff:
             coeff = tableau[-1, j]

From ba0228406651e326ca50311594c046bf9e5df881 Mon Sep 17 00:00:00 2001
From: Daisuke Oyama <oyama@e.u-tokyo.ac.jp>
Date: Fri, 30 Apr 2021 14:15:48 +0900
Subject: [PATCH 6/6] Fixes in docstring

Co-authored-by: Quentin Batista <quentin.batista@mail.mcgill.ca>
---
 quantecon/optimize/linprog_simplex.py | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/quantecon/optimize/linprog_simplex.py b/quantecon/optimize/linprog_simplex.py
index 6c02892e2..deb8e6570 100644
--- a/quantecon/optimize/linprog_simplex.py
+++ b/quantecon/optimize/linprog_simplex.py
@@ -120,7 +120,7 @@ def linprog_simplex(c, A_ub=np.empty((0, 0)), b_ub=np.empty((0,)),
 
     References
     ----------
-    * K. C. Border, "The Gauss–Jordan and Simplex Algorithms," 2004.
+    * K. C. Border, "The Gauss–Jordan and Simplex Algorithms" 2004.
 
     """
     n, m, k = c.shape[0], A_ub.shape[0], A_eq.shape[0]
@@ -208,8 +208,8 @@ def _initialize_tableau(A_ub, b_ub, A_eq, b_eq, tableau, basis):
         A_eq @ x     == b_eq
         x, s         >= 0
 
-    Then, let (z1, z2) be a vector of artificial variables for Phase 1:
-    we solve the following LP:
+    Then, let (z1, z2) be a vector of artificial variables for Phase 1.
+    We solve the following LP:
 
     maximize::
 
@@ -453,7 +453,7 @@ def _pivot_col(tableau, skip_aux, piv_options):
     -------
     found : bool
         True iff there is a positive element in the last row of the
-        tableau (and then pivotting should be conducted).
+        tableau (and then pivoting should be conducted).
 
     pivcol : int
         The index of column containing the pivot element. (-1 if `found
@@ -489,8 +489,7 @@ def get_solution(tableau, basis, x, lambd, b_signs):
         where L=m+k.
 
     basis : ndarray(int, ndim=1)
-        Empty ndarray of shape (L,) to store the basic variables.
-        Modified in place.
+        ndarray of shape (L,) containing the optimal basis.
 
     x : ndarray(float, ndim=1)
         Empty ndarray of shape (n,) to store the primal solution.