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