function [centers, radii] = mykmeans(k, data, err, minimal)
% Learning and neural computation project by - Igor Bogodlov and Ron Sidi
% Function: mykmeans(k, data, err, minimal)
%
% K-Means clustering algorithm
%
% Output:
%    centers - the centers of the clusters
%    radii   - the radiuses of the clusters (this is a vector of radiuses equal in all directions)
%   empty clusters are not returned (the number of clusters could be less then requested)
%
% Input:
%    k       - requested number of clusters
%    data    - dataset
%    err     - maximal centers change for termination criterion 
%    minimal - minimal number of points in cluster
%

% Getting data size
dsize = size(data, 1);
ddim = size(data, 2);

% if minimal is not given, set it
if (nargin < 4)
  % if a cluster has less than 'minimal' vectors in it we try to
  % find a new random cluster center. This threshold is set to 
  % 30% of the average cluster size
  minimal = ceil(.3*dsize/k + .5);
end

% if maximal error is not given then set it to eps
if (nargin < 3)
  err = eps;
end

% We choose initial cluster centers set to be as much 'scattered'
% as possible. The set is built iteratively, and each time a vector
% with largest minimum distance to the current set is added 

alldist = dist(data,data');
centeridx = [unidrnd(dsize)];
for i = 1:k
	[x,y] = min(alldist(:,centeridx)',[],1);
	[r,t] = max(x);
	centeridx = [centeridx, t];
end

newcenters = data(centeridx',:); 

% Initialization of the loop
iterno = 1;
maxdiff = inf;
prevdiff = inf;

% Only clusters with points are considered active and will be returned
active = repmat(1==1,1,k);

% This is main k-means clustering loop. In each iteration points are
% assigned to nearest clusters, and then cluster centers are recalculated
% The process stops when the difference between centers before and after
% the iteration gets smaller than the requested error, or a certain amount
% of iterations were performed.

while ((maxdiff > err) & (iterno < 20))

   prevdiff = maxdiff;   
   centers = newcenters;
   % disp(['Iteration ', num2str(iterno), ', max change is ', num2str(maxdiff)]);
  
   % Find the distance of all of the points from all the cluster centers
   distances = dist(data,centers');
   
   % For each point, find the closest cluster center
   [mins, places] = min(distances');

   % Find cluster with most points
   [p, jumps] = unique(sort(places));
   [b, t] = max(diff([0,jumps]));
   fatcluster = p(t);

   % Assign points to cluster centers
   for i = 1:k
      
      % Find all the points in the current cluster
      b = places==i;
      
      % Count them
      s = sum(b);
      
      if (s < minimal)
         % There are less than minimum points in cluster.
         % Select another cluster center.
         
         % Choose a point from the "fattest" cluster, which is
         % a maximal distance away from the cluster center. If
         % the cluster has only one point, we choose a point, 
         % which is a maximal distance from all
         % the other cluster centers, among randomly selected subset
         % of data points. Introduced randomness help escape local
         % minima.
         % Remark: Since this shouldn't happen much the "fattest" cluster
         %         is calculated outside of the loop in order to make it
         %         run faster.
         
         % Find the distance of all points from all the points that were
         % ever taken as a center
         [x,y] = min(alldist(:,centeridx)',[],1);
         range = 1:dsize;
         
         % Take all the points that are in the biggest cluster
         excerpt = range(places==fatcluster);
         
         % If the biggest cluster is too small, randomly take half of the input set
         if size(excerpt,2) <= 2
            excerpt = unidrnd(dsize,1,ceil(dsize/2));
         end
         
         % Find a point from the "fattest" cluster which is the frathest from all the 
         % points that where ever selected as centers 
         [r,t] = max(x(excerpt));
         t = excerpt(t);
         
         % add the selected points to the list of points that were selected as centers
         % in order not to select a point twice
         centeridx = [centeridx, t];
         
         % Update the new centers
         newcenters(i,:)=data(t,:);
         
      else
         % If there are enough points in the cluster,
         % set the new center to be mean of all its points
         newcenters(i,:)=mean(data(b,:),1);
      end
   end
   iterno = iterno + 1;
%  maxdiff = max(sqrt(sum((newcenters'-centers').^2)));
   maxdiff = max(mean((newcenters'-centers').^2));
end

% Calculate the radiuses of the selected clusters
radii = [];

% For each cluster, if it is active calculate the radius, otherwise mark it as not active
for i = 1:k
   % Calc the size of the cluster
   s = sum(places==i);
   
   % Clusters with 0 points are not active
   if s > 0
      if s == 1
         % If there is only one point in the cluster then the radius is the distance to 
         % the closest point in all the other clusters
         
         % rad = abs(data(places==i,:)-centers(i,:));
         rad = distances(places==i,i);
         rad = rad(1);
         if rad <= eps
            % others = 1:dsize;
            % others = others(:,places~=i);
            % [mn, place] = min(distances(others,i));
            % rad = abs(data(others(place),:)-centers(i,:));
            rad = min(distances(places~=i,i));
	 end

	 % Uncomment out if different radii in all directions wanted
	 % rad = rad(1,ones(1,ddim)) ./ sqrt(ddim);
      else
         % c = centers(i,:);
         % rad = sqrt(mean((data(places==i,:)-c(ones(1,s),:)).^2));
         rad = mean(distances(places==i,i));

	 % Uncomment out if different radii in all directions wanted
	 % rad = std(data(places==i,:));
      end
      
      radii = [radii; rad];
   else
      % Mark the cluster as not active
      active(i) = logical(0);
    end;
end

% zeroplaces = radii==0;
% radii(zeroplaces) = 0.01;

% Remove the centers of all non-active clusters
centers = centers(active,:);

% plot(data(:,1),data(:,2),'r.');
% hold on;
% plot(centers(:,1),centers(:,2),'b*');
% hold off;
% refresh;
