Homec4science
Diffusion ISFC clustering (master) Browse

ISFCClustering/2019_03_03_BCT/core_periphery_dir.m
3209cc2b0b1dmaster

Last Change
Raw File

core_periphery_dir.m
View Options

function [C, q]=core_periphery_dir(W,gamm,C)
%CORE_PERIPHERY_DIR Core/periphery structure and core-ness statistic
%
% C = core_periphery_dir(W)
% [C,q] = core_periphery_dir(W,gamm,C0)
%
% The optimal core/periphery subdivision is a partition of the network
% into two non-overlapping groups of nodes, a core group and a periphery
% group, in a way that maximizes the number/weight of within core-group
% edges, and minimizes the number/weight of within periphery-group edges.
%
% The core-ness is a statistic which quantifies the goodness of the
% optimal core/periphery subdivision.
%
% Input: W directed (weighted or binary) connection matrix.
% gamma, core-ness resolution parameter (optional)
% gamma>1 detects small core/large periphery
% 0<=gamma<1 detects large core/small periphery
% default is gamma=1
%
% Outputs: C, binary vector of optimal core structure
% C = 1 represents nodes in the core
% C = 0 represents nodes in the periphery
% q, maximized core-ness statistic
%
% Algorithm: A version of Kernighan-Lin algorithm for graph partitioning
% used in community detection (Newman, 2006) applied to optimize a
% core-structure objective described in Borgatti and Everett (2000).
%
% Reference: Borgatti and Everett (2000) Soc Networks 21:375–395.
% Newman (2006) Phys Rev E 74:036104, PNAS 23:8577-8582.
% Rubinov, Ypma et al. (2015) PNAS 112:10032-7
%
% 2015, Mika Rubinov, U Cambridge
n = length(W); % number of nodes
W = double(W); % convert from logical
W(1:n+1:end) = 0; % clear diagonal
if ~exist('gamm','var')
gamm = 1;
end
if ~exist('C','var')
C = (rand(1,n)<0.5);
else
C = logical(reshape(C,1,n));
end
% Methodological note: cf. community detection, the core-detection
% null model is not corrected for degree (to enable detection of hubs).
s = sum(W(:));
p = mean(W(:));
b = W - gamm*p;
B = (b+b.')/(2*s); % directed core-ness matrix
q = sum(sum(B(C,C))) - sum(sum(B(~C,~C))); % core-ness statistic
f=1; % loop flag
while f
f=0;
Idx = 1:n; % initial node indices
Ct = C;
while any(Idx)
Qt = zeros(1,n); % check swaps of node indices
q0 = sum(sum(B(Ct,Ct))) - sum(sum(B(~Ct,~Ct)));
Qt( Ct) = q0 - 2*sum(B( Ct, :),2);
Qt(~Ct) = q0 + 2*sum(B(~Ct, :),2);
%%% verification that the above update is equivalent to:
% for u=Idx
% Ct(u) = ~Ct(u);
% Qt(u) = sum(sum(B(Ct,Ct))) - sum(sum(B(~Ct,~Ct)));
% Ct(u) = ~Ct(u);
% end
max_Qt = max(Qt(Idx)); % make swap with maximal
u = find(abs(Qt(Idx)-max_Qt)<1e-10);% increase in core-ness
u = u(randi(numel(u)));
Ct(Idx(u)) = ~Ct(Idx(u));
Idx(u)=[]; % remove index from consideration
if max_Qt-q>1e-10 % recompute core-ness statistic
f = 1;
C = Ct;
q = sum(sum(B(C,C))) - sum(sum(B(~C,~C)));
end
end
end
q = sum(sum(B(C,C))) - sum(sum(B(~C,~C))); % return core-ness statistic
c4science · Help