def hamiltonian_path(G):
    """Returns a Hamiltonian path in the given tournament graph.

    Each tournament has a Hamiltonian path. If furthermore, the
    tournament is strongly connected, then the returned Hamiltonian path
    is a Hamiltonian cycle (by joining the endpoints of the path).

    Parameters
    ----------
    G : NetworkX graph
        A directed graph representing a tournament.

    Returns
    -------
    bool
        Whether the given graph is a tournament graph.

    Notes
    -----
    This is a recursive implementation with an asymptotic running time
    of $O(n^2)$, ignoring multiplicative polylogarithmic factors, where
    $n$ is the number of nodes in the graph.

    """
    if len(G) == 0:
        return []
    if len(G) == 1:
        return [arbitrary_element(G)]
    v = arbitrary_element(G)
    hampath = hamiltonian_path(G.subgraph(set(G) - {v}))
    # Get the index of the first node in the path that does *not* have
    # an edge to `v`, then insert `v` before that node.
    index = index_satisfying(hampath, lambda u: v not in G[u])
    hampath.insert(index, v)
    return hampath

def test_edge_cutset_random_graphs():
    for flow_func in flow_funcs:
        for i in range(3):
            G = nx.fast_gnp_random_graph(50, 0.25)
            if not nx.is_connected(G):
                ccs = iter(nx.connected_components(G))
                start = arbitrary_element(next(ccs))
                G.add_edges_from((start, arbitrary_element(c)) for c in ccs)
            cutset = nx.minimum_edge_cut(G, flow_func=flow_func)
            assert_equal(nx.edge_connectivity(G), len(cutset), msg=msg.format(flow_func.__name__))
            G.remove_edges_from(cutset)
            assert_false(nx.is_connected(G), msg=msg.format(flow_func.__name__))

    def test_quotient_graph_edge_relation(self):
        """Tests for specifying an alternate edge relation for the quotient
        graph.

        """
        G = nx.path_graph(5)
        identity = lambda u, v: u == v
        same_parity = lambda b, c: (arbitrary_element(b) % 2
                                    == arbitrary_element(c) % 2)
        actual = nx.quotient_graph(G, identity, same_parity)
        expected = nx.Graph()
        expected.add_edges_from([(0, 2), (0, 4), (2, 4)])
        expected.add_edge(1, 3)
        assert_true(nx.is_isomorphic(actual, expected))

def equivalence_classes(iterable, relation):
    """Returns the set of equivalence classes of the given ``iterable`` under
    the specified equivalence relation.

    ``relation`` must be a Boolean-valued function that takes two argument. It
    must represent an equivalence relation (that is, the relation induced by
    the function must be reflexive, symmetric, and transitive).

    The return value is a set of sets. It is a partition of the elements of
    ``iterable``; duplicate elements will be ignored so it makes the most sense
    for ``iterable`` to be a :class:`set`.

    """
    # For simplicity of implementation, we initialize the return value as a
    # list of lists, then convert it to a set of sets at the end of the
    # function.
    blocks = []
    # Determine the equivalence class for each element of the iterable.
    for y in iterable:
        # Each element y must be in *exactly one* equivalence class.
        #
        # Each block is guaranteed to be non-empty
        for block in blocks:
            x = arbitrary_element(block)
            if relation(x, y):
                block.append(y)
                break
        else:
            # If the element y is not part of any known equivalence class, it
            # must be in its own, so we create a new singleton equivalence
            # class for it.
            blocks.append([y])
    return {frozenset(block) for block in blocks}

def _connected_chordal_graph_cliques(G):
    """Returns the set of maximal cliques of a connected chordal graph."""
    if G.number_of_nodes() == 1:
        x = frozenset(G.nodes())
        return set([x])
    else:
        cliques = set()
        unnumbered = set(G.nodes())
        v = arbitrary_element(G)
        unnumbered.remove(v)
        numbered = set([v])
        clique_wanna_be = set([v])
        while unnumbered:
            v = _max_cardinality_node(G, unnumbered, numbered)
            unnumbered.remove(v)
            numbered.add(v)
            new_clique_wanna_be = set(G.neighbors(v)) & numbered
            sg = G.subgraph(clique_wanna_be)
            if _is_complete_graph(sg):
                new_clique_wanna_be.add(v)
                if not new_clique_wanna_be >= clique_wanna_be:
                    cliques.add(frozenset(clique_wanna_be))
                clique_wanna_be = new_clique_wanna_be
            else:
                raise nx.NetworkXError("Input graph is not chordal.")
        cliques.add(frozenset(clique_wanna_be))
        return cliques

