function [ comm_radius_array ] = comm_radius(partitions,locations)

% COMM_RADIUS - function to calculate the radius 
%         of detected communities
%
% Inputs: 
%         'partitions'    --- an Nx1 array where N is the number of nodes
%         in the network.  Each entry contains the community index for
%         node i (where communities are sequentially indexed as 1:M and M
%         is the total number of detected communities).
% 
%         'locations'   ---  a Nxdim array where N is the number of nodes and
%         dim is the spatial dimensionality of the network (1,2,or 3).  The
%         columns contain the (x,y,z) coordinates of the nodes in euclidean
%         space
% 
% Outputs: 
%         comm_radius_array   --- an (M+1)x2 array where M is the number
%         of communities.  The first column contains the community index and
%         the second column contains radius of that community. 
%         The M+1th entry contains the average community radius of the
%         network (denoted as community 0).
%  
% Written by Sarah Feldt Muldoon

%Reference:  Lohse C, Bassett DS, Lim KO, Carlson JM (2014) Resolving 
%         anatomical and functional structure in human brain organization: 
%         identifying mesoscale organization in weighted network representations. 
%         PLoS Comput Biol 10:e1003712.

number_nodes = length(partitions);
number_communities = max(partitions);
comm_radius_array = zeros(number_communities+1, 2);
comm_radius_array(:,1) = [1:number_communities 0];
number_nodes_in_communities = zeros(number_communities, 1);

% calculate the radius for each community
for i = 1:number_communities
    community_i_nodes = find(partitions==i);
    number_nodes_in_i = length(community_i_nodes);
    number_nodes_in_communities(i) = number_nodes_in_i;
    if number_nodes_in_i >= 2
        position_vectors_i = locations(community_i_nodes,:);
        std_position_vectors_i = std(position_vectors_i,1);
        radius_i = sqrt(sum(std_position_vectors_i.^2,2));
        comm_radius_array(i,2) = radius_i;
    end
end


%calculate the average community radius for the partition
std_all_position_vectors = std(locations,1);
radius_all = sqrt(sum(std_all_position_vectors.^2));
average_community_radius = 1 / (number_nodes*radius_all) * sum(number_nodes_in_communities .* comm_radius_array(1:number_communities, 2));
comm_radius_array(number_communities + 1,2) = average_community_radius;


end

