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