def strategy_connected_sequential(G, colors, traversal='bfs'):
    """Returns an iterable over nodes in ``G`` in the order given by a
    breadth-first or depth-first traversal.

    ``traversal`` must be one of the strings ``'dfs'`` or ``'bfs'``,
    representing depth-first traversal or breadth-first traversal,
    respectively.

    The generated sequence has the property that for each node except
    the first, at least one neighbor appeared earlier in the sequence.

    ``G`` is a NetworkX graph. ``colors`` is ignored.

    """
    if traversal == 'bfs':
        traverse = nx.bfs_edges
    elif traversal == 'dfs':
        traverse = nx.dfs_edges
    else:
        raise nx.NetworkXError("Please specify one of the strings 'bfs' or"
                               " 'dfs' for connected sequential ordering")
    for component in nx.connected_component_subgraphs(G):
        source = arbitrary_element(component)
        # Yield the source node, then all the nodes in the specified
        # traversal order.
        yield source
        for (_, end) in traverse(component, source):
            yield end

def _find_chordality_breaker(G, s=None, treewidth_bound=sys.maxsize):
    """ Given a graph G, starts a max cardinality search
    (starting from s if s is given and from an arbitrary node otherwise)
    trying to find a non-chordal cycle.

    If it does find one, it returns (u,v,w) where u,v,w are the three
    nodes that together with s are involved in the cycle.
    """

    unnumbered = set(G)
    if s is None:
        s = arbitrary_element(G)
    unnumbered.remove(s)
    numbered = set([s])
#    current_treewidth = None
    current_treewidth = -1
    while unnumbered:  # and current_treewidth <= treewidth_bound:
        v = _max_cardinality_node(G, unnumbered, numbered)
        unnumbered.remove(v)
        numbered.add(v)
        clique_wanna_be = set(G[v]) & numbered
        sg = G.subgraph(clique_wanna_be)
        if _is_complete_graph(sg):
            # The graph seems to be chordal by now. We update the treewidth
            current_treewidth = max(current_treewidth, len(clique_wanna_be))
            if current_treewidth > treewidth_bound:
                raise nx.NetworkXTreewidthBoundExceeded(
                    "treewidth_bound exceeded: %s" % current_treewidth)
        else:
            # sg is not a clique,
            # look for an edge that is not included in sg
            (u, w) = _find_missing_edge(sg)
            return (u, v, w)
    return ()

def hamiltonian_path(G, source):
    source = arbitrary_element(G)
    neighbors = set(G[source]) - set([source])
    n = len(G)
    for target in neighbors:
        for path in nx.all_simple_paths(G, source, target):
            if len(path) == n:
                yield path

def is_aperiodic(G):
    """Return True if G is aperiodic.

    A directed graph is aperiodic if there is no integer k > 1 that
    divides the length of every cycle in the graph.

    Parameters
    ----------
    G : NetworkX DiGraph
        Graph

    Returns
    -------
    aperiodic : boolean
        True if the graph is aperiodic False otherwise

    Raises
    ------
    NetworkXError
        If G is not directed

    Notes
    -----
    This uses the method outlined in [1]_, which runs in O(m) time
    given m edges in G. Note that a graph is not aperiodic if it is
    acyclic as every integer trivial divides length 0 cycles.

    References
    ----------
    .. [1] Jarvis, J. P.; Shier, D. R. (1996),
       Graph-theoretic analysis of finite Markov chains,
       in Shier, D. R.; Wallenius, K. T., Applied Mathematical Modeling:
       A Multidisciplinary Approach, CRC Press.
    """
    if not G.is_directed():
        raise nx.NetworkXError(
            "is_aperiodic not defined for undirected graphs")

    s = arbitrary_element(G)
    levels = {s: 0}
    this_level = [s]
    g = 0
    l = 1
    while this_level:
        next_level = []
        for u in this_level:
            for v in G[u]:
                if v in levels:  # Non-Tree Edge
                    g = gcd(g, levels[u] - levels[v] + 1)
                else:  # Tree Edge
                    next_level.append(v)
                    levels[v] = l
        this_level = next_level
        l += 1
    if len(levels) == len(G):  # All nodes in tree
        return g == 1
    else:
        return g == 1 and nx.is_aperiodic(G.subgraph(set(G) - set(levels)))

