function [ comm_laterality_array ] = comm_laterality(partitions,categories)

% COMM_LATERALITY - function to calculate the laterality 
% of detected communities

% Inputs: 
%         'partitions'    --- an Nx1 array where N is the number of nodes
%         in the network.  Each entry contains the communitiy index for
%         node i (where communities are sequentially indexed as 1:M and M
%         is the total number of detected communities).
% 
%         'categories'   ---  a Nx1 array where N is the number of nodes and
%         each entry is either a '0' or '1' denoting the assignment of each 
%         node to one of two communities.
%
% Outputs: 
%         comm_laterality_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 the laterality of that community. 
%         The M+1th entry contains the  laterality of the partition of the
%         network (denoted as community 0).
%  
% Written by Sarah Feldt Muldoon

%Reference: Doron KW, Bassett DS, Gazzaniga MS (2012) Dynamic network 
%         structure of interhemispheric coordination. Proc Natl Acad Sci USA 
%         109:18661?18668.



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

%calcuate the laterality 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,1) = number_nodes_in_i;
    categories_i = categories(community_i_nodes);
    nodes_in_comm_0 = find(categories_i == 0);
    number_nodes_in_comm_0 = length(nodes_in_comm_0);
    number_nodes_in_communities(i,2) = number_nodes_in_comm_0;
    nodes_in_comm_1 = find(categories_i == 1);
    number_nodes_in_comm_1 = length(nodes_in_comm_1);
    number_nodes_in_communities(i,3) = number_nodes_in_comm_1;
    laterality_i = abs(number_nodes_in_comm_0 - number_nodes_in_comm_1) / number_nodes_in_i;
    comm_laterality_array(i,2) = laterality_i;
end


%calcuate the laterality for the network partition

%need to calcuated "expected" laterality from each community from surrogate
%data
total_nodes_assignments = sum(number_nodes_in_communities, 1)
total_number_in_comm_1 = total_nodes_assignments(3);
number_surrogates = 1000;
surrogate_laterality = zeros(number_communities,number_surrogates);
for j = 1:number_surrogates
    randomized_categories = zeros(1,number_nodes);
    rand_perm_assignment = randperm(number_nodes);
    place_in_comm_1 = rand_perm_assignment(1:total_number_in_comm_1);
    randomized_categories(place_in_comm_1) = 1;
    %now calculate the community laterality for the randomized data
    for i = 1:number_communities
        community_i_nodes = find(partitions == i);
        number_nodes_in_i = length(community_i_nodes);
        randomized_categories_i = randomized_categories(community_i_nodes);
        rand_nodes_in_comm_0 = find(randomized_categories_i == 0);
        rand_number_nodes_in_comm_0 = length(rand_nodes_in_comm_0);
        rand_nodes_in_comm_1 = find(randomized_categories_i == 1);
        rand_number_nodes_in_comm_1 = length(rand_nodes_in_comm_1);
        rand_laterality_i = abs(rand_number_nodes_in_comm_0 - rand_number_nodes_in_comm_1) / number_nodes_in_i;
        surrogate_laterality(i,j) = rand_laterality_i;
    end
end
expected_comm_laterality = sum(surrogate_laterality, 2) / number_surrogates;

network_laterality = (1 / number_nodes) * sum(number_nodes_in_communities(:, 1) .* comm_laterality_array(1:number_communities, 2) - expected_comm_laterality);
comm_laterality_array(number_communities + 1, 2) = network_laterality;

end
