Repository: networkx/networkx

Issue:
greedy modularity communities fails
I am trying to find communities of a directed graph (Number of Nodes: 53663 and Number of Edges: 953380) with `greedy_modularity_communities`. 

This is the following error I receive:

```
 line 127, in greedy_modularity_communities
    if dq_heap[j].h[0] == (-dq, j, i):
IndexError: list index out of range
```

However, running greedy_modularity_communities on the same graph with undirected edges is successful.

Discussion:
elplatt: CNM is Clauset-Newman-Moore, the authors of the algorithm.
ThanasisMattas: Thanks @dschult and @elplatt ! Indeed it is not so clear now that I moved it to a new function, so I made it more verbose.

As a last comment, although it seemed that much are going on, I tried to simplify it a bit, but everything seems to be in place.

Patch:
diff --git a/networkx/algorithms/centrality/betweenness.py b/networkx/algorithms/centrality/betweenness.py
index dd478dcd2..aab55367f 100644
--- a/networkx/algorithms/centrality/betweenness.py
+++ b/networkx/algorithms/centrality/betweenness.py
@@ -5,12 +5,12 @@ import warnings
 
 from networkx.utils import py_random_state
 from networkx.utils.decorators import not_implemented_for
+from networkx.algorithms.shortest_paths.weighted import _weight_function
 
 __all__ = ["betweenness_centrality", "edge_betweenness_centrality", "edge_betweenness"]
 
 
 @py_random_state(5)
-@not_implemented_for("multigraph")
 def betweenness_centrality(
     G, k=None, normalized=True, weight=None, endpoints=False, seed=None
 ):
