function section2_hard(k_in)
%% Hard EM (k-means) example
% CS224N NLP Section 2 EM  2004.8.18
% Author: Paul Baumstarck
global P C k ix lw
if nargin < 1
    k = 2;
else
    k = k_in;
end
c1 = [1;2]+1.5;
c2 = [5;4]+1.5;
c3 = [6.5;1]+1.5;
pre = 3;
ms = 20;
ms1 = 16;
lw = 5;
n1 = pre*15;
n2 = pre*20;
n3 = pre*20;
n = n1+n2+n3;
ix = zeros(1,n);
s = 1;
P = [bsxfun(@plus,s*randn(2,n1),c1) ...
     bsxfun(@plus,s*randn(2,n2),c2) ...
     bsxfun(@plus,s*randn(2,n3),c3)];
figure(1), clf, hold on, grid on, axis equal;
plot(P(1,:),P(2,:),'.','MarkerSize',ms,'Color',0.5*ones(1,3));
title('Point cloud we wish to segment into 2 clusters');
z = axis;
z([1 3]) = 0;
axis(z);
pause;
C = reshape(P(:,ceil(n*rand(1,k))),[2 1 k]);
if k == 2
    kcolors = [1 0 0; 0 0 1];
elseif k == 3
    kcolors = eye(3);
else
    kcolors = lines(k);
end
for i=1:10
    figure(1), clf, hold on, grid on, axis equal;
    plot(P(1,:),P(2,:),'.','MarkerSize',ms,'Color',0.5*ones(1,3));
    for j=1:k
        plot(C(1,1,j),C(2,1,j),'^','MarkerSize',ms1,...
            'MarkerFaceColor',kcolors(j,:),'Color',kcolors(j,:));
    end
    title(sprintf('Cluster centroids, itertion %d',i-1));
    axis(z);
    pause;
    
    E(); % E step: compute expected labels ix
    
    figure(1), clf, hold on, grid on, axis equal;
    plot(P(1,:),P(2,:),'.','MarkerSize',ms,'Color',0.5*ones(1,3));
    for j=1:k
        plot(C(1,1,j),C(2,1,j),'^','MarkerSize',ms1,...
            'MarkerFaceColor',kcolors(j,:),'Color',kcolors(j,:));
    end
    plotline();
    axis(z);
    title('Minimum distance to centroid used to decide labels');
    pause;
	
    figure(1), clf, hold on, grid on, axis equal;
    plot(P(1,:),P(2,:),'.','MarkerSize',ms,'Color',0.5*ones(1,3));
    for j=1:k
        plot(C(1,1,j),C(2,1,j),'^','MarkerSize',ms1,...
            'MarkerFaceColor',kcolors(j,:),'Color',kcolors(j,:));
    end
    for j=1:k
        plot(P(1,ix==j),P(2,ix==j),'o','MarkerSize',ms,...
            'LineWidth',2,'Color',kcolors(j,:));
    end
    plotline();
    axis(z);
    title('Labeled points');
    pause;
    
    Cp = C;
    M();
    
    figure(1), clf, hold on, grid on, axis equal;
    plot(P(1,:),P(2,:),'.','MarkerSize',ms,'Color',0.5*ones(1,3));
    for j=1:k
        plot(Cp(1,1,j),Cp(2,1,j),'^','MarkerSize',ms1,...
            'MarkerFaceColor',kcolors(j,:),'Color',kcolors(j,:));
    end
    for j=1:k
        plot(P(1,ix==j),P(2,ix==j),'o','MarkerSize',ms,...
            'LineWidth',2,'Color',kcolors(j,:));
    end
%     plotline();
    axis(z);
%     plot(C(1,:),C(2,:),'rp'); % new centroids
    for i=1:n
        plot([P(1,i) C(1,1,ix(i))],...
             [P(2,i) C(2,1,ix(i))],':','Color',kcolors(ix(i),:));
    end
    title('New cluster centroids');
    pause;    
end


function E()
% Compute expected labels from current centroids
global P C k ix
d = sum(bsxfun(@minus,P,C).^2,1);
ix(:) = 0;
for i=1:k
    ix(min(d,[],3)==d(:,:,i)) = i;
end


function M()
% Compute best centroids from current labels
global P C k ix
for i=1:k
    C(:,:,i) = mean(P(:,ix==i),2);
end


function plotline()
global C k lw
if k == 2
    t = [-15 15];
    for i=1:k
        for j=i+1:k
            c = mean(C(:,:,[i j]),3);
            n = C(:,:,i)-C(:,:,j);
            n = n/norm(n);
            v = [-n(2);n(1)];
            plot(c(1)+t*v(1),c(2)+t*v(2),'-','LineWidth',lw,'Color',0.5*ones(1,3));
        end
    end
else
    voronoi(C(1,:,:),C(2,:,:));
end