function [E,V] = lapbin(DATA,lodim,nn,step,stepsbtwnreport)
% LAPBIN Laplacian Eigenmaps Algorithm with Binary weights.
%
% [E,V] = lapbin(DATA,lodim,[nn,step,stepsbtwnreport])
%
% DATA is a N x hidim matrix representing N points in R^hidim
% lodim is a positive integer; you want to map the points in DATA
% to R^lodim.
% nn is the number of nearest neighbors used to create the graph
% (10 by default)
% step is a value that can be modified to optimize this code's
% running on your machine. L2_distance is called on a step x N matrix.
%
% E is a N x lodim matrix represnting the same N points in R^lodim.
% V is a vector with the lodim smallest positive eigenvectors.
%
% E(:,i) and V(i) are the i-th smallest positive eigenvectors/values.
%
% Original code by Misha Belkin (misha@cs.uchicago.edu)
% Modified by Dinoj Surendran (dinoj@cs.uchicago.edu), April/May 2004
if nargin < 3, nn = 10;end
if nargin < 4, step=500;end
if nargin < 5, stepsbtwnreport=5;end
n=size(DATA,1);
hidim = size(DATA,2);
W=sparse(n,n);
sprintf ('Finding a mapping of %d points in R^%d to R^%d\n using the Laplacian Algorithm with binary weights and %d nearest neighbors\n\n', n, hidim, lodim, nn)
% calculate the adjacency matrix for DATA
count=0;
for i1=1:step:n
count=count+1;
i2 = min(i1+step-1,n);
XX= DATA(i1:i2,:);
dt = L2_distance(XX',DATA');
[Z,I] = sort ( dt,2);
for i=i1:i2
W(i,I(i-i1+1,2:nn+1)) = ones(1,nn);
W(I(i-i1+1,2:nn+1),i) = ones(nn,1);
end
if (mod(count,stepsbtwnreport) ==0)
disp(sprintf('%d points processed.', i1+step-1));
end;
end;
% D = sum(W,2);
L = spdiags(sum(W,2),0,sparse(n))-W;
opts.tol = 1e-9;
opts.maxit = 200;
opts.issym=1;
opts.disp = 2;
[EE,VV] = eigs(L,1+lodim,'sm',opts);
% the smallest eigenvector is in the (1+lodim)-th column (it's not positive)
% we want the first column to correspond to the second smallest evalue
% we want the second column to correspond to the third smallest evalue, etc
%
V = zeros (1,lodim);
for i = 1 : lodim
E(:,i) = EE(:,1+lodim-i);
V(i) = VV(1+lodim-i,1+lodim-i);
end