docs: def comb_indices(n, k):
      """``n``-dimensional version of itertools.combinations.
  
      Args:
          a (np.ndarray): The array from which to get combinations.
          k (int): The desired length of the combinations.
  
      Returns:
          np.ndarray: Indices that give the ``k``-combinations of ``n`` elements.
  
      Example:
          >>> n, k = 3, 2
          >>> data = np.arange(6).reshape(2, 3)
          >>> data[:, comb_indices(n, k)]
          array([[[0, 1],
                  [0, 2],
                  [1, 2]],
          <BLANKLINE>
                 [[3, 4],
                  [3, 5],
                  [4, 5]]])
      """
      # Count the number of combinations for preallocation
      count = comb(n, k, exact=True)
      # Get numpy iterable from ``itertools.combinations``
      indices = np.fromiter(
          chain.from_iterable(combinations(range(n), k)),
          int,
          count=(count * k))
      # Reshape output into the array of combination indicies
      return indices.reshape(-1, k)