@@ -236,6 +236,8 @@ def edge_betweenness_centrality(G, k=None, normalized=True, weight=None, seed=No
     betweenness = _rescale_e(
         betweenness, len(G), normalized=normalized, directed=G.is_directed()
     )
+    if G.is_multigraph():
+        betweenness = _add_edge_keys(G, betweenness, weight=weight)
     return betweenness
 
 
@@ -276,6 +278,7 @@ def _single_source_shortest_path_basic(G, s):
 
 
 def _single_source_dijkstra_path_basic(G, s, weight):
+    weight = _weight_function(G, weight)
     # modified from Eppstein
     S = []
     P = {}
@@ -298,7 +301,7 @@ def _single_source_dijkstra_path_basic(G, s, weight):
         S.append(v)
         D[v] = dist
         for w, edgedata in G[v].items():
-            vw_dist = dist + edgedata.get(weight, 1)
+            vw_dist = dist + weight(v, w, edgedata)
             if w not in D and (w not in seen or vw_dist < seen[w]):
                 seen[w] = vw_dist
                 push(Q, (vw_dist, next(c), v, w))
@@ -394,3 +397,39 @@ def _rescale_e(betweenness, n, normalized, directed=False, k=None):
         for v in betweenness:
             betweenness[v] *= scale
     return betweenness
+
+
+@not_implemented_for("graph")
+def _add_edge_keys(G, betweenness, weight=None):
+    r"""Adds the corrected betweenness centrality (BC) values for multigraphs.
+
+    Parameters
+    ----------
+    G : NetworkX graph.
+
+    betweenness : dictionary
+        Dictionary mapping adjacent node tuples to betweenness centrality values.
+
+    weight : string or function
+        See `_weight_function´ for details. Defaults to `None`.
+
+    Returns
+    -------
+    edges : dictionary
+        The parameter `betweenness` including edges with keys and their
+        betweenness centrality values.
+
+    The BC value is divided among edges of equal weight.
+    """
+    _weight = _weight_function(G, weight)
+
+    edge_bc = dict.fromkeys(G.edges, 0.0)
+    for u, v in betweenness:
+        d = G[u][v]
+        wt = _weight(u, v, d)
+        keys = [k for k in d if _weight(u, v, {k: d[k]}) == wt]
+        bc = betweenness[(u, v)] / len(keys)
+        for k in keys:
+            edge_bc[(u, v, k)] = bc
+
+    return edge_bc
diff --git a/networkx/algorithms/centrality/betweenness_subset.py b/networkx/algorithms/centrality/betweenness_subset.py
index ed073855c..473b42ca6 100644
--- a/networkx/algorithms/centrality/betweenness_subset.py
+++ b/networkx/algorithms/centrality/betweenness_subset.py
@@ -3,9 +3,8 @@ import warnings
 
 from networkx.algorithms.centrality.betweenness import (
     _single_source_dijkstra_path_basic as dijkstra,
-)
-from networkx.algorithms.centrality.betweenness import (
     _single_source_shortest_path_basic as shortest_path,
+    _add_edge_keys,
 )
 
 __all__ = [
@@ -193,6 +192,8 @@ def edge_betweenness_centrality_subset(
     for n in G:  # remove nodes to only return edges
         del b[n]
     b = _rescale_e(b, len(G), normalized=normalized, directed=G.is_directed())
+    if G.is_multigraph():
+        b = _add_edge_keys(G, b, weight=weight)
     return b
 
 
diff --git a/networkx/algorithms/community/modularity_max.py b/networkx/algorithms/community/modularity_max.py
index debf5fb3a..10c5d01ce 100644
--- a/networkx/algorithms/community/modularity_max.py
+++ b/networkx/algorithms/community/modularity_max.py
@@ -1,8 +1,11 @@
 """Functions for detecting communities based on modularity."""
 
-from networkx.algorithms.community.quality import modularity
+from collections import Counter
 
+import networkx as nx
+from networkx.algorithms.community.quality import modularity
 from networkx.utils.mapped_queue import MappedQueue
+from networkx.utils import not_implemented_for
 
 __all__ = [
     "greedy_modularity_communities",
@@ -11,11 +14,114 @@ __all__ = [
 ]
 
 
+def _greedy_modularity_communities_init(G, weight=None, resolution=1):
+    r"""Initializes the data structures for greedy_modularity_communities().
+
+    Clauset-Newman-Moore Eq 8-9. Eq 8 was missing a factor of 2 (from A_ij + A_ji).
+    See [2]_ at :func:`greedy_modularity_communities`.
+
+    Parameters
+    ----------
+    G : NetworkX graph
+
+    weight : string or None, optional (default=None)
+        The name of an edge attribute that holds the numerical value used
+        as a weight.  If None, then each edge has weight 1.
+        The degree is the sum of the edge weights adjacent to the node.
+
+    resolution : float (default=1)
+        If resolution is less than 1, modularity favors larger communities.
+        Greater than 1 favors smaller communities.
+
+    Returns
+    -------
+    dq_dict : dict of dict's
+        dq_dict[i][j]: dQ for merging community i, j
+
+    dq_heap : dict of MappedQueue's
+        dq_heap[i][n] : (-dq, i, j) for communitiy i nth largest dQ
+
+    H : MappedQueue
+        (-dq, i, j) for community with nth largest max_j(dQ_ij)
+
+    a, b : dict
+        undirected:
+            a[i]: fraction of (total weight of) edges within community i
+            b : None
+        directed:
+            a[i]: fraction of (total weight of) edges with tails within community i
+            b[i]: fraction of (total weight of) edges with heads within community i
+
+    See Also
+    --------
+    :func:`greedy_modularity_communities`
+    :func:`~networkx.algorithms.community.quality.modularity`
+    """
+    # Count nodes and edges (or the sum of edge-weights for weighted graphs)
+    N = G.number_of_nodes()
+    m = G.size(weight)
+
+    # Calculate degrees
+    if G.is_directed():
+        k_in = dict(G.in_degree(weight=weight))
+        k_out = dict(G.out_degree(weight=weight))
+        q0 = 1.0 / m
+    else:
+        k_in = k_out = dict(G.degree(weight=weight))
+        q0 = 1.0 / (2.0 * m)
+
+    a = {node: kout * q0 for node, kout in k_out.items()}
+    if G.is_directed():
+        b = {node: kin * q0 for node, kin in k_in.items()}
+    else:
+        b = None
+
+    dq_dict = {
+        i: {
+            j: q0
+            * (
+                G.get_edge_data(i, j, default={weight: 0}).get(weight, 1.0)
+                + G.get_edge_data(j, i, default={weight: 0}).get(weight, 1.0)
+                - resolution * q0 * (k_out[i] * k_in[j] + k_in[i] * k_out[j])
+            )
+            for j in nx.all_neighbors(G, i)
+            if j != i
+        }
+        for i in G.nodes()
+    }
+
+    # dq correction for multi-edges
+    # In case of multi-edges, get_edge_data(i, j) returns the key: data dict of the i, j
+    # edges, which does not have a 'weight' key. Therefore, when calculating dq for i, j
+    # Aij is always 1.0 and a correction is required.
+    if G.is_multigraph():
+        edges_count = dict(Counter(G.edges()))
+        multi_edges = [edge for edge, count in edges_count.items() if count > 1]
+        for edge in multi_edges:
+            total_wt = sum(d.get(weight, 1) for d in G.get_edge_data(*edge).values())
+            if G.is_directed():
+                # The correction applies only to the direction of the edge. The edge at
+                # the other direction is either not a multiedge (where the weight is
+                # added correctly), non-existent or it is also a multiedge, in which
+                # case it will be handled singly when its turn in the loop comes.
+                q00 = q0
+            else:
+                q00 = 2 * q0
+            dq_dict[edge[0]][edge[1]] += q00 * (total_wt - 1)
+            dq_dict[edge[1]][edge[0]] += q00 * (total_wt - 1)
+
+    dq_heap = {
+        i: MappedQueue([(-dq, i, j) for j, dq in dq_dict[i].items()]) for i in G.nodes()
+    }
+    H = MappedQueue([dq_heap[i].h[0] for i in G.nodes() if len(dq_heap[i]) > 0])
+
+    return dq_dict, dq_heap, H, a, b
+
+
 def greedy_modularity_communities(G, weight=None, resolution=1, n_communities=1):
     r"""Find communities in G using greedy modularity maximization.
 
     This function uses Clauset-Newman-Moore greedy modularity maximization [2]_.
-    This method currently supports the Graph class.
 
     Greedy modularity maximization begins with each node in its own community
     and joins the pair of communities that most increases modularity until no
@@ -48,8 +154,8 @@ def greedy_modularity_communities(G, weight=None, resolution=1, n_communities=1)
 
     Returns
     -------
-    list
-        A list of sets of nodes, one for each community.
+    partition: list
+        A list of frozensets of nodes, one for each community.
         Sorted by length with largest communities first.
 
     Examples
@@ -66,63 +172,30 @@ def greedy_modularity_communities(G, weight=None, resolution=1, n_communities=1)
 
     References
     ----------
-    .. [1] M. E. J Newman "Networks: An Introduction", page 224
+    .. [1] Newman, M. E. J. "Networks: An Introduction", page 224
        Oxford University Press 2011.
     .. [2] Clauset, A., Newman, M. E., & Moore, C.
        "Finding community structure in very large networks."
        Physical Review E 70(6), 2004.
     .. [3] Reichardt and Bornholdt "Statistical Mechanics of Community
        Detection" Phys. Rev. E74, 2006.
+    .. [4] Newman, M. E. J."Analysis of weighted networks"
+       Physical Review E 70(5 Pt 2):056131, 2004.
     """
-
-    # Count nodes and edges
-    N = len(G.nodes())
-    m = sum([d.get(weight, 1) for u, v, d in G.edges(data=True)])
-    q0 = 1.0 / (2.0 * m)
-
+    N = G.number_of_nodes()
     if (n_communities < 1) or (n_communities > N):
         raise ValueError(
-            f"n_communities must be between 1 and {len(G.nodes())}. Got {n_communities}"
+            f"n_communities must be between 1 and {N}. Got {n_communities}"
         )
 
-    # Map node labels to contiguous integers
-    label_for_node = {i: v for i, v in enumerate(G.nodes())}
-    node_for_label = {label_for_node[i]: i for i in range(N)}
-
-    # Calculate degrees
-    k_for_label = G.degree(G.nodes(), weight=weight)
-    k = [k_for_label[label_for_node[i]] for i in range(N)]
-
-    # Initialize community and merge lists
-    communities = {i: frozenset([i]) for i in range(N)}
-    merges = []
-
-    # Initial modularity
-    partition = [[label_for_node[x] for x in c] for c in communities.values()]
-    q_cnm = modularity(G, partition, resolution=resolution)
-
     # Initialize data structures
-    # CNM Eq 8-9 (Eq 8 was missing a factor of 2 (from A_ij + A_ji)
-    # a[i]: fraction of edges within community i
-    # dq_dict[i][j]: dQ for merging community i, j
-    # dq_heap[i][n] : (-dq, i, j) for communitiy i nth largest dQ
-    # H[n]: (-dq, i, j) for community with nth largest max_j(dQ_ij)
-    a = [k[i] * q0 for i in range(N)]
-    dq_dict = {
-        i: {
-            j: 2
-            * q0
-            * G.get_edge_data(label_for_node[i], label_for_node[j]).get(weight, 1.0)
-            - 2 * resolution * k[i] * k[j] * q0 * q0
-            for j in [node_for_label[u] for u in G.neighbors(label_for_node[i])]
-            if j != i
-        }
-        for i in range(N)
-    }
-    dq_heap = [
-        MappedQueue([(-dq, i, j) for j, dq in dq_dict[i].items()]) for i in range(N)
-    ]
-    H = MappedQueue([dq_heap[i].h[0] for i in range(N) if len(dq_heap[i]) > 0])
+    dq_dict, dq_heap, H, a, b = _greedy_modularity_communities_init(
+        G, weight, resolution
+    )
+    # Initialize single-node communities
+    communities = {i: frozenset([i]) for i in G.nodes()}
+    # Initial modularity
+    q_cnm = modularity(G, communities.values(), resolution=resolution)
 
     # Merge communities until we can't improve modularity or until desired number of
     # communities (n_communities) is reached.
@@ -159,7 +232,6 @@ def greedy_modularity_communities(G, weight=None, resolution=1, n_communities=1)
         # Perform merge
         communities[j] = frozenset(communities[i] | communities[j])
         del communities[i]
-        merges.append((i, j, dq))
         # New modularity
         q_cnm += dq
         # Get list of communities connected to merged communities
@@ -173,10 +245,16 @@ def greedy_modularity_communities(G, weight=None, resolution=1, n_communities=1)
             if k in both_set:
                 dq_jk = dq_dict[j][k] + dq_dict[i][k]
             elif k in j_set:
-                dq_jk = dq_dict[j][k] - 2.0 * resolution * a[i] * a[k]
+                if G.is_directed():
+                    dq_jk = dq_dict[j][k] - resolution * (a[i] * b[k] + a[k] * b[i])
+                else:
+                    dq_jk = dq_dict[j][k] - 2.0 * resolution * a[i] * a[k]
             else:
                 # k in i_set
-                dq_jk = dq_dict[i][k] - 2.0 * resolution * a[j] * a[k]
+                if G.is_directed():
+                    dq_jk = dq_dict[i][k] - resolution * (a[j] * b[k] + a[k] * b[j])
+                else:
+                    dq_jk = dq_dict[i][k] - 2.0 * resolution * a[j] * a[k]
             # Update rows j and k
             for row, col in [(j, k), (k, j)]:
                 # Save old value for finding heap index
@@ -237,13 +315,16 @@ def greedy_modularity_communities(G, weight=None, resolution=1, n_communities=1)
         # Merge i into j and update a
         a[j] += a[i]
         a[i] = 0
+        if G.is_directed():
+            b[j] += b[i]
+            b[i] = 0
 
-    communities = [
-        frozenset([label_for_node[i] for i in c]) for c in communities.values()
-    ]
-    return sorted(communities, key=len, reverse=True)
+    partition = sorted(communities.values(), key=len, reverse=True)
+    return partition
 
 
+@not_implemented_for("directed")
+@not_implemented_for("multigraph")
 def naive_greedy_modularity_communities(G, resolution=1):
     r"""Find communities in G using greedy modularity maximization.
 
diff --git a/networkx/algorithms/community/quality.py b/networkx/algorithms/community/quality.py
index 9f374435f..ce79727bc 100644
--- a/networkx/algorithms/community/quality.py
+++ b/networkx/algorithms/community/quality.py
@@ -283,9 +283,9 @@ def modularity(G, communities, weight="weight", resolution=1):
         These node sets must represent a partition of G's nodes.
 
     weight : string or None, optional (default="weight")
-            The edge attribute that holds the numerical value used
-            as a weight. If None or an edge does not have that attribute,
-            then that edge has weight 1.
+        The edge attribute that holds the numerical value used
+        as a weight. If None or an edge does not have that attribute,
+        then that edge has weight 1.
 
     resolution : float (default=1)
         If resolution is less than 1, modularity favors larger communities.
diff --git a/networkx/algorithms/connectivity/edge_augmentation.py b/networkx/algorithms/connectivity/edge_augmentation.py
index dc8b86075..2939d3bc2 100644
--- a/networkx/algorithms/connectivity/edge_augmentation.py
+++ b/networkx/algorithms/connectivity/edge_augmentation.py
@@ -1124,15 +1124,16 @@ def complement_edges(G):
     >>> sorted(complement_edges(G))
     []
     """
+    G_adj = G._adj  # Store as a variable to eliminate attribute lookup
     if G.is_directed():
         for u, v in it.combinations(G.nodes(), 2):
-            if v not in G.adj[u]:
+            if v not in G_adj[u]:
                 yield (u, v)
-            if u not in G.adj[v]:
+            if u not in G_adj[v]:
                 yield (v, u)
     else:
         for u, v in it.combinations(G.nodes(), 2):
-            if v not in G.adj[u]:
+            if v not in G_adj[u]:
                 yield (u, v)
 
 
diff --git a/networkx/algorithms/link_prediction.py b/networkx/algorithms/link_prediction.py
index 5922360dc..8c065f492 100644
--- a/networkx/algorithms/link_prediction.py
+++ b/networkx/algorithms/link_prediction.py
@@ -246,7 +246,8 @@ def common_neighbor_centrality(G, ebunch=None, alpha=0.8):
         Default value: None.
 
     alpha : Parameter defined for participation of Common Neighbor
-            and Centrality Algorithm share. Default value set to 0.8
+            and Centrality Algorithm share. Values for alpha should
+            normally be between 0 and 1. Default value set to 0.8
             because author found better performance at 0.8 for all the
             dataset.
             Default value: 0.8
@@ -275,11 +276,17 @@ def common_neighbor_centrality(G, ebunch=None, alpha=0.8):
            Sci Rep 10, 364 (2020).
            https://doi.org/10.1038/s41598-019-57304-y
     """
-    shortest_path = nx.shortest_path(G)
+
+    spl = dict(nx.shortest_path_length(G))
+    inf = float("inf")
 
     def predict(u, v):
-        return alpha * len(list(nx.common_neighbors(G, u, v))) + (1 - alpha) * (
-            G.number_of_nodes() / (len(shortest_path[u][v]) - 1)
+        if u == v:
+            raise nx.NetworkXAlgorithmError("Self links are not supported")
+        path_len = spl[u].get(v, inf)
+
+        return alpha * sum(1 for _ in nx.common_neighbors(G, u, v)) + (1 - alpha) * (
+            G.number_of_nodes() / path_len
         )
 
     return _apply_prediction(G, predict, ebunch)
diff --git a/networkx/drawing/nx_pylab.py b/networkx/drawing/nx_pylab.py
index e66d0c3a4..bbffab224 100644
--- a/networkx/drawing/nx_pylab.py
+++ b/networkx/drawing/nx_pylab.py
@@ -852,10 +852,10 @@ def draw_networkx_edges(
     # Draw the edges
     if use_linecollection:
         edge_viz_obj = _draw_networkx_edges_line_collection()
-        # Make sure selfloop edges are also drawn.
-        edgelist = list(nx.selfloop_edges(G))
-        if edgelist:
-            edge_pos = np.asarray([(pos[e[0]], pos[e[1]]) for e in edgelist])
+        # Make sure selfloop edges are also drawn
+        selfloops_to_draw = [loop for loop in nx.selfloop_edges(G) if loop in edgelist]
+        if selfloops_to_draw:
+            edge_pos = np.asarray([(pos[e[0]], pos[e[1]]) for e in selfloops_to_draw])
             arrowstyle = "-"
             _draw_networkx_edges_fancy_arrow_patch()
     else:


Test patch:
diff --git a/networkx/algorithms/approximation/tests/test_traveling_salesman.py b/networkx/algorithms/approximation/tests/test_traveling_salesman.py
index db86d7b38..5b8ff39ec 100644
--- a/networkx/algorithms/approximation/tests/test_traveling_salesman.py
+++ b/networkx/algorithms/approximation/tests/test_traveling_salesman.py
@@ -732,6 +732,7 @@ def test_sample_spanning_tree():
     assert nx.utils.edges_equal(solution.edges, sampled_tree.edges)
 
 
+@pytest.mark.slow
 def test_sample_spanning_tree_large_sample():
     """
     Sample a single spanning tree from the distribution created in the last test
@@ -1018,6 +1019,7 @@ def test_asadpour_empty_graph():
     pytest.raises(nx.NetworkXError, nx_app.asadpour_atsp, G)
 
 
+@pytest.mark.slow
 def test_asadpour_integral_held_karp():
     """
     This test uses an integral held karp solution and the held karp function
diff --git a/networkx/algorithms/centrality/tests/test_betweenness_centrality.py b/networkx/algorithms/centrality/tests/test_betweenness_centrality.py
index d563db917..d66d4a310 100644
--- a/networkx/algorithms/centrality/tests/test_betweenness_centrality.py
+++ b/networkx/algorithms/centrality/tests/test_betweenness_centrality.py
@@ -56,6 +56,7 @@ class TestBetweennessCentrality:
             assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
 
     def test_sample_from_P3(self):
+        """Betweenness centrality: P3 sample"""
         G = nx.path_graph(3)
         b_answer = {0: 0.0, 1: 1.0, 2: 0.0}
         b = nx.betweenness_centrality(G, k=3, weight=None, normalized=False, seed=1)
@@ -513,6 +514,44 @@ class TestWeightedBetweennessCentrality:
         for n in sorted(G):
             assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
 
+    def test_G3(self):
+        """Weighted betweenness centrality: G3"""
+        G = nx.MultiGraph(weighted_G())
+        es = list(G.edges(data=True))[::2]  # duplicate every other edge
+        G.add_edges_from(es)
+        b_answer = {0: 2.0, 1: 0.0, 2: 4.0, 3: 3.0, 4: 4.0, 5: 0.0}
+        b = nx.betweenness_centrality(G, weight="weight", normalized=False)
+        for n in sorted(G):
+            assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
+
+    def test_G4(self):
+        """Weighted betweenness centrality: G4"""
+        G = nx.MultiDiGraph()
+        G.add_weighted_edges_from(
+            [
+                ("s", "u", 10),
+                ("s", "x", 5),
+                ("s", "x", 6),
+                ("u", "v", 1),
+                ("u", "x", 2),
+                ("v", "y", 1),
+                ("v", "y", 1),
+                ("x", "u", 3),
+                ("x", "v", 5),
+                ("x", "y", 2),
+                ("x", "y", 3),
+                ("y", "s", 7),
+                ("y", "v", 6),
+                ("y", "v", 6),
+            ]
+        )
+
+        b_answer = {"y": 5.0, "x": 5.0, "s": 4.0, "u": 2.0, "v": 2.0}
+
+        b = nx.betweenness_centrality(G, weight="weight", normalized=False)
+        for n in sorted(G):
+            assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
+
 
 class TestEdgeBetweennessCentrality:
     def test_K5(self):
@@ -598,6 +637,7 @@ class TestWeightedEdgeBetweennessCentrality:
             assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
 
     def test_weighted_graph(self):
+        """Edge betweenness centrality: weighted"""
         eList = [
             (0, 1, 5),
             (0, 2, 4),
@@ -627,6 +667,7 @@ class TestWeightedEdgeBetweennessCentrality:
             assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
 
     def test_normalized_weighted_graph(self):
+        """Edge betweenness centrality: normalized weighted"""
         eList = [
             (0, 1, 5),
             (0, 2, 4),
@@ -655,3 +696,84 @@ class TestWeightedEdgeBetweennessCentrality:
         norm = len(G) * (len(G) - 1) / 2
         for n in sorted(G.edges()):
             assert b[n] == pytest.approx(b_answer[n] / norm, abs=1e-7)
+
+    def test_weighted_multigraph(self):
+        """Edge betweenness centrality: weighted multigraph"""
+        eList = [
+            (0, 1, 5),
+            (0, 1, 4),
+            (0, 2, 4),
+            (0, 3, 3),
+            (0, 3, 3),
+            (0, 4, 2),
+            (1, 2, 4),
+            (1, 3, 1),
+            (1, 3, 2),
+            (1, 4, 3),
+            (1, 4, 4),
+            (2, 4, 5),
+            (3, 4, 4),
+            (3, 4, 4),
+        ]
+        G = nx.MultiGraph()
+        G.add_weighted_edges_from(eList)
+        b = nx.edge_betweenness_centrality(G, weight="weight", normalized=False)
+        b_answer = {
+            (0, 1, 0): 0.0,
+            (0, 1, 1): 0.5,
+            (0, 2, 0): 1.0,
+            (0, 3, 0): 0.75,
+            (0, 3, 1): 0.75,
+            (0, 4, 0): 1.0,
+            (1, 2, 0): 2.0,
+            (1, 3, 0): 3.0,
+            (1, 3, 1): 0.0,
+            (1, 4, 0): 1.5,
+            (1, 4, 1): 0.0,
+            (2, 4, 0): 1.0,
+            (3, 4, 0): 0.25,
+            (3, 4, 1): 0.25,
+        }
+        for n in sorted(G.edges(keys=True)):
+            assert b[n] == pytest.approx(b_answer[n], abs=1e-7)
+
+    def test_normalized_weighted_multigraph(self):
+        """Edge betweenness centrality: normalized weighted multigraph"""
+        eList = [
+            (0, 1, 5),
+            (0, 1, 4),
+            (0, 2, 4),
+            (0, 3, 3),
+            (0, 3, 3),
+            (0, 4, 2),
+            (1, 2, 4),
+            (1, 3, 1),
+            (1, 3, 2),
+            (1, 4, 3),
+            (1, 4, 4),
+            (2, 4, 5),
+            (3, 4, 4),
+            (3, 4, 4),
+        ]
+        G = nx.MultiGraph()
+        G.add_weighted_edges_from(eList)
+        b = nx.edge_betweenness_centrality(G, weight="weight", normalized=True)
+        b_answer = {
+            (0, 1, 0): 0.0,
+            (0, 1, 1): 0.5,
+            (0, 2, 0): 1.0,
+            (0, 3, 0): 0.75,
+            (0, 3, 1): 0.75,
+            (0, 4, 0): 1.0,
+            (1, 2, 0): 2.0,
+            (1, 3, 0): 3.0,
+            (1, 3, 1): 0.0,
+            (1, 4, 0): 1.5,
+            (1, 4, 1): 0.0,
+            (2, 4, 0): 1.0,
+            (3, 4, 0): 0.25,
+            (3, 4, 1): 0.25,
+        }
+        norm = len(G) * (len(G) - 1) / 2
+        for n in sorted(G.edges(keys=True)):
+            assert b[n] == pytest.approx(b_answer[n] / norm, abs=1e-7)
diff --git a/networkx/algorithms/community/tests/test_modularity_max.py b/networkx/algorithms/community/tests/test_modularity_max.py
index 8a0cddbc9..3bae85285 100644
--- a/networkx/algorithms/community/tests/test_modularity_max.py
+++ b/networkx/algorithms/community/tests/test_modularity_max.py
@@ -53,6 +53,33 @@ def test_greedy_modularity_communities_relabeled():
     assert greedy_modularity_communities(G) == expected
 
 
+def test_greedy_modularity_communities_directed():
+    G = nx.DiGraph(
+        [
+            ("a", "b"),
+            ("a", "c"),
+            ("b", "c"),
+            ("b", "d"),  # inter-community edge
+            ("d", "e"),
+            ("d", "f"),
+            ("d", "g"),
+            ("f", "g"),
+            ("d", "e"),
+            ("f", "e"),
+        ]
+    )
+    expected = [frozenset({"f", "g", "e", "d"}), frozenset({"a", "b", "c"})]
+    assert greedy_modularity_communities(G) == expected
+
+    # with loops
+    G = nx.DiGraph()
+    G.add_edges_from(
+        [(1, 1), (1, 2), (1, 3), (2, 3), (1, 4), (4, 4), (5, 5), (4, 5), (4, 6), (5, 6)]
+    )
+    expected = [frozenset({1, 2, 3}), frozenset({4, 5, 6})]
+    assert greedy_modularity_communities(G) == expected
+
+
 def test_modularity_communities_weighted():
     G = nx.balanced_tree(2, 3)
     for (a, b) in G.edges:
@@ -66,6 +93,161 @@ def test_modularity_communities_weighted():
     assert greedy_modularity_communities(G, weight="weight") == expected
 
 
+def test_modularity_communities_directed_weighted():
+    G = nx.DiGraph()
+    G.add_weighted_edges_from(
+        [
+            (1, 2, 5),
+            (1, 3, 3),
+            (2, 3, 6),
+            (2, 6, 1),
+            (1, 4, 1),
+            (4, 5, 3),
+            (4, 6, 7),
+            (5, 6, 2),
+            (5, 7, 5),
+            (5, 8, 4),
+            (6, 8, 3),
+        ]
+    )
+    expected = [frozenset({4, 5, 6, 7, 8}), frozenset({1, 2, 3})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+    # A large weight of the edge (2, 6) causes 6 to change group, even if it shares
+    # only one connection with the new group and 3 with the old one.
+    G[2][6]["weight"] = 20
+    expected = [frozenset({1, 2, 3, 6}), frozenset({4, 5, 7, 8})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+
+def test_greedy_modularity_communities_multigraph():
+    G = nx.MultiGraph()
+    G.add_edges_from(
+        [
+            (1, 2),
+            (1, 2),
+            (1, 3),
+            (2, 3),
+            (1, 4),
+            (2, 4),
+            (4, 5),
+            (5, 6),
+            (5, 7),
+            (5, 7),
+            (6, 7),
+            (7, 8),
+            (5, 8),
+        ]
+    )
+    expected = [frozenset({1, 2, 3, 4}), frozenset({5, 6, 7, 8})]
+    assert greedy_modularity_communities(G) == expected
+
+    # Converting (4, 5) into a multi-edge causes node 4 to change group.
+    G.add_edge(4, 5)
+    expected = [frozenset({4, 5, 6, 7, 8}), frozenset({1, 2, 3})]
+    assert greedy_modularity_communities(G) == expected
+
+
+def test_greedy_modularity_communities_multigraph_weighted():
+    G = nx.MultiGraph()
+    G.add_weighted_edges_from(
+        [
+            (1, 2, 5),
+            (1, 2, 3),
+            (1, 3, 6),
+            (1, 3, 6),
+            (2, 3, 4),
+            (1, 4, 1),
+            (1, 4, 1),
+            (2, 4, 3),
+            (2, 4, 3),
+            (4, 5, 1),
+            (5, 6, 3),
+            (5, 6, 7),
+            (5, 6, 4),
+            (5, 7, 9),
+            (5, 7, 9),
+            (6, 7, 8),
+            (7, 8, 2),
+            (7, 8, 2),
+            (5, 8, 6),
+            (5, 8, 6),
+        ]
+    )
+    expected = [frozenset({1, 2, 3, 4}), frozenset({5, 6, 7, 8})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+    # Adding multi-edge (4, 5, 16) causes node 4 to change group.
+    G.add_edge(4, 5, weight=16)
+    expected = [frozenset({4, 5, 6, 7, 8}), frozenset({1, 2, 3})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+    # Increasing the weight of edge (1, 4) causes node 4 to return to the former group.
+    G[1][4][1]["weight"] = 3
+    expected = [frozenset({1, 2, 3, 4}), frozenset({5, 6, 7, 8})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+
+def test_greed_modularity_communities_multidigraph():
+    G = nx.MultiDiGraph()
+    G.add_edges_from(
+        [
+            (1, 2),
+            (1, 2),
+            (3, 1),
+            (2, 3),
+            (2, 3),
+            (3, 2),
+            (1, 4),
+            (2, 4),
+            (4, 2),
+            (4, 5),
+            (5, 6),
+            (5, 6),
+            (6, 5),
+            (5, 7),
+            (6, 7),
+            (7, 8),
+            (5, 8),
+            (8, 4),
+        ]
+    )
+    expected = [frozenset({1, 2, 3, 4}), frozenset({5, 6, 7, 8})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+
+def test_greed_modularity_communities_multidigraph_weighted():
+    G = nx.MultiDiGraph()
+    G.add_weighted_edges_from(
+        [
+            (1, 2, 5),
+            (1, 2, 3),
+            (3, 1, 6),
+            (1, 3, 6),
+            (3, 2, 4),
+            (1, 4, 2),
+            (1, 4, 5),
+            (2, 4, 3),
+            (3, 2, 8),
+            (4, 2, 3),
+            (4, 3, 5),
+            (4, 5, 2),
+            (5, 6, 3),
+            (5, 6, 7),
+            (6, 5, 4),
+            (5, 7, 9),
+            (5, 7, 9),
+            (7, 6, 8),
+            (7, 8, 2),
+            (8, 7, 2),
+            (5, 8, 6),
+            (5, 8, 6),
+        ]
+    )
+    expected = [frozenset({1, 2, 3, 4}), frozenset({5, 6, 7, 8})]
+    assert greedy_modularity_communities(G, weight="weight") == expected
+
+
 def test_resolution_parameter_impact():
     G = nx.barbell_graph(5, 3)
 
diff --git a/networkx/algorithms/tests/test_link_prediction.py b/networkx/algorithms/tests/test_link_prediction.py
index bbb58299a..5ce738e06 100644
--- a/networkx/algorithms/tests/test_link_prediction.py
+++ b/networkx/algorithms/tests/test_link_prediction.py
@@ -165,6 +165,46 @@ class TestAdamicAdarIndex:
         )
 
 
+class TestCommonNeighborCentrality:
+    @classmethod
+    def setup_class(cls):
+        cls.func = staticmethod(nx.common_neighbor_centrality)
+        cls.test = partial(_test_func, predict_func=cls.func)
+
+    def test_K5(self):
+        G = nx.complete_graph(5)
+        self.test(G, [(0, 1)], [(0, 1, 3.0)], alpha=1)
+        self.test(G, [(0, 1)], [(0, 1, 5.0)], alpha=0)
+
+    def test_P3(self):
+        G = nx.path_graph(3)
+        self.test(G, [(0, 2)], [(0, 2, 1.25)], alpha=0.5)
+
+    def test_S4(self):
+        G = nx.star_graph(4)
+        self.test(G, [(1, 2)], [(1, 2, 1.75)], alpha=0.5)
+
+    @pytest.mark.parametrize("graph_type", (nx.DiGraph, nx.MultiGraph, nx.MultiDiGraph))
+    def test_notimplemented(self, graph_type):
+        assert pytest.raises(
+            nx.NetworkXNotImplemented, self.func, graph_type([(0, 1), (1, 2)]), [(0, 2)]
+        )
+
+    def test_no_common_neighbor(self):
+        G = nx.Graph()
+        G.add_nodes_from([0, 1])
+        self.test(G, [(0, 1)], [(0, 1, 0)])
+
+    def test_equal_nodes(self):
+        G = nx.complete_graph(4)
+        assert pytest.raises(nx.NetworkXAlgorithmError, self.test, G, [(0, 0)], [])
+
+    def test_all_nonexistent_edges(self):
+        G = nx.Graph()
+        G.add_edges_from([(0, 1), (0, 2), (2, 3)])
+        self.test(G, None, [(0, 3, 1.5), (1, 2, 1.5), (1, 3, 2 / 3)], alpha=0.5)
+
+
 class TestPreferentialAttachment:
     @classmethod
     def setup_class(cls):
diff --git a/networkx/drawing/tests/test_pylab.py b/networkx/drawing/tests/test_pylab.py
index d533cebd1..cb79708b7 100644
--- a/networkx/drawing/tests/test_pylab.py
+++ b/networkx/drawing/tests/test_pylab.py
@@ -423,3 +423,13 @@ def test_draw_networkx_arrows_default_directed(drawing_func):
     )
     assert ax.patches
     plt.delaxes(ax)
+
+
+def test_edgelist_kwarg_not_ignored():
+    # See gh-4994
+    G = nx.path_graph(3)
+    G.add_edge(0, 0)
+    fig, ax = plt.subplots()
+    nx.draw(G, edgelist=[(0, 1), (1, 2)], ax=ax)  # Exclude self-loop from edgelist
+    assert not ax.patches
+    plt.delaxes(ax)