wolffd@0: function [t,r,Cd,S] = som_gapindex(sM, base, between)
wolffd@0:  
wolffd@0: % SOM_GAPINDEX Gap clustering evaluation index.
wolffd@0: %
wolffd@0: % [t,r] = som_gapindex(sM, base, [between])
wolffd@0: %
wolffd@0: %  Input and output arguments ([]'s are optional):  
wolffd@0: %    sM        (struct) map struct
wolffd@0: %    base      (vector) clusters indeces for each map unit, map units
wolffd@0: %                       with index<=0 or NaN are not taken into account
wolffd@0: %    [between] (vector) indices of prototypes which are "between" clusters:
wolffd@0: %                       the associated distances are doubled
wolffd@0: % 
wolffd@0: %    t         (scalar) Gap index index for the clustering (=mean(r))
wolffd@0: %    r         (vector) maximum Gap index for each cluster (size max(base) x 1)    
wolffd@0: % 
wolffd@0: % See also  KMEANS, KMEANS_CLUSTERS, SOM_GAPINDEX.
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0: if nargin<3, between = find(isnan(base)); end
wolffd@0: 
wolffd@0: nc = max(base); 
wolffd@0: cinds = cell(nc,1); 
wolffd@0: for i=1:nc, cinds{i} = find(base==i); end
wolffd@0:  
wolffd@0: % distances between neighboring prototypes
wolffd@0: Ne = som_neighbors(sM,'N1'); 
wolffd@0: Md = som_mdist(sM.codebook,2,[],Ne);
wolffd@0: Md(Ne==0) = NaN;
wolffd@0: 
wolffd@0: Md(between,:) = Md(between,:)*2;
wolffd@0: Md(:,between) = Md(:,between)*2;
wolffd@0: Md(between,between) = Md(between,between)/2;
wolffd@0:  
wolffd@0: % dispersion in each cluster 
wolffd@0: S = zeros(nc,1);
wolffd@0: for i=1:nc, 
wolffd@0:     inds = setdiff(cinds{i},between);    
wolffd@0:     if any(inds), 
wolffd@0:         indist = Md(inds,inds); 
wolffd@0:         for j=1:size(indist,1), indist(j,j) = NaN; end
wolffd@0:         indist = indist(isfinite(indist(:))); 
wolffd@0:         if any(indist), S(i) = mean(indist); end
wolffd@0:     end
wolffd@0: end 
wolffd@0:  
wolffd@0: % distances between clusters
wolffd@0: Cd = zeros(nc,nc) + NaN;
wolffd@0: for i=1:nc,
wolffd@0:     inds1 = cinds{i}; 
wolffd@0:     for j=1:nc, 
wolffd@0:         inds2 = cinds{j}; 
wolffd@0: 	od = Md(inds1,inds2); 
wolffd@0: 	od = od(isfinite(od(:)));
wolffd@0: 	if any(od), Cd(i,j) = mean(od(:)); end        
wolffd@0:     end    
wolffd@0: end
wolffd@0: 
wolffd@0: % Gap index
wolffd@0: R = NaN * zeros(nc);
wolffd@0: for i = 1:nc
wolffd@0:   for j = i+1:nc
wolffd@0:     R(i,j) = (S(i) + S(j))/Cd(i,j);
wolffd@0:     R(j,i) = R(i,j); 
wolffd@0:   end
wolffd@0: end
wolffd@0: r = max(R,[],2);
wolffd@0:  
wolffd@0: t = mean(r(isfinite(r)));
wolffd@0:  
wolffd@0: return;