Find one mixed-action Nash equilibrium of a 2-player normal form
    game by the Lemke-Howson algorithm [2]_, implemented with
    "complementary pivoting" (see, e.g., von Stengel [3]_ for details).

    Parameters
    ----------
    g : NormalFormGame
        NormalFormGame instance with 2 players.

    init_pivot : scalar(int), optional(default=0)
        Initial pivot, an integer k such that 0 <= k < m+n, where
        integers 0, ..., m-1 and m, ..., m+n-1 correspond to the actions
        of players 0 and 1, respectively.

    max_iter : scalar(int), optional(default=10**6)
        Maximum number of pivoting steps.

    capping : scalar(int), optional(default=None)
        If supplied, the routine is executed with the heuristics
        proposed by Codenotti et al. [1]_; see Notes below for details.

    full_output : bool, optional(default=False)
        If False, only the computed Nash equilibrium is returned. If
        True, the return value is `(NE, res)`, where `NE` is the Nash
        equilibrium and `res` is a `NashResult` object.

    Returns
    -------
    NE : tuple(ndarray(float, ndim=1))
        Tuple of computed Nash equilibrium mixed actions.

    res : NashResult
        Object containing information about the computation. Returned
        only when `full_output` is True. See `NashResult` for details.

    Examples
    --------
    Consider the following game from von Stengel [3]_:

    >>> np.set_printoptions(precision=4)  # Reduce the digits printed
    >>> bimatrix = [[(3, 3), (3, 2)],
    ...             [(2, 2), (5, 6)],
    ...             [(0, 3), (6, 1)]]
    >>> g = NormalFormGame(bimatrix)

    Obtain a Nash equilibrium of this game by `lemke_howson` with player
    0's action 1 (out of the three actions 0, 1, and 2) as the initial
    pivot:

    >>> lemke_howson(g, init_pivot=1)
    (array([ 0.    ,  0.3333,  0.6667]), array([ 0.3333,  0.6667]))
    >>> g.is_nash(_)
    True

    Additional information is returned if `full_output` is set True:

    >>> NE, res = lemke_howson(g, init_pivot=1, full_output=True)
    >>> res.converged  # Whether the routine has converged
    True
    >>> res.num_iter  # Number of pivoting steps performed
    4

    Notes
    -----
    * This routine is implemented with floating point arithmetic and
      thus is subject to numerical instability.

    * If `capping` is set to a positive integer, the routine is executed
      with the heuristics proposed by [1]_:

      * For k = `init_pivot`, `init_pivot` + 1, ..., `init_pivot` +
        (m+n-2), (modulo m+n), the Lemke-Howson algorithm is executed
        with k as the initial pivot and `capping` as the maximum number
        of pivoting steps. If the algorithm converges during this loop,
        then the Nash equilibrium found is returned.

      * Otherwise, the Lemke-Howson algorithm is executed with
        `init_pivot` + (m+n-1) (modulo m+n) as the initial pivot, with a
        limit `max_iter` on the total number of pivoting steps.

      Accoding to the simulation results for *uniformly random games*,
      for medium- to large-size games this heuristics outperforms the
      basic Lemke-Howson algorithm with a fixed initial pivot, where
      [1]_ suggests that `capping` be set to 10.

    References
    ----------
    .. [1] B. Codenotti, S. De Rossi, and M. Pagan, "An Experimental
       Analysis of Lemke-Howson Algorithm," arXiv:0811.3247, 2008.

    .. [2] C. E. Lemke and J. T. Howson, "Equilibrium Points of Bimatrix
       Games," Journal of the Society for Industrial and Applied
       Mathematics (1964), 413-423.

    .. [3] B. von Stengel, "Equilibrium Computation for Two-Player Games
       in Strategic and Extensive Form," Chapter 3, N. Nisan, T.
       Roughgarden, E. Tardos, and V. Vazirani eds., Algorithmic Game
       Theory, 2007.

