-
Notifications
You must be signed in to change notification settings - Fork 0
/
community_louvain.m
198 lines (183 loc) · 7.6 KB
/
community_louvain.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
function [M,Q]=community_louvain(W,gamma,M0,B)
%COMMUNITY_LOUVAIN Optimal community structure
%
% M = community_louvain(W);
% [M,Q] = community_louvain(W,gamma);
% [M,Q] = community_louvain(W,gamma,M0);
% [M,Q] = community_louvain(W,gamma,M0,'potts');
% [M,Q] = community_louvain(W,gamma,M0,'negative_asym');
% [M,Q] = community_louvain(W,[],[],B);
%
% The optimal community structure is a subdivision of the network into
% nonoverlapping groups of nodes which maximizes the number of within-
% group edges, and minimizes the number of between-group edges.
%
% This function is a fast and accurate multi-iterative generalization of
% the Louvain community detection algorithm. This function subsumes and
% improves upon,
% modularity_louvain_und.m, modularity_finetune_und.m,
% modularity_louvain_dir.m, modularity_finetune_dir.m,
% modularity_louvain_und_sign.m
% and additionally allows to optimize other objective functions (includes
% built-in Potts-model Hamiltonian, allows for custom objective-function
% matrices).
%
% Inputs:
% W,
% directed/undirected weighted/binary connection matrix with
% positive and possibly negative weights.
% gamma,
% resolution parameter (optional)
% gamma>1, detects smaller modules
% 0<=gamma<1, detects larger modules
% gamma=1, classic modularity (default)
% M0,
% initial community affiliation vector (optional)
% B,
% objective-function type or custom objective matrix (optional)
% 'modularity', modularity (default)
% 'potts', Potts-model Hamiltonian (for binary networks)
% 'negative_sym', symmetric treatment of negative weights
% 'negative_asym', asymmetric treatment of negative weights
% B, custom objective-function matrix
%
% Note: see Rubinov and Sporns (2011) for a discussion of
% symmetric vs. asymmetric treatment of negative weights.
%
% Outputs:
% M,
% community affiliation vector
% Q,
% optimized community-structure statistic (modularity by default)
%
% Example:
% % Iterative community finetuning.
% % W is the input connection matrix.
% n = size(W,1); % number of nodes
% M = 1:n; % initial community affiliations
% Q0 = -1; Q1 = 0; % initialize modularity values
% while Q1-Q0>1e-5; % while modularity increases
% Q0 = Q1; % perform community detection
% [M, Q1] = community_louvain(W, [], M);
% end
%
% References:
% Blondel et al. (2008) J. Stat. Mech. P10008.
% Reichardt and Bornholdt (2006) Phys. Rev. E 74, 016110.
% Ronhovde and Nussinov (2008) Phys. Rev. E 80, 016109
% Sun et al. (2008) Europhysics Lett 86, 28004.
% Rubinov and Sporns (2011) Neuroimage 56:2068-79.
%
% Mika Rubinov, U Cambridge 2015-2016
% Modification history
% 2015: Original
% 2016: Included generalization for negative weights.
% Enforced binary network input for Potts-model Hamiltonian.
% Streamlined code and expanded documentation.
W=double(W); % convert to double format
n=length(W); % get number of nodes
s=sum(sum(W)); % get sum of edges
if ~exist('B','var') || isempty(B)
type_B = 'modularity';
elseif ischar(B)
type_B = B;
else
type_B = 0;
if exist('gamma','var') && ~isempty(gamma)
warning('Value of gamma is ignored in generalized mode.')
end
end
if ~exist('gamma','var') || isempty(gamma)
gamma = 1;
end
if strcmp(type_B,'negative_sym') || strcmp(type_B,'negative_asym')
W0 = W.*(W>0); %positive weights matrix
s0 = sum(sum(W0)); %weight of positive links
B0 = W0-gamma*(sum(W0,2)*sum(W0,1))/s0; %positive modularity
W1 =-W.*(W<0); %negative weights matrix
s1 = sum(sum(W1)); %weight of negative links
if s1 %negative modularity
B1 = W1-gamma*(sum(W1,2)*sum(W1,1))/s1;
else
B1 = 0;
end
elseif min(min(W))<-1e-10
err_string = [
'The input connection matrix contains negative weights.\nSpecify ' ...
'''negative_sym'' or ''negative_asym'' objective-function types.'];
error(sprintf(err_string)) %#ok<SPERR>
end
if strcmp(type_B,'potts') && any(any(W ~= logical(W)))
error('Potts-model Hamiltonian requires a binary W.')
end
if type_B
switch type_B
case 'modularity'; B = (W-gamma*(sum(W,2)*sum(W,1))/s)/s;
case 'potts'; B = W-gamma*(~W);
case 'negative_sym'; B = B0/(s0+s1) - B1/(s0+s1);
case 'negative_asym'; B = B0/s0 - B1/(s0+s1);
otherwise; error('Unknown objective function.');
end
else % custom objective function matrix as input
B = double(B);
if ~isequal(size(W),size(B))
error('W and B must have the same size.')
end
end
if ~exist('M0','var') || isempty(M0)
M0=1:n;
elseif numel(M0)~=n
error('M0 must contain n elements.')
end
[~,~,Mb] = unique(M0);
M = Mb;
B = (B+B.')/2; % symmetrize modularity matrix
Hnm=zeros(n,n); % node-to-module degree
for m=1:max(Mb) % loop over modules
Hnm(:,m)=sum(B(:,Mb==m),2);
end
Q0 = -inf;
Q = sum(B(bsxfun(@eq,M0,M0.'))); % compute modularity
first_iteration = true;
while Q-Q0>1e-10
flag = true; % flag for within-hierarchy search
while flag
flag = false;
for u=randperm(n) % loop over all nodes in random order
ma = Mb(u); % current module of u
dQ = Hnm(u,:) - Hnm(u,ma) + B(u,u);
dQ(ma) = 0; % (line above) algorithm condition
[max_dQ,mb] = max(dQ); % maximal increase in modularity and corresponding module
if max_dQ>1e-10 % if maximal increase is positive
flag = true;
Mb(u) = mb; % reassign module
Hnm(:,mb) = Hnm(:,mb)+B(:,u); % change node-to-module strengths
Hnm(:,ma) = Hnm(:,ma)-B(:,u);
end
end
end
[~,~,Mb] = unique(Mb); % new module assignments
M0 = M;
if first_iteration
M=Mb;
first_iteration=false;
else
for u=1:n % loop through initial module assignments
M(M0==u)=Mb(u); % assign new modules
end
end
n=max(Mb); % new number of modules
B1=zeros(n); % new weighted matrix
for u=1:n
for v=u:n
bm=sum(sum(B(Mb==u,Mb==v))); % pool weights of nodes in same module
B1(u,v)=bm;
B1(v,u)=bm;
end
end
B=B1;
Mb=1:n; % initial module assignments
Hnm=B; % node-to-module strength
Q0=Q;
Q=trace(B); % compute modularity
end