function [ClusterA,ClusterB,anc]=sepclust(dist,factor)
% SEPCLUST  Make clusters which make a separation problem solvable
%	[ClusterA,ClusterB]=sepclust(dist) based on a similarity matrix,
%	dist, where the elements correspond to the distance in similarity
%	between vertices, and the row and column numbers correspond to the
%	vertex numbers, cluster pairs are formed such that
%	 - all vertices within a cluster have less distance to each other
%	   than to any vertices in the paired cluster
%	 - given any unknown vertex and the clusters of the pairings,
%	   it is possible to determine which vertex it is
%	Vertices may correspond to classes in a classification problem.
%	
%	The problem (which is assumed to be NP) is solved in an non-optimal
%	manner. The solution starts with the two most distant vertices,
%	grows each of them into the maximal clusters, takes the remaining
%	two most distant vertices that are not yet separated (may still be
%	the same), grows into maximum clusters etc.
%	
%	[ClusterA,ClusterB]=sepclust(dist,factor) increases the distance 
%	requirement, so that some "buffer" is ensured. E.g. factor=0.9
%	requires that (distance within ) < 0.9*(distance between)
%	
%	[ClusterA,ClusterB,anc]=sepclust(dist) also returns the average number
%	of clusters, anc, each vertex occurs in.

if exist('factor')~=1
  factor=1;
end

N_T=size(dist,1);
ClusterA=[]; ClusterB=[];
NotYetSeparated=ones(N_T)-diag(ones(N_T,1));

while sum(NotYetSeparated(:))
  [I,J]=find(dist==max(max(dist.*NotYetSeparated)));
  cau=I(1);
  cbu=J(1);
  
  Prospects=notin(find(sum(NotYetSeparated([cau cbu],:))),[cau cbu]);
  for p=Prospects
    da2a=nonzero(dist(cau,cau)); if isempty(da2a), da2a=0; end
    db2b=nonzero(dist(cbu,cbu)); if isempty(db2b), db2b=0; end
    dap2b=nonzero(dist([cau p],cbu));
    dbp2a=nonzero(dist([cbu p],cau));
    dap2a=nonzero(dist([cau p],cau));
    dbp2b=nonzero(dist([cbu p],cbu));
    if (max(dap2a) < factor*min(dap2b)) & (max(db2b) < factor*min(dap2b))
      cau=[cau p];
    elseif (max(dbp2b) < factor*min(dbp2a)) & (max(da2a) < factor*min(dbp2a))
      cbu=[cbu p];
    end
  end
  % Set elements that will now be separated to zero in NotYetSeparated
  NotYetSeparated(cau,cbu)=zeros(length(cau),length(cbu));
  NotYetSeparated(cbu,cau)=zeros(length(cbu),length(cau)); 
  % Get sizes
  [M_CA,N_CA]=size(ClusterA); [M_cau,N_cau]=size(cau); 
  [M_CB,N_CB]=size(ClusterB); [M_cbu,N_cbu]=size(cbu);
  % Update cluster databases
  ClusterA=[[ClusterA zeros(M_CA,N_cau-N_CA)];[cau zeros(M_cau,N_CA-N_cau)]];
  ClusterB=[[ClusterB zeros(M_CB,N_cbu-N_CB)];[cbu zeros(M_cbu,N_CB-N_cbu)]];
end
anc=sum(sum([ClusterA ClusterB]>0))/N_T;






 