function [center, U, distortion] = kmeans(dataSet, clusterNum, plotOpt)
%KMEANS K-means clustering using Forgy's batch-mode method
%	Usage: [center, U, distortion] = KMEANS(dataSet, clusterNum)
%		dataSet: data set to be clustered; where each row is a sample data
%		clusterNum: number of clusters (greater than one)
%		center: final cluster centers, where each row is a center
%		U: final fuzzy partition matrix (or MF matrix)
%		distortion: values of the objective function during iterations 
%
%	Type "kmeans" for a self demo.

%	Roger Jang, 20030330

if nargin==0, selfdemo; return; end
if nargin<3, plotOpt=0; end

maxLoopCount = 100;			% Max. iteration
distortion = zeros(maxLoopCount, 1);	% Array for objective function
center = initCenter(clusterNum, dataSet)	% Initial cluster centers

if plotOpt
	plot(dataSet(:,1), dataSet(:,2), 'b.');
	centerH=line(center(:,1), center(:,2), 'color', 'r', 'marker', 'o', 'linestyle', 'none', 'linewidth', 2);
	axis image
end;

% Main loop
distMat=vecdist(center, dataSet);
for i = 1:maxLoopCount,
	[center, distortion(i), distMat, U] = updateCenter(center, dataSet, distMat);
	center
	fprintf('Iteration count = %d, distortion = %f\n', i, distortion(i));
	if plotOpt, 
		set(centerH, 'xdata', center(:,1), 'ydata', center(:,2));
		drawnow;
	end
	% check termination condition
	if i > 1,
		if abs(distortion(i-1) - distortion(i))/distortion(i-1) < eps, break; end,
	end
end
loopCount = i;	% Actual number of iterations 
distortion(loopCount+1:maxLoopCount) = [];

if plotOpt
	color = {'r', 'g', 'c', 'y', 'm', 'b', 'k'};
	figure;
	plot(dataSet(:, 1), dataSet(:, 2), 'o');
	maxU = max(U);
	clusterNum = size(center,1);
	for i=1:clusterNum,
		index = find(U(i, :) == maxU);
		colorIndex = rem(i, length(color))+1;  
		line(dataSet(index, 1), dataSet(index, 2), 'linestyle', 'none', 'marker', '*', 'color', color{colorIndex});
		line(center(:,1), center(:,2), 'color', 'r', 'marker', 'o', 'linestyle', 'none', 'linewidth', 2);
	end
	axis image;
end


% ========== subfunctions ==========
% ====== Find the initial centers
function center = initCenter(clusterNum, dataSet, method)
if nargin<3; method=4; end
switch method
	case 1
		% ====== Method 1: Randomly pick clusterNum data points as cluster centers
		dataNum = size(dataSet, 1);
		tmp = randperm(dataNum);
		index = tmp(1:clusterNum);
		center = dataSet(index, :);
	case 2
		% ====== Method 2: Choose clusterNum data points closest to mean vector
		meanVec = mean(dataSet);
		distMat = vecdist(meanVec, dataSet);
		[a,b] = sort(distMat);
		center = dataSet(b(1:clusterNum), :);
	case 3
		% ====== Method 3: Choose clusterNum data points furthest to the mean vector
		meanVec = mean(dataSet);
		distMat = vecdist(meanVec, dataSet);
		[a,b] = sort(distMat);
		b = fliplr(b);
		center = dataSet(b(1:clusterNum), :);
	case 4
		% ====== Method 4: ㄏノ戈