% INPUT variables

% X: N-by-D data matrix;
% K: number of clusters
% sigma: tuning parameter

% OUTPUT variable

% U: N-by-K matrix, its rows are new coordinates in the embedded space

load twocircles.mat

K = 2;
sigma = 0.1;
% in general many sigmas need to be tested

N = size(X,1);

G = X*X';
norms = diag(G);
dists = repmat(norms,1,N) + repmat(norms',N,1) - 2*G;

isigma = 1/(2*sigma^2);

W = exp(-dists.*isigma);
W(1:(N+1):N^2) = 0;

degrees = sum(W,2);
degrees(degrees==0) = 1;
D_sqrt_mat = repmat(1./sqrt(degrees),1,N);
Z = D_sqrt_mat.*W.*D_sqrt_mat';

Z = real(Z);
Z = (Z+Z')/2;
opt.isreal = true;
opt.issym = true;
opt.disp = 0;

[Ui,eigvals] = eigs(Z,2*K,'LM',opt);
eigvals = diag(eigvals);

[sort_eigvals, ind_sort] = sort(abs(eigvals),'descend');

U = Ui(:,ind_sort(1:K));
% normalize U along the rows to have unit length
U = U./repmat(sqrt(sum(U.^2,2)),1,K);