"""
Implements a linear programming solver using the "Simplex" method.
"""
import numpy as np
from numba import jit, generated_jit, types
infeasible_err_msg = "The problem appears to be infeasible"
unbounded_obj = "The problem appears to be unbounded."
max_iter_reached = "The maximum number of iterations has been reached."
@jit(nopython=True, cache=True)
def linprog_simplex(tableau, N, M_ub=0, M_eq=0, maxiter=10000, tol=1e-10):
"""
Solve the following LP problem using the original Simplex method with
Bland's rule:
Minimize: c.T @ x
Subject to: A_ub @ x ≤ b_ub
A_eq @ x = b_eq
x ≥ 0
Jit-compiled in `nopython` mode.
Parameters
----------
tableau : ndarray(float, ndim=2)
2-D array containing the standardized LP problem in detached
coefficients form augmented with the artificial variables, the
infeasibility form and with nonnegative constant terms. The
infeasibility form is assumed to be placed in the last row, and the
objective function in the second last row.
N : scalar(int)
Number of control variables.
M_ub : scalar(int)
Number of inequality constraints.
M_eq : scalar(int)
Number of equality constraints.
maxiter : scalar(int), optional(default=10000)
Maximum number of pivot operation for each phase.
tol : scalar(float), optional(default=1e-10)
Tolerance for treating an element as nonpositive.
Return
----------
sol : ndarray(float, ndim=1)
A basic solution to the linear programming problem.
References
----------
[1] Dantzig, George B., Linear programming and extensions. Rand Corporation
Research Study Princeton Univ. Press, Princeton, NJ, 1963
[2] Bland, Robert G. (May 1977). "New finite pivoting rules for the
simplex method". Mathematics of Operations Research. 2 (2): 103–107.
[2] http://mat.gsia.cmu.edu/classes/QUANT/NOTES/chap7.pdf
"""
M = M_ub + M_eq
# Phase I
num_vars = N + M_ub + M
simplex_algorithm(tableau, M, num_vars, maxiter=maxiter)
if abs(tableau[-1, -1]) > tol:
raise ValueError(infeasible_err_msg)
nonpos_coeff = tableau[-1, :] <= tol
tableau = tableau[:-1, nonpos_coeff]
# Phase II
num_vars = 0
for coef in nonpos_coeff[:-1]:
if coef:
num_vars += 1
simplex_algorithm(tableau, M, num_vars, maxiter=maxiter)
sol = _find_basic_solution(tableau, nonpos_coeff, N, M, tol=tol)
return sol
@jit(nopython=True, cache=True)
def simplex_algorithm(tableau, M, num_vars, maxiter=10000, tol=1e-10):
"""
Execute the simplex algorithm on `tableau` using Bland's rule. Jit-compiled
in `nopython` mode.
Parameters
----------
tableau : ndarray(float, ndim=2)
2-D array to be modified inplace which contains the LP problem in
detached coefficients form. The objective is assumed to be placed in
the last row.
M : scalar(int)
Total number of constraints in the LP problem.
numvars : scalar(int)
Number of variables in the objective function.
maxiter : scalar(int), optional(default=10000)
Maximum number of pivot operation for each phase.
tol : scalar(float), optional(default=1e-10)
Tolerance for treating an element as nonpositive.
"""
pivot_col = _search_negative_cost_factor(tableau, num_vars, tol=tol)
argmins = np.arange(0, M)
for num_iter in range(maxiter):
if pivot_col == -1:
return
num_argmins = min_ratio_test(tableau, pivot_col, -1, argmins, M)
# Check if there is no lower bound
if num_argmins == 0:
raise ValueError(unbounded_obj)
pivot_row = argmins[:num_argmins].min() # Bland's rule
pivot_operation(tableau, (pivot_row, pivot_col))
# Update `pivot_col`
pivot_col = _search_negative_cost_factor(tableau, num_vars, tol=tol)
# Restore `argmins`
for i in range(M):
argmins[i] = i
raise ValueError(max_iter_reached)
@jit(nopython=True, cache=True)
def _search_negative_cost_factor(tableau, num_vars, tol=1e-10):
"""
Search for the first negative relative cost factor in the last row of
`tableau`. Jit-compiled in `nopython` mode.
Parameters
----------
tableau : ndarray(float, ndim=2)
2-D array to be modified inplace which contains the LP problem in
detached coefficients form.
num_vars : scalar(int)
Number of variables in the objective function.
tol : scalar(float), optional(default=1e-10)
Tolerance for treating an element as nonpositive.
Return
----------
idx : scalar(int)
The index of the variable with a negative coefficient and the lowest
column index. If all variables have nonnegative coefficients, return
-1.
"""
for idx in range(num_vars):
if tableau[-1, idx] < -tol:
return idx
return -1
@jit(nopython=True, cache=True)
def _find_basic_solution(tableau, nonpos_coeff, N, M, tol=1e-10):
"""
Find a basic solution to the LP problem in `tableau`. Jit-compiled in
`nopython` mode.
Parameters
----------
tableau : ndarray(float, ndim=2)
2-D array to be modified inplace which contains the LP problem in
detached coefficients form.
nonpos_coeff : ndarray(bool, ndim=1)
1-D array indicating which coefficients of the infeasibility form are
nonpositive at the end of Phase I.
N : scalar(int)
Number of control variables.
M : scalar(int)
Total number of constraints in the LP problem.
Return
----------
sol : ndarray(float, ndim=1)
A basic solution to the LP problem.
"""
sol = 1. * nonpos_coeff[:N]
idx = 0
for i in range(N):
if sol[i] == 1.: # Variable has not been dropped
if abs(tableau[M, idx]) > tol:
sol[i] = 0.
idx += 1
continue
not_found_one = True
for j in range(M):
if abs(tableau[j, idx]) <= tol:
continue
elif abs(tableau[j, idx]-1.) <= tol and not_found_one:
sol[i] = tableau[j, -1]
not_found_one = False
else:
sol[i] = 0.
break
idx += 1
return sol
"""
Useful routines for manipulating linear equation systems.
"""
import numpy as np
from numba import jit, generated_jit, types
cons_err_msg = "At least one tpye of constraints must be specified."
ub_err_msg = "Inequality constraints are not properly specified."
eq_err_msg = "Equality constraints are not properly specified."
@jit(nopython=True, cache=True)
def make_tableau(c, A_ub=np.array([[]]).T, b_ub=np.array([[]]),
A_eq=np.array([[]]).T, b_eq=np.array([[]])):
"""
Create a tableau for an LP problem given an objective function and
constraints by transforming the problem to its standard form, making,
constants nonnegative, augmenting it with artificial variables and with an
infeasibility form. Jit-compiled in `nopython` mode.
Parameters
----------
c : ndarray(float, ndim=1)
1-D array containing the coefficients of the `N` variables of the
linear objective function to be minimized.
A_ub : ndarray(float, ndim=2)
2-D array containing the coefficients of the left hand side of the
`M_ub` inequality constraints in `N` variables.
b_ub : ndarray(float, ndim=1)
1-D array containing the values of the right hand side of the `M_ub`
inequality constraints.
A_eq : ndarray(float, ndim=2)
2-D array containing the coefficients of the left hand side of the
`M_eq` equality constraints in `N` variables.
b_eq : ndarray(float, ndim=1)
1-D array containing the values of the right hand side of the `M_eq`
equality constraints.
Return
----------
tableau : ndarray(float, ndim=2)
2-D array containing the standardized LP problem in detached
coefficients form augmented with the artificial variables, the
infeasibility form and with nonnegative constant terms.
"""
M_ub, N_ub = A_ub.shape
M_eq, N_eq = A_eq.shape
M = M_ub + M_eq
N = max(N_ub, N_eq)
tableau = np.zeros((M+2, N+M_ub+M+1))
standardize_lp_problem(c, A_ub, b_ub, A_eq, b_eq, tableau)
# Make constraints nonnegative
for i in range(M):
if tableau[i, -1] < 0:
tableau[i, :] = -tableau[i, :]
# Add artificial variables
for (row_idx, col_idx) in zip(range(M), range(N+M_ub, N+M_ub+M)):
tableau[row_idx, col_idx] = 1.
# Add infeasability form
for col_idx in range(N+M_ub):
for row_idx in range(M):
tableau[-1, col_idx] -= tableau[row_idx, col_idx]
for row_idx in range(M):
tableau[-1, -1] -= tableau[row_idx, -1]
return tableau
@jit(nopython=True, cache=True)
def standardize_lp_problem(c, A_ub=np.array([[]]).T, b_ub=np.array([[]]),
A_eq=np.array([[]]).T, b_eq=np.array([[]]),
tableau=None):
"""
Standardize an LP problem of the following form:
Objective: c.T @ x
Subject to: A_ub @ x ≤ b_ub
A_eq @ x = b_eq
Jit-compiled in `nopython` mode.
Parameters
----------
c : ndarray(float, ndim=1)
1-D array containing the coefficients of the `N` variables of the
linear objective function to be minimized.
A_ub : ndarray(float, ndim=2)
2-D array containing the coefficients of the left hand side of the
`M_ub` inequality constraints in `N` variables.
b_ub : ndarray(float, ndim=1)
1-D array containing the values of the right hand side of the `M_ub`
inequality constraints.
A_eq : ndarray(float, ndim=2)
2-D array containing the coefficients of the left hand side of the
`M_eq` equality constraints in `N` variables.
b_eq : ndarray(float, ndim=1)
1-D array containing the values of the right hand side of the `M_eq`
equality constraints.
tableau : ndarray(float, ndim=2), optional(default=None)
2-D array to be modified inplace which will contain the standardized
LP problem in detached coefficients form. If there are any, the
inequality constrains are stacked on top of the equality constraints.
The constant terms are placed in the last column. The objective is
placed in row `M_ub+M_eq`.
Return
----------
tableau : ndarray(float, ndim=2)
View of `tableau`.
"""
M_ub, N_ub = A_ub.shape
M_eq, N_eq = A_eq.shape
M = M_ub + M_eq
N = max(N_ub, N_eq)
if M_ub != b_ub.size:
raise ValueError(ub_err_msg)
if M_eq != b_eq.size:
raise ValueError(eq_err_msg)
if M_ub > 0 and N > 0: # At least inequality constraints
if tableau is None:
tableau = np.zeros((M+1, N+M_ub+1))
# Place the inequality contraints in the tableau
tableau[0:M_ub, -1] = b_ub
tableau[0:M_ub, 0:N] = A_ub
if M > M_ub: # Both type of constraints
# Place equality constraints in the tableau
tableau[M_ub:M, -1] = b_eq
tableau[M_ub:M, 0:N] = A_eq
# Add the slack variables
for (row_idx, col_idx) in zip(range(M_ub), range(N, N+M_ub)):
# Make diagonal elements equal to one for slack variables part
tableau[row_idx, col_idx] = 1.
# Standardize the objective function
tableau[M, 0:N] = c
return tableau
elif M_eq > 0 and N > 0: # Only equality constraints
if tableau is None:
tableau = np.zeros((M+1, N+1))
# Place the equality constraints in the tableau
tableau[0:M, -1] = b_eq
tableau[0:M, 0:N] = A_eq
# Standardize the objective function
tableau[M, 0:N] = c
return tableau
else:
raise ValueError(cons_err_msg)
@jit(nopython=True, cache=True)
def pivot_operation(tableau, pivot):
"""
Perform a pivoting step. Modify `tableau` in place.
Parameters
----------
tableau : ndarray(float, ndim=2)
Array containing the tableau.
pivot : tuple(int)
Tuple containing the row and column index of the pivot element.
Returns
-------
tableau : ndarray(float, ndim=2)
View of `tableau`.
"""
nrows, ncols = tableau.shape
pivot_row, pivot_col = pivot
pivot_elt = tableau[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_col]
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(tableau, pivot_col, test_col, argmins, num_candidates,
tol_piv=1e-10, tol_ratio_diff=1e-15):
"""
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_col : scalar(int)
Index of the column of the pivot element.
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`.
tol_piv : scalar(float), optional(default=1e-10)
Tolerance for treating an element of the pivot column as nonpositive.
tol_ratio_diff : scalar(float), optional(default=1e-15)
Tolerance for comparing candidate minimum ratios.
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_col] <= tol_piv: # Treated as nonpositive
continue
ratio = tableau[i, test_col] / tableau[i, pivot_col]
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
import numpy as np
from scipy.optimize import linprog
linprog_test = np.load('linprog_benchmark_files/ADLITTLE.npz')
M_ub, N = linprog_test['A_ub'].shape
M_eq, N = linprog_test['A_eq'].shape
tableau = make_tableau(linprog_test['c'], linprog_test['A_ub'], linprog_test['b_ub'], linprog_test['A_eq'], linprog_test['b_eq'])
x = linprog_simplex(tableau, N, M_ub, M_eq)
%timeit linprog_simplex(tableau, N, M_ub, M_eq)
%timeit linprog(linprog_test['c'], linprog_test['A_ub'], linprog_test['b_ub'], linprog_test['A_eq'], linprog_test['b_eq'], method='simplex')
%timeit linprog(linprog_test['c'], linprog_test['A_ub'], linprog_test['b_ub'], linprog_test['A_eq'], linprog_test['b_eq'], method='interior-point')
1000 loops, best of 3: 240 µs per loop 10 loops, best of 3: 21.3 ms per loop The slowest run took 7.53 times longer than the fastest. This could mean that an intermediate result is being cached. 100 loops, best of 3: 13.2 ms per loop
def gridmake(*arrays):
"""
TODO: finish this docstring
Notes
-----
Based of original function ``gridmake`` in CompEcon toolbox by
Miranda and Fackler
References
----------
Miranda, Mario J, and Paul L Fackler. Applied Computational Economics
and Finance, MIT Press, 2002.
"""
if all([i.ndim == 1 for i in arrays]):
d = len(arrays)
if d == 2:
out = _gridmake2(*arrays)
else:
out = _gridmake2(arrays[0], arrays[1])
for arr in arrays[2:]:
out = _gridmake2(out, arr)
return out
else:
raise NotImplementedError("Come back here")
def _gridmake2(x1, x2):
"""
TODO: finish this docstring
Notes
-----
Based of original function ``gridmake2`` in CompEcon toolbox by
Miranda and Fackler
References
----------
Miranda, Mario J, and Paul L Fackler. Applied Computational Economics
and Finance, MIT Press, 2002.
"""
if x1.ndim == 1 and x2.ndim == 1:
return np.column_stack([np.tile(x1, x2.shape[0]),
np.repeat(x2, x1.shape[0])])
elif x1.ndim > 1 and x2.ndim == 1:
first = np.tile(x1, (x2.shape[0], 1))
second = np.repeat(x2, x1.shape[0])
return np.column_stack([first, second])
else:
raise NotImplementedError("Come back here")
class RepeatedGame:
"""
Class representing an N-player repeated game.
Parameters
----------
stage_game : NormalFormGame
The stage game used to create the repeated game.
delta : scalar(float)
The common discount rate at which all players discount the future.
"""
def __init__(self, stage_game, delta):
self.sg = stage_game
self.delta = delta
self.N = stage_game.N
self.nums_actions = stage_game.nums_actions
def outerapproximation(self, nH=32, tol=1e-8, maxiter=500,
check_pure_nash=True, verbose=False, nskipprint=50,
linprog_method='simplex', tol_linprog=1e-10):
"""
Approximates the set of equilibrium value set for a repeated game with
the outer hyperplane approximation described by Judd, Yeltekin,
Conklin 2002.
Parameters
----------
rpd : RepeatedGame
Two player repeated game.
nH : scalar(int), optional(default=32)
Number of subgradients used for the approximation.
tol : scalar(float), optional(default=1e-8)
Tolerance in differences of set.
maxiter : scalar(int), optional(default=500)
Maximum number of iterations.
check_pure_nash : bool, optional(default=True)
Whether to perform a check about whether a pure Nash equilibrium
exists.
verbose : bool, optional(default=False)
Whether to display updates about iterations and distance.
nskipprint : scalar(int), optional(default=50)
Number of iterations between printing information (assuming
verbose=true).
linprog_method : ,optional(default='simplex')
The name of the scipy linprog method used to solve linear
programming problems.
Returns
-------
vertices : ndarray(float, ndim=2)
Vertices of the outer approximation of the value set.
"""
sg, delta = self.sg, self.delta
p0, p1 = sg.players
flat_payoff_array0 = p0.payoff_array.flatten()
flat_payoff_array1 = p1.payoff_array.flatten()
try:
next(pure_nash_brute_gen(sg))
except StopIteration:
raise ValueError("No pure action Nash equilibrium exists in" +
" stage game")
# Get number of actions for each player and create action space
nA0, nA1 = p0.num_actions, p1.num_actions
nAS = nA0 * nA1
AS = gridmake(np.array(range(nA0)), np.array(range(nA1)))
# Create the unit circle, points, and hyperplane levels
C, H, Z = initialize_sg_hpl(flat_payoff_array0, flat_payoff_array1, nH)
Cnew = C.copy()
# Create matrices for linear programming
c, A, b = initialize_LP_matrices(delta, H)
# bounds on w are [-Inf, Inf] while bounds on slack are [0, Inf]
lb = -np.inf
ub = np.inf
# Allocate space to store all solutions
Cia = np.zeros(nAS)
Wia = np.zeros([2, nAS])
# Set iterative parameters and iterate until converged
itr, dist = 0, 10.0
while (itr < maxiter) and (dist > tol):
# Compute the current worst values for each agent
_w1 = worst_value_1(self, H, C, tol_linprog=tol_linprog)
_w2 = worst_value_2(self, H, C, tol_linprog=tol_linprog)
# Iterate over all subgradients
for ih in range(nH):
#
# Subgradient specific instructions
#
h1, h2 = H[ih, :]
# Update all set constraints - Copies elements 1:nH of C into b
b[:nH] = C
# Put the right objective into c (negative because want
# maximize)
c[0] = -h1
c[1] = -h2
for ia in range(nAS):
#
# Action specific instruction
#
a1, a2 = AS[ia, :]
# Update incentive constraints
b[nH] = (1-delta)*flow_u_1(self, a1, a2) - \
(1-delta)*best_dev_payoff_1(self, a2) - delta*_w1
b[nH+1] = (1-delta)*flow_u_2(self, a1, a2) - \
(1-delta)*best_dev_payoff_2(self, a1) - delta*_w2
# Pull out optimal value and compute
M_ub, N = A.shape
tableau = make_tableau(c, A, b)
w_sol = linprog_simplex(tableau, N, M_ub, tol=tol_linprog)
value = (1-delta)*flow_u(self, a1, a2) + delta*w_sol
# Save hyperplane level and continuation promises
Cia[ia] = h1*value[0] + h2*value[1]
Wia[:, ia] = value
# Action which pushes furthest in direction h_i
astar = np.argmax(Cia)
a1star, a2star = AS[astar, :]
# Get hyperplane level and continuation value
Cstar = Cia[astar]
Wstar = Wia[:, astar]
if Cstar > -1e10:
Cnew[ih] = Cstar
else:
raise Error("Failed to find feasible action/continuation" +
" pair")
# Update the points
Z[:, ih] = \
(1-delta)*flow_u(self, a1star, a2star) + \
delta*np.array([Wstar[0], Wstar[1]])
# Update distance and iteration counter
dist = np.max(abs(C - Cnew))
itr += 1
if verbose and (nskipprint % itr == 0):
#TODO: Fix this
print("")
if itr >= maxiter:
warn("Maximum Iteration Reached")
# Update hyperplane levels
C[:] = Cnew
# Given the H-representation `(H, C)` of the computed polytope of
# equilibrium payoff profiles, we obtain its V-representation
# `vertices` using scipy
p = HalfspaceIntersection(np.column_stack((H, -C)), np.mean(Z, axis=1))
vertices = p.intersections
# Reduce the number of vertices by rounding points to the tolerance
tol_int = int(round(abs(np.log10(tol))) - 1)
# Find vertices that are unique within tolerance level
vertices = np.vstack({tuple(row) for row in
np.round(vertices, tol_int)})
return vertices
def unitcircle(npts):
"""
Places `npts` equally spaced points along the 2 dimensional circle and
returns the points with x coordinates in first column and y coordinates
in second column.
Parameters
----------
npts : scalar(float)
Number of points.
Returns
-------
pts : ndarray(float, dim=2)
The coordinates of equally spaced points.
"""
degrees = np.linspace(0, 2*np.pi, npts+1)
pts = np.zeros((npts, 2))
for i in range(npts):
x = degrees[i]
pts[i, 0] = np.cos(x)
pts[i, 1] = np.sin(x)
return pts
def initialize_hpl(nH, o, r):
"""
Initializes subgradients, extreme points and hyperplane levels for the
approximation of the convex value set of a 2 player repeated game.
Parameters
----------
nH : scalar(int)
Number of subgradients used for the approximation.
o : ndarray(float, ndim=2)
Origin for the approximation.
r: scalar(float)
Radius for the approximation.
Returns
-------
C : ndarray(float, ndim=1)
The array containing the hyperplane levels.
H : ndarray(float, ndim=2)
The array containing the subgradients.
Z : ndarray(float, ndim=2)
The array containing the extreme points of the value set.
"""
# Create unit circle
H = unitcircle(nH)
HT = H.T
# Choose origin and radius for big approximation
Z = np.zeros((2, nH))
C = np.zeros(nH)
for i in range(nH):
temp_val0 = o[0] + r*HT[0, i]
temp_val1 = o[1] + r*HT[1, i]
Z[0, i] = temp_val0
Z[1, i] = temp_val1
C[i] = temp_val0*HT[0, i] + temp_val1*HT[1, i]
return C, H, Z
def initialize_sg_hpl(flat_payoff_array0, flat_payoff_array1, nH):
"""
Initializes subgradients, extreme points and hyperplane levels for the
approximation of the convex value set of a 2 player repeated game by
choosing an appropriate origin and radius.
Parameters
----------
flat_payoff_array0 : ndarray(float, ndim=2)
Flattened payoff array for player 0.
flat_payoff_array1 : ndarray(float, ndim=2)
Flattened payoff array for player 1.
nH : scalar(int)
Number of subgradients used for the approximation.
Returns
-------
C : ndarray(float, ndim=1)
The array containing the hyperplane levels.
H : ndarray(float, ndim=2)
The array containing the subgradients.
Z : ndarray(float, ndim=2)
The array containing the extreme points of the value set.
"""
# Choose the origin to be mean of max and min payoffs
p0_min, p0_max = minmax(flat_payoff_array0)
p1_min, p1_max = minmax(flat_payoff_array1)
o = np.array([(p0_min + p0_max)/2.0, (p1_min + p1_max)/2.0])
r0 = np.max(((p0_max - o[0])**2, (o[0] - p0_min)**2))
r1 = np.max(((p1_max - o[1])**2, (o[1] - p1_min)**2))
r = np.sqrt(r0 + r1)
return initialize_hpl(nH, o, r)
def initialize_LP_matrices(delta, H):
"""
Initializes matrices for linear programming problems.
Parameters
----------
delta : scalar(float)
The common discount rate at which all players discount the future.
H : ndarray(float, ndim=2)
Subgradients used to approximate value set.
Returns
-------
c : ndarray(float, ndim=1)
Vector used to determine which subgradient is being used.
A : ndarray(float, ndim=2)
Matrix with nH set constraints and to be filled with 2 additional
incentive compatibility constraints.
b : ndarray(float, ndim=1)
Vector to be filled with the values for the constraints.
"""
# Total number of subgradients
nH = H.shape[0]
# Create the c vector (objective)
c = np.zeros(2)
# Create the A matrix (constraints)
A = np.zeros((nH+2, 2))
A[0:nH, :] = H
A[nH, :] = -delta, 0.
A[nH+1, :] = 0., -delta
# Create the b vector (constraints)
b = np.zeros(nH + 2)
return c, A, b
# Flow utility in terms of the players actions
def flow_u_1(rpd, a1, a2): return rpd.sg.players[0].payoff_array[a1, a2]
def flow_u_2(rpd, a1, a2): return rpd.sg.players[1].payoff_array[a2, a1]
def flow_u(rpd, a1, a2):
return np.array([flow_u_1(rpd, a1, a2), flow_u_2(rpd, a1, a2)])
# Computes each players best deviation given an opponent's action
def best_dev_i(rpd, i, aj):
return np.argmax(rpd.sg.players[i].payoff_array[:, aj])
def best_dev_1(rpd, a2): return best_dev_i(rpd, 0, a2)
def best_dev_2(rpd, a1): return best_dev_i(rpd, 1, a1)
# Computes the payoff of the best deviation
def best_dev_payoff_i(rpd, i, aj):
return max(rpd.sg.players[i].payoff_array[:, aj])
def best_dev_payoff_1(rpd, a2):
return max(rpd.sg.players[0].payoff_array[:, a2])
def best_dev_payoff_2(rpd, a1):
return max(rpd.sg.players[1].payoff_array[:, a1])
def worst_value_i(rpd, H, C, i, linprog_method='simplex', tol_linprog=1e-10):
"""
Returns the worst possible payoff for player i.
Parameters
----------
rpd : RepeatedGame
Two player repeated game instance.
H : ndarray(float, ndim=2)
Subgradients used to approximate value set.
C : ndarray(float, ndim=1)
Hyperplane levels used to approximate the value set.
i : scalar(int)
The player of interest.
Returns
-------
out : scalar(float)
Worst possible payoff of player i.
"""
# Objective depends on which player we are minimizing
c = np.zeros(2)
c[i] = 1.0
M_ub, N = H.shape
tableau = make_tableau(c, A_ub=H, b_ub=C)
out = linprog_simplex(tableau, N, M_ub, tol=tol_linprog)
return out[i]
def worst_value_1(rpd, H, C, linprog_method='simplex', tol_linprog=1e-10):
return worst_value_i(rpd, H, C, 0, linprog_method, tol_linprog=tol_linprog)
def worst_value_2(rpd, H, C, linprog_method='simplex', tol_linprog=1e-10):
return worst_value_i(rpd, H, C, 1, linprog_method, tol_linprog=tol_linprog)
def worst_values(rpd, H, C, linprog_method='simplex'):
return (worst_value_1(rpd, H, C, linprog_method),
worst_value_2(rpd, H, C, linprog_method))
def minmax(x):
maximum = x[0]
minimum = x[0]
for i in x[1:]:
if i > maximum:
maximum = i
elif i < minimum:
minimum = i
return (minimum, maximum)
import numpy as np
from quantecon.game_theory import Player, NormalFormGame, pure_nash_brute_gen
from scipy.spatial import HalfspaceIntersection
pd_payoff = np.array([[9.0, 1.0], [10.0, 3.0]])
A = Player(pd_payoff)
B = Player(pd_payoff)
nfg = NormalFormGame((A, B))
rpd = RepeatedGame(nfg, 0.75)
rpd.outerapproximation(tol_linprog=1e-10)
array([[10. , 3.9726605], [ 3. , 3. ], [ 3. , 10. ], [ 3.9726605, 10. ], [ 9. , 9. ], [10. , 3. ]])