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

% function [vGraph,vNNInfo]= FastGraphFromData( cX, cOpts, cType, cF )
%
% IN:
%   cOpts       : structure containing:
%                   kNN     : number of nearest neighbors if cType=='nn'
%                   NNsymm  : how to symmetrize the weight matrix if cType=='nn'
%                             'ave' for 1/2(W+W*)
%                             'leftadj' for W*W
%                             'rightadj' for WW*
%                   Delta   : width of gaussian if cType=='gauss'
%                   Precision : determines the range_search radius when cType=='gauss'
%   cType       : 'nn' or 'gauss'
%   cF          : the function whose gradient determines the weights.
%
% OUT:
%   vGraph      : sparse graph obtained from the data
%   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]
%

% Pre-process points
vNNInfo.Atria = nn_prepare( cX );

if strcmpi(cType,'nn') == 1
    % Perform the nn search
    [vNNInfo.NNIdxs vNNInfo.NNDist] = nn_search( cX, vNNInfo.Atria, cX, cOpts.kNN );
    lN = size(vNNInfo.NNIdxs,1);
    % Build the weight matrix
    vGraph.W = sparse(reshape(ones(cOpts.kNN,1)*(1:lN),[cOpts.kNN*lN 1]),reshape(vNNInfo.NNIdxs(:,1:cOpts.kNN)',[cOpts.kNN*lN 1]),ones(cOpts.kNN*lN,1),lN,lN,cOpts.kNN*lN);
    % Symmetrize it
    try
        if strcmpi(cOpts.NNsymm,'ave')==1,
            vGraph.W = (vGraph.W'+vGraph.W)/2;
        elseif strcmpi(cOpts.NNsymm,'leftadj')==1,
            vGraph.W = vGraph.W'*vGraph.W;
        elseif strcmpi(cOpts.NNsymm,'rightadj')==1,
            vGraph.W = vGraph.W*vGraph.W';
        end;
    catch
        vGraph.W = (vGraph.W'+vGraph.W)/2;
    end;
else
    lNumberOfPoints = size(cX,1);
    % Compute the range for requested precision            
    lRange = cOpts.Delta*(log(1/cOpts.Precision))^(0.5);
    % Perform the range search (this could be sped up, we don't need all pairwise computation
    [vNNInfo.NNCounts lTemp] = range_search( cX, vNNInfo.Atria, 1:lNumberOfPoints, lRange, -1);
    % Matlab is stupid in handling arrays of cell, so waste some memory and do this trick! The performance is WAY better (factors of 100s...!!!!) in this way (sigh)...
    lMaxSize = max(vNNInfo.NNCounts);
    lTotalSize = sum(vNNInfo.NNCounts);
    lIdx = 1;
    lIdxsI = zeros(1,lTotalSize);
    lIdxsJ = zeros(1,lTotalSize);
    lEntries = zeros(1,lTotalSize);
    % Check if F is specified
    if (nargin<4) | (isempty(cF))
        lFSpecified = false;
    else
        lFSpecified = true;
        lDeltaSquared = cOpts.Delta^2;
    end;
    % Re-organize the stuff returned from the range search.
    for lk = 1:lNumberOfPoints
        vNNInfo.NNIdxs{lk} = lTemp{lk,1};
        vNNInfo.NNDist{lk} = lTemp{lk,2};
        lIdxsI(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = lk;
        lIdxsJ(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = vNNInfo.NNIdxs{lk};
        if ~lFSpecified
            lEntries(lIdx:lIdx+vNNInfo.NNCounts(lk)-1) = exp(-(vNNInfo.NNDist{lk}/cOpts.Delta).^2);
        else
            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;
    clear lTemp;    
    % Fill-in the matrix.
    vGraph.W = sparse(lIdxsI,lIdxsJ,lEntries,lNumberOfPoints,lNumberOfPoints);
    %        vGraph.W = spalloc(lNumberOfPoints,lNumberOfPoints,sum(vNNInfo.NNCounts));
    %        for lk = 1:lNumberOfPoints                     
    %           vGraph.W(lk,vNNInfo.NNIdxs{lk}) = exp(-(vNNInfo.NNDist{lk}/cOpts.Delta).^2);            % This is unbelievably slow!! It is NOT the loop, but copying things into a SPARSE matrix (it'd be fast if the matrix is full)
    %        end;
end;

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

return;