def dominating_set(G, start_with=None):
    r"""Finds a dominating set for the graph G.

    A *dominating set* for a graph with node set *V* is a subset *D* of
    *V* such that every node not in *D* is adjacent to at least one
    member of *D* [1]_.

    Parameters
    ----------
    G : NetworkX graph

    start_with : node (default=None)
        Node to use as a starting point for the algorithm.

    Returns
    -------
    D : set
        A dominating set for G.

    Notes
    -----
    This function is an implementation of algorithm 7 in [2]_ which
    finds some dominating set, not necessarily the smallest one.

    See also
    --------
    is_dominating_set

    References
    ----------
    .. [1] http://en.wikipedia.org/wiki/Dominating_set

    .. [2] Abdol-Hossein Esfahanian. Connectivity Algorithms.
        http://www.cse.msu.edu/~cse835/Papers/Graph_connectivity_revised.pdf

    """
    all_nodes = set(G)
    if start_with is None:
        start_with = arbitrary_element(all_nodes)
    if start_with not in G:
        raise nx.NetworkXError('node {} is not in G'.format(start_with))
    dominating_set = {start_with}
    dominated_nodes = set(G[start_with])
    remaining_nodes = all_nodes - dominated_nodes - dominating_set
    while remaining_nodes:
        # Choose an arbitrary node and determine its undominated neighbors.
        v = remaining_nodes.pop()
        undominated_neighbors = set(G[v]) - dominating_set
        # Add the node to the dominating set and the neighbors to the
        # dominated set. Finally, remove all of those nodes from the set
        # of remaining nodes.
        dominating_set.add(v)
        dominated_nodes |= undominated_neighbors
        remaining_nodes -= undominated_neighbors
    return dominating_set

    def test_basic(self):
        # This example is from the Wikipedia article "Trie"
        # <https://en.wikipedia.org/wiki/Trie>.
        strings = ['a', 'to', 'tea', 'ted', 'ten', 'i', 'in', 'inn']
        T, root = nx.prefix_tree(strings)

        def source_label(v): return T.node[v]['source']

        # First, we check that the tree has the expected
        # structure. Recall that each node that corresponds to one of
        # the input strings has an edge to the NIL node.
        #
        # Consider the three children at level 1 in the trie.
        a, i, t = sorted(T[root], key=source_label)
        # Check the 'a' branch.
        assert_equal(len(T[a]), 1)
        nil = arbitrary_element(T[a])
        assert_equal(len(T[nil]), 0)
        # Check the 'i' branch.
        assert_equal(len(T[i]), 2)
        nil, in_ = sorted(T[i], key=source_label)
        assert_equal(len(T[nil]), 0)
        assert_equal(len(T[in_]), 2)
        nil, inn = sorted(T[in_], key=source_label)
        assert_equal(len(T[nil]), 0)
        assert_equal(len(T[inn]), 1)
        nil = arbitrary_element(T[inn])
        assert_equal(len(T[nil]), 0)
        # Check the 't' branch.
        te, to = sorted(T[t], key=source_label)
        assert_equal(len(T[to]), 1)
        nil = arbitrary_element(T[to])
        assert_equal(len(T[nil]), 0)
        tea, ted, ten = sorted(T[te], key=source_label)
        assert_equal(len(T[tea]), 1)
        assert_equal(len(T[ted]), 1)
        assert_equal(len(T[ten]), 1)
        nil = arbitrary_element(T[tea])
        assert_equal(len(T[nil]), 0)
        nil = arbitrary_element(T[ted])
        assert_equal(len(T[nil]), 0)
        nil = arbitrary_element(T[ten])
        assert_equal(len(T[nil]), 0)

        # Next, we check that the "sources" of each of the nodes is the
        # rightmost letter in the string corresponding to the path to
        # that node.
        assert_equal(source_label(root), None)
        assert_equal(source_label(a), 'a')
        assert_equal(source_label(i), 'i')
        assert_equal(source_label(t), 't')
        assert_equal(source_label(in_), 'n')
        assert_equal(source_label(inn), 'n')
        assert_equal(source_label(to), 'o')
        assert_equal(source_label(te), 'e')
        assert_equal(source_label(tea), 'a')
        assert_equal(source_label(ted), 'd')
        assert_equal(source_label(ten), 'n')
        assert_equal(source_label(NIL), NIL)

 def _recursive_build(H, A, source, avail):
     # Terminate once the flow has been compute to every node.
     if {source} == avail:
         return
     # pick an arbitrary node as the sink
     sink = arbitrary_element(avail - {source})
     # find the minimum cut and its weight
     value, (S, T) = nx.minimum_cut(H, source, sink)
     if H.is_directed():
         # check if the reverse direction has a smaller cut
         value_, (T_, S_) = nx.minimum_cut(H, sink, source)
         if value_ < value:
             value, S, T = value_, S_, T_
     # add edge with weight of cut to the aux graph
     A.add_edge(source, sink, weight=value)
     # recursively call until all but one node is used
     _recursive_build(H, A, source, avail.intersection(S))
     _recursive_build(H, A, sink, avail.intersection(T))

 def _to_string(formula, root):
     # If there are no children, this is a variable node.
     label = formula.node[root]['label']
     if not formula[root]:
         return label
     # Otherwise, this is an operator.
     children = formula[root]
     # If one child, the label must be a NOT operator.
     if len(children) == 1:
         child = arbitrary_element(children)
         return '{}({})'.format(label, _to_string(formula, child))
     # NB "left" and "right" here are a little misleading: there is
     # no order on the children of a node. That's okay because the
     # Boolean AND and OR operators are symmetric. It just means that
     # the order of the operands cannot be predicted and hence the
     # function does not necessarily behave the same way on every
     # invocation.
     left, right = formula[root]
     left_subformula = _to_string(formula, left)
     right_subformula = _to_string(formula, right)
     return '({} {} {})'.format(left_subformula, label, right_subformula)

    def construct(EdgeComponentAuxGraph, G):
        """Builds an auxiliary graph encoding edge-connectivity between nodes.

        Notes
        -----
        Given G=(V, E), initialize an empty auxiliary graph A.
        Choose an arbitrary source node s.  Initialize a set N of available
        nodes (that can be used as the sink). The algorithm picks an
        arbitrary node t from N - {s}, and then computes the minimum st-cut
        (S, T) with value w. If G is directed the the minimum of the st-cut or
        the ts-cut is used instead. Then, the edge (s, t) is added to the
        auxiliary graph with weight w. The algorithm is called recursively
        first using S as the available nodes and s as the source, and then
        using T and t. Recursion stops when the source is the only available
        node.

        Parameters
        ----------
        G : NetworkX graph
        """
        # workaround for classmethod decorator
        not_implemented_for('multigraph')(lambda G: G)(G)

        def _recursive_build(H, A, source, avail):
            # Terminate once the flow has been compute to every node.
            if {source} == avail:
                return
            # pick an arbitrary node as the sink
            sink = arbitrary_element(avail - {source})
            # find the minimum cut and its weight
            value, (S, T) = nx.minimum_cut(H, source, sink)
            if H.is_directed():
                # check if the reverse direction has a smaller cut
                value_, (T_, S_) = nx.minimum_cut(H, sink, source)
                if value_ < value:
                    value, S, T = value_, S_, T_
            # add edge with weight of cut to the aux graph
            A.add_edge(source, sink, weight=value)
            # recursively call until all but one node is used
            _recursive_build(H, A, source, avail.intersection(S))
            _recursive_build(H, A, sink, avail.intersection(T))

        # Copy input to ensure all edges have unit capacity
        H = G.__class__()
        H.add_nodes_from(G.nodes())
        H.add_edges_from(G.edges(), capacity=1)

        # A is the auxiliary graph to be constructed
        # It is a weighted undirected tree
        A = nx.Graph()

        # Pick an arbitrary node as the source
        if H.number_of_nodes() > 0:
            source = arbitrary_element(H.nodes())
            # Initialize a set of elements that can be chosen as the sink
            avail = set(H.nodes())

            # This constructs A
            _recursive_build(H, A, source, avail)

        # This class is a container the holds the auxiliary graph A and
        # provides access the the k_edge_components function.
        self = EdgeComponentAuxGraph()
        self.A = A
        self.H = H
        return self

