diff --git a/REQUIRE b/REQUIRE index feb3921..0247e8f 100644 --- a/REQUIRE +++ b/REQUIRE @@ -2,3 +2,4 @@ julia 0.7 LightGraphs 1.1.0 Clustering ArnoldiMethod +SimpleWeightedGraphs diff --git a/src/CommunityDetection.jl b/src/CommunityDetection.jl index 38b6d5f..afb83fd 100644 --- a/src/CommunityDetection.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index f884041..7e1ee0d 100644 --- a/test/runtests.jl +++ b/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