function [centres, options, post, errlog] = kmeans(centres, data, options)%KMEANS	Trains a k means cluster model.%%	Description%	 CENTRES = KMEANS(CENTRES, DATA, OPTIONS) uses the batch K-means%	algorithm to set the centres of a cluster model. The matrix DATA%	represents the data which is being clustered, with each row%	corresponding to a vector. The sum of squares error function is used.%	The point at which a local minimum is achieved is returned as%	CENTRES.  The error value at that point is returned in OPTIONS(8).%%	[CENTRES, OPTIONS, POST, ERRLOG] = KMEANS(CENTRES, DATA, OPTIONS)%	also returns the cluster number (in a one-of-N encoding) for each%	data point in POST and a log of the error values after each cycle in%	ERRLOG.    The optional parameters have the following%	interpretations.%%	OPTIONS(1) is set to 1 to display error values; also logs error%	values in the return argument ERRLOG. If OPTIONS(1) is set to 0, then%	only warning messages are displayed.  If OPTIONS(1) is -1, then%	nothing is displayed.%%	OPTIONS(2) is a measure of the absolute precision required for the%	value of CENTRES at the solution.  If the absolute difference between%	the values of CENTRES between two successive steps is less than%	OPTIONS(2), then this condition is satisfied.%%	OPTIONS(3) is a measure of the precision required of the error%	function at the solution.  If the absolute difference between the%	error functions between two successive steps is less than OPTIONS(3),%	then this condition is satisfied. Both this and the previous%	condition must be satisfied for termination.%%	OPTIONS(14) is the maximum number of iterations; default 100.%%	See also%	GMMINIT, GMMEM%%	Copyright (c) Ian T Nabney (1996-2001)[ndata, data_dim] = size(data);[ncentres, dim] = size(centres);if dim ~= data_dim  error('Data dimension does not match dimension of centres')endif (ncentres > ndata)  error('More centres than data')end% Sort out the optionsif (options(14))  niters = options(14);else  niters = 100;endstore = 0;if (nargout > 3)  store = 1;  errlog = zeros(1, niters);end% Check if centres and posteriors need to be initialised from dataif (options(5) == 1)  % Do the initialisation  perm = randperm(ndata);  perm = perm(1:ncentres);  % Assign first ncentres (permuted) data points as centres  centres = data(perm, :);end% Matrix to make unit vectors easy to constructid = eye(ncentres);% Main loop of algorithmfor n = 1:niters  % Save old centres to check for termination  old_centres = centres;    % Calculate posteriors based on existing centres  d2 = dist2(data, centres);  % Assign each point to nearest centre  [minvals, index] = min(d2', [], 1);  post = id(index,:);  num_points = sum(post, 1);  % Adjust the centres based on new posteriors  for j = 1:ncentres    if (num_points(j) > 0)      centres(j,:) = sum(data(find(post(:,j)),:), 1)/num_points(j);    end  end  % Error value is total squared distance from cluster centres  e = sum(minvals);  if store    errlog(n) = e;  end  if options(1) > 0    fprintf(1, 'Cycle %4d  Error %11.6f\n', n, e);  end  if n > 1    % Test for termination    if max(max(abs(centres - old_centres))) < options(2) & ...        abs(old_e - e) < options(3)      options(8) = e;      return;    end  end  old_e = e;end% If we get here, then we haven't terminated in the given number of % iterations.options(8) = e;if (options(1) >= 0)  disp('Warning: Maximum number of iterations has been exceeded');end