def min_edge_cover(G, matching_algorithm=None):
    """Returns a set of edges which constitutes
    the minimum edge cover of the graph.

    A smallest edge cover can be found in polynomial time by finding
    a maximum matching and extending it greedily so that all nodes
    are covered.

    Parameters
    ----------
    G : NetworkX graph
        An undirected bipartite graph.

    matching_algorithm : function
        A function that returns a maximum cardinality matching in a
        given bipartite graph. The function must take one input, the
        graph ``G``, and return a dictionary mapping each node to its
        mate. If not specified,
        :func:`~networkx.algorithms.bipartite.matching.hopcroft_karp_matching`
        will be used. Other possibilities include
        :func:`~networkx.algorithms.bipartite.matching.eppstein_matching`,
        or matching algorithms in the
        :mod:`networkx.algorithms.matching` module.

    Returns
    -------
    min_cover : set

        It contains all the edges of minimum edge cover
        in form of tuples. It contains both the edges `(u, v)` and `(v, u)`
        for given nodes `u` and `v` among the edges of minimum edge cover.

    Notes
    -----
    An edge cover of a graph is a set of edges such that every node of
    the graph is incident to at least one edge of the set.
    The minimum edge cover is an edge covering of smallest cardinality.

    Due to its implementation, the worst-case running time of this algorithm
    is bounded by the worst-case running time of the function
    ``matching_algorithm``.

    Minimum edge cover for bipartite graph can also be found using the
    function present in :mod:`networkx.algorithms.bipartite.covering`
    """
    if nx.number_of_isolates(G) > 0:
        # ``min_cover`` does not exist as there is an isolated node
        raise nx.NetworkXException(
            "Graph has a node with no edge incident on it, "
            "so no edge cover exists.")
    if matching_algorithm is None:
        matching_algorithm = partial(nx.max_weight_matching,
                                     maxcardinality=True)
    maximum_matching = matching_algorithm(G)
    # ``min_cover`` is superset of ``maximum_matching``
    min_cover = set(maximum_matching.items())
    # iterate for uncovered nodes
    uncovered_nodes = set(G) - {v for u, v in min_cover}
    for v in uncovered_nodes:
        # Since `v` is uncovered, each edge incident to `v` will join it
        # with a covered node (otherwise, if there were an edge joining
        # uncovered nodes `u` and `v`, the maximum matching algorithm
        # would have found it), so we can choose an arbitrary edge
        # incident to `v`. (This applies only in a simple graph, not a
        # multigraph.)
        u = arbitrary_element(G[v])
        min_cover.add((u, v))
        min_cover.add((v, u))
    return min_cover

 def same_parity(b, c):
     return (arbitrary_element(b) % 2 == arbitrary_element(c) % 2)

