Repository: networkx/networkx
Pull request: Fix graph_hash iteration counts and DiGraph handling
Fixes #7806 
Includes bugfixes for the issues discussed.
Does not cover the potential extension to ```MultiGraph```s.

Extended testing, but note that a lot of the test in this module in depend on randomly generated graphs. We continued that here. However, that does not necessarily mean all edge-cases are covered (as the original bug showed).
Note that subgraph hashing requires hashing of initial labels for digest size testing. This is not necessarily (and thus omitted) for general graph hashing.

~~A small optimisation is also made to the original graph hashing method. It used to collect _all_ labels calculated during execution in a counter, and hashed that to get the final graph hash. However, the WL algorithm only relies on the colouring (labels) calculated at the last iteration, so the current versionl only hashes those to obtain the final hash.~~



Problem statement:
Bugs WL (sub)graph hashing in relation to use degree at initialisation
### Current Behavior

At initialisation in both `weisfeiler_lehman_graph_hash` and `weisfeiler_lehman_subgraph_hashes` for graphs that do not have node or edge attributes, node weights are initialised using the degree of the nodes. This already aggregates information about neighbours at a 1-hop distance (the direct neighbours). 
This has three main consequences.
- The explanation and examples in the documentation are incorrect. It gives an example that requires 4 neighbourhood aggregations for certain nodes to converge to the same label, but incorrectly finds only 3 (the aggregation during node initialisation is not visible).
- The behaviour of graphs with and without attributes (even if those attributes are trivial/irrelevant) is not in line with each other. This is because graphs with node attributes skip this step.
- It hides a bug in the implementation for directed nodes: certain graphs are now correctly distinguished because of this initial step, but fail to be distinguished by graphs with attributes.


### Expected Behavior

Using a node's degree can be an efficient way to aggregate information about the 1-hop neighbourhood. However, this should then be used as the first iteration, instead of during initialisation. It should not be used for directed graphs (at least not in this way) because apparently it treats outgoing and ingoing edges the same.

Directed graphs should consider both in- and outgoing edges at each iteration of the algorithm.

### Steps to Reproduce

The iteration problem on undirected graphs:
```python
# This is the example used in the documentation.
G1 = nx.Graph()
G1.add_edges_from([(1, 2), (2, 3), (2, 4), (3, 5), (4, 6), (5, 7), (6, 7)])
G2 = nx.Graph()
G2.add_edges_from([(1, 3), (2, 3), (1, 6), (1, 5), (4, 6)])
g1_hashes = nx.weisfeiler_lehman_subgraph_hashes(G1, iterations=3, digest_size=8)
g2_hashes = nx.weisfeiler_lehman_subgraph_hashes(G2, iterations=3, digest_size=8)
# diverge at iteration 3"
print(g2_hashes[5])
print(g1_hashes[1])
# ['a93b64973cfc8897', 'db1b43ae35a1878f', '1716d2a4012fa4bc']
# ['a93b64973cfc8897', 'db1b43ae35a1878f', '57872a7d2059c1c0']


# now give both G1 and G2 node weights, all the same weight
for i in range(1, 8):
    G1.nodes[i]['weight'] = 1
for i in range(1, 7):
    G2.nodes[i]['weight'] = 1
# run the WL algorithm again, but this time with node weights
g1_hashes = nx.weisfeiler_lehman_subgraph_hashes(G1, node_attr='weight', iterations=4, digest_size=8)
g2_hashes = nx.weisfeiler_lehman_subgraph_hashes(G2, node_attr='weight', iterations=4, digest_size=8)
# diverge at iteration 4"
print(g2_hashes[5])
print(g1_hashes[1])
# ['d9614b235de46ab6', 'c2e8da8e81cb1547', 'b04e4278628353d2', 'ad7c8ff452869678']
# ['d9614b235de46ab6', 'c2e8da8e81cb1547', 'b04e4278628353d2', '70167e752a1990c2']
```
The (hidden) problem with directed graphs:
```python
# Digraphs without weights use neighbourhood aggregation, which gives them information about both in- and outgoing edges.
Ga = nx.DiGraph()
Gb = nx.DiGraph()
Ga.add_nodes_from([1, 2, 3, 4])
Gb.add_nodes_from([1, 2, 3, 4])
Ga.add_edges_from([(1,2), (3,2)])
Gb.add_edges_from([(1,2), (3,4)])
Ga_hash = nx.weisfeiler_lehman_graph_hash(Ga)
Gb_hash = nx.weisfeiler_lehman_graph_hash(Gb)
print(Ga_hash == Gb_hash) # False

# Now with trivial weights
for i in range(1, 5):
    Ga.nodes[i]['weight'] = 1
    Gb.nodes[i]['weight'] = 1
Ga_hash = nx.weisfeiler_lehman_graph_hash(Ga, node_attr='weight')
Gb_hash = nx.weisfeiler_lehman_graph_hash(Gb, node_attr='weight')
print(Ga_hash == Gb_hash) # True
```

### Environment

<!--- Please provide details about your local environment -->

Python version: 3.12.3
NetworkX version: 3.4.2

### Additional context

Apologies for posting what are effectively two bugs in one issue, since they both at least partially arise from the same underlying issue I had trouble separating them.

