function [G] = gsp_community(n, param)
%GSP_COMMUNITY Create a community graph
% Usage: G = gsp_community(N);
% G = gsp_community();
% G = gsp_community(N, param );
%
% Input parameters
% N : Number of nodes (default 256)
% param : Structure of optional parameters
% Output parameters
% G : Graph
%
% This function create a 2 dimentional random sensor graph. All the
% coordonates are between 0 and 1.
%
% param is an optional structure with the following fields
%
% param.Nc : Number of communities (default round(sqrt(N)/2) )
% param.verbose*: display parameter - 0 no log - 1 display the errors
% (default 1)
% param.com_sizes : size of the communities. The sum of the sizes has
% to be equal to N. Leave this field empty if you want random sizes.
% param.min_comm : Minimum size of the community
% (default: round(N / param.Nc / 3) )
% param.min_deg*: Minimum degree of each nodes (default:
% round(param.min_comm / 2)) (NOT WORKING YET!)
% param.size_ratio*: ratio between radius of world and radius of
% communities (default 1)
% param.world_density probability of a random edge between any pair
% of nodes (default 1/N)
%
%
% Example:
%
% G = gsp_community();
% paramplot.show_edges = 1;
% gsp_plot_graph(G,paramplot);
%
%
%
% Url: https://epfl-lts2.github.io/gspbox-html/doc/graphs/gsp_community.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: Vassilis Kalofolias, Nathanael Perraudin, Johan Paratte
% Date: March 2014
% Testing: test_graphs
% TODO: This function has to be revisited
if nargin < 2
param = struct;
end
if nargin < 1
n = 256;
end
if isfield(param, 'com_sizes')
if sum(param.com_sizes) ~= n
error(['GSP_COMMUNITY: The sum of the community sizes has ',...
'to be equal to N']);
end
param.Nc = numel(param.com_sizes);
else
if ~isfield(param, 'Nc'), param.Nc = round(sqrt(n)/2); end
param.com_sizes = [];
end
if ~isfield(param, 'min_comm'), param.min_comm = round(n / param.Nc / 3); end
if ~isfield(param, 'min_deg'), param.min_deg = round(param.min_comm/2); end
if ~isfield(param, 'verbose'), param.verbose = 1; end
if ~isfield(param, 'param.size_ratio'), param.size_ratio = 1; end
if ~isfield(param, 'world_density'), param.world_density = 1/n; end
if isempty(param.com_sizes)
com_lims = sort(randperm(n - (param.min_comm-1) * param.Nc - 1, param.Nc-1), 'ascend');
com_lims = com_lims + cumsum((param.min_comm-1) * ones(size(com_lims)));
com_lims = [0, com_lims, n];
param.com_sizes = diff(com_lims);
else
com_lims = [0, cumsum(param.com_sizes)];
end
if param.verbose > 2
X = zeros(10000, param.Nc + 1);
%pick randomly param.Nc-1 points to cut the rows in communtities:
for ii=1:10000
com_lims_temp = sort(randperm(n - (param.min_comm-1) * param.Nc - 1, param.Nc-1), 'ascend');
com_lims_temp = com_lims_temp + cumsum((param.min_comm-1) * ones(size(com_lims_temp)));
X(ii,:) = [0, com_lims_temp, n];
end
dX = diff(X')';
for ii=1:param.Nc; figure;hist(dX(:,ii), 100); title('histogram of row community size'); end
clear X com_lims_temp
end
rad_world = param.size_ratio * sqrt(n);
com_coords = rad_world * [-cos(2*pi*(1:param.Nc)/param.Nc)', sin(2*pi*(1:param.Nc)/param.Nc)'];
G.coords = ones(n, 2);
% create uniformly random points in the unit disc
for ii = 1:n
% use rejection sampling to sample from a unit disc (probability = pi/4)
while norm(G.coords(ii, :)) >= 1/2
% sample from the square and reject anything outside the circle
G.coords(ii, :) = [rand-.5, rand-.5];
end
end
% add the offset for each node depending on which community it belongs to
info.node_com = zeros(n, 1);
for ii = 1:(param.Nc)
com_size = param.com_sizes(ii);
rad_com = sqrt(com_size);
node_ind = (com_lims(ii) + 1) : com_lims(ii+1);
G.coords(node_ind, :) = bsxfun(@plus, rad_com * G.coords(node_ind, :), com_coords(ii, :));
info.node_com(node_ind) = ii;
end
% TODO: this can (and should to prevent overlap) be done for each community separately!
% D = gsp_distanz(G.coords');
% %R = rmse_mv(G.coords')*sqrt(2);
% %W = graph_k_NN(exp(-D.^2), param.min_deg);
% W = exp(-D.^2);
% W(W<1e-3) = 0;
%TODO: this could be more sophisticated (e.g. one sigma for communities,
%one sigma for inter-community connections, exp(-d^2/sigma) weights!
% W = W + abs(sprandsym(N, param.world_density));
% W = double(abs(W) > 0);
% G.W = sparse(W);
% % Fast (and scalable) implementation of the above
% kdt = KDTreeSearcher(G.coords, 'distance', 'euclidean');
% epsilon = sqrt(-log(1e-3));
% [NN, D] = rangesearch(kdt, G.coords, epsilon, 'distance', 'euclidean' );
%
% %Counting non-zero elements
% count = 0;
% for ii = 1:N
% count = count + length(NN{ii}) - 1;
% end
%
% spi = zeros(count,1);
% spj = zeros(count,1);
% spv = ones(count,1);
% start = 1;
%
%
% % Fill the 3-col values with [i, j, exp(-d(i,j)^2 / sigma)]
% for ii = 1:N
% len = length(NN{ii}) - 1;
% spi(start:start+len-1) = repmat(ii, len, 1);
% spj(start:start+len-1) = NN{ii}(2:end)';
% % spv(start:start+len-1) = exp(-D{ii}(2:end).^2);
% start = start + len;
% end
%
% W = sparse(spi, spj, spv, N, N);
paramnn.rescale = 0;
paramnn.center = 0;
paramnn.type = 'radius';
epsilon = sqrt(-log(1e-3));
paramnn.epsilon = epsilon;
[spi, spj] = gsp_nn_distanz(G.coords',G.coords',paramnn);
W = sparse(spi, spj, ones(size(spi)), n, n);
W = (W+W')/2;
% Adding the sparse rand connections
W = W + abs(sprandsym(n, param.world_density));
W = double(abs(W) > 0);
% get rid of self loops:
W(1:n+1:end) = 0;
G.W = W;
G.type = 'community';
% return additional info about the communities
info.com_lims = com_lims;
info.com_coords = com_coords;
info.com_sizes = param.com_sizes;
G.info = info;
G = gsp_graph_default_parameters( G );
end