function sheet08

% Generate data
centers = [0, 0; 7, 3; -2, 4; 0, 10; -5, -5];
% (2): You centers here
% centers = [...]
X = generate_data(centers, 50);

MAXK = 8; % maximal number of clusters
ITERS = 10; % how many restarts.

% compute (in-)stability for k-means and a random clustering
DKMEANS = evaluate_clusterer(@k_means_clustering, X, MAXK, ITERS);
DRAND = evaluate_clusterer(@random_clustering, X, MAXK, ITERS);

% plot the two stability values
figure(1)
errorbar(2:8, mean(DKMEANS), std(DKMEANS));
hold on
errorbar(2:8, mean(DRAND), std(DRAND), 'r');
hold off
legend('k-means stability', 'random clusterer')

% normalize the stability by the random clusterer
D = normalize(DKMEANS, DRAND);

figure(2)
errorbar(2:8, mean(D), std(D));
title('normalized k-means')

% show clustering for most stable number of clusters
[dummy, OPTK] = min(mean(D)); OPTK = OPTK + 1;

Y = k_means_clustering(OPTK, X);
figure(3);
plot_clustering(X, Y);

% generate data from some centers
function X = generate_data(C, N)
X = [];
for I = 1:size(C, 1)
  X = [ X ; randn(N, 2) + repmat(C(I, :), N, 1) ];
end

% A variant of the k-means clustering algorithm which just iterates 100
% times.
function Y = k_means_clustering(K, X)
[N, D] = size(X);

% randomly select K points as centers
P = randperm(N);
MEANS = X(P(1:K), :);

% set up variables
for I = 1:100
  % compute nearest neighbor assignments
  dist = pwdist(MEANS, X);
  [dummy, Y] = min(dist);
  
  % compute new means
  for J = 1:K
    NJ = sum(Y == J);
    if NJ > 0
      MEANS(J, :) = sum(X(find(Y == J), :))/NJ;
    end
  end
end

% Plot the clustering and the centers.
function plot_clustering(X, Y)
gscatter(X(:, 1), X(:, 2), Y);

% Compute all pairwise distances quickly.
function D = pwdist(X, Y)
D = size(X, 2);
N = size(X, 1);
M = size(Y, 1);
  
XX = sum(X.*X, 2);
YY = sum(Y.*Y, 2);
D = repmat(XX, 1, M) + repmat(YY', N, 1) - 2*X*Y';

% evaluate a clustering algorithm for clusters 2:MAXK
% for ITERS restarts.
function D = evaluate_clusterer(method, X, MAXK, ITERS)
D = zeros(ITERS, MAXK-1);
for K = 2:MAXK
  fprintf('K = %d\n', K);
  for I = 1:ITERS
    Y1 = feval(method, K, X);
    Y2 = feval(method, K, X);
    D(I, K-1) = cocluster_distance(Y1, Y2);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fill in your solutions below
%

% 1 (a). The random clusterer. It should return random labels from 1..K
% for each data point in X (stored in the rows.
function Y = random_clustering(K, X)
% ...

% 1 (b). Compute the co-cluster distance for two labelings Y1 and Y2.
function D = cocluster_distance(Y1, Y2)
% ...

% 1 (c). Normalize DKMEANS by dividing the result for a fixed
% K by the mean value of the corrsponding values in DRAND.
function D = normalize(DKMEANS, DRAND)
% ...