If convenient, I'd be willing to make a pull request to fix this. Especially the fix for undirected graphs should be fairly straightforward. I'd need some more time to figure out how to fix the bug on directed graphs.

Hints from issue discussion:
If I understand correctly you are saying: setting initial node labels to the node's degree:
- shifts the iteration count by one compared to the same graph with trivial node attributes
- gives hashes for `DiGraph`s that differ from what is produced with trivial node attributes

So, a short-term fix for users is to always use node attributes. 

And tests of the function should include comparing hashes between the case with no node attributes and the case with trivial node attributes. And these tests will be needed for both `nx.Graph` and `nx.DiGraph`.

A long term fix would involve adjusting the iteration count for the case of no node attribute, and to initialize the DiGraph node labels to be a 2-tuple of `(in_degree, out_degree)` values.

Have I missed (or incorrectly summarized) anything?
How do edge attributes fit with all of this?
For (undirected) `Graph`s  this is all correct!
I believe the version with edge attributes should work similar to the one with node attributes. I haven't tried it, but in the code I see it is initialised in a similar way (without using the degree of the nodes explicitly).

The concern with `DiGraph`s is a bit more complex. As shown in the example above, the method with node attributes cannot distinguish between two clearly isomorphic graphs. So that would suggest that there instead we should rely on the version *without* node attributes. That feels untrustworthy as well, however, because that relies on G.degree, which (if I understand correctly) does not make any distinction between in/outgoing nodes. I have, however, not yet been able to construct an example where this makes the algorithm fail.  
In fixing that version, your suggestion on initialisation seems great. However, to guarantee that it always works, the following iterations should be handled similarly. Currently, it uses G.neighbours, which only considers the outgoing neighbours of each node, whereas it should consider both in  and outgoing neighbours. (To be absolutely fair with you, I could imagine changing only the initialisation would work too. However, in the literature, WL for directed graphs is defined as considering both in/outgoing nodes at every iteration, so only doing that at initialisation is not *guaranteed* to work).

Thus a type distinction for DiGraph with a neighborhood aggregation step using `DiGraph.successor` and `DiGraph.predecessor` would be needed, I think. Also, we would need to account for the separation of successor and predecessor labels in each aggregation step. Thus we aggregate them in a tuple to differentiate between successor node labels and predecessor node labels as @dschult suggests. 

```
def _neighborhood_aggregate(G, node, node_labels, edge_attr=None):
    """
    Compute new labels for a given node by aggregating
    the labels of each node's neighbors.
    """
    label_list = []
    for nbr in G.neighbors(node):
        prefix = "" if edge_attr is None else str(G[node][nbr][edge_attr])
        label_list.append(prefix + node_labels[nbr])
    return node_labels[node] + "".join(sorted(label_list)) 
``` 

would then be changed to differentiate between Graph and DiGraph and in the case of a DiGraph be in line like: 

```
def _digraph_neighborhood_aggregate(G, node, node_labels, edge_attr=None):
    """
    Compute new labels for a given node in a DiGraph by aggregating
    the labels of each successor and predecessor node separately.
    """
    succ_list = []
    for succ in G.successor(node):
        prefix = "" if edge_attr is None else str(G[node][succ][edge_attr])
        succ_list.append(prefix + node_labels[succ_list])

    pred_list = []
    for pred in G.predecessors(node):
        prefix = "" if edge_attr is None else str(G[node][pred][edge_attr])
        pred_list.append(prefix + node_labels[pred])

    return node_labels[node] + ("".join(sorted(succ_list), "".join(sorted(pred_list)) 
``` 

An alternative to using the tuple could be to add a "succ" / "pred" prefix via the `prefix` that is used for the edge attribute 

```
def _digraph_neighborhood_aggregate(G, node, node_labels, edge_attr=None):
    """
    Compute new labels for a given node in a DiGraph by aggregating
    the labels of each successor and predecessor node separately.
    """
    label_list = []
    for succ in G.successor(node):
        prefix =  "_succ_" + "" if edge_attr is None else str(G[node][succ][edge_attr])
        label_list.append(prefix + node_labels[succ_list])

    pred_list = []
    for pred in G.predecessors(node):
        prefix = "_pred_" + "" if edge_attr is None else str(G[node][pred][edge_attr])
        label_list.append(prefix + node_labels[pred])

    return node_labels[node] + "".join(sorted(label_list) 
``` 





It is not clear to me that for DiGraph we need to include both predecessors and successors in the construction. If we cover all the successors for the entire network, then we know where all the predecessors are.  Said another way, specifying all successors of all nodes is sufficient for determining the DiGraph completely. The predecessor information is encoded in the successors. So the hash shouldn't need to explicitly compute the predecessors. 

Also, to unify/simplify the code for `Graph` and `DiGraph` cases, you can use `G[node]` in place of `G.neighbors(node)`. Something like:
```python
def _neighborhood_aggregate(G, node, node_labels, edge_attr=None):
    """
    Compute new labels for a given node by aggregating
    the labels of each node's neighbors.
    """
    label_list = []
    for nbr in G[node]:
        prefix = "" if edge_attr is None else str(G[node][nbr][edge_attr])
        label_list.append(prefix + node_labels[nbr])
    return node_labels[node] + "".join(sorted(label_list)) 
```
I believe this should work for both `DiGraph` and `Graph`.

