function [U] = kkMeans2(K, nc, epsilon)
% Kernel k-means. 
% Returns the (nc x size(X,1)) matrix of cluster memberships.
%
% Input:
%   K -- kernel matrix
%   nc -- number of clusters
%   epsilon -- threshold for termination
%   varargin -- kernel parameters according to 'kernel' (see kernel.m)
% Output: 
%   U -- membership matrix (nc x size(X,1))
%
% (c) 2008-2016 Zalan Bodo

N = size(K, 1);

iC = randperm(N); % Indices of initial centers 
iC = iC(1:nc);

t1 = ones(nc, N)*diag(diag(K));

D = t1 + diag(diag(K(iC,iC)))*ones(nc, N) - 2*K(iC,:);
[val, ind] = min(D);
U = zeros(N, nc);
U(sub2ind([N, nc], 1:N, ind)) = 1;
U = U'; % Initial cluster assignment
%sum(U')

t = 1;
while (1)
    % DISTANCES, F:
    s = diag(sum(U'));
    mU = (s^-1)*U;
    dd = diag(mU*K*mU');
    t2 = repmat(dd, [1,N]);
    t3 = -2*mU*K;
    D = t1 + t2 + t3;
    F(t) = sum(diag(U*D'));
    % Convergence:
    if (t > 1) && (abs(F(t-1)-F(t)) < epsilon)
        break;
    end
    % U:
    U = zeros(nc, N);
    for j = 1:N
        [value, index] = min(D(:,j));
        U(index,j) = 1;
    end
%     disp(t);
%     fflush(stdout);
    t = t+1;
end
