function [vEigenVectors,vEigenValues,vLaplacian] = FastLaplacianEigs( cOpts, cDegrees, cWeights, cNumberOfEigenValues, cPrecision )

% function [vEigenVectors,vEigenValues,vLaplacian] = FastLaplacianEigs( cOpts, cDegrees, cWeights, cNumberOfEigenValues, cPrecision )
%
% IN:
%   cOpts              : structure containing
%       NormalizationType :  'graph':   D^{-1/2}(D-W)D^{-1/2}
%                            'graphmarkov' : D^{-1/2}(D-W)D^{-1/2} but the returned eigenvectors are rescaled by D^{-1/2}, so they are the same as 'ave'
%                            'beltrami' : D^{-1}(D-W)D^{-1}
%                            'ave':     D^{-1}(D-W)
%                            'noave':   D-W
%                            'none':    W
%       WhichEigenVals : 'small' or 'large'
%   cDegrees    : 1xN vector: d_{i}=order of the vertex i
%   cWeights    : NxN matrix : w_{ij}=weight between vertex i and vertex j, a symmetric matrix.
%   cPrecision  : Desired precision
%
% OUT:
%   vEigenVectors, vEigenValues: eigenvectors (row-wise) and eigenvalues of the Laplacian.
%   vLaplacian : the Laplacian sparse matrix
%
% USES:
%
% Original code by MM. Modified by AL, May 2004

if ( strcmpi(cOpts.NormalizationType,'graph')==1 ) | ( strcmpi(cOpts.NormalizationType,'graphmarkov')==1 )                  % D^{-1/2}(D-W)D^{-1/2}
    %DD = speye(length(cDegrees));
    % Normalized Laplacian
    %for i=1:length(cDegrees)
    %    if (cDegrees(i) ~= 0)
    %        DD(i,i) = 1/sqrt(cDegrees(i));
    %    else 
    %         DD(i,i) = 0;
    %    end
    % end
    N=length(cDegrees); temp=zeros(N,1); idx=find(cDegrees>eps);
    temp(idx)=1 ./sqrt(cDegrees);
    DD=spdiags(temp,0,N,N);
    
    % Normalized laplacian
    vLaplacian = DD*(spdiags(cDegrees,0,speye(size(cWeights,1)))-cWeights)*DD;
    % Force simmetry to precision
    vLaplacian = (vLaplacian+vLaplacian')/2;
else
    if strcmpi(cOpts.NormalizationType,'ave')==1                 % D^{-1}(D-W)
        %DD = speye(length(cDegrees));
        % Normalized Laplacian
        %for i=1:length(cDegrees)
        %    if (cDegrees(i) ~= 0)
        %        DD(i,i) = 1/cDegrees(i);
        %    else 
        %         DD(i,i) = 0;
        %    end
        % end
        N=length(cDegrees); temp=zeros(N,1); idx=find(cDegrees>eps);
        temp(idx)=1 ./cDegrees;
        DD=spdiags(temp,0,N,N);
    
        % beltrami laplacian
        vLaplacian = speye(size(cWeights,1))-DD*cWeights;%DD*(spdiags(cDegrees,0,speye(size(cWeights,1)))-cWeights);
    else
        if strcmpi(cOpts.NormalizationType,'beltrami')==1                 % D^{-1}(D-W)D^{-1}
            %DD = speye(length(cDegrees));
            % Normalized Laplacian
            %for i=1:length(cDegrees)
            %    if (cDegrees(i) ~= 0)
            %        DD(i,i) = 1/cDegrees(i);
            %    else 
            %        DD(i,i) = 0;
            %    end
            %end
            N=length(cDegrees); temp=zeros(N,1); idx=find(cDegrees>eps);
            temp(idx)=1 ./cDegrees;
            DD=spdiags(temp,0,N,N);
        
            % Normalized laplacian
            %vLaplacian = DD*(spdiags(cDegrees,0,speye(size(cWeights,1)))-cWeights)*DD;
            lKtilde=DD*cWeights*DD;
            lDegrees = sum(lKtilde,2);
            % Normalized Beltrami Laplacian
            %for i=1:length(lDegrees)
            %    if (lDegrees(i) ~= 0)
            %        DD(i,i) = 1/lDegrees(i);
            %    else 
            %         DD(i,i) = 0;
            %     end
            % end
            N=length(lDegrees); temp=zeros(N,1); idx=find(lDegrees>eps);
            temp(idx)=1 ./sqrt(lDegrees);  %???
            DD=spdiags(temp,0,N,N);
           
            vLaplacian = DD*(spdiags(lDegrees,0,speye(size(cWeights,1)))-lKtilde)*DD;
            % Force simmetry to precision
            vLaplacian = (vLaplacian+vLaplacian')/2;
            clear lKtilde lDegrees;
        else
            if strcmpi(cOpts.NormalizationType,'noave')==1           % D-W
                % Normalized laplacian
                vLaplacian = spdiags(cDegrees,0,speye(size(cWeights,1)))-cWeights;
            else                                                % W
                vLaplacian = cWeights;
            end;
        end;
    end;
end;
    
% Then we want to find eigenvalues and eigenvectors
opts.tol = cPrecision;
opts.issym=1;                           % No matter what the Laplacian, it's always symmetric since we assume W symmetric 
opts.disp = 0; 
opts.isreal = 1;
if strcmpi(cOpts.WhichEigenVals,'small')==1
    [vEigenVectors, vEigenValues]=eigs(speye(size(cWeights,1))-vLaplacian,cNumberOfEigenValues,'la',opts);
    vEigenValues = 1-diag(vEigenValues);
else
    [vEigenVectors, vEigenValues]=eigs(vLaplacian,cNumberOfEigenValues,'la',opts);
    vEigenValues = diag(vEigenValues);
end;
vEigenVectors = vEigenVectors';

% Sort the eigenvalues from small to big, do the same on the eigenvectors
[lTemp lIdx] = sort(vEigenValues);
vEigenVectors = vEigenVectors(lIdx,:);

if (strcmpi(cOpts.NormalizationType,'graphmarkov')==1) | (strcmpi(cOpts.NormalizationType,'beltrami')==1)
    for lk = 1:size(vEigenVectors,1)
        vEigenVectors(lk,:) = vEigenVectors(lk,:)*DD;
        vEigenVectors(lk,:) = vEigenVectors(lk,:)/norm(vEigenVectors(lk,:));
    end;
end;

return;