%  [preproc, processed] = prockmeans(data, preproc)

%

%    Perform processing of data using k-means algorithm for clustering.

%

%    Input:   data      =  vector of coordinates

%             preproc   =  previously found k-means preprocessing (optional); if supplied,

%                          function will use this preprocessing to process the data rather

%                          finding new proprocessing

%             inputdim  =  number of desired dimensions

%

%    Output:  preproc   =  calculated (or supplied) preprocessing, found using k-means algorithm

%             processed =  processed data, according to the revealed preprocessing



function [preproc, processed] = prockmeans(data, preproc, inputdim)



% check number of arguments

if ((nargin < 1) | (nargin > 3))

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

end



% preliminary -- find data dimension and count, conclude the size of k (num of clusters)

[dataNum, dataDim] = size(data);



if (nargin < 3 | isempty(inputdim))

   inputdim = floor(sqrt(dataNum));

end



if ((nargin == 1) | isempty(preproc))

   % preprocessing not supplied, therefore we should find it...

   

   % estimate k (the expected number of clusters)

   k = inputdim;



   % find k centers, using two consecutive runs of k-means algorithm:

   if (k > 2)

      % one run to find k*log(k) centers...

      mid_centers = kmeans(ceil(k * log(k)), data, 0, 1e-5, 30);

   else

      mid_centers = data;

   end

   % second run to find k centers over those k*log(k) found...

   centers = kmeans(k, mid_centers, 0, 1e-5, 30);

   

   % divide the vectors into clusters, by minimal MSE between each vector 

   % and the found / supplied centers:

   % first, concatenate cluster centers and data vectors, for better performance

   tmpCenters = ncat(centers, dataNum);

   tmpData    = ndup(data, k);

   % second, calculate MSE of distance for each center and data vector

   distVec = mserr(tmpData, tmpCenters, 2);

   % third, put distances vector into matrix

   distMtx = reshape(distVec, k, dataNum);

   % fourth, find location of minimum value within distances matrix

   [distMinVal, distMinIndx] = min(distMtx);

   

   % find the variance of each dimension, per each cluster separately

   clustVar = zeros(k, dataDim);

   for i = 1 : k

      clustVar(i, :) = var(data((distMinIndx == i), :));

   end

   % add an "epsilon" element to the var, to avoid normalizing (division) by 

   % zero variance later on

   clustVar = clustVar + 0.1;



   preproc(:, :, 1) = centers;

   preproc(:, :, 2) = clustVar;

else

   % retrieve preprocessing parameters

   centers  = preproc(:, :, 1);

   clustVar = preproc(:, :, 2);

   k = size(centers, 1);

   tmpCenters = ncat(centers, dataNum);

   tmpData    = ndup(data, k);

end



% now, use the preprocessing to process the data and produce new coordinates

processedVec = sum(((tmpData - tmpCenters) .^ 2) ./ ncat(clustVar, dataNum), 2);

processed = reshape(processedVec, k, dataNum)';

