function maxcliques = findmaxcliques(EorA)
%FINDMAXCLIQUES Finds returns the set of maximal cliques in a graph.
% maxcliques = findmaxcliques(E or A)
%     E = |E(G)|-by-2 matrix enumerating the edges of graph G
%     A = |V(G)|-by-|V(G)| binary adjacency matrix for graph G
%
% THIS ALGORITHM HAS EXPONENTIAL RUNNING TIME; IT IS BASED ON THE ALGORITHM
% BY BRON/KERBOSCH (1971) AND REFINED BY TOMITA/TANAKA/TAKAHASHI (2006).

% (by Delbert Dueck, delbert@cs.columbia.edu, 2009)

if size(EorA,2)==2, E=EorA; end;
if size(EorA,1)==size(EorA,2), A=EorA; end;

if exist('E','var'), V=max(max(E));
    for v=1:V, neighbors{v,1}=unique([E(find(E(:,1)==v),2)',E(find(E(:,2)==v),1)']); end; % enumerate neighbors
elseif exist('A','var'), V=size(A,1);
    for v=1:V, neighbors{v,1}=find(A(v,:)); neighbors{v}(find(neighbors{v}==v))=[]; end; % enumerate neighbors
end;

col=1; newline=1; max_cliques=0;

expand(1:V,1:V);

while (all(max_cliques(end,:)==0)), max_cliques(end,:)=[]; end;
for i=1:size(max_cliques,1),
	temp = find(max_cliques(i,:)~=0,1,'first')-1;
    if temp>0, max_cliques(i,1:temp)=max_cliques(i-1,1:temp); end;
end;
maxcliques = cell(0,0);
for i=1:size(max_cliques,1),
    temp = find(max_cliques(i,:)==-1,1,'first');
    if ~isempty(temp),
        if V<=intmax('uint8'), maxcliques{end+1,1} = uint8(max_cliques(i,1:temp-1));
        elseif V<=intmax('uint16'), maxcliques{end+1,1} = uint16(max_cliques(i,1:temp-1));
        elseif V<=intmax('uint32'), maxcliques{end+1,1} = uint32(max_cliques(i,1:temp-1));
        else maxcliques{end+1,1} = max_cliques(i,1:temp-1); end;
    end;
end;

    
    function expand(SUBG, CAND)
        if isempty(SUBG), max_cliques(end,col)=-1; newline=1; return; end; % fprintf('clique,');

        % u is a vertex in SUBG that maximizes |CAND ^ N(u)|
        u = SUBG(1); card=numel(intersect(CAND,neighbors{u}));
        for uu=SUBG(2:end),
            temp = numel(intersect(CAND,neighbors{uu}));
            if temp>card, u=uu; card=temp; end;
        end;

        while ~isempty(setdiff(CAND,neighbors{u})),
            q = setdiff(CAND,neighbors{u}); q=q(1);
            max_cliques(end,col)=q; col=col+1; newline=1; % fprintf('%d,',q);
            expand(intersect(SUBG,neighbors{q}),intersect(CAND,neighbors{q}));
            CAND(find(CAND==q))=[];
            col=col-1; if newline, max_cliques(end+1,1)=0; newline=0; end; % fprintf('back,');
        end;
    end

i=1; while i<=numel(maxcliques), if numel(maxcliques{i})<=1, maxcliques(i)=[]; else i=i+1; end; end;

end




% A = [0 1 0 0 0 0 0 0 1; ...
%      1 0 1 0 0 0 0 0 1; ...
%      0 1 0 1 0 0 0 1 1; ...
%      0 0 1 0 1 1 1 1 0; ...
%      0 0 0 1 0 1 0 0 0; ...
%      0 0 0 1 1 0 1 1 0; ...
%      0 0 0 1 0 1 0 1 0; ...
%      0 0 1 1 0 1 1 0 0; ...
%      1 1 1 0 0 0 0 0 0];
% for i=1:length(A), A(i,i)=1; end;
%      
% A = rand(20); A=.5*A+.5*A';
% A = (A>.5);
% for i=1:length(A), A(i,i)=1; end;
% maxcliques = findmaxcliques(A);
% 
