Repository: networkx/networkx
Pull request: Correct louvain formula
Closes #5175 

After playing around with the modularity formula I figured we weren't accounting for the modularity cost when we remove a node from its original community.

The following example can prove the correctness of the new formula:

Initially we have a partition and we move node 2 from the 0th community into the 1st community. We use NetworkX's `modularity` function to measure the actual gain.

```Python
>>> import networkx as nx
>>> from networkx.algorithms.community import louvain_communities, quality
>>> G = nx.karate_club_graph()
>>> partition = louvain_communities(G, weight="weight", seed=123)
>>> print(partition)
[{1, 2, 3, 7, 12, 13}, {0, 4, 5, 6, 10, 11, 16, 17, 19, 21}, {23, 24, 25, 27, 28, 31}, {32, 33, 8, 9, 14, 15, 18, 20, 22, 26, 29, 30}]
>>> before = quality.modularity(G, partition, weight="weight")
>>> partition[0].remove(2)
>>> partition[1].add(2)
>>> print(partition)
[{1, 3, 7, 12, 13}, {0, 2, 4, 5, 6, 10, 11, 16, 17, 19, 21}, {23, 24, 25, 27, 28, 31}, {32, 33, 8, 9, 14, 15, 18, 20, 22, 26, 29, 30}]
>>> after = quality.modularity(G, partition, weight="weight")
>>> print(after - before)
-0.0553494124922696
```
The below snippet computes both the remove and insert parameters to calculate the gain.
``` Python3
Stot = 0
RStot = 0
for n in partition[1]:
    if n == 2:
        continue
    Stot += G.degree(n, weight="weight")
for n in partition[0]:
    RStot += G.degree(n, weight="weight")

r_wt = 0
wt = 0
for u, w in G[2].items():
    if u in partition[0]:
        r_wt += w["weight"]
    if u in partition[1]:
        wt += w["weight"]

m = G.size(weight="weight")
ki = G.degree(2, weight="weight")
```

Finally:
```Python3
>>> gain = (wt) / (m) - (Stot * ki) / (2 * m ** 2) # Original formula didn't take into consideration the removal of node from initial community
>>> remove_cost = -r_wt / m + (RStot * ki) / (2 * m ** 2)
>>> print(gain) # Incorrect Result
-0.011131725417439706
>>> print(remove_cost)
-0.04421768707482993
>>> print(remove_cost + gain) # Actual modularity gain from whole movement process
-0.05534941249226964
```

PS. I have added a test for the issue's case which doesn't loop forever with the new formula. Feel free to propose any better test cases that I might miss.

Problem statement:
louvain for large directed graph cannot stop
<!-- If you have a general question about NetworkX, please use the discussions tab to create a new discussion -->

<!--- Provide a general summary of the issue in the Title above -->
the nb_moves cannot decrease

### Current Behavior
<!--- Tell us what happens instead of the expected behavior -->

### Expected Behavior
<!--- Tell us what should happen -->

### Steps to Reproduce
<!--- Provide a minimal example that reproduces the bug -->

### Environment
<!--- Please provide details about your local environment -->
Python version:
NetworkX version:


### Additional context
<!--- Add any other context about the problem here, screenshots, etc. -->

Hints from issue discussion:
Try starting with a small graph and increasing the size to see how the time varies.
It might be that you graph is so big it will just take a long time.
my graph has 30000 nodes and the nb_moves is always around 27000 and cannot drop
I also try a graph with only 100 edges, but the same thing happen

@dschult 
my graph seemingly  has very bad structure, and i also test a graph  with 10 edges, For some bad structures the program seems to get stuck in a dead loop, I tried to shuffle  the neighbor_com list, which seems to work for me
> and i also test a graph with 10 edges, For some bad structures the program seems to get stuck in a dead loop, 

It would be very useful if you could share a minimal reproducing example that results in the bad behavior: for example, the graph with only 10 edges that results in the dead loop. This can help us determine whether there is something in the code that can be improved.
import networkx as nx

G = nx.DiGraph()
for i in range(10):
    G.add_node(i)
G.add_edges_from([(0, 2), (0, 1), (1, 0), (2, 1), (2, 0), (3, 4), (4, 3), (7, 8), (8, 7), (9, 10), (10, 9)])

pa=nx.algorithms.community.louvain_communities(G,seed=123)
@rossbar  
The above  dead loop situation could probably be solved by changing the code in louvain.py 

`for nbr_com, wt in weights2com.items():`

 to
```python
            if (nbrs[u] == {}):
                continue
            r_list = (list(weights2com.items()))
            random.shuffle(r_list)
            for nbr_com, wt in r_list:
                if is_directed:
```
However, this function still takes a lot of time when the input is a large directed graph, and the result seems not  good
So it seems that we also need to check the neighbor communities of a node in a random order to fix the issue.

