linprog_simplex [WIP].ipynb

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