# Source code for quantecon.optimize.linprog_simplex

"""
Contain a Numba-jitted linear programming solver 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

SimplexResult = namedtuple(
'SimplexResult', ['x', 'lambd', 'fun', 'success', 'status', 'num_iter']
)
SimplexResult.__doc__ = \
'namedtuple containing the result from linprog_simplex.'

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)
PivOptions.__doc__ = 'namedtuple to hold tolerance values for pivoting'

# Delete useless docstring for fields of namedtuple
def _del_field_docstring(nt):
for field in nt._fields:
getattr(nt, field).__doc__ = ''

for nt in [SimplexResult, PivOptions]:
_del_field_docstring(nt)

[docs]@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, piv_options=PivOptions(),
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.

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.

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

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,
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
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,
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):
"""
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)

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):
"""
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]

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

[docs]@jit(nopython=True, cache=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.

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 : 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.

piv_options : PivOptions, optional
PivOptions namedtuple to set the tolerance values.

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

# 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, piv_options)

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,
piv_options.tol_piv, piv_options.tol_ratio_diff
)

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, piv_options):
"""
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.

piv_options : PivOptions, optional
PivOptions namedtuple to set the tolerance values.

Returns
-------
found : bool
True iff there is a positive element in the last row of the
tableau (and then pivoting 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:
criterion_row_stop -= L

found = False
pivcol = -1
coeff = piv_options.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

[docs]@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)
ndarray of shape (L,) containing the optimal basis.

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

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