function Diffuser_GlobalEig( cFcn, cPoints, cExtraParam )

%
% function Diffuser_GlobalEig( cFcn, cPoints, cExtraParam )
%
% IN:
%   gDiffuseObj : global variable containing the state of the diffusion object. Created by calling this with cFcn='create'
%                 It's a structure containing at least the following fields:
%                 Heat : current distribution on cPoints                                                             % NumClasses-by-M matrix
%                 A : current adjacency matrix for cPoints (non-normalized)
%                 D : normalizing measure on the points, i.e. D(x,x)=\int_{y\in cPoints}W(x,y)
%                 P           : normalized transition matrix
%                 RangeRadius : radius for range searches
%                 Atria       : preparation for fast range search
%                 Heat        : current heat (A^t x);                                                               % NumClasses-by-M matrix
%                 Delta       : scale parameter in the distance kernel
%                 BoundaryIdxs: indices of the points at the boundary of diffusion (points directly to data matrix) 
%                 AFlags      : M vector where comp i is 1 if node i is included in A; 0 otherwise                 
%                 M           : number of data points                                                              
%                 NumClasses  : number of classes
%                 GetWeightsFcn : function for weight computation, of the form GetWeights_Template
%                 
%   cFcn : string specifying the action to be taken. Supports:
%          'create'     : creates and initializes gDiffuseObj. Expects cExtraParam.P to be an Numclass-by-M matrix representing the initial distribution to be diffused.
%          'diffuse'    : diffuses one step
%          'plug'       : plugs the points specified by cExtraParam, which is a vector of indices into cPoints (and hence into gDiffuseObj.A), so that no extra diffusion occurs.
%   cPoints  : M points in R^N in a M by N matrix
%   cExtraParam : extra parameters. 
%                 When cFcn=='create' is should contain at least the following fields:
%                 Heat : initial heat distribution on the points
%                 Delta : scale of the kernel
%                 kNN: maximum number of nearest neighbors to consider
%                 Precision : precision for the calculations
%                 MaxEigenValues: Maximum number of eigenvalues to keep
%                 EigEnergyThres: fraction of energy of the first MaxEigenValues to preserve, it's a number in (0,1]
%                 EigFlag: 1 if direct computation of matrix power, 0 if step-wise global diffusion
%                 t: time (needed if EigFlag=1)
%
% OUT:
%   gDiffuseObj : updated global variable gDiffuseObj:
%                 If cFcn = 'create', gDiffuseObj gets initialized, possibly based on cExtraParam
%                 If cFcn = 'diffuse', gDiffuseObj.A,.D,.P are updated after a step in the diffusion of gDiffuseObj.P is done.
%   gEigObj:      If cFcn = 'create' and cExtraParam.EigFlag=1, gEigObj is intialized
%                 EigenVecs, EigenVals: eigenvectors and eigenvalues of symmetrized Laplacian of graph (columnwise)
%                 EigenVecsL, EigenVecsR: eigenvectors normalized with multiplicative factors 1/sqrt(D) and sqrt(D), resp. 

% USES:
%   nn_prepare,range_search, FastLaplacianDetEigs(+FastLaplacianEigs)
%
% SOURCE CONTROL:
%   Ann Lee, 5/27/04
%   Modified,   AL 3/19/04 [Calculate number of eigencomps; set small negative heat to zero]
%               AL 3/24/04 [Normalize P_t]
%               AL 5/27/04 [Added: gDiffuseObj.Metric and gDiffuseObj.MetricFlags for points with fixed warped metric]

global gDiffuseObj;         % Thanks to Matlab which hasn't got pointers...
global gEigObj;   % variables specific to global diffusion using eigenvectors

