Daniel@0: function Cd = som_cldist(D,clinds1,clinds2,cldist,q,mask)
Daniel@0: 
Daniel@0: % SOM_CLDIST Distances between two clusters.
Daniel@0: % 
Daniel@0: %   Cd = som_cldist(Md,c1,c2,'single')
Daniel@0: %   Cd = som_cldist(Md,c1,c2,'average')
Daniel@0: %   Cd = som_cldist(Md,c1,c2,'complete')
Daniel@0: %   Cd = som_cldist(Md,c1,c2,'neighf',H)
Daniel@0: %   Cd = som_cldist(Md,c1,[],...)
Daniel@0: %   Cd = som_cldist(D,c1,c2,'centroid',q,mask)
Daniel@0: %   Cd = som_cldist(D,c1,c2,'ward',q,mask)
Daniel@0: %   Cd = som_cldist(D,c1,[],...)
Daniel@0: %
Daniel@0: %  Input and output arguments ([]'s are optional):
Daniel@0: %   D        (matrix) size dlen x dim, the data set
Daniel@0: %            (struct) map or data struct
Daniel@0: %   Md       (matrix) size dlen x dlen, mutual distance matrix, see SOM_MDIST
Daniel@0: %   c1       (cell array) size n1 x 1, indices of clusters from which 
Daniel@0: %                     the distances should be calculated, each cell
Daniel@0: %                     contains indices of vectors that belong to that
Daniel@0: %                     cluster (indices are between 1...dlen)
Daniel@0: %   c2       (cell array) size n2 x 1, same as c1 but have the clusters
Daniel@0: %                     to which the distances should be calculated
Daniel@0: %            (empty)  c1 is used in place of c2
Daniel@0: %   [q]      (scalar) distance norm, default = 2
Daniel@0: %   [mask]   (vector) size dim x 1, the weighting mask, a vector of ones
Daniel@0: %                     by default
Daniel@0: %   H        (matrix) size dlen x dlen, neighborhood function values
Daniel@0: %
Daniel@0: %   Cd       (matrix) size n1 x n2, distances between the clusters
Daniel@0: %
Daniel@0: % See also SOM_MDIST. 
Daniel@0: 
Daniel@0: % Copyright (c) 2000 by Juha Vesanto
Daniel@0: % Contributed to SOM Toolbox on XXX by Juha Vesanto
Daniel@0: % http://www.cis.hut.fi/projects/somtoolbox/
Daniel@0:  
Daniel@0: % Version 2.0beta juuso 250800
Daniel@0: 
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: 
Daniel@0: [dlen dim] = size(D); 
Daniel@0: if nargin<5, q = 2; end
Daniel@0: if nargin<6, mask = ones(dim,1); end
Daniel@0: if ~iscell(clinds1), clinds1 = {clinds1}; end
Daniel@0: if ~isempty(clinds2) & ~iscell(clinds2), clinds2 = {clinds2}; end
Daniel@0: 
Daniel@0: n1 = length(clinds1); 
Daniel@0: n2 = length(clinds2); 
Daniel@0: if n2>0, Cd = zeros(n1,n2); else Cd = zeros(n1); end
Daniel@0: if n1==0, return; end
Daniel@0: 
Daniel@0: switch cldist, 
Daniel@0:   
Daniel@0:  % centroid distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Daniel@0:  case 'centroid',  
Daniel@0: 
Daniel@0:   C1 = zeros(n1,dim); for i=1:n1, C1(i,:) = mean(D(clinds1{i},:),1); end
Daniel@0:   C2 = zeros(n2,dim); for i=1:n2, C2(i,:) = mean(D(clinds2{i},:),1); end
Daniel@0:   if n2==0, 
Daniel@0:     for i=1:n1-1, 
Daniel@0:       for j=i+1:n1, 
Daniel@0: 	diff = C1(i,:)-C1(j,:); 
Daniel@0: 	switch q, 
Daniel@0: 	 case 1,    Cd(i,j)=abs(diff)*mask;
Daniel@0: 	 case 2,    Cd(i,j)=sqrt((diff.^2)*mask);  
Daniel@0: 	 case Inf,  Cd(i,j)=max(diag(mask)*abs(diff),[],2);
Daniel@0: 	 otherwise, Cd(i,j)=((abs(diff).^q)*mask).^(1/q);
Daniel@0: 	end   
Daniel@0:       end
Daniel@0:       Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
Daniel@0:     end
Daniel@0:   else
Daniel@0:     for i=1:n1, 
Daniel@0:       for j=1:n2, 
Daniel@0: 	diff = C1(i,:)-C2(j,:); 
Daniel@0: 	switch q, 
Daniel@0: 	 case 1,    Cd(i,j)=abs(diff)*mask;
Daniel@0: 	 case 2,    Cd(i,j)=sqrt((diff.^2)*mask);  
Daniel@0: 	 case Inf,  Cd(i,j)=max(diag(mask)*abs(diff),[],2);
Daniel@0: 	 otherwise, Cd(i,j)=((abs(diff).^q)*mask).^(1/q);
Daniel@0: 	end   
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0: 
Daniel@0:  % ward distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:  case 'ward',
Daniel@0: 
Daniel@0:   C1 = zeros(n1,dim); nn1 = zeros(n1,dim); 
Daniel@0:   for i=1:n1, C1(i,:) = mean(D(clinds1{i},:),1); nn1(i) = length(clinds1{i}); end
Daniel@0:   C2 = zeros(n2,dim); nn2 = zeros(n2,dim); 
Daniel@0:   for i=1:n2, C2(i,:) = mean(D(clinds2{i},:),1); nn2(i) = length(clinds2{i}); end
Daniel@0:   if n2==0, 
Daniel@0:     for i=1:n1-1, 
Daniel@0:       for j=i+1:n1, 
Daniel@0: 	diff = C1(i,:) - C1(j,:); 
Daniel@0: 	f = 2*nn1(i)*nn1(j) / (nn1(i)+nn1(j)); 
Daniel@0: 	switch q, 
Daniel@0: 	 case 1,    Cd(i,j)=f*abs(diff)*mask;
Daniel@0: 	 case 2,    Cd(i,j)=f*sqrt((diff.^2)*mask);  
Daniel@0: 	 case Inf,  Cd(i,j)=f*max(diag(mask)*abs(diff),[],2);
Daniel@0: 	 otherwise, Cd(i,j)=f*((abs(diff).^q)*mask).^(1/q);
Daniel@0: 	end   
Daniel@0:       end
Daniel@0:       Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
Daniel@0:     end
Daniel@0:   else
Daniel@0:     for i=1:n1, 
Daniel@0:       for j=1:n2, 
Daniel@0: 	diff = C1(i,:) - C2(j,:); 
Daniel@0: 	f = 2*nn1(i)*nn2(j) / (nn1(i)+nn2(j)); 
Daniel@0: 	switch q, 
Daniel@0: 	 case 1,    Cd(i,j)=f*abs(diff)*mask;
Daniel@0: 	 case 2,    Cd(i,j)=f*sqrt((diff.^2)*mask);  
Daniel@0: 	 case Inf,  Cd(i,j)=f*max(diag(mask)*abs(diff),[],2);
Daniel@0: 	 otherwise, Cd(i,j)=f*((abs(diff).^q)*mask).^(1/q);
Daniel@0: 	end   
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end  
Daniel@0: 
Daniel@0:  % single linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:  case 'single',
Daniel@0: 
Daniel@0:   if n2==0, 
Daniel@0:     for i=1:n1-1, 
Daniel@0:       for j=i+1:n1, 
Daniel@0: 	vd = D(clinds1{i},clinds1{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), Cd(i,j) = min(vd(fi)); else Cd(i,j) = Inf; end
Daniel@0:       end
Daniel@0:       Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
Daniel@0:     end
Daniel@0:   else
Daniel@0:     for i=1:n1, 
Daniel@0:       for j=1:n2, 
Daniel@0: 	vd = D(clinds1{i},clinds2{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), Cd(i,j) = min(vd(fi)); else Cd(i,j) = Inf; end
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0: 
Daniel@0:  % average linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0:  case 'average',
Daniel@0:   
Daniel@0:   if n2==0, 
Daniel@0:     for i=1:n1-1, 
Daniel@0:       for j=i+1:n1, 
Daniel@0: 	vd = D(clinds1{i},clinds1{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), Cd(i,j) = mean(vd(fi)); else Cd(i,j) = Inf; end
Daniel@0:       end
Daniel@0:       Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
Daniel@0:     end
Daniel@0:   else
Daniel@0:     for i=1:n1, 
Daniel@0:       for j=1:n2, 
Daniel@0: 	vd = D(clinds1{i},clinds2{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), Cd(i,j) = mean(vd(fi)); else Cd(i,j) = Inf; end
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0:     
Daniel@0:  % complete linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Daniel@0:  case 'complete',
Daniel@0:  
Daniel@0:    if n2==0, 
Daniel@0:     for i=1:n1-1, 
Daniel@0:       for j=i+1:n1, 
Daniel@0: 	vd = D(clinds1{i},clinds1{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), Cd(i,j) = max(vd(fi)); else Cd(i,j) = Inf; end
Daniel@0:       end
Daniel@0:       Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
Daniel@0:     end
Daniel@0:   else
Daniel@0:     for i=1:n1, 
Daniel@0:       for j=1:n2, 
Daniel@0: 	vd = D(clinds1{i},clinds2{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), Cd(i,j) = max(vd(fi)); else Cd(i,j) = Inf; end
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0:  
Daniel@0:  % neighborhood function linkage distance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Daniel@0:  case 'neighf',
Daniel@0:   
Daniel@0:   if n2==0, 
Daniel@0:     for i=1:n1-1, 
Daniel@0:       for j=i+1:n1, 
Daniel@0: 	vd = D(clinds1{i},clinds1{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), 
Daniel@0: 	  hd = q(clinds1{i},clinds1{j}); 
Daniel@0: 	  hd = hd(fi); 
Daniel@0: 	  Cd(i,j) = sum(hd.*vd(fi))/sum(hd); 	  
Daniel@0: 	else Cd(i,j) = Inf; 
Daniel@0: 	end
Daniel@0:       end
Daniel@0:       Cd([(i+1):n1],i) = Cd(i,[(i+1):n1])';
Daniel@0:     end
Daniel@0:   else
Daniel@0:     for i=1:n1, 
Daniel@0:       for j=1:n2, 
Daniel@0: 	vd = D(clinds1{i},clinds2{j}); 
Daniel@0: 	fi = isfinite(vd(:));
Daniel@0: 	if any(fi), 
Daniel@0: 	  hd = q(clinds1{i},clinds2{j}); 
Daniel@0: 	  hd = hd(fi); 
Daniel@0: 	  Cd(i,j) = sum(hd.*vd(fi))/sum(hd); 	  
Daniel@0: 	else Cd(i,j) = Inf; 
Daniel@0: 	end
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0: 
Daniel@0:  otherwise, error(['Unknown cluster distance metric: ' cldist]); 
Daniel@0: end
Daniel@0:   
Daniel@0: return;
Daniel@0:   
Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniel@0: