Updating the payoff convex hull by iterating all action pairs.
    Using the R operator proposed by Abreu and Sannikov 2014.

    Parameters
    ----------
    delta : scalar(float)
            The common discount rate at which all players discount
            the future.

    nums_actions : tuple(int)
        Tuple of the numbers of actions, one for each player.

    payoff_arrays : tuple(ndarray(float, ndim=2))
        Tuple of the payoff arrays, one for each player.

    best_dev_gains : tuple(ndarray(float, ndim=2))
        Tuple of the normalized best deviation payoff gain arrays.
        best_dev_gains[i][ai, a-i] is payoff gain player i
        can get if originally players are choosing ai and a-i,
        and player i deviates to the best response action.

    points : ndarray(float, ndim=2)
        Coordinates of the points in the W, which construct a
        feasible payoff convex hull.

    vertices : ndarray(float, ndim=1)
        Indices of points forming the vertices of the convex hull,
        which are in counterclockwise order.

    equations : ndarray(float, ndim=2)
        [normal, offset] forming the hyperplane equation of the facet

    u : ndarray(float, ndim=1)
        The threat points.

    IC : ndarray(float, ndim=1)
        The minimum IC continuation values.

    action_profile_payoff : ndarray(float, ndim=1)
        Array of payoff for one action profile.

    extended_payoff : ndarray(float, ndim=2)
        The array [payoff0, payoff1, 1] for checking if
        [payoff0, payoff1] is in the feasible payoff convex hull.

    new_pts : ndarray(float, ndim=1)
        The 4 by 2 array for storing the generated potential
        extreme points of one action profile. One action profile
        can only generate at most 4 points.

    W_new : ndarray(float, ndim=2)
        Array for storing the coordinates of the generated potential
        extreme points that construct a new feasible payoff convex hull.

    tol: scalar(float), optional(default=1e-10)
        The tolerance for checking if two values are equal.

    Returns
    -------
    W_new : ndarray(float, ndim=2)
        The coordinates of the generated potential extreme points
        that construct a new feasible payoff convex hull.

    n_new_pt : scalar(int)
        The number of points in W_new that construct the feasible
        payoff convex hull.

def _R(delta, nums_actions, payoff_arrays, best_dev_gains, points,
       vertices, equations, u, IC, action_profile_payoff,
       extended_payoff, new_pts, W_new, tol=1e-10):
    """
    Updating the payoff convex hull by iterating all action pairs.
    Using the R operator proposed by Abreu and Sannikov 2014.

    Parameters
    ----------
    delta : scalar(float)
            The common discount rate at which all players discount
            the future.

    nums_actions : tuple(int)
        Tuple of the numbers of actions, one for each player.

    payoff_arrays : tuple(ndarray(float, ndim=2))
        Tuple of the payoff arrays, one for each player.

    best_dev_gains : tuple(ndarray(float, ndim=2))
        Tuple of the normalized best deviation payoff gain arrays.
        best_dev_gains[i][ai, a-i] is payoff gain player i
        can get if originally players are choosing ai and a-i,
        and player i deviates to the best response action.

    points : ndarray(float, ndim=2)
        Coordinates of the points in the W, which construct a
        feasible payoff convex hull.

    vertices : ndarray(float, ndim=1)
        Indices of points forming the vertices of the convex hull,
        which are in counterclockwise order.

    equations : ndarray(float, ndim=2)
        [normal, offset] forming the hyperplane equation of the facet

    u : ndarray(float, ndim=1)
        The threat points.

    IC : ndarray(float, ndim=1)
        The minimum IC continuation values.

    action_profile_payoff : ndarray(float, ndim=1)
        Array of payoff for one action profile.

    extended_payoff : ndarray(float, ndim=2)
        The array [payoff0, payoff1, 1] for checking if
        [payoff0, payoff1] is in the feasible payoff convex hull.

    new_pts : ndarray(float, ndim=1)
        The 4 by 2 array for storing the generated potential
        extreme points of one action profile. One action profile
        can only generate at most 4 points.

    W_new : ndarray(float, ndim=2)
        Array for storing the coordinates of the generated potential
        extreme points that construct a new feasible payoff convex hull.

    tol: scalar(float), optional(default=1e-10)
        The tolerance for checking if two values are equal.

    Returns
    -------
    W_new : ndarray(float, ndim=2)
        The coordinates of the generated potential extreme points
        that construct a new feasible payoff convex hull.

    n_new_pt : scalar(int)
        The number of points in W_new that construct the feasible
        payoff convex hull.
    """
    n_new_pt = 0
    for a0 in range(nums_actions[0]):
        for a1 in range(nums_actions[1]):
            action_profile_payoff[0] = payoff_arrays[0][a0, a1]
            action_profile_payoff[1] = payoff_arrays[1][a1, a0]
            IC[0] = u[0] + best_dev_gains[0][a0, a1]
            IC[1] = u[1] + best_dev_gains[1][a1, a0]

            # check if payoff is larger than IC
            if (action_profile_payoff >= IC).all():
                # check if payoff is inside the convex hull
                extended_payoff[:2] = action_profile_payoff
                if (np.dot(equations, extended_payoff) <= tol).all():
                    W_new[n_new_pt] = action_profile_payoff
                    n_new_pt += 1
                    continue

            new_pts, n = _find_C(new_pts, points, vertices, equations,
                                 extended_payoff, IC, tol)

            for i in range(n):
                W_new[n_new_pt] = \
                    delta * new_pts[i] + (1-delta) * action_profile_payoff
                n_new_pt += 1

    return W_new, n_new_pt