%  curCenters = kmeans(k, data, reduce, accuracy, maxIterations)

%

%    Run (enhanced) K-means clustering algorithm over a given set of data.

%

%    Input:   k             =  num of cluster centers

%             data          =  dataset

%             reduce        =  numeral indicating the ratio required to reduce "unjustified" centers 

%                              (zero for no reduction)

%             accuracy      =  required accuracy of equality between to subsequent iterations 

%                              (zero for equal)

%             maxIterations =  maximum count of iterations performed by algorithm (zero of no-limit)

%

%    Output:  curCenters    =  k (or less) centers found by the k-means algorithm

%



function curCenters = kmeans(k, data, reduce, accuracy, maxIterations)



% check arguments

if (nargin ~= 5)

   error('kmeans: invalid number of arguments');

end



% utility: init indices vector (for future use)

for i = 1 : k

   indices(i) = i;

end



% find data vectors dimension and count, values boundaries within data

[vecsNum, vecsDim] = size(data);

dataValMax   = max(data);

dataValMin   = min(data);

dataValRange = dataValMax - dataValMin;



% initiate a random matrix of k centers for the clusters

newCenters = (ncat(dataValRange, k) .* rand(k, vecsDim)) + ncat(dataValMin, k);



curIteration  = 0;

isDone = 0;

% iterate to improve the centers

while (~isDone)

   % acquire new centers to be current centers

   curCenters = newCenters;

   

   % calculate the MSE between each center and each vector of data

   % step-1: concatenate cluster centers and data vectors, for better performance

   tmpCenters = ncat(curCenters, vecsNum);

   tmpData    = ndup(data, k);

   % step-2: calculate MSE of distance for each center and data vector

   distVec = mserr(tmpData, tmpCenters, 2);

   % step-3: put distances vector into matrix

   distMtx = reshape(distVec, k, vecsNum);

   

   % find location of minimum value within distances matrix

   [distMinVal, distMinIndx] = min(distMtx);

   

   % calculate new centers

   for i = 1 : k

      % get the number of vectors in each cluster

      numVecsPerClust(i) = sum(distMinIndx == i);

      % if cluster contains vectors, recalculate its center

      if (numVecsPerClust(i) > 0)

         % improve the center, to be the mean of vectors belonging to it

         newCenters(i, :) = mean(data((distMinIndx == i), :), 1);

         % calculate square root of MSE of distances from cluster center

         clustMeanDist(i, :) = mserr(data((distMinIndx == i), :), ncat(newCenters(i, :), numVecsPerClust(i))) .^ 0.5;

      end

   end

   

   % wishful assumption: no more iterations required...

   isDone = 1;

   

   % remove empty clusters, if any

   numEmptyClusters = sum(numVecsPerClust == 0);

   if (numEmptyClusters > 0)

      newCenters = newCenters((numVecsPerClust > 0), :);

      clustMeanDist = clustMeanDist((numVecsPerClust > 0), :);

      numVecsPerClust = numVecsPerClust((numVecsPerClust > 0));

      k = k - numEmptyClusters;

      

      % another iteration is required...

      isDone = 0;

   end

      

   if ((curIteration == 0) | (reduce == 0))

      % case-1: first iteration, or reduction not allowed

      

      % for each empty cluster, random choose a data vector to replace its previous center

      if (numEmptyClusters > 0)

         newCenters = [newCenters; data(ceil(rand(1, numEmptyClusters) * vecsNum), :)];

         k = k + numEmptyClusters;

      end

      

      % another iteration required...

      isDone = 0;

      

   elseif ((reduce > 0) & (k > 1))

      % case-2: advanced iterations, reduction factor defined

      

      % update utility indices

      indices = indices(1 : k);

   

      % calculate "justification ratio" for each cluster center

      justRatio = zeros(k, vecsDim);

      for i = 1 : k

         otherCenters = newCenters((indices ~= i), :);

         clustOthersDist = otherCenters - ncat(newCenters(i, :), (k - 1));

         [clustMinOthersDistVal, clustMinOthersDistIndx] = min(mserr(otherCenters, ncat(newCenters(i, :), (k - 1)), 2));

         clustMinOthersDist = clustOthersDist(clustMinOthersDistIndx, :);

         

         if (numVecsPerClust(i) > 1)

            justRatio(i, :) = clustMinOthersDist ./ clustMeanDist(i, :);

         else

            justRatio(i, :) = clustMinOthersDist;

         end

      end

      

      % eliminate zero justification ratios for single-vector clusters

      numSingleVecClusts = sum((numVecsPerClust == 1));

      if (numSingleVecClusts > 0)

         [maxJustCenterVal, maxJustCenterIndx] = max(mean((justRatio .^ 2), 2));

         justRatio((numVecsPerClust == 1), :) = ncat(justRatio(maxJustCenterIndx, :), numSingleVecClusts);

      end

      

      % now, choose center with smallest justification

      [minJustCenterVal, minJustCenterIndx] = min(mean((justRatio .^ 2), 2));

      minJustCenter = newCenters(minJustCenterIndx, :);

      

      % if center is "not justified", reduce it

      i = minJustCenterIndx;

      otherCenters = newCenters((indices ~= i), :);

      clustOthersDist = otherCenters - ncat(minJustCenter, (k - 1));

      [clustMinOthersDistVal, clustMinOthersDistIndx] = min(mserr(otherCenters, ncat(minJustCenter, (k - 1)), 2));

      clustMinOthersDist = clustOthersDist(clustMinOthersDistIndx, :);

      

      if (abs(clustMinOthersDist) < (clustMeanDist(i, :) .* reduce))

         isDone = 0;

         % cluster does not satisfy requirements, reduce it as required

         newCenters = otherCenters;

         clustMeanDist = clustMeanDist((indices ~= i), :);

         numVecsPerClust = numVecsPerClust((indices ~= i));

         k = k - 1;

      end

   end

   

   if (isDone & (max(abs(curCenters - newCenters)) > accuracy))

      isDone = 0;

   end

   

   % advance iteration count, check against maximum count

   % (NOTE: in case of zero supplied max-iteration argument, this check will never satisfy)

   curIteration = curIteration + 1;

   if (curIteration == maxIterations)

      isDone = 1;

   end

end

