function [vGraph,vNNInfo]= FastGraphFromData( cX, cOpts, cType, cF )

% function [vGraph,vNNInfo]= FastGraphFromData( cX, cOpts, cType, cF )
%
% IN:
%   cType         : 'nn' or 'gauss'
%   cOpts         : structure containing:
%                     NNsymm  : how to symmetrize the weight matrix (empty otherwise)
%                               'ave' for 1/2(W+W*)
%                               'leftadj' for W*W
%                               'rightadj' for WW*
%                     if cType == 'nn':
%                               kNN     : number of nearest neighbors if cType=='nn'
%                               Delta         : width of gaussian, if Delta == Inf, then uses 1 for all edge weights
%                     if cType == 'gauss'
%                               Delta         : width of gaussian
%                               Precision     : determines the range_search radius when cType=='gauss'
%                     Epsilon   : no diffusion outside any local metric ball of radius Epsilon (used only by flexible metrics)
%   GetWeightsFcn : pointer to function to compute weights, as in GetWeightsFcn_Template.m
%   cF            : the function whose gradient determines the weights.
%
% OUT:
%   vGraph      : sparse graph obtained from the data. Contains the fields D and W. W is the adjancency matrix (with weights), D is just the vector of row sums of W.
%   vNNInfo     : structure containing the following.
%                   Atria  : tree for fast nn and range searches
%                 If cType == 'nn':
%                   NNIdxs  : M by cOpts.kNN matrix whose (i,j) entry is the j-th nn to the i-th point
%                   NNDist  : M by cOpts.kNN matrix whose (i,j) entry is the distance of the j-th nn to the i-th point (these are sorted in increasing order
%                 If cType == 'gauss':
%                   NNCounts : M column vector whose i-th entry is the number of \epsilon neighbors of the i-th point
%                   NNIdxs   : cell array whose i-th cell contains the indices of the \epsilon neighbors of the i-th point
%                   NNDist   :  cell array whose i-th cell contains the distances of the \epsilon neighbors of the i-th point from the i-th point, NOT sorted
%
% USES:
%   nn_prepare, nn_search, range_search
%
%
% SOURCE CONTROL:
%   MM      03/13/04    [it's pretty good now, fast and includes range_searches]
%   MM      05/31/04    [added generic weight function also for nn, reorganized code, added epsilon for calls to flexible metrics, updated doc]
%   MM      07/06/04    [added some parameter checking]


% Check if a weight function is specified
if isfield(cOpts,'GetWeightsFcn')
    lWeightsFcnSpecified = true;
    if ~isfield(cOpts,'Epsilon')
        cOpts.Epsilon = 0;
    end;
else
    lWeightsFcnSpecified = false;
end;
% Number of points
lNumberOfPoints = size(cX,1);   
% Pre-process points
vNNInfo.Atria = nn_prepare( cX );
% Check if vector field F is specified
if (nargin<4) | (isempty(cF))
    lFSpecified = false;
else
    lFSpecified = true;
    lDeltaSquared = cOpts.Delta^2;
end;

if strcmpi(cType,'nn') == 1
    % Perform the nn search
    [vNNInfo.NNIdxs vNNInfo.NNDist] = nn_search( cX, vNNInfo.Atria, cX, cOpts.kNN );
    % Total number of entries in the W matrix
    lTotalSize = cOpts.kNN*lNumberOfPoints;
    % We need to create the sparse W matrix all at once (Matlab is VERY slow in adding one element to a sparse matrix)
    % Allocate memory for needed indices
    lIdx = 1;lIdxsI=[];lIdxsJ=[];lIdxsI(lTotalSize) = int16(0);lIdxsJ(lTotalSize) = int16(0);lEntries = zeros(1,lTotalSize);
    % Re-organize the stuff returned from the nn search.
    for lk = 1:lNumberOfPoints
        % Build the sparse row and column index sets
        lIdxsI(lIdx:lIdx+cOpts.kNN-1) = lk;
        lIdxsJ(lIdx:lIdx+cOpts.kNN-1) = vNNInfo.NNIdxs(lk,:);
        if ~lFSpecified 
            % No vector field specified
            if lWeightsFcnSpecified
                % A weight function is specified
                lExtraParam.Delta = cOpts.Delta;
                lExtraParam.Distances = vNNInfo.NNDist(lk,:);
                lExtraParam.Epsilon = cOpts.Epsilon;
                % Call the weight function
                try
                    lW = feval( cOpts.GetWeightsFcn,cX,vNNInfo.NNIdxs(lk,:),cX(lk,:), lExtraParam );
                    lEntries(lIdx:lIdx+cOpts.kNN-1) = lW;
                catch
                    warning('FastGraphFromData:','Could not call the GetWeightsFcn, using default Gaussian weights');
                    % No weight function is specified: use standard Gaussian weight
                    lEntries(lIdx:lIdx+cOpts.kNN-1) = exp(-(vNNInfo.NNDist(lk,:)/cOpts.Delta).^2);
                end;
            else
                % No weight function is specified: use standard Gaussian weight
                lEntries(lIdx:lIdx+cOpts.kNN-1) = exp(-(vNNInfo.NNDist(lk,:)/cOpts.Delta).^2);
            end;
        else
            % A vector field is specified: flow along it
            lEntries(lIdx:lIdx+cOpts.kNN-1) = exp(-((vNNInfo.NNDist(lk,:)).^2+abs(cF(vNNInfo.NNIdxs(lk,:))-cF(lk)))/lDeltaSquared);
        end;
        lIdx = lIdx + cOpts.kNN;
    end;    
else
    % Compute the range for requested precision            
    if isfield(cOpts,'Precision')
        if cOpts.Precision <= 0, cOpts.Precision = 1e-4;end;
    else
        cOpts.Precision = 1e-4;
    end;
    lRange = cOpts.Delta*(log(1/cOpts.Precision))^(0.5);
    % Perform the range search
    [vNNInfo.NNCounts lTemp] = range_search( cX, vNNInfo.Atria, 1:lNumberOfPoints, lRange, -1);
    % Total number of entries in the W matrix
    lTotalSize = sum(vNNInfo.NNCounts);
    % We need to create the sparse W matrix all at once (Matlab is VERY slow in adding one element to a sparse matrix)
    % Allocate memory for needed indices
    lIdx = 1;lIdxsI=[];lIdxsJ=[];lIdxsI(lTotalSize) = int16(0);lIdxsJ(lTotalSize) = int16(0);lEntries = zeros(1,lTotalSize);    
    % Prepare the vNNInfo structure as nn-search does
    for lk = 1:lNumberOfPoints
        vNNInfo.NNIdxs{lk} = lTemp{lk,1};
        vNNInfo.NNDist{lk} = lTemp{lk,2};
    end;
    % Re-organize the stuff returned from the range search.
    for lk = 1:lNumberOfPoints
        % Build the sparse row and column index sets
        lIdxsI(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = lk;
        lIdxsJ(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = vNNInfo.NNIdxs{lk};
        if ~lFSpecified 
            % No vector field specified
            if lWeightsFcnSpecified
                % A weight function is specified
                lExtraParam.Delta = cOpts.Delta;
                lExtraParam.Distances = vNNInfo.NNDist{lk};
                lExtraParam.Epsilon = cOpts.Epsilon;
                % Call the weight function
                lW = feval( cOpts.GetWeightsFcn,cX,vNNInfo.NNIdxs{lk},cX(lk,:), lExtraParam );
                lEntries(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = lW;
            else
                % No weight function is specified: use standard Gaussian weight
                if cOpts.Delta == Inf
                    lEntries(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = 1;
                else
                    lEntries(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = exp(-(vNNInfo.NNDist{lk}/cOpts.Delta).^2);
                end;
            end;
        else
            % A vector field is specified: flow along it
            lEntries(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = exp(-((vNNInfo.NNDist{lk}).^2+abs(cF(vNNInfo.NNIdxs{lk})-cF(lk)))/lDeltaSquared);
        end;
        lIdx = lIdx + vNNInfo.NNCounts(lk);
    end;
end;
clear lTemp;    
% Fill-in the matrix.
vGraph.W = sparse(lIdxsI,lIdxsJ,lEntries,lNumberOfPoints,lNumberOfPoints);
% Symmetrize if requested    
if isfield(cOpts,'NNsymm'),   
    % Symmetrize it
    [lI,lJ,lV] = find(vGraph.W);
    lGraphWT = sparse(lJ,lI,lV,size(vGraph.W,1),size(vGraph.W,2),length(lI));
    try
        if strcmpi(cOpts.NNsymm,'ave')==1,
            vGraph.W = (lGraphWT+vGraph.W)/2;
        elseif strcmpi(cOpts.NNsymm,'leftadj')==1,
            vGraph.W = lGraphWT*vGraph.W;
        elseif strcmpi(cOpts.NNsymm,'rightadj')==1,
            vGraph.W = vGraph.W*lGraphWT;
        end;
    catch
        % Default is to average with adjoint
        vGraph.W = (lGraphWT+vGraph.W)/2;
    end;
end;

vGraph.D = sum(vGraph.W,2);

return;