function y = spectral(X,k,r,beta)
%     Arguments: X - the data
%                k - number of resulting clusters
%                r - number of neighbors (default: 4)
%                BETA - the weight falloff parameter (default: 1)
%
%     Returns in labels(i) the cluster index associated with X(i,:)

if (nargin <2)
  k    = 2;  % number of clusters
end
if (nargin < 3)
  r    = 4;  % nearest neighbors
end
if (nargin < 4)
  beta = 1;  % weight exponent
end

W = weights(X,r,beta); % Unnormalized transfer weights
Dinv    = diag(1./sum(W,2));  % Required normalization
M = sqrt(Dinv)*W*sqrt(Dinv); % Stochastic (roms sum to one) matrix of transition probs

[V,E] = eig(M); % eigenvectors and eigenvalues of M
[sortedE,I] = sort(diag(E)); % sort eigenvalues
V = V(:,I(end-1)); % eigenvector corresp. to second largest eigenvalue

y = (3+sign(V))./2;