function Diffuser_ToNeighbors( cFcn, cPoints, cExtraParam )

%
% function Diffuser_ToNeighbors( 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
%                 Precision : precision for the calculations
%                   
%
% 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.
% USES:
%   nn_prepare,range_search
%
% SOURCE CONTROL:
%   Original code by Mauro Maggioni, Feb 2004         [growing matrix, slow because of multiple memory allocations]
%   Modified by Ann Lee, 2/27/04                      [global matrix, simplified indexing, still slow]
%   Modified by Mauro Maggioni, 3/1/04                [fixed bugs, speed up by incremental filling, reshaped according to generic diffusor template, standardized code]
%   Modified, AL, 3/12/04                             [multiple classes; bugs with Aflags and gDiffuseObj.Stop fixed; AIdxs reomoved]
%   Modified, 3/23/04                                 [new variables: gDiffuseObj.Metric and gDiffuseObj.MetricFlags]

global gDiffuseObj;         % Thanks to Matlab which hasn't got pointers...

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.Metric{size(cPoints,1)}=struct([]); % MM, 2/26/04
    %gDiffuseObj.ModFlags      = uint8(zeros(gDiffuseObj.M,1)); % AL 3/23/04 (=1 if metric changed)
    gDiffuseObj.Metric = cExtraParam.Metric; gDiffuseObj.MetricFlags=cExtraParam.MetricFlags;    % NEW!
    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.GetWeightsFcn = cExtraParam.GetWeightsFcn;
    gDiffuseObj.Epsilon = cExtraParam.Epsilon;  % MM, 3/26/04
    
elseif strcmpi(cFcn,'diffuse'),
    
    % Diffusion    
    
    lNewBoundary=[]; lFlags=uint8(zeros(gDiffuseObj.M,1));  lFlags(gDiffuseObj.BoundaryIdxs)=uint8(1); % lFlags prevents duplicates
    %lDidDiffuse = false;           
   
    
    % Find the new points to diffuse to from the current boundary points
    %if length(gDiffuseObj.BoundaryIdxs) == 0, return;end;   
    if ~isempty(gDiffuseObj.BoundaryIdxs)                            % if new boundary pts, then update A and P, else only diffuse
        % Find neighbors of current boundary points => NNIdxs, Count
        lBdryPoints = cPoints(gDiffuseObj.BoundaryIdxs,:); 
        % NNIdxs are indices into cPoints
        %[Count NNIdxs] = range_search( cPoints, gDiffuseObj.Atria, lBdryPoints, 2*gDiffuseObj.RangeRadius, 0 );
        [Count NNIdxs] = range_search( cPoints, gDiffuseObj.Atria, lBdryPoints, gDiffuseObj.RangeRadius, 0 );
        % If no neighbors to boundary pts, return saying the diffusion has stopped
        if Count==0,    % If no new neighbors, diffuse one step, then exit
            gDiffuseObj.Heat=gDiffuseObj.Heat*gDiffuseObj.P; 
            gDiffuseObj.BoundaryIdxs = [];
            return; 
        end;    
        
        % Go through the list of boundary points
        for lk = 1:size(lBdryPoints,1)
            
            % Find which points close to the boundary points are not already (completely) done   
            lNewIdxsI=find(gDiffuseObj.AFlags(NNIdxs{lk,1})==0);
            if length(lNewIdxsI) == 0, gDiffuseObj.AFlags(gDiffuseObj.BoundaryIdxs(lk))=1;continue;end;
            %lDidDiffuse = true;                      
            lNewIdxs=NNIdxs{lk,1}(lNewIdxsI);
            
            % Add the indices of these points to the matrix indices
            %gDiffuseObj.AIdxs = [gDiffuseObj.AIdxs lNewIdxs]; 
            %gDiffuseObj.AFlags(lNewIdxs)=1;  
            
            % These new indices in the matrix are going to be the next boundary points 
            AddIdxs =  lNewIdxs(find(lFlags(lNewIdxs)==0));  
            lNewBoundary=[lNewBoundary, AddIdxs]; lFlags(AddIdxs)=uint8(1);
            
            % Compute the weights of the links from the current boundary points to the new boundary points
            lWeightsExtraParam.Delta = gDiffuseObj.Delta;
            lWeightsExtraParam.Epsilon = gDiffuseObj.Epsilon;
               
            if 1 % NEW!
                % Isotropic diffusion (does not recalculate weights)
                lWeightsExtraParam.Distances = NNIdxs{lk,2}(lNewIdxsI);
                %lW = feval(gDiffuseObj.GetWeightsFcn,cPoints,lNewIdxs,lBdryPoints(lk),lWeightsExtraParam);
                if gDiffuseObj.MetricFlags(gDiffuseObj.BoundaryIdxs(lk))==1 % fixed metric
                    lW = feval(gDiffuseObj.GetWeightsFcn,cPoints,lNewIdxs,lBdryPoints(lk,:),lWeightsExtraParam,gDiffuseObj.Metric{gDiffuseObj.BoundaryIdxs(lk)});
                else % identitity matrix defines metric       
                    lW = feval(gDiffuseObj.GetWeightsFcn,cPoints,lNewIdxs,lBdryPoints(lk,:),lWeightsExtraParam,[]); % weights between boundary point k and neighbors        
                end
                gDiffuseObj.A(gDiffuseObj.BoundaryIdxs(lk),lNewIdxs) = lW;
                gDiffuseObj.A(lNewIdxs,gDiffuseObj.BoundaryIdxs(lk)) = lW';
            else
                % Anisotropic diffusion or isotropic diffusion
                lWeightsExtraParam.Distances = NNIdxs{lk,2}; % AL 4/23/04
                if gDiffuseObj.MetricFlags(gDiffuseObj.BoundaryIdxs(lk))==1 % fixed metric
                    lW = feval(gDiffuseObj.GetWeightsFcn,cPoints,lNewIdxs,lBdryPoints(lk,:),lWeightsExtraParam,gDiffuseObj.Metric{gDiffuseObj.BoundaryIdxs(lk)});
                else  % estimate metric       
                    lW = feval(gDiffuseObj.GetWeightsFcn,cPoints,lNewIdxs,lBdryPoints(lk,:),lWeightsExtraParam,[]); % weights between boundary point k and neighbors        
                end
                gDiffuseObj.A(gDiffuseObj.BoundaryIdxs(lk),lNewIdxs) = lW(lNewIdxsI);   % stuff the weights in the matrix
                gDiffuseObj.A(lNewIdxs,gDiffuseObj.BoundaryIdxs(lk)) = lW(lNewIdxsI)';
            end   
            gDiffuseObj.AFlags(gDiffuseObj.BoundaryIdxs(lk))=1;
        end;
        
        % If didn't diffuse, turn on the stop flag and return
        %if lDidDiffuse==false, gDiffuseObj.Stop=true;return;end;   
        
        % Update the boundary of the diffusion
        gDiffuseObj.BoundaryIdxs = lNewBoundary;
        
        % Compute the new normalized transition matrix
        lD = sum(gDiffuseObj.A,2);
        lTemp = spdiags(lD.^(-1),0,speye(size(lD,1)));
        lTemp(find(lTemp==inf))=0;
        gDiffuseObj.P = lTemp*gDiffuseObj.A;    
    end
    % Apply the transition matrix to the current heat distribution.
    gDiffuseObj.Heat = gDiffuseObj.Heat * gDiffuseObj.P;  
    %tmp=sort(gDiffuseObj.Heat,2); size(tmp),
    %figure(11), plot(1:gDiffuseObj.M, tmp(1,:)); drawnow,
    %size(gDiffuseObj.AIdxs),
    
elseif strcmpi(cFcn,'plug')==1,
    
    % Plug the points connected to the specified points
    for lk = 1:length(cExtraParam)
       lIdxs = find( gDiffuseObj.A(cExtraParam(lk),:)~=0 );
       for lj = 1:length(lIdxs)
           lIdxs2 = find(gDiffuseObj.A(lIdxs(lj),:)~=0);
           lIdxs2 = intersect(lIdxs,gDiffuseObj.BoundaryIdxs);
           figure(1);hold on;plot(cPoints(lIdxs2,1),cPoints(lIdxs2,2),'rx');
           gDiffuseObj.Heat(lIdxs2) = 0;
           gDiffuseObj.BoundaryIdxs = setdiff(gDiffuseObj.BoundaryIdxs,lIdxs2);
       end;
    end;

    pause;
    
    return;
    
    figure(1);hold on;plot(cPoints(cExtraParam,1),cPoints(cExtraParam,2),'rx');pause(0.1);
    gDiffuseObj.Heat(cExtraParam)=0;
    gDiffuseObj.BoundaryIdxs = setdiff(gDiffuseObj.BoundaryIdxs,cExtraParam);
    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;