%HCLUST hierarchical clustering
% 
%  [LABELS, DENDROGRAM] = HCLUST(D,TYPE,K)
%   DENDROGRAM = HCLUST(D,TYPE)
% 
% INPUT
%  D     dissimilarity matrix
%  TYPE  string name of clustering criterion (optional)
%        - 's' or 'single'   : single linkage (default)
%        - 'c' or 'complete' : complete linkage
%        - 'a' or 'average'  : average linkage (weighted over cluster sizes)
%  K     number of clusters (optional)
%
% OUTPUT
%  LABELS       vector with labels
%  DENDROGRAM   matrix with dendrogram
%
% DESCRIPTION 
% Computation of cluster labels and a dendrogram between the clusters for
% objects with a given distance matrix D. K is the desired number of
% clusters. The dendrogram is a 2*K matrix. The first row yields all
% cluster sizes. The second row is the cluster level on which the set of
% clusters starting at that position is merged with the set of clusters
% just above it in the dendrogram. A dendrogram may be plotted by PLOTDG.
%
%   DENDROGRAM = HCLUST(D,TYPE)
%
% As in this case no clustering level is supplied, just the entire
% dendrogram is returned. The first row now contains the object indices.
%
% EXAMPLE
% a = gendats([25 25],20,5);    % 50 points in 20 dimensional feature space
% d = sqrt(distm(a));           % Euclidean distances
% dendg = hclust(d,'complete'); % dendrogram
% plotdg(dendg)
% lab = hclust(d,'complete',2); % labels
% confmat(lab,getlabels(a));    % confusion matrix
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% PLOTDG, PRKMEANS, KCENTRES, MODESEEK, EMCLUST

% 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: hclust.m,v 1.3 2008/02/14 11:54:43 duin Exp $

function [labels, dendrogram] = hclust(D,type,k)   
		
	if nargin < 3, k = []; end
	if nargin < 2, type = 's'; end
	[m,m1] = size(D);
	if m ~= m1
		error('Input matrix should be square')
	end
	D = D + diag(inf*ones(1,m));     % set diagonal at infinity.
	W = [1:m+1];               % starting points of clusters in linear object set.
	V = [1:m+2];               % positions of objects in final linear object set.
	F = inf * ones(1,m+1);     % distance of next cluster to previous cluster to be stored at first point of second cluster
	Z = ones(1,m); % number of samples in a cluster (only for average linkage)

  t = sprintf('Analysing %i cluster levels: ',m);
  prwaitbar(m,t);
	for n = 1:m-1
    prwaitbar(m,n,[t int2str(n)]);
		% find minimum distance D(i,j) i<j
		[di,I] = min(D); [dj,j] = min(di); i = I(j);
		if i > j, j1 = j; j = i; i = j1; end
		% combine clusters i,j
		switch type
		case {'s','single'}
			D(i,:) = min(D(i,:),D(j,:));
		case {'c','complete'}
			D(i,:) = max(D(i,:),D(j,:));
		case {'a','average'}
			D(i,:) = (Z(i)*D(i,:) + Z(j)*D(j,:))/(Z(i)+Z(j));
			Z(i:j-1) = [Z(i)+Z(j),Z(i+1:j-1)]; Z(j) = [];
		otherwise
			error('Unknown clustertype desired')
		end
		D(:,i) = D(i,:)';
		D(i,i) = inf;
		D(j,:) = []; D(:,j) = [];
		% store cluster distance
		F(V(j)) = dj;
		% move second cluster in linear ordering right after first cluster
		IV = [1:V(i+1)-1,V(j):V(j+1)-1,V(i+1):V(j)-1,V(j+1):m+1];
		W = W(IV); F = F(IV);   
		% keep track of object positions and cluster distances
		V = [V(1:i),V(i+1:j) + V(j+1) - V(j),V(j+2:m-n+3)];
  end
  prwaitbar(0);
	if ~isempty(k) | nargout == 2
		if isempty(k), k = m; end
		labels = zeros(1,m); 
		[S,J] = sort(-F);      % find cluster level
		I = sort(J(1:k+1));    % find all indices where cluster starts
		for i = 1:k            % find for all objects cluster labels
			labels(W(I(i):I(i+1)-1)) = i * ones(1,I(i+1)-I(i));
		end                    % compute dendrogram
		dendrogram = [I(2:k+1) - I(1:k); F(I(1:k))];
		labels = labels';
	else
		labels = [W(1:m);F(1:m)]; % full dendrogram
	end
return