function [ G ] = gsp_nn_graph(Xin, param)
%GSP_NN_GRAPH Create a nearest neighbors graph from a point cloud
%   Usage :  G = gsp_nn_graph( Xin );
%            G = gsp_nn_graph( Xin, param );
%
%   Input parameters:
%       Xin         : Input points
%       param       : Structure of optional parameters
%
%   Output parameters:
%       G           : Resulting graph
%
%   'gsp_nn_graph( Xin, param )' creates a graph from positional data. The points are 
%   connected to their neighbors (either belonging to the k nearest 
%   neighbors or to the epsilon-closest neighbors. 
%
%   Example:
%
%           P = gsp_pointcloud('bunny');
%           param.type = 'knn';
%           G = gsp_nn_graph(double(P), param);
%           gsp_plot_graph(G);
%
%   Additional parameters
%   ---------------------
%
%    param.type      : ['knn', 'radius']   the type of graph (default 'knn')
%    param.use_flann : [0, 1]              use the FLANN library (default 0)
%    param.use_full  : [0, 1] - Compute the full distance matrix and then
%     sparsify it (default 0) 
%    param.center    : [0, 1]              center the data (default 0)
%    param.rescale   : [0, 1]              rescale the data on a 1-ball (def 0)
%    param.sigma     : float               the variance of the distance kernel
%    param.k         : int                 number of neighbors for knn (def 10)
%    param.epsilon   : float               the radius for the range search
%    param.use_l1    : [0, 1]              use the l1 distance (def 1)
%    param.symmetrize_type*: ['average','full'] symmetrization type (default 'full')
%    param.min_weight*: float               constant additive weight for each edge (default 0) 
%    param.zero_diagonal*: [0, 1]           zero out the diagonal (default 1)
%    param.weight_kernel*: function         edge-weighting kernel (default gaussian)
%
%
%   See also: gsp_nn_distanz gsp_pointcloud
%
%
%   Url: https://epfl-lts2.github.io/gspbox-html/doc/graphs/gsp_nn_graph.html

% Copyright (C) 2013-2016 Nathanael Perraudin, Johan Paratte, David I Shuman.
% This file is part of GSPbox version 0.7.5
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% If you use this toolbox please kindly cite
%     N. Perraudin, J. Paratte, D. Shuman, V. Kalofolias, P. Vandergheynst,
%     and D. K. Hammond. GSPBOX: A toolbox for signal processing on graphs.
%     ArXiv e-prints, Aug. 2014.
% http://arxiv.org/abs/1408.5781

% Author: Johan Paratte, Nathanael Perraudin
% Date: 16 June 2014
% Testing: test_rmse

    if nargin < 2
    % Define parameters
        param = {};
    end
    
    %Parameters
    if ~isfield(param, 'type'), param.type = 'knn'; end
    if ~isfield(param, 'use_flann'), param.use_flann = 0; end
    if ~isfield(param, 'center'), param.center = 0; end
    if ~isfield(param, 'zero_diagonal'), param.zero_diagonal = 1; end
    if ~isfield(param, 'rescale'), param.rescale = 0; end
    if ~isfield(param, 'k'), param.k = 10; end
    if ~isfield(param, 'epsilon'), param.epsilon = 0.01; end
    if ~isfield(param, 'use_l1'), param.use_l1 = 0; end
    if ~isfield(param, 'target_degree'), param.target_degree = 0; end;
    if ~isfield(param, 'symmetrize_type'), param.symmetrize_type = 'average'; end
    if ~isfield(param, 'min_weight'), param.min_weight = 0; end
    if ~isfield(param, 'light'); param.light = 0; end
    if ~isfield(param, 'weight_kernel'); param.weight_kernel = @(x, sigma) exp(-x/sigma); end
    
    paramnn = param;
    paramnn.k = param.k +1;
    kdist = @(x1,x2) gsp_nn_distanz(x1',x2',paramnn);
    [indx, indy, dist, Xout, ~, epsilon] = kdist(Xin,Xin);
    Xout = transpose(Xout);
    switch param.type
        case 'knn'
            if param.use_l1
                if ~isfield(param, 'sigma'), param.sigma = mean(dist); end
            else
                if ~isfield(param, 'sigma'), param.sigma = mean(dist)^2; end
            end
        case 'radius'
            if param.use_l1
                if ~isfield(param, 'sigma'), param.sigma = epsilon/2; end
            else
                if ~isfield(param, 'sigma'), param.sigma = epsilon.^2/2; end
            end
        otherwise
            error('Unknown graph type')
    end
    
    
    if param.use_l1
        Wmat = @(indx,indy,dist,n,m) sparse(indx, indy, param.min_weight+double(param.weight_kernel(dist, param.sigma)), n, m);
    else
        Wmat = @(indx,indy,dist,n,m) sparse(indx, indy, param.min_weight+double(param.weight_kernel(dist.^2, param.sigma)), n, m);
    end
    
    n = size(Xin,1);
    W = Wmat(indx,indy,dist,n,n);
    
    Wkernel = @(x) create_kernel(Xin,x,kdist,Wmat);
    
    if param.zero_diagonal
        % We need zero diagonal
        W(1:(n+1):end) = 0;     % W = W-diag(diag(W));
    end
    
    % Computes the average degree when using the epsilon-based neighborhood
    if (strcmp(param.type,'radius'))
        text = sprintf('Average number of connection = %d', nnz(W)/size(W, 1));
        disp(text);
    end

    % Sanity check
    if size(W,1) ~= size(W,2), error('Weight matrix W is not square'); end
    
    % Symmetry checks
    %if issymmetric(W)
    if (norm(W - W', 'fro') == 0)
        disp('The matrix W is symmetric');
    else
         W = gsp_symmetrize(W,param.symmetrize_type);
    end
    
    %Fill in the graph structure
    G.N = n;
    G.W = W;
    G.coords = Xout;
    G.Wkernel = Wkernel;
    %G.limits=[-1e-4,1.01*max(x),-1e-4,1.01*max(y)];
    if param.use_l1
        G.type = 'nearest neighbors l1';
    else
        G.type = 'nearest neighbors';
    end
    %G.vertex_size=30;
    G.sigma = param.sigma;
    if param.light
        G = gsp_graph_lightweight_parameters(G);
    else
        G = gsp_graph_default_parameters(G);
    end
end

function W = create_kernel(Xin,x,kdist,Wmat)
[indx, indy, dist] = kdist(Xin,x);
n = size(Xin,1);
m = size(x,1);
W = Wmat(indx,indy,dist,n,m)';

end

