%MODECLUST KNN mode-seeking clustering
% 
% 	LAB = MODECLUST(A,K)
% 
% INPUT
%   A       Dataset
%   K       Number of neighbours to search for local mode
%
% OUTPUT
  %   LAB   Indices of modal samples
%
% DESCRIPTION
% A K-NN modeseeking method is used to assign objects to their nearest mode.
% Object densities are defined by one over the distance to the K-th nearest 
% neighbour. Clusters are defined by recursively jumping for every object 
% to object with the highest density in the local neighborhood.
%
% K can also be a vector of neighborhood sizes, which is much faster.
% Default K: a set of values determined by a geometric series.
%
% EXAMPLE
% delfigs
% a = gendatm(1000);      % generate 1000 objects in 8 
% lab = modeclust(a); 
% for j=1:size(lab,2)
%   nclust = numel(unique(lab(:,j)));
%   if nclust < 20 & nclust > 1
%     figure; scatterd(prdataset(a,lab(:,j)));
%     title(['Number of clusters: ' int2str(nclust)]);
%   end
% end
% showfigs
% 
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% MAPPINGS, DATASETS, PRKMEANS, HCLUST, KCENTRES, PROXM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function labels = modeclust(a,K)

	% prepare data
  a = +a;                 % get rid of possible PRTools dataset structure
	m = size(a,1);          % number of objects
	if (nargin < 2  | isempty(K))
    K = intgeos(2,1.2,m); % default set of nb-sizes (neighborhood-sizes)
  end
  nk = numel(K);          % number of neighborhood sizes (clusterings)
	f = zeros(m,nk);        % space for storing densities
  batchsize = 100;        % batch sizes (avoid large distance matrices)
  nbatch = ceil(m/batchsize); % number of batches
  
  % first loop over the data: compute densities for every nb-size
	t = sprintf('Computing %i densities: ',m);
	prwaitbar(m,t);
  for j=1:nbatch
		prwaitbar(m,j*batchsize,[t num2str(j*batchsize)]);
    batch = (j-1)*batchsize+1:j*batchsize;  % take a subset of the objects
    if any(batch>m)              % arrange last batch
      batch = batch(batch <= m);
    end
		prwaitbar(m,(j-1)*batchsize,[t num2str((j-1)*batchsize)]);
    dj = distm(a(batch,:),a);    % distances between batch and all data
    [dd,L] = sort(dj,2);         % sort them per object
    f(batch,:) = 1./dd(:,K);     % density estimates for all nb-sizes
  end
  prwaitbar(0);
  
  % second loop: find indices of maximum density point in every nb
  NN = zeros(m,nk);              % space for storing these indices
	t = sprintf('Checking %i neighborhoods: ',m);
	prwaitbar(m,t);
  for j=1:nbatch
		prwaitbar(m,j*batchsize,[t num2str(j*batchsize)]);
    batch = (j-1)*batchsize+1:j*batchsize;   % take a subset of the objects
    if any(batch>m)              % arrange last batch
      batch = batch(batch <= m)';
    end
    mj = numel(batch);           % number of data points in batch
    dj = distm(a(batch,:),a);    % recomputes distances
    [dd,J] = sort(dj,2);         % sort them per object
    for n=1:nk                   % run over all nb's
      ff = f(:,n);               % densities for this nb-size
      [dummy,I] = max(reshape(ff(J(:,1:K(n))),mj,K(n)),[],2);
      NN(batch,n)  = J([1:mj]'+(I-1)*mj); % for all points the index of
    end                          %  point with maximum local density 
  end
  prwaitbar(0);
  
  % third loop: follow gradient to density modes
	t = sprintf('Computing %i clusterings: ',nk);
	prwaitbar(nk,t);
  for n=1:nk
		prwaitbar(nk,n,[t num2str(n)]);
		N = NN(:,n);                 % initial nb points with maximum density
		% Re-assign samples to the sample their nearest neighbour is assigned to.
		% Iterate until assignments don't change anymore. Samples that then point 
		% to themselves are modes; all other samples point to the closest mode.
		M = N(N);
    while (any(M~=N))
			N = M; M = N(N);
    end
    labels(:,n) = M;  % labels are the indices of the modes to which data 
		                  % points are assigned to. 
    % Data points assigned to the same mode constitute a cluster
	end
	prwaitbar(0);

return