function [grps,maxEV] = multileadevsplit(A,varargin)

% MULTILEADEVSPLIT partition graph using leading eigenvector of modularity matrix
%   [S,U] = MULTILEADEVSPLIT(A) splits the vertices of the graph in adjacency matrix
%   A into n groups, defined in the vector S (each element taking value
%   {1,...,n} to indicate group membership of the vertex at that index). It
%   also returns U, the leading eigenvector of the modularity matrix.
%
%   MULTILEADEVSPLIT(...,FLAG) where FLAG is a string featuring any
%   combination of:
%       'r': adds the refinement step after each sub-graph,
%               so that each split maximises Q.
%       'a': keeps all splits found, even if Q < 0
%
%
%   Notes: 
%   (1) this repeats the leading eigenvector split on each sub-graph
%   until the change in modularity is not positive...
%
%   (2) or could use singular value decomposition to find eigenvalues,
%   as standard MatLab routine EIG() does not converge for some sample
%   directed graphs (e.g. mac95) even though they are made symmetric.
%   (Note though that using SVD does not split the dolphin graph into 3)
%
%   References: 
%   (1) Newman, M. E. J. (2006a) "Finding community structure in
%       networks using the eigenvectors of matrices". Phys Rev E, in press.
%   (2) Newman, M. E. J. (2006b) "Modularity and community structure in
%       networks". PNAS, in press.
%   (3) Reichardt & Bornholdt (2006) "Statistical mechanics of community detection" 
%       Phys Rev E. 74, 016110
%   
%   Mark Humphries 27/1/2009

blnRefine = 0;
blnAll = 0;
if nargin >= 2 
    if strfind(varargin{1},'r') blnRefine = 1; end
    if strfind(varargin{1},'a') blnAll = 1; end
end

% apply correction to A for directed graphs. If graph undirected then this
% has no effect
Abar = (A + A') / 2;    

[n c] = size(A);
grps = zeros(n,1);
ks = sum(Abar);        
m = sum(ks)/2;         

[S,maxEV,B] = splitsubgraph(A,ks,m);

Q = S' * B * S;

% check modularity has increased
if ~blnAll & Q <= 1e-5    % very close or less than zero (deals with rounding errors)
    % if not increased then return nothing
    S = zeros(n,1);
    return
end

% make subgraphs
n1 = sum(S>0);
n2 = sum(S<=0);
subG1 = A(S>0,S>0);
subG2 = A(S<=0,S<=0);
grps = real(S>0) + 1;     % group 1 and group 2

subqueue = {subG1, subG2};
grpvertices = {find(S>0), find(S<=0)};

while ~isempty(subqueue)
    newA = subqueue{1};
    [S,maxEV,B] = splitsubgraph(newA,ks(grpvertices{1}),m); 
  
    if blnRefine
        S = refinesplit(S,B);
    end
      
    % compute Q
    % keyboard
    try
        Q = S' * B * S; 
    catch
        keyboard
    end
    
    %keyboard
    
    % split this group if:
    % (a) all splits are kept, and there is a split OR
    % (b) not all splits are kept, but Q > 0 (Newman 2006a suggestion)  
    if (blnAll & any(S>0) & any(S<0)) | Q > 0
        subG1 = newA(S>0,S>0);
        subG2 = newA(S<=0,S<=0);
        grp1 = grpvertices{1}(S>0);     % the ids of the original vertices
        grp2 = grpvertices{1}(S<=0);    % the ids of the original vertices
        grps(grp1) = max(grps)+1;       % re-number first group, second group retains numbers
        %grps(grp2) = max(grps)+1;
        % delete split graph from stack
        subqueue(1) = [];
        grpvertices(1) = [];
        % add new ones
        grpvertices = [{grp1,grp2} grpvertices];        
        subqueue = [{subG1,subG2} subqueue];
    else
        % don't split any further, delete top of stack
        subqueue(1) = [];
        grpvertices(1) = [];
    end

    % keyboard
end



function [S,maxEV,B] = splitsubgraph(A,ks,m) 
    % Note: if graph is directed, the correction to A should be
    % applied
    Abar = (A + A') / 2;   
    
    S = zeros(length(A),1);
    
    % generate expected graph
    P = expectedA(A);
    
    %keyboard
    
    % correct P
    Pbar = (P + P') / 2;
    
    kg = sum(Abar);        % in-degrees of sub-graph
    dg = sum(kg);       % total edges of sub-graph
    
    % this correction should apply without problem to a directed graph
    % because by this stage the graph will be undirected. The key is pass
    % the correct value of ks to this function.
    correction = kg - ks.*dg/(2*m);     % definitely 2*m here: is divided by 2 at start :-)
    
    % create modularity matrix - correction applied only to self-self
    % connections
    B = Abar - Pbar - diag(correction);   
        
    % find eigenvalues (D) and eigenvectors (V) of B
    [V,D] = eig(B);         % assumes symmetry
    eg = diag(D);
    if ~isreal(eg)
        warning('Complex eigenvalues have occurred')
    end
    
    % keyboard
    
     % largest positive eigenvalue is last element
    maxeig = max(eg); %% note - can have the same eigenvalue if nodes have no links in subgraph!
    ix = find(eg == maxeig(1));
    maxEV = V(:,ix(1)); % if more than one with the same max eigenvalue, just pick first...
    
        
    %%%% use singular value decomposition %%%%%%%%%%%%
%     [u,s,v] = svd(B);
%     sv = diag(s);
%     maxEV = v(:,find(sv == max(sv)));
     
    
    S(maxEV > 0) = 1;       % group 1
    S(maxEV <= 0) = -1;     % group 2