def commuting_sets_by_indices(pauli_sums, commutation_check):
    """
    For a list of pauli sums, find commuting sets and keep track of which pauli sum they came from.

    :param pauli_sums: A list of PauliSum
    :param commutation_check: a function that checks if all elements of a list
                              and a single pauli term commute.
    :return: A list of commuting sets. Each set is a list of tuples (i, j) to find the particular
        commuting term. i is the index of the pauli sum from whence the term came. j is the
        index within the set.
    """
    assert isinstance(pauli_sums, list)

    group_inds = []
    group_terms = []
    for i, pauli_sum in enumerate(pauli_sums):
        for j, term in enumerate(pauli_sum):
            if len(group_inds) == 0:
                # Initialization
                group_inds.append([(i, j)])
                group_terms.append([term])
                continue

            for k, group in enumerate(group_terms):
                if commutation_check(group, term):
                    group_inds[k] += [(i, j)]
                    group_terms[k] += [term]
                    break
            else:
                # for ... else means loop completed without a `break`
                # Which means this needs to start its own group.
                group_inds.append([(i, j)])
                group_terms.append([term])

    return group_inds