wolffd@0: function [centers,clusters,errors,ind] = kmeans_clusters(sD, n_max, c_max, verbose)
wolffd@0: 
wolffd@0: % KMEANS_CLUSTERS Clustering with k-means with different values for k.
wolffd@0: %
wolffd@0: % [c, p, err, ind] = kmeans_clusters(sD, [n_max], [c_max], [verbose])
wolffd@0: %
wolffd@0: %   [c, p, err, ind] = kmeans_clusters(sD);
wolffd@0: %  
wolffd@0: %  Input and output arguments ([]'s are optional):
wolffd@0: %   D         (struct) map or data struct
wolffd@0: %             (matrix) size dlen x dim, the data 
wolffd@0: %   [n_max]   (scalar) maximum number of clusters, default is sqrt(dlen)
wolffd@0: %   [c_max]   (scalar) maximum number of k-means runs, default is 5
wolffd@0: %   [verbose] (scalar) verbose level, 0 by default
wolffd@0: %
wolffd@0: %   c         (cell array) c{i} contains cluster centroids for k=i
wolffd@0: %   p         (cell array) p{i} contains cluster indeces for k=i
wolffd@0: %   err       (vector) squared sum of errors for each value of k
wolffd@0: %   ind       (vector) Davies-Bouldin index value for each clustering
wolffd@0: %
wolffd@0: % Makes a k-means to the given data set with different values of
wolffd@0: % k. The k-means is run multiple times for each k, and the best of
wolffd@0: % these is selected based on sum of squared errors. Finally, the
wolffd@0: % Davies-Bouldin index is calculated for each clustering. 
wolffd@0: %
wolffd@0: % For example to cluster a SOM: 
wolffd@0: %    [c, p, err, ind] = kmeans_clusters(sM); % find clusterings
wolffd@0: %    [dummy,i] = min(ind); % select the one with smallest index
wolffd@0: %    som_show(sM,'color',{p{i},sprintf('%d clusters',i)}); % visualize
wolffd@0: %    colormap(jet(i)), som_recolorbar % change colormap
wolffd@0: %  
wolffd@0: % See also SOM_KMEANS.
wolffd@0: 
wolffd@0: % References: 
wolffd@0: %   Jain, A.K., Dubes, R.C., "Algorithms for Clustering Data", 
wolffd@0: %   Prentice Hall, 1988, pp. 96-101.
wolffd@0: %
wolffd@0: %   Davies, D.L., Bouldin, D.W., "A Cluster Separation Measure", 
wolffd@0: %   IEEE Transactions on Pattern Analysis and Machine Intelligence, 
wolffd@0: %   vol. PAMI-1, no. 2, 1979, pp. 224-227.
wolffd@0: %
wolffd@0: %   Vesanto, J., Alhoniemi, E., "Clustering of the Self-Organizing
wolffd@0: %   Map", IEEE Transactions on Neural Networks, 2000.
wolffd@0: 
wolffd@0: % Contributed to SOM Toolbox vs2, February 2nd, 2000 by Esa Alhoniemi
wolffd@0: % Copyright (c) by Esa Alhoniemi
wolffd@0: % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0: 
wolffd@0: % ecco 301299 juuso 020200 211201
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% input arguments and initialization
wolffd@0: 
wolffd@0: if isstruct(sD), 
wolffd@0:   if isfield(sD,'data'), D = sD.data; 
wolffd@0:   else D = sD.codebook; 
wolffd@0:   end
wolffd@0: else D = sD; 
wolffd@0: end
wolffd@0: [dlen dim] = size(D);
wolffd@0: 
wolffd@0: if nargin < 2 | isempty(n_max) | isnan(n_max), n_max = ceil(sqrt(dlen)); end
wolffd@0: if nargin < 3 | isempty(c_max) | isnan(c_max), c_max = 5; end
wolffd@0: if nargin < 4 | isempty(verbose) | isnan(verbose), verbose = 0; end
wolffd@0: 
wolffd@0: centers   = cell(n_max,1); 
wolffd@0: clusters  = cell(n_max,1);
wolffd@0: ind       = zeros(1,n_max)+NaN;
wolffd@0: errors    = zeros(1,n_max)+NaN;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% action
wolffd@0: 
wolffd@0: % the case k=1 is trivial, but Davies-Boulding index cannot be evaluated
wolffd@0: m = zeros(1,dim);
wolffd@0: for i=1:dim, m(i)=mean(D(isfinite(D(:,i)),i)); end
wolffd@0: centers{1} = m;
wolffd@0: clusters{1} = ones(dlen,1);
wolffd@0: [dummy qerr] = som_bmus(m,D);
wolffd@0: errors(1) = sum(qerr.^2);
wolffd@0: ind(1) = NaN; 
wolffd@0: 
wolffd@0: if verbose, fprintf(2,'Doing k-means for 2-%d clusters\n',n_max); end
wolffd@0: 
wolffd@0: for i = 2:n_max, % number of clusters
wolffd@0: 
wolffd@0:   % make k-means with k=i for c_max times and select the best based
wolffd@0:   % on sum-of-squared errors (SSE)
wolffd@0:   best = realmax;  
wolffd@0:   for j = 1:c_max     % run number j for cluster i      
wolffd@0:     if verbose,
wolffd@0:       fprintf('%d/%d clusters, k-means run %d/%d\r', i, n_max,j, c_max);
wolffd@0:     end      
wolffd@0:     [c, k, err] = som_kmeans('batch', D, i, 100, 0);
wolffd@0:     if err < best, k_best = k'; c_best = c; best = err; end
wolffd@0:     % ' added in k_best = k'; by kr 1.10.02
wolffd@0:   end
wolffd@0:   if verbose, fprintf(1, '\n');  end
wolffd@0: 
wolffd@0:   % store the results  
wolffd@0:   centers{i}  = c_best;
wolffd@0:   clusters{i} = k_best;
wolffd@0:   errors(i)   = best;
wolffd@0: %  ind(i)      = db_index(D, c_best, k_best, 2); wrong version in somtbx ??
wolffd@0:   ind(i)      = db_index(D, k_best, c_best, 2); % modified by kr 1.10.02
wolffd@0: 
wolffd@0:   % if verbose mode, plot the index & SSE
wolffd@0:   if verbose
wolffd@0:     subplot(2,1,1), plot(ind), grid
wolffd@0:     title('Davies-Bouldin''s index')
wolffd@0:     subplot(2,1,2), plot(errors), grid
wolffd@0:     title('SSE')
wolffd@0:     drawnow
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: return; 
wolffd@0: 
wolffd@0: 
wolffd@0: