function [Y,Its,evs] = ltp(X,KNN,L,flag,max_it,eta)
% 
% Syntax:
%    [Y,Its,evs] = ltp(X,KNN,L,flag,max_it,eta)
% 

[N,D] = size(X);

if ~exist('L','var') || isempty(L)
    L = 0; 
end; 

if ~exist('flag','var') || isempty(flag)
    flag = 0;
end; 

if ~exist('max_it','var') || isempty(max_it)
    max_it = 20;
end

if ~exist('eta','var') || isempty(eta)
    eta = 1.0;
end

Its = 0; 
convcrit = 1;

while convcrit

    [Ngbs,md] = ltpknn(X,X,KNN,flag);
    
    fprintf('Average distance of each data points to its %d nearest neighbors: %g\n',KNN,mean(md));
    
    ev = zeros(N,min(D,KNN+1));
 
    tX = X;
    if L>0
      for i=1:N
        Xngbs = X(Ngbs(i,:),:); 
        % Find neighborhood for X(i,:).
        [R,lambda,xc] = pca(Xngbs,L);
        tX(i,:) = ((X(i,:)-xc)*R)*R' + xc;
        % Remove the motion in tangent space.
        ev(i,:) = lambda;
      end
    else
      for i=1:N
        Xngbs = X(Ngbs(i,:),:); 
        tX(i,:) = mean(Xngbs,1);
      end
      ev = [];
    end
    
    X = X+(tX-X)*eta;
    Its = Its + 1;
    Y{Its} = X;
    evs{Its} = ev;
 
    convcrit = (Its < max_it);   % We need some convergence criterion.  
end
