function [W] = weights(X,r,beta)

% n data points in d dimensions
[n,d] = size(X);


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

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