function [conf_mat, clust_class_distr, accuracy] = k_means_iris(num_clust)

% function [conf_mat, clust_class_distr, accuracy] = k_means_iris(num_clust)
% 
% Use k-means clustering on Iris dataset. Algorithm trained of 70% of data,
% test results reported on 30% of data.
% 
% Inputs : num_clusters      - number of clusters to be found.
% Outputs: conf_mat          - confusion matrix (rows represent true 
%                              labels,columns represent predicted labels).
%          clust_class_distr - class distribution over clusters, for every
%                              cluster found, reports what fraction of
%                              points belongs to classes 1, 2, and 3.
%          accuracy          - accuracy of k-means method.
%
% Author: Nemanja Djuric
% Date  : Sept. 25th, 2010

% check input
if (num_clust < 2 || num_clust > 20)
    disp('Error: Invalid input, num_clust is an integer between 2 and 20.');
    return;
end;

% preprocessing, nothing special, simply move label to the end
data = load('iris.txt');
data = [data(:, 2 : 5), data(:, 1)];

[len label] = size(data);
num_features = label - 1;

% use 70% for training, and the remaining 30% for testing
[testdata traindata] = divideset(data, 70);

% normalization, however, there is no need for that, so I commented it
% [meanv, stdv, traindata(:, 1 : num_features)] = normalize(traindata(:, 1 : num_features), [], []);
% [meanvv, stdvv, testdata(:, 1 : num_features)] = normalize(testdata(:, 1 : num_features), meanv, stdv);
    
% run k-means on the training data
[idx, cent] = kmeans(traindata(:, 1 : num_features), num_clust);

% find the class distribution over all clusters
clust_class_distr = [];
for clust = 1 : num_clust
    cluster_points = find(idx == clust);
    len = length(cluster_points);
    cluster_class1(clust) = sum(traindata(cluster_points, label) == 1) / len;
    cluster_class2(clust) = sum(traindata(cluster_points, label) == 2) / len;
    cluster_class3(clust) = sum(traindata(cluster_points, label) == 3) / len;
    [trash cluster_label(clust)] = max([cluster_class1(clust) cluster_class2(clust) cluster_class3(clust)]);
    
    clust_class_distr = [clust_class_distr; cluster_class1(clust) cluster_class2(clust) cluster_class3(clust)];
end;

% create nice strings to put as legend on the plot of Iris dataset
class_cluster = {'', '', ''};
class_names = {' - Iris Setosa', ' - Iris Versicolour', ' - Iris Virginica'};
for i = 1 : 3
    [trash ind] = find(cluster_label == i);
    for j = 1 : length(ind)
        if (j > 1)
            class_cluster{i} = [class_cluster{i}, ', '];
        end;        
        class_cluster{i} = [class_cluster{i}, num2str(ind(j))];
    end;
    
    % if class not correctly found, print Not/Found
    if (strcmp('', class_cluster{i}))
        class_cluster{i} = 'N/F';
    end;
    class_cluster{i} = [class_cluster{i}, class_names{i}];
end;

% test phase, assign each point to its nearest cluster, prediction will be
%   the majority class from that cluster found during training
% result are reported in confusion matrix, where rows are true labels, and
%   columns represent predicted values
conf_mat = zeros(3);
for i = 1 : length(testdata)
    new_point = testdata(i, 1 : num_features);
    true_label = testdata(i, label);
    
    for j = 1 : num_clust
        dist(j) = eucl_dist(new_point, cent(j, :));
    end;
    
    [trash predicted_cluster] = min(dist);
    conf_mat(true_label, cluster_label(predicted_cluster)) = conf_mat(true_label, cluster_label(predicted_cluster)) + 1;
end;

% outputs for all to see
clust_class_distr = 100 * clust_class_distr
conf_mat
accuracy = trace(conf_mat) / sum(sum(conf_mat))

try
    % first plot what the Iris data actually represents
    Image = imresize(imread('petal-sepal.jpg'), 2);
    figure(2), clf;
    imshow(Image);
    title('Features in Iris dataset (lengths and widths)');
    pause;
    close;
catch
end;

% plot centroids in all possible pairs of features
for i = 1 : num_features
    for j = i + 1 : num_features
        plot_iris([testdata; traindata], i, j, cent, cluster_label);
        legend(class_cluster{1}, class_cluster{2}, class_cluster{3}, 'Location', 'Best');
        pause;
    end;
end;
close;

% function that calculates Euclidean distance between points
function e_dist = eucl_dist(point1, point2)
e_dist = sqrt(sum((point1 - point2).^2));