Re degree for initial node labels: I think using a 2-tuple `(out_degree, in_degree)` works for DiGraph in the same way that `degree` works for Graph. Am I missing something there?  
@dschult I agree with the approach for initial node labels.

With regards to your question about including both predecessors and successors, I believe there are edge-cases for which ignoring one or the other will not work:

![Image](https://github.com/user-attachments/assets/e8aad8f0-9bc8-42a9-97b1-c3d0bbd9fa27)

Note that this is writing from my bachelor's thesis, not from a peer-reviewed paper or something. Of course, if the nodes are initialised using their degree, this _would_ work. However, currently this does not happen for graphs with node and/or edge attributes.

So I think the 'safe' approach would be to consider both successors and predecessors in directed graphs as suggested by @MichaelMartini-Celonis . Alternatively, we could figure out a way that they are considered at least once in initialisation for all types of graphs. I think that should also work and make the algorithm slightly faster, but that approach would be based more on inituition than how WL for digraphs is (to my knowledge) previously described.
I tried to find an example, and at least for the case of **labeled** DiGraphs, it appears that without using both successor and predecessor nodes in the aggregation, we will not discover hash-mismatches for certain non-isomorphic graphs that are "quite simple" after a sufficient number (#Nodes) of aggregations. 

I created two trees, that are not isomorphic given their node labels. Using the current implementation we get the same hash value. 

Using an adapted `init_function` as previously discussed, but not changing the neighborhood aggregation step, we still get the same graph hash: 

```
def _init_node_labels_with_succ_and_pred(G, edge_attr, node_attr):
    if node_attr:
        return {
            u: str(dd[node_attr]) + str((
              "".join(sorted(G.nodes[succ][node_attr] for succ in G.successors(u))), 
			  "".join(sorted(G.nodes[pred][node_attr] for pred in G.predecessors(u))) 
            )) for u, dd in G.nodes(data=True)}

<!--- ignoring other cases ---> 
```

Only using the adapted aggregation that both uses successors and predecessors, do we get a mismatch in hash. 

```
def _digraph_neighborhood_aggregate(G, node, node_labels, edge_attr=None):
    """
    Compute new labels for a given node in a DiGraph by aggregating
    the labels of each successor and predecessor node separately.
    """
    label_list = []
    for succ in G.successors(node):
        prefix =  "_succ_" + "" if edge_attr is None else str(G[node][succ][edge_attr])
        label_list.append(prefix + node_labels[succ_list])

    pred_list = []
    for pred in G.predecessors(node):
        prefix = "_pred_" + "" if edge_attr is None else str(G[node][pred][edge_attr])
        label_list.append(prefix + node_labels[pred])

    return node_labels[node] + "".join(sorted(label_list) 
```

![Image](https://github.com/user-attachments/assets/b3e87df5-d3ed-4009-badf-fa6346ced6b6)

My logic here is that the `blue` nodes in both trees will get a splitting color after 1 iteration using the adapted aggregation function, but not on initialization. Thus using only the adapted initialization function they would remain with the same color. This is because the "information" of the mismatching leaf-node colors (`orange` and `grey`) only reaches the blue nodes with "delay" induced by the `green` nodes. As we can add further layers of `green` nodes we can increase the number of iterations that are needed to pass the information up to the `blue` nodes, as is shown here: 

![Image](https://github.com/user-attachments/assets/600e17ae-73b4-41da-b53a-c00f79afa2e7)
 




Thanks!!   @Anemoonvis for the simple example and @MichaelMartini-Celonis for the nice explanation for why it is needed and how the information filters down through the nodes.

Nice pictures and descriptions!

I think the summary is we need a long term fix involving:
- adjusting the iteration count for the case of no node attribute, 
- initialize DiGraph node labels to be a 2-tuple of (in_degree, out_degree) values,
- aggregate using label identifiers for successors and predecessors when handling a DiGraph.

Do we need to do anything to handle "no edge attribute" better than we currently do?

The current solution for "no edge attribute" seems fine to me. Feels a bit sloppy perhaps to repeatedly just append "empty" prefixes, but I don't expect that would have a big impact on performance and much less on 'correctness'.

The fix seems good to me, though it is important to note that the iteration count only needs to be adjusted if there is no node attribute _and_ no edge attribute.

PS: I think there is no theoretical reason why multigraphs shouldn't work, but for it you need to know the number of edges between nodes. Am I correct that in a multigraph there isn't a very straightforward way to find how many times a node is connected to each neighbour? 
The official way is:  `G.number_of_edges(n, nbr)`

Also for a multgraph, `G.degree` counts how many multiedges are incident to each node. 

The data structure is built on: `G[node][nbr]` is a dictionary keyed by edge-key (any hashable, but by default an integer) indicating which edge, and the value is an attribute dict for that edge.  So to find the number of edges between `u` and `v` you could do `len(G[u][v])`. That's often cleaner if you already have `G[u][v]` for some other reason. `G.number_of_edges(u, v)` is doing that.
Aha, I'd missed that. With that information the implementation could be extended to include multigraphs. That would be fairly straightforward in an undirected multigraph w/o edge attributes (then the number of edges between neighbours can effectively just be used as an edge attribute), but might be a bit more complex if there are edge attributes/directions. 
Maybe it would be better to include that in an issue/fix separate from the one we've discussed here, to keep this fix a bit simpler?

Commit message: only test (un)directed once [skip CI]

Patch to review:
diff --git a/networkx/algorithms/graph_hashing.py b/networkx/algorithms/graph_hashing.py
index 7ded847f057..9f88d012ad1 100644
--- a/networkx/algorithms/graph_hashing.py
+++ b/networkx/algorithms/graph_hashing.py
@@ -22,12 +22,15 @@ def _init_node_labels(G, edge_attr, node_attr):
     elif edge_attr:
         return {u: "" for u in G}
     else:
-        return {u: str(deg) for u, deg in G.degree()}
+        if nx.is_directed(G):
+            return {u: str(G.in_degree(u)) + "_" + str(G.out_degree(u)) for u in G}
+        else:
+            return {u: str(deg) for u, deg in G.degree()}
 
 
-def _neighborhood_aggregate(G, node, node_labels, edge_attr=None):
+def _neighborhood_aggregate_undirected(G, node, node_labels, edge_attr=None):
     """
-    Compute new labels for given node by aggregating
+    Compute new labels for given node in an undirected graph by aggregating
     the labels of each node's neighbors.
     """
     label_list = []
@@ -37,6 +40,27 @@ def _neighborhood_aggregate(G, node, node_labels, edge_attr=None):
     return node_labels[node] + "".join(sorted(label_list))
 
 
+def _neighborhood_aggregate_directed(G, node, node_labels, edge_attr=None):
+    """
+    Compute new labels for given node in a directed graph by aggregating
+    the labels of each node's neighbors.
+    """
+    successor_labels = []
+    for nbr in G.successors(node):
+        prefix = "s_" + "" if edge_attr is None else str(G[node][nbr][edge_attr])
+        successor_labels.append(prefix + node_labels[nbr])
+
+    predecessor_labels = []
+    for nbr in G.predecessors(node):
+        prefix = "p_" + "" if edge_attr is None else str(G[nbr][node][edge_attr])
+        predecessor_labels.append(prefix + node_labels[nbr])
+    return (
+        node_labels[node]
+        + "".join(sorted(successor_labels))
+        + "".join(sorted(predecessor_labels))
+    )
+
+
 @nx.utils.not_implemented_for("multigraph")
 @nx._dispatchable(edge_attrs={"edge_attr": None}, node_attrs="node_attr")
 def weisfeiler_lehman_graph_hash(
@@ -77,6 +101,11 @@ def weisfeiler_lehman_graph_hash(
     h : string
         Hexadecimal string corresponding to hash of the input graph.
 
+    Raises
+    ------
+    ValueError
+        If `iterations` is not a positve number.
+
     Examples
     --------
     Two graphs with edge attributes that are isomorphic, except for
@@ -104,17 +133,17 @@ def weisfeiler_lehman_graph_hash(
     Omitting the `edge_attr` option, results in identical hashes.
 
     >>> nx.weisfeiler_lehman_graph_hash(G1)
-    '7bc4dde9a09d0b94c5097b219891d81a'
+    '6b8542663d7310ddce1732f5f323519a'
     >>> nx.weisfeiler_lehman_graph_hash(G2)
-    '7bc4dde9a09d0b94c5097b219891d81a'
+    '6b8542663d7310ddce1732f5f323519a'
 
     With edge labels, the graphs are no longer assigned
     the same hash digest.
 
     >>> nx.weisfeiler_lehman_graph_hash(G1, edge_attr="label")
-    'c653d85538bcf041d88c011f4f905f10'
+    '4e26b130c32fc56cd84eab0b3d3adebd'
     >>> nx.weisfeiler_lehman_graph_hash(G2, edge_attr="label")
-    '3dcd84af1ca855d0eff3c978d88e7ec7'
+    'b588270cd000e62a38728e4e5317c400'
 
     Notes
     -----
@@ -135,6 +164,11 @@ def weisfeiler_lehman_graph_hash(
     weisfeiler_lehman_subgraph_hashes
     """
 
+    if G.is_directed():
+        _neighborhood_aggregate = _neighborhood_aggregate_directed
+    else:
+        _neighborhood_aggregate = _neighborhood_aggregate_undirected
+
     def weisfeiler_lehman_step(G, labels, edge_attr=None):
         """
         Apply neighborhood aggregation to each node
@@ -147,18 +181,25 @@ def weisfeiler_lehman_step(G, labels, edge_attr=None):
             new_labels[node] = _hash_label(label, digest_size)
         return new_labels
 
+    if iterations <= 0:
+        raise ValueError(
+            "The WL algorithm is only meaningful when run for a positve number of iterations."
+        )
+
     # set initial node labels
     node_labels = _init_node_labels(G, edge_attr, node_attr)
 
-    subgraph_hash_counts = []
+    # If the graph has no attributes, initial labels are the nodes' degrees.
+    # This is equivalent to doing the first iterations of WL.
+    if not edge_attr and not node_attr:
+        iterations -= 1
+
     for _ in range(iterations):
         node_labels = weisfeiler_lehman_step(G, node_labels, edge_attr=edge_attr)
-        counter = Counter(node_labels.values())
-        # sort the counter, extend total counts
-        subgraph_hash_counts.extend(sorted(counter.items(), key=lambda x: x[0]))
 
-    # hash the final counter
-    return _hash_label(str(tuple(subgraph_hash_counts)), digest_size)
+    # hash the final labelling
+    counter = Counter(node_labels.values())
+    return _hash_label(str(sorted(counter.items(), key=lambda x: x[0])), digest_size)
 
 
 @nx.utils.not_implemented_for("multigraph")
@@ -178,7 +219,7 @@ def weisfeiler_lehman_subgraph_hashes(
     Each hash corresponds to a subgraph rooted at a given node u in `G`.
     Lists of subgraph hashes are sorted in increasing order of depth from
     their root node, with the hash at index i corresponding to a subgraph
-    of nodes at most i edges distance from u. Thus, each list will contain
+    of nodes at most i-hops distance from u. Thus, each list will contain
     `iterations` elements - a hash for a subgraph at each depth. If
     `include_initial_labels` is set to `True`, each list will additionally
     have contain a hash of the initial node label (or equivalently a
@@ -241,6 +282,12 @@ def weisfeiler_lehman_subgraph_hashes(
         A dictionary with each key given by a node in G, and each value given
         by the subgraph hashes in order of depth from the key node.
 
+
+    Raises
+    ------
+    ValueError
+        If `iterations` is not a positve number.
+
     Examples
     --------
     Finding similar nodes in different graphs:
@@ -250,25 +297,25 @@ def weisfeiler_lehman_subgraph_hashes(
     >>> G2 = nx.Graph()
     >>> G2.add_edges_from([(1, 3), (2, 3), (1, 6), (1, 5), (4, 6)])
     >>> g1_hashes = nx.weisfeiler_lehman_subgraph_hashes(
-    ...     G1, iterations=3, digest_size=8
+    ...     G1, iterations=4, digest_size=8
     ... )
     >>> g2_hashes = nx.weisfeiler_lehman_subgraph_hashes(
-    ...     G2, iterations=3, digest_size=8
+    ...     G2, iterations=4, digest_size=8
     ... )
 
     Even though G1 and G2 are not isomorphic (they have different numbers of edges),
     the hash sequence of depth 3 for node 1 in G1 and node 5 in G2 are similar:
 
     >>> g1_hashes[1]
-    ['a93b64973cfc8897', 'db1b43ae35a1878f', '57872a7d2059c1c0']
+    ['f6fc42039fba3776', 'a93b64973cfc8897', 'db1b43ae35a1878f', '57872a7d2059c1c0']
     >>> g2_hashes[5]
-    ['a93b64973cfc8897', 'db1b43ae35a1878f', '1716d2a4012fa4bc']
+    ['f6fc42039fba3776', 'a93b64973cfc8897', 'db1b43ae35a1878f', '1716d2a4012fa4bc']
 
-    The first 2 WL subgraph hashes match. From this we can conclude that it's very
-    likely the neighborhood of 2 hops around these nodes are isomorphic.
+    The first 3 WL subgraph hashes match. From this we can conclude that it's very
+    likely the neighborhood of 3 hops around these nodes are isomorphic.
 
-    However the 3-hop neighborhoods of ``G1`` and ``G2`` are not isomorphic since the
-    3rd hashes in the lists above are not equal.
+    However the 4-hop neighborhoods of ``G1`` and ``G2`` are not isomorphic since the
+    4th hashes in the lists above are not equal.
 
     These nodes may be candidates to be classified together since their local topology
     is similar.
@@ -296,6 +343,11 @@ def weisfeiler_lehman_subgraph_hashes(
     weisfeiler_lehman_graph_hash
     """
 
+    if G.is_directed():
+        _neighborhood_aggregate = _neighborhood_aggregate_directed
+    else:
+        _neighborhood_aggregate = _neighborhood_aggregate_undirected
+
     def weisfeiler_lehman_step(G, labels, node_subgraph_hashes, edge_attr=None):
         """
         Apply neighborhood aggregation to each node
@@ -312,7 +364,13 @@ def weisfeiler_lehman_step(G, labels, node_subgraph_hashes, edge_attr=None):
             node_subgraph_hashes[node].append(hashed_label)
         return new_labels
 
+    if iterations <= 0:
+        raise ValueError(
+            "The WL algorithm is only meaningful when run for a positve number of iterations."
+        )
+
     node_labels = _init_node_labels(G, edge_attr, node_attr)
+
     if include_initial_labels:
         node_subgraph_hashes = {
             k: [_hash_label(v, digest_size)] for k, v in node_labels.items()
@@ -320,6 +378,14 @@ def weisfeiler_lehman_step(G, labels, node_subgraph_hashes, edge_attr=None):
     else:
         node_subgraph_hashes = defaultdict(list)
 
+    # If the graph has no attributes, initial labels are the nodes' degrees.
+    # This is equivalent to doing the first iterations of WL.
+    if not edge_attr and not node_attr:
+        iterations -= 1
+        for node in G.nodes():
+            hashed_label = _hash_label(node_labels[node], digest_size)
+            node_subgraph_hashes[node].append(hashed_label)
+
     for _ in range(iterations):
         node_labels = weisfeiler_lehman_step(
             G, node_labels, node_subgraph_hashes, edge_attr
diff --git a/networkx/algorithms/tests/test_graph_hashing.py b/networkx/algorithms/tests/test_graph_hashing.py
index 0828069d1c3..a8128080c67 100644
--- a/networkx/algorithms/tests/test_graph_hashing.py
+++ b/networkx/algorithms/tests/test_graph_hashing.py
@@ -1,7 +1,35 @@
+import copy
+
 import pytest
 
 import networkx as nx
-from networkx.generators import directed
+
+# Unit tests relevant for both functions in this module.
+
+
+def test_positive_iters():
+    G1 = nx.empty_graph()
+    with pytest.raises(
+        ValueError,
+        match="The WL algorithm is only meaningful when run for a positve number of iterations.",
+    ):
+        nx.weisfeiler_lehman_graph_hash(G1, iterations=-3)
+    with pytest.raises(
+        ValueError,
+        match="The WL algorithm is only meaningful when run for a positve number of iterations.",
+    ):
+        nx.weisfeiler_lehman_subgraph_hashes(G1, iterations=-3)
+    with pytest.raises(
+        ValueError,
+        match="The WL algorithm is only meaningful when run for a positve number of iterations.",
+    ):
+        nx.weisfeiler_lehman_graph_hash(G1, iterations=0)
+    with pytest.raises(
+        ValueError,
+        match="The WL algorithm is only meaningful when run for a positve number of iterations.",
+    ):
+        nx.weisfeiler_lehman_subgraph_hashes(G1, iterations=0)
+
 
 # Unit tests for the :func:`~networkx.weisfeiler_lehman_graph_hash` function
 
@@ -48,8 +76,8 @@ def test_directed():
 def test_reversed():
     """
     A directed graph with no bi-directional edges should yield different a graph hash
-    to the same graph taken with edge directions reversed if there are no hash collisions.
-    Here we test a cycle graph which is the minimal counterexample
+    to the same graph taken with edge directions reversed if there are no hash
+    collisions. Here we test a cycle graph which is the minimal counterexample
     """
     G = nx.cycle_graph(5, create_using=nx.DiGraph)
     nx.set_node_attributes(G, {n: str(n) for n in G.nodes()}, name="label")
@@ -250,7 +278,161 @@ def test_digest_size():
         assert len(h32) == 32 * 2
 
 
-# Unit tests for the :func:`~networkx.weisfeiler_lehman_hash_subgraphs` function
+def test_directed_bugs():
+    """
+    These were bugs for directed graphs as discussed in issue #7806
+    """
+    Ga = nx.DiGraph()
+    Gb = nx.DiGraph()
+    Ga.add_nodes_from([1, 2, 3, 4])
+    Gb.add_nodes_from([1, 2, 3, 4])
+    Ga.add_edges_from([(1, 2), (3, 2)])
+    Gb.add_edges_from([(1, 2), (3, 4)])
+    Ga_hash = nx.weisfeiler_lehman_graph_hash(Ga)
+    Gb_hash = nx.weisfeiler_lehman_graph_hash(Gb)
+    assert Ga_hash != Gb_hash
+
+    Tree1 = nx.DiGraph()
+    Tree1.add_edges_from([(0, 4), (1, 5), (2, 6), (3, 7)])
+    Tree1.add_edges_from([(4, 8), (5, 8), (6, 9), (7, 9)])
+    Tree1.add_edges_from([(8, 10), (9, 10)])
+    nx.set_node_attributes(
+        Tree1, {10: "s", 8: "a", 9: "a", 4: "b", 5: "b", 6: "b", 7: "b"}, "weight"
+    )
+    Tree2 = copy.deepcopy(Tree1)
+    nx.set_node_attributes(Tree1, {0: "d", 1: "c", 2: "d", 3: "c"}, "weight")
+    nx.set_node_attributes(Tree2, {0: "d", 1: "d", 2: "c", 3: "c"}, "weight")
+    Tree1_hash_short = nx.weisfeiler_lehman_graph_hash(
+        Tree1, iterations=1, node_attr="weight"
+    )
+    Tree2_hash_short = nx.weisfeiler_lehman_graph_hash(
+        Tree2, iterations=1, node_attr="weight"
+    )
+    assert Tree1_hash_short == Tree2_hash_short
+    Tree1_hash = nx.weisfeiler_lehman_graph_hash(
+        Tree1, node_attr="weight"
+    )  # Default is 3 iterations
+    Tree2_hash = nx.weisfeiler_lehman_graph_hash(Tree2, node_attr="weight")
+    assert Tree1_hash != Tree2_hash
+
+
+def test_trivial_labels_isomorphism():
+    """
+    Trivial labelling of the graph should not change isomorphism verdicts.
+    """
+    n, r = 100, 10
+    p = 1.0 / r
+    for i in range(1, r + 1):
+        G1 = nx.erdos_renyi_graph(n, p * i, seed=500 + i)
+        G2 = nx.erdos_renyi_graph(n, p * i, seed=42 + i)
+        G1_hash = nx.weisfeiler_lehman_graph_hash(G1)
+        G2_hash = nx.weisfeiler_lehman_graph_hash(G2)
+        equal = G1_hash == G2_hash
+
+        nx.set_node_attributes(G1, values=1, name="weight")
+        nx.set_node_attributes(G2, values=1, name="weight")
+        G1_hash_node = nx.weisfeiler_lehman_graph_hash(G1, node_attr="weight")
+        G2_hash_node = nx.weisfeiler_lehman_graph_hash(G2, node_attr="weight")
+        equal_node = G1_hash_node == G2_hash_node
+
+        nx.set_edge_attributes(G1, values="a", name="e_weight")
+        nx.set_edge_attributes(G2, values="a", name="e_weight")
+        G1_hash_edge = nx.weisfeiler_lehman_graph_hash(G1, edge_attr="e_weight")
+        G2_hash_edge = nx.weisfeiler_lehman_graph_hash(G2, edge_attr="e_weight")
+        equal_edge = G1_hash_edge == G2_hash_edge
+
+        G1_hash_both = nx.weisfeiler_lehman_graph_hash(
+            G1, edge_attr="e_weight", node_attr="weight"
+        )
+        G2_hash_both = nx.weisfeiler_lehman_graph_hash(
+            G2, edge_attr="e_weight", node_attr="weight"
+        )
+        equal_both = G1_hash_both == G2_hash_both
+
+        assert equal == equal_node
+        assert equal_node == equal_edge
+        assert equal_edge == equal_both
+
+
+def test_trivial_labels_isomorphism_directed():
+    """
+    Trivial labelling of the graph should not change isomorphism verdicts on digraphs.
+    """
+    n, r = 100, 10
+    p = 1.0 / r
+    for i in range(1, r + 1):
+        G1 = nx.erdos_renyi_graph(n, p * i, directed=True, seed=500 + i)
+        G2 = nx.erdos_renyi_graph(n, p * i, directed=True, seed=42 + i)
+        G1_hash = nx.weisfeiler_lehman_graph_hash(G1)
+        G2_hash = nx.weisfeiler_lehman_graph_hash(G2)
+        equal = G1_hash == G2_hash
+
+        nx.set_node_attributes(G1, values=1, name="weight")
+        nx.set_node_attributes(G2, values=1, name="weight")
+        G1_hash_node = nx.weisfeiler_lehman_graph_hash(G1, node_attr="weight")
+        G2_hash_node = nx.weisfeiler_lehman_graph_hash(G2, node_attr="weight")
+        equal_node = G1_hash_node == G2_hash_node
+
+        nx.set_edge_attributes(G1, values="a", name="e_weight")
+        nx.set_edge_attributes(G2, values="a", name="e_weight")
+        G1_hash_edge = nx.weisfeiler_lehman_graph_hash(G1, edge_attr="e_weight")
+        G2_hash_edge = nx.weisfeiler_lehman_graph_hash(G2, edge_attr="e_weight")
+        equal_edge = G1_hash_edge == G2_hash_edge
+
+        G1_hash_both = nx.weisfeiler_lehman_graph_hash(
+            G1, edge_attr="e_weight", node_attr="weight"
+        )
+        G2_hash_both = nx.weisfeiler_lehman_graph_hash(
+            G2, edge_attr="e_weight", node_attr="weight"
+        )
+        equal_both = G1_hash_both == G2_hash_both
+
+        assert equal == equal_node
+        assert equal_node == equal_edge
+        assert equal_edge == equal_both
+
+    # Specific case that was found to be a bug in issue #7806
+    # Without weights worked
+    Ga = nx.DiGraph()
+    Ga.add_nodes_from([1, 2, 3, 4])
+    Gb = copy.deepcopy(Ga)
+    Ga.add_edges_from([(1, 2), (3, 2)])
+    Gb.add_edges_from([(1, 2), (3, 4)])
+    Ga_hash = nx.weisfeiler_lehman_graph_hash(Ga)
+    Gb_hash = nx.weisfeiler_lehman_graph_hash(Gb)
+    assert Ga_hash != Gb_hash
+
+    # Now with trivial weights
+    nx.set_node_attributes(Ga, values=1, name="weight")
+    nx.set_node_attributes(Gb, values=1, name="weight")
+    Ga_hash = nx.weisfeiler_lehman_graph_hash(Ga, node_attr="weight")
+    Gb_hash = nx.weisfeiler_lehman_graph_hash(Gb, node_attr="weight")
+    assert Ga_hash != Gb_hash
+
+
+def test_trivial_labels_hashes():
+    """
+    Test that 'empty' labelling of nodes or edges shouldn't have a different impact
+    on the calculated hash. Note that we cannot assume it trivial weights have no
+    impact at all. Without (trivial) weights, a node will start with hashing its
+    degree. This step isomitted when there áre weights.
+    """
+    n, r = 100, 10
+    p = 1.0 / r
+    for i in range(1, r + 1):
+        G1 = nx.erdos_renyi_graph(n, p * i, seed=500 + i)
+        nx.set_node_attributes(G1, values="", name="weight")
+        first = nx.weisfeiler_lehman_graph_hash(G1, node_attr="weight")
+        nx.set_edge_attributes(G1, values="", name="e_weight")
+        second = nx.weisfeiler_lehman_graph_hash(G1, edge_attr="e_weight")
+        assert first == second
+        third = nx.weisfeiler_lehman_graph_hash(
+            G1, edge_attr="e_weight", node_attr="weight"
+        )
+        assert second == third
+
+
+# Unit tests for the :func:`~networkx.weisfeiler_lehman_subgraph_hashes` function
 
 
 def is_subiteration(a, b):
@@ -263,9 +445,10 @@ def is_subiteration(a, b):
 
 def hexdigest_sizes_correct(a, digest_size):
     """
-    returns True if all hex digest sizes are the expected length in a node:subgraph-hashes
-    dictionary. Hex digest string length == 2 * bytes digest length since each pair of hex
-    digits encodes 1 byte (https://docs.python.org/3/library/hashlib.html)
+    returns True if all hex digest sizes are the expected length in a
+    node:subgraph-hashes dictionary. Hex digest string length == 2 * bytes digest
+    length since each pair of hex digits encodes 1 byte
+    (https://docs.python.org/3/library/hashlib.html)
     """
     hexdigest_size = digest_size * 2
     list_digest_sizes_correct = lambda l: all(len(x) == hexdigest_size for x in l)
@@ -293,8 +476,8 @@ def test_empty_graph_subgraph_hash():
 
 def test_directed_subgraph_hash():
     """
-    A directed graph with no bi-directional edges should yield different subgraph hashes
-    to the same graph taken as undirected, if all hashes don't collide.
+    A directed graph with no bi-directional edges should yield different subgraph
+    hashes to the same graph taken as undirected, if all hashes don't collide.
     """
     r = 10
     for i in range(r):
@@ -309,9 +492,9 @@ def test_directed_subgraph_hash():
 
 def test_reversed_subgraph_hash():
     """
-    A directed graph with no bi-directional edges should yield different subgraph hashes
-    to the same graph taken with edge directions reversed if there are no hash collisions.
-    Here we test a cycle graph which is the minimal counterexample
+    A directed graph with no bi-directional edges should yield different subgraph
+    hashes to the same graph taken with edge directions reversed if there are no
+    hash collisions. Here we test a cycle graph which is the minimal counterexample
     """
     G = nx.cycle_graph(5, create_using=nx.DiGraph)
     nx.set_node_attributes(G, {n: str(n) for n in G.nodes()}, name="label")
@@ -326,8 +509,8 @@ def test_reversed_subgraph_hash():
 
 def test_isomorphic_subgraph_hash():
     """
-    the subgraph hashes should be invariant to node-relabeling when the output is reindexed
-    by the same mapping and all hashes don't collide.
+    the subgraph hashes should be invariant to node-relabeling when the output is
+    reindexed by the same mapping and all hashes don't collide.
     """
     n, r = 100, 10
     p = 1.0 / r
@@ -513,7 +696,7 @@ def test_isomorphic_edge_attr_and_node_attr_subgraph_hash():
 def test_iteration_depth():
     """
     All nodes should have the correct number of subgraph hashes in the output when
-    using degree as initial node labels
+    using degree as initial node labels.
     Subsequent iteration depths for the same graph should be additive for each node
     """
     n, r = 100, 10


Code review comments:
diff_hunk: @@ -104,17 +133,17 @@ def weisfeiler_lehman_graph_hash(
     Omitting the `edge_attr` option, results in identical hashes.
 
     >>> nx.weisfeiler_lehman_graph_hash(G1)
-    '7bc4dde9a09d0b94c5097b219891d81a'
+    '6b8542663d7310ddce1732f5f323519a'
     >>> nx.weisfeiler_lehman_graph_hash(G2)
-    '7bc4dde9a09d0b94c5097b219891d81a'
+    '6b8542663d7310ddce1732f5f323519a'
 
     With edge labels, the graphs are no longer assigned
     the same hash digest.
 
     >>> nx.weisfeiler_lehman_graph_hash(G1, edge_attr="label")
-    'c653d85538bcf041d88c011f4f905f10'
+    '4e26b130c32fc56cd84eab0b3d3adebd'
path: networkx/algorithms/graph_hashing.py
text: @user1:
Here we use attributes and we don't specify an iteration count. And G1 is not directed.
Why is the hash changing?  Maybe I am not understanding something.... :(

@author:
Your confusion is understandable. Something that we did not discuss in the issue, but I noticed (again) when writing the pull request, is that the previous implementation did something a bit weird in constructing the final hash. It counts *all* hashes/labels that are made in the full iteration process and hashes those to obtain the final graph hash. This is not in line with how the Weisfeiler-Leman algorithm is defined. In WL, we (normally) only consider the final labeling of the graphs when comparing them. 

When I first read the NetworkX implementation, I was pretty confused about this detail. Additionally, it is of course inefficient, because we're storing/hashing (many) more labels than necessary. 
For that reason, I updated that step in  [one of the commits](https://github.com/networkx/networkx/pull/7834/commits/35025e8ad326d5a04dc946b5877a1b721766d439) in this pull request. I mentioned that shortly in the description of the pull request and again later in this discussion about continuity. But because of the brevity in which I did that and the fact that we didn't discuss it in the original issue I can imagine that wasn't clear. Apologies for that.

An important note is that the old approach is not wrong, in the sense that it will never produce wrong isomorphism verdicts (barring any bugs we've missed). However, it is atypical (and thus potentially a bit confusing to users who know the algorithm) and of course inefficient. For those reasons, and because there are already significant continuity breaks due to this bugfix, I would propose changing it now. Potentially doing so later would only lead to more continuity breaks. Of course, then I need to document these changes better than I've done so far. On the other hand, because the old approach already worked, you can also decide you prefer to keep as much continuity as possible. Then I'll revert that commit and only document the changes to directed graphs/graphs with attributes. 

Apologies for being a bit opaque, and thanks a lot for being so careful in the review process :)