def _select_starting_cell(G, starting_edge=None):
    """ Select a cell to initiate _find_partition

    Parameters
    ----------
    G : NetworkX Graph
    starting_edge: an edge to build the starting cell from

    Returns
    -------
    Tuple of vertices in G

    Raises
    ------
    NetworkXError
        If it is determined that G is not a line graph

    Notes
    -----
    If starting edge not specified then pick an arbitrary edge - doesn't
    matter which. However, this function may call itself requiring a
    specific starting edge. Note that the r, s notation for counting
    triangles is the same as in the Roussopoulos paper cited above.
    """
    if starting_edge is None:
        e = arbitrary_element(list(G.edges()))
    else:
        e = starting_edge
        if e[0] not in G[e[1]]:
            msg = 'starting_edge (%s, %s) is not in the Graph'
            raise nx.NetworkXError(msg % e)
    e_triangles = _triangles(G, e)
    r = len(e_triangles)
    if r == 0:
        # there are no triangles containing e, so the starting cell is just e
        starting_cell = e
    elif r == 1:
        # there is exactly one triangle, T, containing e. If other 2 edges
        # of T belong only to this triangle then T is starting cell
        T = e_triangles[0]
        a, b, c = T
        # ab was original edge so check the other 2 edges
        ac_edges = [x for x in _triangles(G, (a, c))]
        bc_edges = [x for x in _triangles(G, (b, c))]
        if len(ac_edges) == 1:
            if len(bc_edges) == 1:
                starting_cell = T
            else:
                return _select_starting_cell(G, starting_edge=(b, c))
        else:
            return _select_starting_cell(G, starting_edge=(a, c))
    else:
        # r >= 2 so we need to count the number of odd triangles, s
        s = 0
        odd_triangles = []
        for T in e_triangles:
            if _odd_triangle(G, T):
                s += 1
                odd_triangles.append(T)
        if r == 2 and s == 0:
            # in this case either triangle works, so just use T
            starting_cell = T
        elif r - 1 <= s <= r:
            # check if odd triangles containing e form complete subgraph
            # there must be exactly s+2 of them
            # and they must all be connected
            triangle_nodes = set([])
            for T in odd_triangles:
                for x in T:
                    triangle_nodes.add(x)
            if len(triangle_nodes) == s + 2:
                for u in triangle_nodes:
                    for v in triangle_nodes:
                        if u != v and (v not in G.neighbors(u)):
                            msg = "G is not a line graph (odd triangles " \
                                  "do not form complete subgraph)"
                            raise nx.NetworkXError(msg)
                # otherwise then we can use this as the starting cell
                starting_cell = tuple(triangle_nodes)
            else:
                msg = "G is not a line graph (odd triangles " \
                      "do not form complete subgraph)"
                raise nx.NetworkXError(msg)
        else:
            msg = "G is not a line graph (incorrect number of " \
                  "odd triangles around starting edge)"
            raise nx.NetworkXError(msg)
    return starting_cell

