%KCENTRES Finds K center objects from a distance matrix
% 
%  [LAB,J,DM] = KCENTRES(D,K,N)
% 
% INPUT
%   D    Distance matrix between, e.g. M objects (may be a dataset)
%   K    Number of center objects to be found (optional; default: 1)
%   N    Number of trials starting from a random initialization
%        (optional; default: 1)
%
% OUTPUT
%   LAB  Integer labels: each object is assigned to its nearest center
%   J    Indices of the center objects
%   DM   A list of distances corresponding to J: for each center in J 
%        the maximum distance of the objects assigned to this center.	
%
% DESCRIPTION  
% Finds K center objects from a symmetric distance matrix D. The center 
% objects are chosen from all M objects such that the maximum of the 
% distances over all objects to the nearest center is minimized. For K > 1,
% the results depend on a random initialization. The procedure is repeated 
% N times and the best result is returned. 
%
% If N = 0, initialisation is not random, but done by a systematic
% selection based on a greedy approach.
% 
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% HCLUST, PRKMEANS, EMCLUST, MODESEEK

% 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

% $Id: kcentres.m,v 1.5 2008/07/03 09:08:43 duin Exp $

function [labels,Jopt,dm] = kcentres(d,k,n)
	
	if (nargin < 3) | isempty(n)
      n = 1;
      prwarning(4,'Number of trials not supplied, assuming one.');
   end

	if (nargin < 2) | isempty(k)
      k = 1;
      prwarning(4,'Number of centers not supplied, assuming one.');
   end
	
	if(isdataset(d))
		d = +d;
		prwarning(4,'Distance matrix is convert to double.');
	end

	[m,m2] = size(d);
	if ( ~issym(d,1e-12) | any(diag(d) > 1e-14) )
		error('Distance matrix should be symmetric and have zero diagonal')
	end

% checking for a zero diagonal 
	t = eye(m) == 1;	
	if(~all(d(t)==0))
		error('Distance matrix should have a zero diagonal.')
	end

	if (k == 1)
		dmax = max(d);
		[dm,Jopt] = min(dmax);
		labels = repmat(1,m,1);
		return;
	end

	if k > m
		error('Number of centres should not exceed number of objects')
	end
	
	% We are here only if K (> 1) centers are to be found.
	% Loop over number of trials.
	dmax = max(max(d));
	dopt = inf;
	s = sprintf('k-centres, %i attempts: ',n);
	prwaitbar(n,s,n>1);
	if n == 0
		nrep = 1; 
	else 
		nrep = n; 
	end
	for tri = 1:nrep
		prwaitbar(n,tri,[s int2str(tri)]);
		if n == 0
			M = kcentresort(d,k);           % systematic initialisation
		else
			M = randperm(m); M = M(1:k);		% Random initializations
		end
 		J = zeros(1,k); 					    	  % Center objects to be found.

		% Iterate until J == M. See below.
		while 1,
			[dm,I] = min(d(M,:));

			% Find K centers. 
			for i = 1:k                   
				JJ = find(I==i); 
				if (isempty(JJ))
					%JJ can be empty if two or more objects are in the same position of 
					% feature space in dataset
					J(i) = 0;									
				else
						% Find objects closest to the object M(i)
					[dummy,j,dm] = kcentres(d(JJ,JJ),1,1);
					J(i) = JJ(j);
				end
			end

			J(find(J==0)) = [];
			k = length(J);
			if k == 0
				error('kcentres fails as some objects are identical: add some noise')
			end
			if (length(M) == k) & (all(M == J))
				% K centers are found.
				[dmin,labs] = min(d(J,:));
				dmin = max(dmin);
				break;
			end
			M = J;
		end

		% Store the optimal centers in JOPT.
	
		if (dmin <= dopt)
			dopt = dmin;
			labels = labs';
			Jopt = J;
		end
	
	end
	prwaitbar(0)
	
	% Determine the best centers over the N trials.
	dm = zeros(1,k);   
	for i=1:k
		L = find(labels==i);
		dm(i) = max(d(Jopt(i),L));
	end

return;

%KCENTRESORT Sort objects given by dissimilarity matrix
%
%               N = KCENTRESORT(D,P,CRIT)
%
% INPUT
%   D     Square dissimilarity matrix, zeros on diagonal
%   P     Number of prototypes to be selected
%   CRIT  'dist' or 'centre'
%
% OUTPUT
%   N     Indices of selected prototypes
%
% DESCRIPTION
% Sort objects given by square dissim matrix D using a greedy approach
% such that the maximum NN distance from all objects (prototypes)
% to the first K: max(min(D(:,N(1:K),[],2)) is minimized.
%
% This routines tries to sample the objects such that they are evenly
% spaced judged from their dissimilarities. This may be used as
% initialisation in KCENTRES. It works reasonably, but not very good.
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% KCENTRES

% 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 N = kcentresort(d,p,crit);

d = +d;
m = size(d,1);
if nargin < 3, crit = 'dist'; end
if nargin < 2 | isempty(p), p = m; end
L = [1:m];
N = zeros(1,p);
[dd,n] = min(max(d,[],2));  % this is the first (central) prototype
e = d(:,n);                 % store here the distances to the nearest prototype (dNNP)
f = min(d,repmat(e,1,m));   % replace distances that are larger than dNNP by dNNP
N(1) = n;                   % ranking of selected prototypes
L(n) = [];                  % candidate prototypes (all not yet selected objects)

for k=2:p                   % extend prototype set
        if strcmp(crit,'centre')
        [dd,n] = min(max(f(L,L),[],1)); % This selects the next prototype out of candidates in L
        e = min([d(:,L(n)) e],[],2);    % update dNNP
        f = min(d,repmat(e,1,m));       % update replacement of distances that are larger
                                              % than dNNP by dNNP
        elseif strcmp(crit,'dist')
                [dd,n] = max(mean(d(N(find(N > 0)),L)));
        else
                error('Illegal crit')
        end
  N(k) = L(n);                    % update list of selected prototypes
  L(n) = [];                      % update list of candidate prototypes
end