function [Y,Its,evs] = mbms(X,sigma,KNN,L,flag,max_it,eta,alpha)
% mbms denoise data set according to its underlying manifold.
% 
% Syntax:
%    [Y,Its,evs] = mbms(X,sigma,KNN,L,flag,max_it,eta,alpha)
% 
% Input:
%         X = data set, each row is a data point.
%     sigma = the kernel width used in Gaussian kernel.
%       KNN = number of nearest neighbors used for estimating tangent space
%             at each data point.
%         L = the underlying dimension of manifold from which noisy data 
%             set is sampled, when L equals 0, no tangent space is computed
%             and the algorithm is running usual bms algorithm.
%      flag = indicates whether make Gaussian affinity matrix sparse,
%             0(default) no, otherwise yes.
%    max_it = maximum number of iterations, making sure the algorithm stops anyway.
%       eta = controls the speed of motion in explicit bms algorithm, that 
%             is, Xnew = eta*P*Xold+(1-eta)*Xold, default value is 1.0 .
%     alpha = normalize gaussian affinity matrix as done in
%             "diffusion maps", valid value: [0,1], default: 0.
% 
% Output:
%         Y = cell array storing the whole evolution of data set.
%       Its = the actual number of iteration used.
%       evs = cell array storing the whole evolution of eigenvalues of
%             local covariance matrix.
% 
% History:
%    Weiran Wang        May 22, 2009  recreated.

[N,D] = size(X);

if ~exist('L','var') || isempty(L)
    L = 0;           % Usual gbms algorithm.
end; 

if ~exist('flag','var') || isempty(flag)
    flag = 0;        % Default: don't use sparse affinity matrix.
end

if ~exist('max_it','var') || isempty(max_it)
    max_it = 20;     % Usually enough for denoising.
end

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

if ~exist('alpha','var') || isempty(alpha)
    alpha = 0;
end
% Default: do not use the normalization in "diffusion maps".

fprintf('\n');
fprintf('Number of data points: %d\n',N);
fprintf('Dimension of data points: %d\n',D);
fprintf('Kernel width: %g\n',sigma);
fprintf('Nearest neighbors, K: %d\n',KNN);
fprintf('Latent dimension: %d\n',L);
if flag==0
    fprintf('Using FULL graph!\n');
else
    fprintf('Using KNN graph!\n');
end
fprintf('Maximum number of iterations: %d\n',max_it);
fprintf('Explicit GBMS, eta: %f\n',eta);
fprintf('Normalization, alpha: %f\n',alpha);
fprintf('\n');

Its = 0; 
convcrit = 1;

while convcrit
    [W,Ngbs,md] = gaussknn(X,X,sigma,KNN,flag);
    fprintf('Average distance of each data points to its %d nearest neighbors: %g\n',KNN,mean(md));

    if alpha>0
	px = sum(W,2).^(-alpha);
        py = sum(W,1).^(-alpha);  
        W = diag(sparse(px))*W*diag(sparse(py));
    end
    
    deg = sum(W,2);   
    tX = diag(sparse(deg.^(-1)))*(W*X);
    dX = tX-X;
    
    if L>0    % Tangent space is of positive dimension.
      ev = zeros(N,min(D,KNN+1));
      for i=1:N
        Xngbs = X(Ngbs(i,:),:); 
        % Find neighborhood for X(i,:).
        [R,lambda] = pca(Xngbs,L);
        dX(i,:) = dX(i,:)-((dX(i,:))*R)*R';
        % Remove the motion in tangent space.
        ev(i,:) = lambda;
      end
    else
      ev = [];
    end
    
    X = X+eta*dX;
    Its = Its + 1;
    Y{Its} = X; evs{Its} = ev;
    
    convcrit = (Its < max_it);   % We need some convergence criterion.  
end

