Louvain method

REQUIRE

@@ -2,3 +2,4 @@ julia 0.7
 LightGraphs 1.1.0
 Clustering
 ArnoldiMethod
+SimpleWeightedGraphs

src/CommunityDetection.jl

@@ -3,8 +3,11 @@ using LightGraphs
 using ArnoldiMethod: LR, SR
 using LinearAlgebra: I, Diagonal
 using Clustering: kmeans
+using SparseArrays
+using SimpleWeightedGraphs
+using Random: shuffle
 
-export community_detection_nback, community_detection_bethe
+export community_detection_nback, community_detection_bethe, community_detection_louvain
 
 """
     community_detection_nback(g::AbstractGraph, k::Int)
@@ -107,4 +110,131 @@ function community_detection_bethe(g::AbstractGraph, k::Int=-1; kmax::Int=15)
     return labels
 end
 
+
+"""
+    community_detection_louvain(g::AbstractGraph; tol = 1e-6)
+
+Perform fast-unfolding to maximize the modularity of a graph. `tol` is the minimum amount of improvement that should be made to consider a pass
+to make progress.  Return a vector containing the vertex assignments.
+
+### References
+- [Blondel et al.](https://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008/meta)
+"""
+function community_detection_louvain(g::AbstractGraph; tol::Real = 1e-6)
+    graph = WGraph(copy(g))   # Create a weighted graph
+    comms = collect(1:nv(g))  # Assume at first that each vertex is a community
+
+    improvement = true
+    levels = Vector{Int}[]    # Hold onto the group identity at each pass through `make_pass`
+    while improvement
+        improvement = louvain_make_pass!(graph, comms, tol)  # one pass of moving nodes and aggregating into communities
+        improvement || break
+
+        # create a new aggregated graph
+        Anew = louvain_aggregate_communities!(graph, comms)
+        push!(levels, comms)
+        graph = WGraph(Anew)
+        comms = collect(1:nv(graph))
+    end
+
+    # unpack and return communities containing the original nodes
+    for i in length(levels):-1:2
+        for node in eachindex(levels[i-1])
+            levels[i-1][node] = levels[i][levels[i-1][node]]
+        end
+    end
+    return levels[1]
+end
+
+
+function louvain_make_pass!(graph::AbstractGraph, comms::Vector{Int}, tol::Real)
+    W = sum(weights(graph)) # total weight of graph
+    node_W = reshape(sum(weights(graph),dims=1), :)
+    comm_W = zeros(maximum(comms))
+    for i in 1:nv(graph)
+        comm_W[comms[i]] += node_W[i]
+    end
+
+    # `node_comm_W[i]` is the total weight between node and community `i`
+    node_comm_W = zeros(maximum(comms))
+
+    made_progress = true
+    improvement = false
+    nb_moves = 0  # number of moves
+    while made_progress
+        best_gain = 0.0
+        for i in shuffle(1:nv(graph))
+            current_comm = comms[i]
+            improvement = false
+            best_gain = 0.0
+            best_comm = current_comm
+
+            # compute weights between node `i` and its neighboring communities
+            for neighbor in neighbors(graph, i)
+                neighbor_comm = comms[neighbor]
+                if neighbor != i
+                    node_comm_W[neighbor_comm] += 2weights(graph)[i, neighbor]
+                else
+                    node_comm_W[neighbor_comm] += weights(graph)[i, neighbor]
+                end
+            end
+
+            # change of modularity if we remove node `i` from old community
+            ΔM1 = (comm_W[current_comm] - 1.0) * node_W[i] / W - node_comm_W[current_comm]
+            node_comm_W[current_comm] = 0
+
+            # using `shuffle` to select one from multiple proposed communities which have the same maximum modularity gain
+            for neighbor in shuffle(neighbors(graph,i))
+                neighbor_comm = comms[neighbor]
+                if node_comm_W[neighbor_comm] > 0 && neighbor_comm != current_comm
+                    # change of modularity if we add node to new community
+                    ΔM2 = node_comm_W[neighbor_comm] - comm_W[neighbor_comm] * node_W[i] / W
+                    modularity_gain = ΔM1 + ΔM2
+                    if modularity_gain > best_gain
+                        best_gain = modularity_gain
+                        best_comm = neighbor_comm
+                    end
+                    # this neighboring community has seen, so set the weight to zero to skip it next time
+                    node_comm_W[neighbor_comm] = 0
+                end
+            end
+            if best_comm != current_comm
+                nb_moves += 1
+                comms[i] = best_comm  # move node to new community
+                comm_W[current_comm] -= node_W[i]  # remove node weight from old community
+                comm_W[best_comm] += node_W[i]  # add node weight to new community
+            end
+        end
+
+        # check whether the algorithm is making progress for this pass
+        improvement = nb_moves > 0
+        made_progress = improvement && 2best_gain/W > tol
+    end
+    return improvement
+end
+
+
+#Aggregate matrix communities.
+function louvain_aggregate_communities!(graph::AbstractGraph,comms::Vector{Int})
+    #Renumber communities starting from 1
+    label_counters = zero(comms)
+    j = 1
+    for i in eachindex(comms)
+        if label_counters[comms[i]] == 0
+            label_counters[comms[i]] = j
+            comms[i] = j
+            j += 1
+        else
+            comms[i] = label_counters[comms[i]]
+        end
+    end
+
+    #Create a group identify matrix and aggregate into larger communities
+    n = nv(graph)
+    A = graph.weights
+    G = sparse(1:n, comms, ones(n))
+    Anew = G'*A*G  #New weighted adjacency matrix
+    return Anew
+end
+
 end #module

test/runtests.jl

@@ -119,4 +119,23 @@ end
     end
 end
 
+@testset "community_detection_louvain(z)" begin
+    n=10
+    g10 = CompleteGraph(n)
+    z = copy(g10)
+    for k=2:5
+        z = blockdiag(z, g10)
+        add_edge!(z, (k-1)*n, k*n)
+
+        c = community_detection_louvain(z)
+        @test sort(union(c)) == [1:k;]
+        a = collect(n:n:k*n)
+        @test length(c[a]) == length(unique(c[a]))
+        for i=1:k
+            cluster_range = (1:n) .+ (i-1)*n
+            @test length(unique(c[cluster_range])) == 1
+        end
+    end
+end
+
 end