When I do that on the example graph provided by @ginandsherry I get a result.
Interesting - many thanks @ginandsherry for the reproducer!

It does indeed seem that node order has an effect here. A slight modification of the initial example that results in a different insertion order has no problem:

```python
>>> import networkx as nx
>>> G = nx.DiGraph()
>>> G.add_edges_from([
...     (0, 2), (0, 1), (1, 0), (2, 1), (2, 0), (3, 4), (4, 3), (7, 8), (8, 7), (9, 10), (10, 9)
... ])
>>> G.add_node(5)  # Add node that was not included when adding edges
>>> G.nodes()
NodeView((0, 2, 1, 3, 4, 7, 8, 9, 10, 5))
>>> pa = nx.algorithms.community.louvain_communities(G, seed=123)  # returns immediately
```

However, when running the reproducer as presented (i.e. adding nodes prior to edges) I get the stated behavior:

```python
>>> import networkx as nx
>>> G = nx.DiGraph()
>>> G.add_nodes_from(n for n in range(10))
>>> G.add_edges_from([
...     (0, 2), (0, 1), (1, 0), (2, 1), (2, 0), (3, 4), (4, 3), (7, 8), (8, 7), (9, 10), (10, 9)
... ])
>>> G.nodes()
NodeView((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
>>> pa = nx.algorithms.community.louvain_communities(G, seed=123)  # hangs
```
The order in which the nodes are considered definitely affects the outcome as stated by the authors of the algorithm [here](https://arxiv.org/pdf/0803.0476.pdf)

> One should also note that the output of the algorithm depends on the order in which
the nodes are considered. Preliminary results on several test cases seem to indicate that
the ordering of the nodes does not have a significant influence on the modularity that
is obtained. However the ordering can influence the computation time. The problem of
choosing an order is thus worth studying since it could give good heuristics to enhance
the computation time.


Here is the original implementation in C++ for the Directed Louvain https://github.com/nicolasdugue/DirectedLouvain

This is the function that computes the modularity gain on each step
https://github.com/nicolasdugue/DirectedLouvain/blob/00398d5659ee49973171f67db4f6dacb4230bfbf/include/community.h#L121-L134

It looks similar to the one we have in NetworkX and they don't shuffle the nodes after every pass. Maybe I am missing something in the formula, but this would be a good place to investigate before implementing a shuffling after every pass which could reduce the algorihtm's performance

Commit message: We need both the in and out weights to nodes in community!

Patch to review:
diff --git a/networkx/algorithms/community/louvain.py b/networkx/algorithms/community/louvain.py
index e279f8c3c14..c98a6a1c4ed 100644
--- a/networkx/algorithms/community/louvain.py
+++ b/networkx/algorithms/community/louvain.py
@@ -219,10 +219,18 @@ def _one_level(G, m, partition, resolution=1, is_directed=False, seed=None):
         out_degrees = dict(G.out_degree(weight="weight"))
         Stot_in = [deg for deg in in_degrees.values()]
         Stot_out = [deg for deg in out_degrees.values()]
+        # Calculate weights for both in and out neighbours
+        nbrs = {}
+        for u in G:
+            nbrs[u] = defaultdict(float)
+            for _, n, wt in G.out_edges(u, data=True):
+                nbrs[u][n] += wt["weight"]
+            for n, _, wt in G.in_edges(u, data=True):
+                nbrs[u][n] += wt["weight"]
     else:
         degrees = dict(G.degree(weight="weight"))
         Stot = [deg for deg in degrees.values()]
-    nbrs = {u: {v: data["weight"] for v, data in G[u].items() if v != u} for u in G}
+        nbrs = {u: {v: data["weight"] for v, data in G[u].items() if v != u} for u in G}
     rand_nodes = list(G.nodes)
     seed.shuffle(rand_nodes)
     nb_moves = 1
@@ -238,22 +246,36 @@ def _one_level(G, m, partition, resolution=1, is_directed=False, seed=None):
                 out_degree = out_degrees[u]
                 Stot_in[best_com] -= in_degree
                 Stot_out[best_com] -= out_degree
+                remove_cost = (
+                    -weights2com[best_com] / m
+                    + resolution
+                    * (out_degree * Stot_in[best_com] + in_degree * Stot_out[best_com])
+                    / m**2
+                )
             else:
                 degree = degrees[u]
                 Stot[best_com] -= degree
+                remove_cost = -weights2com[best_com] / m + resolution * (
+                    Stot[best_com] * degree
+                ) / (2 * m**2)
             for nbr_com, wt in weights2com.items():
                 if is_directed:
                     gain = (
-                        wt
+                        remove_cost
+                        + wt / m
                         - resolution
                         * (
                             out_degree * Stot_in[nbr_com]
                             + in_degree * Stot_out[nbr_com]
                         )
-                        / m
+                        / m**2
                     )
                 else:
-                    gain = 2 * wt - resolution * (Stot[nbr_com] * degree) / m
+                    gain = (
+                        remove_cost
+                        + wt / m
+                        - resolution * (Stot[nbr_com] * degree) / (2 * m**2)
+                    )
                 if gain > best_mod:
                     best_mod = gain
                     best_com = nbr_com
diff --git a/networkx/algorithms/community/tests/test_louvain.py b/networkx/algorithms/community/tests/test_louvain.py
index 350784f2aba..7aa9ccbb4d4 100644
--- a/networkx/algorithms/community/tests/test_louvain.py
+++ b/networkx/algorithms/community/tests/test_louvain.py
@@ -43,6 +43,58 @@ def test_partition():
     assert part == partition
 
 
+def test_directed_partition():
+    """
+    Test 2 cases that were looping infinitely
+    from issues #5175 and #5704
+    """
+    G = nx.DiGraph()
+    H = nx.DiGraph()
+    H.add_nodes_from([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
+    G_edges = [
+        (0, 2),
+        (0, 1),
+        (1, 0),
+        (2, 1),
+        (2, 0),
+        (3, 4),
+        (4, 3),
+        (7, 8),
+        (8, 7),
+        (9, 10),
+        (10, 9),
+    ]
+    H_edges = [
+        (1, 2),
+        (1, 6),
+        (1, 9),
+        (2, 3),
+        (2, 4),
+        (2, 5),
+        (3, 4),
+        (4, 3),
+        (4, 5),
+        (5, 4),
+        (6, 7),
+        (6, 8),
+        (9, 10),
+        (9, 11),
+        (10, 11),
+        (11, 10),
+    ]
+    G.add_edges_from(G_edges)
+    H.add_edges_from(H_edges)
+
+    G_expected_partition = [{0, 1, 2}, {3, 4}, {8, 7}, {9, 10}]
+    G_partition = louvain_communities(G, seed=123, weight=None)
+
+    H_expected_partition = [{2, 3, 4, 5}, {8, 1, 6, 7}, {9, 10, 11}]
+    H_partition = louvain_communities(H, seed=123, weight=None)
+
+    assert G_partition == G_expected_partition
+    assert H_partition == H_expected_partition
+
+
 def test_none_weight_param():
     G = nx.karate_club_graph()
     nx.set_edge_attributes(
@@ -84,9 +136,9 @@ def test_quality():
     quality4 = partition_quality(J, partition4)[0]
 
     assert quality >= 0.65
-    assert quality2 >= 0.85
+    assert quality2 >= 0.65
     assert quality3 >= 0.65
-    assert quality4 >= 0.85
+    assert quality4 >= 0.65
 
 
 def test_multigraph():
@@ -120,7 +172,7 @@ def test_threshold():
     G = nx.LFR_benchmark_graph(
         250, 3, 1.5, 0.009, average_degree=5, min_community=20, seed=10
     )
-    partition1 = louvain_communities(G, threshold=0.2, seed=2)
+    partition1 = louvain_communities(G, threshold=0.3, seed=2)
     partition2 = louvain_communities(G, seed=2)
     mod1 = modularity(G, partition1)
     mod2 = modularity(G, partition2)


Code review comments:
diff_hunk: @@ -219,10 +219,18 @@ def _one_level(G, m, partition, resolution=1, is_directed=False, seed=None):
         out_degrees = dict(G.out_degree(weight="weight"))
         Stot_in = [deg for deg in in_degrees.values()]
         Stot_out = [deg for deg in out_degrees.values()]
+        # Calculate weights for both in and out neighbours
+        nbrs = {}
+        for u in G:
+            nbrs[u] = defaultdict(float)
+            for _, n, wt in G.out_edges(u, data=True):
+                nbrs[u][n] += wt["weight"]
+            for n, _, wt in G.in_edges(u, data=True):
+                nbrs[u][n] += wt["weight"]
path: networkx/algorithms/community/louvain.py
text: @user1:
Just a nit, but you can pass the edge attr name directly to the `data` kwarg to make this a little cleaner:

```suggestion
            for _, n, wt in G.out_edges(u, data="weight"):
                nbrs[u][n] += wt
            for n, _, wt in G.in_edges(u, data="weight"):
                nbrs[u][n] += wt
```