function S = refinesplit(S,B,varargin)
% REFINESPLIT improve modularity maximisation
% N = REFINESPLIT(S,B) attempts to further maximise the modularity value
% of the network, described in modularity matrix B, given the current 2 group
% partition in vector S. It returns the vector N of new groups that
% maximise modularity (values in {1,+1}).
%
% REFINESPLIT(...,'i') returns the groups in {1,2} instead of {1,1},
% which is more useful for postprocessing.
%
% Notes:
% #1: This works by taking each vertex in turn, assigning it to the
% other group, and recomputing modularity Q. The group assignment that
% maximises Q is kept. This is repeated until no further increases in Q
% are obtained (with each vertex only moving once to a new group  i.e.
% once moved, they cannot move back!).
%
% #2: Requires that group membership be +1 and 1, so converts the
% membership vector S if required.
%
% References:
% (1) Newman, M. E. J. (2006) "Finding community structure in
% networks using the eigenvectors of matrices". Phys Rev E, in press.
% (2) Newman, M. E. J. (2006) "Modularity and community structure in
% networks". PNAS, in press.
%
% Mark Humphries 16/10/2006
% keyboard
Q = S' * B * S;
newQ = Q;
prevQ = Q;
blnSwap = 1;
G = sort(S);
g1 = 1;
g2 = 1;
if G(1) == G(end)
% then no split to refine!
return
end
S(S==G(1)) = g1;
S(S==G(end)) = g2;
[n c] = size(B);
moved = 0;
while blnSwap
tempQ = zeros(n,1);
for i = 1:n
if i ~= moved
Sp = S;
if Sp(i) == g1 Sp(i) = g2;
else Sp(i) = g1;
end
% recompute modularity
tempQ(i) = Sp' * B * Sp;
else
% compute modularity of unmoved vertex
Sp = S;
tempQ(i) = Sp' * B * Sp;
end
end
bestnode = find(tempQ == max(tempQ) & tempQ > newQ);
if ~isempty(bestnode)
% there could be more than one with exact same change in modularity
% so take the first
moved = [moved bestnode(1)];
end
if max(tempQ) > prevQ
newQ = max(tempQ);
prevQ = newQ;
% swap permanently
if S(bestnode) == g1 S(bestnode) = g2;
else S(bestnode) = g1;
end
else
% no increase in Q
blnSwap = 0;
end
end
if nargin >= 3 & varargin{1} == 'i'
S(S==g2) = 2;
S(S==g1) = 1;
end