if strcmpi(cFcn,'create')
    % Constructor
    
    %lM = size(cPoints,1);                                                 
    gDiffuseObj = [];
    gDiffuseObj.Heat        = cExtraParam.Heat;                                 % NumClasses-by-M matrix
    gDiffuseObj.M=size(cPoints,1); gDiffuseObj.NumClasses=size(gDiffuseObj.Heat,1);  
    %%gDiffuseObj.A = speye(gDiffuseObj.M);  gDiffuseObj.P=gDiffuseObj.A;  gDiffuseObj.D=sum(gDiffuseObj.A,2);  
    gDiffuseObj.Atria       = nn_prepare(cPoints);
    %gDiffuseObj.AIdxs       = find(gDiffuseObj.Heat~=0); 
    tmp=max(abs(gDiffuseObj.Heat),[],1); %gDiffuseObj.AIdxs=find(tmp>eps);   
    gDiffuseObj.BoundaryIdxs= find(tmp>eps);
    gDiffuseObj.Delta       = cExtraParam.Delta; 
    gDiffuseObj.Precision   = cExtraParam.Precision;
    gDiffuseObj.RangeRadius = gDiffuseObj.Delta*(log(1/gDiffuseObj.Precision))^(0.5);
    gDiffuseObj.AFlags      = uint8(zeros(gDiffuseObj.M,1)); 
    gDiffuseObj.Stop        = false;
    gDiffuseObj.Metric=cExtraParam.Metric; gDiffuseObj.MetricFlags=cExtraParam.MetricFlags; % NEW (points with fixed metric) 
    
    lOpts.Type='gauss';  lOpts.NormalizationType=cExtraParam.NormalizationType; lOpts.WhichEigenVals='small';   % range search
    lOpts.Delta=cExtraParam.Delta; lOpts.Precision=cExtraParam.Precision; lOpts.kNN=cExtraParam.kNN; 
    lOpts.MaxEigenVals=cExtraParam.MaxEigenVals; lOpts.EigEnergyThres=cExtraParam.EigEnergyThres;    
    lOpts.NNsymm=cExtraParam.NNsymm; %'ave';
    
    % Prepare graph:
    %[lNNInfo.NNIdxs lNNInfo.NNDist] = nn_search(cPoints,gDiffuseObj.Atria, cPoints,lOpts.kNN);
    [lGraph, lEigenVecs, lEigenVals] = FastLaplacianDetEigs( cPoints, lOpts,[]); NEigenVals=length(lEigenVals),
    gDiffuseObj.A=lGraph.W; gDiffuseObj.D=lGraph.D;
    % Compute the new normalized transition matrix
    lTemp = spdiags(lGraph.D.^(-1),0,speye(size(lGraph.D,1)));
    lTemp(find(lTemp==inf))=0;
    gDiffuseObj.P = lTemp*gDiffuseObj.A;    

    % global variables from eigenvector calculations
    gEigObj.EigenVecs=lEigenVecs; gEigObj.EigenVals=lEigenVals;
    if  (cExtraParam.EigFlag==1) % store normalized eigenvectors for faster computation
        dsqrt= sqrt(lGraph.D)*ones(1,NEigenVals); 
        gEigObj.EigenVecsL=lEigenVecs./dsqrt;  % eigv/sqrt(d)
        gEigObj.EigenVecsR=lEigenVecs.*dsqrt; clear dsqrt;       %eigv*sqrt(d)
    end
    