def tree_all_pairs_lowest_common_ancestor(G, root=None, pairs=None):
    r"""Yield the lowest common ancestor for sets of pairs in a tree.

    Parameters
    ----------
    G : NetworkX directed graph (must be a tree)

    root : node, optional (default: None)
        The root of the subtree to operate on.
        If None, assume the entire graph has exactly one source and use that.

    pairs : iterable or iterator of pairs of nodes, optional (default: None)
        The pairs of interest. If None, Defaults to all pairs of nodes
        under `root` that have a lowest common ancestor.

    Returns
    -------
    lcas : generator of tuples `((u, v), lca)` where `u` and `v` are nodes
        in `pairs` and `lca` is their lowest common ancestor.

    Notes
    -----
    Only defined on non-null trees represented with directed edges from
    parents to children. Uses Tarjan's off-line lowest-common-ancestors
    algorithm. Runs in time $O(4 \times (V + E + P))$ time, where 4 is the largest
    value of the inverse Ackermann function likely to ever come up in actual
    use, and $P$ is the number of pairs requested (or $V^2$ if all are needed).

    Tarjan, R. E. (1979), "Applications of path compression on balanced trees",
    Journal of the ACM 26 (4): 690-715, doi:10.1145/322154.322161.

    See Also
    --------
    all_pairs_lowest_common_ancestor (similar routine for general DAGs)
    lowest_common_ancestor           (just a single pair for general DAGs)
    """
    if len(G) == 0:
        raise nx.NetworkXPointlessConcept("LCA meaningless on null graphs.")
    elif None in G:
        raise nx.NetworkXError("None is not a valid node.")

    # Index pairs of interest for efficient lookup from either side.
    if pairs is not None:
        pair_dict = defaultdict(set)
        # See note on all_pairs_lowest_common_ancestor.
        if not isinstance(pairs, (Mapping, Set)):
            pairs = set(pairs)
        for u, v in pairs:
            for n in (u, v):
                if n not in G:
                    msg = "The node %s is not in the digraph." % str(n)
                    raise nx.NodeNotFound(msg)
            pair_dict[u].add(v)
            pair_dict[v].add(u)

    # If root is not specified, find the exactly one node with in degree 0 and
    # use it. Raise an error if none are found, or more than one is. Also check
    # for any nodes with in degree larger than 1, which would imply G is not a
    # tree.
    if root is None:
        for n, deg in G.in_degree:
            if deg == 0:
                if root is not None:
                    msg = "No root specified and tree has multiple sources."
                    raise nx.NetworkXError(msg)
                root = n
            elif deg > 1:
                msg = "Tree LCA only defined on trees; use DAG routine."
                raise nx.NetworkXError(msg)
    if root is None:
        raise nx.NetworkXError("Graph contains a cycle.")

    # Iterative implementation of Tarjan's offline lca algorithm
    # as described in CLRS on page 521 (2nd edition)/page 584 (3rd edition)
    uf = UnionFind()
    ancestors = {}
    for node in G:
        ancestors[node] = uf[node]

    colors = defaultdict(bool)
    for node in nx.dfs_postorder_nodes(G, root):
        colors[node] = True
        for v in (pair_dict[node] if pairs is not None else G):
            if colors[v]:
                # If the user requested both directions of a pair, give it.
                # Otherwise, just give one.
                if pairs is not None and (node, v) in pairs:
                    yield (node, v), ancestors[uf[v]]
                if pairs is None or (v, node) in pairs:
                    yield (v, node), ancestors[uf[v]]
        if node != root:
            parent = arbitrary_element(G.pred[node])
            uf.union(parent, node)
            ancestors[uf[parent]] = parent