function labels = spectral(X,k,beta)
%  labels = SPECTRAL(X,K,BETA)
%
%     runs spectral clustering on X and displays the result as a graph,
%     with points in X shown in blue or red (depending on the cluster to
%     which they were assigned). Also shown are the edges of the
%     neighborhood graph constructed by the algorithm.
%     Arguments: X - the data 
%                K - 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    = 3;  % nearest neighbors
end
if (nargin < 3)
  beta = 1;  % weight exponent
end

labels = zeros(size(X,1),1);

W = weights(X,k,beta);
V = eigvector(W);

I1 = find(V>0);
I0 = find(V<=0);
labels(I0)=0;
labels(I1)=1;


plot(X(I1,1),X(I1,2),'bo',X(I0,1),X(I0,2),'ro'); hold on;
plotgraph(X,W); hold off;

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

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;

% --------------------------------------
function [V] = eigvector(W)

w    = sum(W,2);
D    = diag(w); 
Dinv = diag(1./w); 

M = sqrt(Dinv)*W*sqrt(Dinv);

[V,E] = eig(M);
[et,I] = sort(diag(E)); 

V = V(:,I(end-1)); % eigenvector corresp. to second largest eigenvalue


function [W] = weights(X,k,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:(k+1));

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