function y = spectral(X,k,r,beta)
%  labels = SPECTRAL(X,K,BETA)
%
%     runs spectral clustering on X and displays the resulting clustering, 
%     overlayed on the neighbourhood graph used.
%     Arguments: X - the data
%                k - number of resulting clusters
%                r - number of neighbors (default: 3)
%                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    = 3;  % 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

% evaluate the k-dimensional representations of the points and
% run the k-means algorithm to get y

V = ...
y = ...

% Also plot the result

plotgraph(X,W);
hold on; plotclust(X,y); hold off;

% ----- function to calculate the weights
function [W] = weights(X,r,beta)

[n,d] = size(X);

W = zeros(n,n);
for i=1:n,
  d = sqrt( sum((X-repmat(X(i,:),n,1)).^2,2) );
  [dt,ind] = sort(d);
  ind = ind(1:(r+1));

  W(i,ind) = exp(-beta*d(ind)');
  W(ind,i) = exp(-beta*d(ind));
end;

% --------------------------------------
function [] = plotgraph(X,W)

clf; hold on;
for i=1:size(W,1),
  ind = find(W(i,:));
  ind = ind(find(ind~=i)); 

  for j=1:length(ind),
    plot([X(i,1);X(ind(j),1)],[X(i,2);X(ind(j),2)],'k');
  end;
end;
hold off;