def lemke_howson(g, init_pivot=0, max_iter=10**6, capping=None,
                 full_output=False):
    """
    Find one mixed-action Nash equilibrium of a 2-player normal form
    game by the Lemke-Howson algorithm [2]_, implemented with
    "complementary pivoting" (see, e.g., von Stengel [3]_ for details).

    Parameters
    ----------
    g : NormalFormGame
        NormalFormGame instance with 2 players.

    init_pivot : scalar(int), optional(default=0)
        Initial pivot, an integer k such that 0 <= k < m+n, where
        integers 0, ..., m-1 and m, ..., m+n-1 correspond to the actions
        of players 0 and 1, respectively.

    max_iter : scalar(int), optional(default=10**6)
        Maximum number of pivoting steps.

    capping : scalar(int), optional(default=None)
        If supplied, the routine is executed with the heuristics
        proposed by Codenotti et al. [1]_; see Notes below for details.

    full_output : bool, optional(default=False)
        If False, only the computed Nash equilibrium is returned. If
        True, the return value is `(NE, res)`, where `NE` is the Nash
        equilibrium and `res` is a `NashResult` object.

    Returns
    -------
    NE : tuple(ndarray(float, ndim=1))
        Tuple of computed Nash equilibrium mixed actions.

    res : NashResult
        Object containing information about the computation. Returned
        only when `full_output` is True. See `NashResult` for details.

    Examples
    --------
    Consider the following game from von Stengel [3]_:

    >>> np.set_printoptions(precision=4)  # Reduce the digits printed
    >>> bimatrix = [[(3, 3), (3, 2)],
    ...             [(2, 2), (5, 6)],
    ...             [(0, 3), (6, 1)]]
    >>> g = NormalFormGame(bimatrix)

    Obtain a Nash equilibrium of this game by `lemke_howson` with player
    0's action 1 (out of the three actions 0, 1, and 2) as the initial
    pivot:

    >>> lemke_howson(g, init_pivot=1)
    (array([ 0.    ,  0.3333,  0.6667]), array([ 0.3333,  0.6667]))
    >>> g.is_nash(_)
    True

    Additional information is returned if `full_output` is set True:

    >>> NE, res = lemke_howson(g, init_pivot=1, full_output=True)
    >>> res.converged  # Whether the routine has converged
    True
    >>> res.num_iter  # Number of pivoting steps performed
    4

    Notes
    -----
    * This routine is implemented with floating point arithmetic and
      thus is subject to numerical instability.

    * If `capping` is set to a positive integer, the routine is executed
      with the heuristics proposed by [1]_:

      * For k = `init_pivot`, `init_pivot` + 1, ..., `init_pivot` +
        (m+n-2), (modulo m+n), the Lemke-Howson algorithm is executed
        with k as the initial pivot and `capping` as the maximum number
        of pivoting steps. If the algorithm converges during this loop,
        then the Nash equilibrium found is returned.

      * Otherwise, the Lemke-Howson algorithm is executed with
        `init_pivot` + (m+n-1) (modulo m+n) as the initial pivot, with a
        limit `max_iter` on the total number of pivoting steps.

      Accoding to the simulation results for *uniformly random games*,
      for medium- to large-size games this heuristics outperforms the
      basic Lemke-Howson algorithm with a fixed initial pivot, where
      [1]_ suggests that `capping` be set to 10.

    References
    ----------
    .. [1] B. Codenotti, S. De Rossi, and M. Pagan, "An Experimental
       Analysis of Lemke-Howson Algorithm," arXiv:0811.3247, 2008.

    .. [2] C. E. Lemke and J. T. Howson, "Equilibrium Points of Bimatrix
       Games," Journal of the Society for Industrial and Applied
       Mathematics (1964), 413-423.

    .. [3] B. von Stengel, "Equilibrium Computation for Two-Player Games
       in Strategic and Extensive Form," Chapter 3, N. Nisan, T.
       Roughgarden, E. Tardos, and V. Vazirani eds., Algorithmic Game
       Theory, 2007.

    """
    try:
        N = g.N
    except:
        raise TypeError('g must be a 2-player NormalFormGame')
    if N != 2:
        raise NotImplementedError('Implemented only for 2-player games')

    payoff_matrices = g.payoff_arrays
    nums_actions = g.nums_actions
    total_num = sum(nums_actions)

    msg = '`init_pivot` must be an integer k' + \
          'such that 0 <= k < {0}'.format(total_num)

    if not isinstance(init_pivot, numbers.Integral):
        raise TypeError(msg)

    if not (0 <= init_pivot < total_num):
        raise ValueError(msg)

    if capping is None:
        capping = max_iter

    tableaux = tuple(
        np.empty((nums_actions[1-i], total_num+1)) for i in range(N)
    )
    bases = tuple(np.empty(nums_actions[1-i], dtype=int) for i in range(N))

    converged, num_iter, init_pivot_used = \
        _lemke_howson_capping(payoff_matrices, tableaux, bases, init_pivot,
                              max_iter, capping)
    NE = _get_mixed_actions(tableaux, bases)

    if not full_output:
        return NE

    res = NashResult(NE=NE,
                     converged=converged,
                     num_iter=num_iter,
                     max_iter=max_iter,
                     init=init_pivot_used)

    return NE, res