% K-means Clustering Algorithm 
% Two dimensional clustering demonstration
% Karl Dyer - October 2010

%% HOUSE CLEAN
% *************
clear all; close all; clc;

%% DEFINE CONSTANTS
% ******************
prompt = {'Enter number of points (0 for daisies):',...
    'Enter number of clusters:',...
    'Enter axis bounds:'};
name = 'K - Cluster Algorithm Demo'; numlines=1;
defaultanswer = {'500', '5', '10'};
design_parameters = inputdlg(prompt,name,numlines,defaultanswer, 'on');

N = str2double(design_parameters{1});            % Number of data points
B = str2double(design_parameters{3});            % Axis bounds

if (N == 0)
    K = 5;                                       % Number of data points
    pause_time = 0.75;                           % Delay between iterations
else
    K = str2double(design_parameters{2});        % Number of clusters
    pause_time = 0.3;                            % Delay between iterations
end

% Create names for legend
for i = 1:K
    names{i} = (['Cluster ', num2str(i)]);
end

% Create color schemes for plotting
color = {'r.', 'b.', 'g.', 'k.',...
         'kx', 'rx', 'bx', 'gx',...
         'go', 'ko', 'ro', 'bo',...
         'b+', 'g+', 'k+', 'r+',};

%% ERROR CHECK
% *************
if (N ~= 0 && N < K)
    display('Error - Number of points must exceed number of classes');
    break
end

if (K > 16)
    display('WARNING - Not enough unique marker combinations defined.')
    display('Add more values to color matrix at Line 36 or reduce K')
    break
end

%% CREATE DATA
% *************
if (N == 0)                 % If 0 entered create 'daisy dataset'
    x(:,1:10) = B/8*ones(2,10) + (rand(2, 10) - 0.5);
    x(:,11:20) = B/2*ones(2,10) + (rand(2, 10) - 0.5);
    x(:,21:30) = B*7/8*ones(2,10) + (rand(2, 10) - 0.5);
    x(:,31:40) = B/8*ones(2,10) + (rand(2, 10) - 0.5); 
        x(2,31:40) = x(2,31:40) + B*3/4;
    x(:,41:50) = B/8*ones(2,10) + (rand(2, 10) - 0.5); 
        x(1,41:50) = x(1,41:50) + B*3/4;
else                        % Otherwise create specified random points
    x = B*rand(2, N);       % Coordinates of data points Row 1: X, Row 2: Y
end

%% INITIALIZATION
% ****************
t = 1;                  % Time index
dummy = 1;              % Dummy variable to keep loop running

I = randperm(length(x));
c = x(:, I(1:K));       % Pick random point as the start value for centroid

%% CALCULATE DISTANCES AND GROUPS
% ********************************
d = dist(c', x);               % Calculate distance from each centroid to each point 
d(d==0) = B^2;                 % Remove 0's from distance matrix so min calculation is not skewed
[min_val groups] = min(d,[],1);% Calculate minimum distance and find closest centroid
temp = groups;

while (dummy == 1)
    %% UPDATE CENTROID
    % *****************
    for i = 1:K
        c(:,i) = sum(x(:,groups==i),2)/sum(groups==i,2);
    end
    
    %% CALCULATE DISTANCES AND GROUPS
    % ********************************
    d = dist(c', x);               % Calculate distance from each centroid to each point 
    d(d==0) = B^2;                 % Remove 0's from distance matrix so min calculation is not skewed
    [min_val groups] = min(d,[],1);% Calculate minimum distance and find closest centroid
    
    %% PLOT CURRENT DATA
    % *******************
    %Plot data points by cluster
    figure(1)
    hold off;
    for i = 1:K
        plot(x(1,groups==i), x(2,groups==i), color{i})
        hold on
    end

    %Plot centroids
    plot(c(1,:),c(2,:),'kx','MarkerSize',10,'LineWidth',2)
    plot(c(1,:),c(2,:),'ko','MarkerSize',10,'LineWidth',2)

    % Title from with model parameters
    title({'K - Means Clustered Data';...
        ['N = ' num2str(length(x)), ', K = ', num2str(K)];...
        ['Iteration: ', num2str(t)];})
    legend(names,'Location','EastOutside')
    axis([0 B 0 B]);
    
    %Capture frames for movie
    M(t) = getframe;            % Record plot for movie

    %% CHECKS AND UPDATES
    % ********************
    if (temp == groups)
        break                   % Exit loop if optimal centroid found
    else
        temp = groups;          % Update group comparison
        t = t + 1;              % Increment time
        pause(pause_time);      % Delay plot update
    end
end

%% SUMMARIZE RESULTS
% *******************
clc;
display(['After ', num2str(t), ' iterations:'])
for i = 1:K
    display(['  Cluster ', num2str(i), ' has ', num2str(sum(groups==i)),...
        ' points with a centroid at (', num2str(c(1,i)), ', ', ...
        num2str(c(2,i)), ')']);
end
display(' ');
display('To see animated sequence again type: movie(M)')
    