elseif strcmpi(cFcn,'diffuse'),
    
    if cExtraParam.EigFlag==1  % Global diffusion using eigenvectors (at time t)
        mu_t=(1-gEigObj.EigenVals).^cExtraParam.t; % (eigenvalues of weight matrix)^t
        if cExtraParam.EigEnergyThres < 1
            mu_2t=mu_t.*mu_t; tmp_perc = cumsum(mu_2t)/sum(mu_2t); % percentage of total energy
            %figure(20), plot(mu_2t,'o'),
            tmp_idx=find(tmp_perc >= cExtraParam.EigEnergyThres); Ncomps=tmp_idx(1), % number of comps to include
        else
            Ncomps=cExtraParam.MaxEigenVals;
        end
        EigenVecsL=gEigObj.EigenVecsL(:,1:Ncomps); EigenVecsR=gEigObj.EigenVecsR(:,1:Ncomps);
        EigenVals=1-gEigObj.EigenVals(1:Ncomps);  % eigenvalues of weight matrix
        
        % Commented out 5/27/04
        % Mark points where the heat is zero at time t-1
        %if (cExtraParam.tstep~=1 & cExtraParam.t>1) % perform calculations at time t-1
        %    Sigma_t=sparse(diag(EigenVals).^(cExtraParam.t-1));
        %    P_t = EigenVecsL*Sigma_t*EigenVecsR'; % Matrix power of transition matrix after normalization; M-by-M matrix
        %    temp_norm=sum(P_t,2); temp_N=diag(1./temp_norm); P_t=temp_N*P_t;  % Normalize the rows of P_t
        %    gDiffuseObj.Heat = cExtraParam.Heat * P_t;  % apply to initial heat distribution
        %end
        %tmp1=max(abs(full(gDiffuseObj.Heat)),[],1)<1E-10; % unreliable calc?? mark points where the heat is zero 
        
        % Calculations at time t 
        Sigma_t=sparse(diag(mu_t(1:Ncomps)));
        %%A_t = EigenVecs*Sigma_t*EigenVecs';  % Matrix power of transition matrix before normalization, A^t=G*Sigma^t*G
        P_t = EigenVecsL*Sigma_t*EigenVecsR'; % Matrix power of transition matrix after normalization
        temp_norm=sum(P_t,2); temp_N=diag(1./temp_norm); P_t=temp_N*P_t;  % Normalize the rows of P_t
        gDiffuseObj.Heat = cExtraParam.Heat * P_t;  % apply to initial heat distribution
        %min(gDiffuseObj.Heat(:))
        gDiffuseObj.Heat(find(gDiffuseObj.Heat<0))=0; % Added, AL 3/19/04 (negative heat because of numerical error)
        check_norm=sum(gDiffuseObj.Heat,2);  gDiffuseObj.Heat=diag(1./check_norm)*gDiffuseObj.Heat; % normalize heat
        tmp2=max(abs(full(gDiffuseObj.Heat)),[],1)>1E-10; % unreliable calc? mark points where the heat is non-zero at time t
        tmp2idx=find(tmp2==1); gDiffuseObj.AFlags(tmp2idx)=uint8(1);
        %gDiffuseObj.BoundaryIdxs= find(tmp1.*tmp2==1);  % find boundary points % commented out 5/27/04
        
    else % step-wise global diffusion; OK!
        tmp1=max(abs(gDiffuseObj.Heat),[],1)<1E-10; % mark points where the heat is zero at time t-1 
        % Apply the transition matrix to the initial heat distribution.
        gDiffuseObj.Heat = gDiffuseObj.Heat * gDiffuseObj.P;  %check_norm=sum(gDiffuseObj.Heat,2),
        tmp2=max(abs(gDiffuseObj.Heat),[],1)>1E-10; % mark points where the heat is non-zero at time t
        tmp2idx=find(tmp2==1); gDiffuseObj.AFlags(tmp2idx)=uint8(1);
        gDiffuseObj.BoundaryIdxs= find(tmp1.*tmp2==1);
    end           
    % Update gDiffuseObj.BoundaryIdxs
elseif strcmpi(cFcn,'plug')==1,
    
    % Plug some points
    figure(1);hold on;plot(cPoints(cExtraParam,1),cPoints(cExtraParam,2),'rx');pause(0.1);
    for lk = 1:length(cExtraParam)
        lIdxs = find( gDiffuseObj.A(cExtraParam(lk),:)~=0 );
        gDiffuseObj.A(cExtraParam(lk),lIdxs) = 0;
        lIdxs = find( gDiffuseObj.A(:,cExtraParam(lk))~=0 );
        gDiffuseObj.A(lIdxs,cExtraParam(lk)) = 0;            
        gDiffuseObj.A(cExtraParam(lk),cExtraParam(lk)) = 1;
    end;
    [lInt lIntIdx] = intersect(gDiffuseObj.BoundaryIdxs,cExtraParam);
    gDiffuseObj.BoundaryIdxs(lIntIdx)=[];
    
end;

return;