function [SWP,delta_C,delta_L] = small_world_propensity(A, varargin)

% a function for calculating the small world propensity of
% a given network - assumes that matrix is undirected (symmeteric) and if
% not, creates a symmetric matrix which is used for the calculations

%NOTE:  This code requires the Bioinformatics Toolbox to be installed
%        (uses graphallshortestpaths.m)

%Inputs:
%   A           the connectivity matrix, weighted or binary
%   varargin    a string corresponding to the method of clustering
%               to be used, where 'O' is Onnela, 'Z' is Zhang, 
%               'B' is Barrat, 'bin' is binary (default is Onnela).
%               If the user specifies binary analysis, a 
%               weighted matrix will be converted to a binary matrix 
%               before proceeding.
%        

%Outputs:
%   SWP         the small world propensity of the matrix
%   delta_C     the fractional deviation from the expected culstering coefficient of a
%                   random network
%   delta_L     the fractional deviation from the expected path length of a
%                   random network

%written by Eric Bridgeford and modified by Sarah F. Muldoon

% Reference: Muldoon, Bridgeford, and Bassett (2015) "Small-World Propensity in Weighted, 
%               Real-World Networks" http://arxiv.org/abs/1505.02194

if isempty(varargin)
    varargin{1} = 'O';
end

if sum(sum(A)) > 0
    
bin_matrix = 0;
if strcmp(varargin{1},'bin') == 1
   bin_matrix = 1;
   A = A > 0;
end

%check to see if matrix is symmeteric
symcheck=abs(A-A');
if sum(sum(symcheck)) > 0
    % adjust the input matrix to symmeterize
    disp('Input matrix is not symmetric. Symmetrizing.')
    W = symm_matrix(A, bin_matrix);
else
    W=A;
end

%calculate the number of nodes
n = length(W);  
%compute the weighted density of the network
dens_net = sum(sum(W))/(max(max(W))*n*(n-1));

%compute the average degree of the unweighted network, to give
%the approximate radius
numb_connections = length(find(W>0));
avg_deg_unw = numb_connections/n;
avg_rad_unw = avg_deg_unw/2;
avg_rad_eff = ceil(avg_rad_unw);


%compute the regular and random matrix for the network W
W_reg = regular_matrix_generator(W, avg_rad_eff);
W_rand = randomize_matrix(W);

%compute all path length calculations for the network
reg_path = avg_path_matrix(1./W_reg);      %path of the regular network
rand_path = avg_path_matrix(1./W_rand);    %path of the random netowork
net_path = avg_path_matrix(1./W);          %path of the network

A = (net_path - rand_path);
if A < 0
    A = 0;
end
diff_path =  A/ (reg_path - rand_path);
if net_path == Inf || rand_path == Inf || reg_path == Inf
    diff_path = 1;
end
if diff_path > 1
    diff_path = 1;
end


%compute all clustering calculations for the network
reg_clus = avg_clus_matrix(W_reg,varargin{1});
rand_clus = avg_clus_matrix(W_rand,varargin{1});
net_clus = avg_clus_matrix(W,varargin{1});

B = (reg_clus - net_clus);
if B < 0
    B = 0;
end
    
diff_clus = B / (reg_clus - rand_clus);
if isnan(reg_clus) || isnan(rand_clus) || isnan(net_clus)
    diff_clus = 1;
end
if diff_clus > 1
    diff_clus = 1;
end

%calculate small world value, the root sum of the squares of
%diff path and diff clus
SWP = 1 - (sqrt(diff_clus^2 + diff_path^2)/sqrt(2));
delta_C=diff_clus;
delta_L=diff_path;
end
end


%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%the above code requires the following funcitons

function [Clus] = avg_clus_matrix(W, met)
%a function to compute the average clusteirng coefficient for a 
%input matrix M

%Inputs:
%   W     a matrix, weighted or unweighted
%   met   a string, to represent the method to be used for computing
%         the clustering coefficient
%         possible strings: 'O' (Onnela), 'Z' (Zhang), 'B' (Barrat),
%         'bin' (binary)
%         default if none is chosen is Onnela

%Outputs:
%   Clus  the average clustering coefficient

%written by Eric Bridgeford

n = length(W);
[C] = clustering_coef_matrix(W, met);
Clus = nanmean(C);
end


function [Len] = avg_path_matrix(M)

%a function to compute the average path length of a given matrix
%using the graphallshortestpaths built-in matlab function

%written by Eric Bridgeford

n = length(M);
M = sparse(M);
D = graphallshortestpaths(M);

%checks if a node is disconnected from the system, and replaces
%its value with 0
for i = 1:n
    for j = 1:n
        if isinf(D(i,j)) == 1
            D(i,j) = 0;
        end
    end
end

Len = mean(mean(D));
end


function [C] = clustering_coef_matrix(W, met)

%a modification of the clustering coefficient function provided
%in the brain connectivity toolbox

%improved definition of Onnela Clustering Coefficient, as well as
%implementation of function for Zhang and Barrat clustering values

%Reference: 
%   Onnela et al., Phys. Rev. E71, 065103(R)(2005)
%   B.Zhang and S. Horvath, Stat. App. Genet. Mol. Biol.4, 17(2005)
%   Barrat et al., Proc. Natl. Acad. Sci. U.S.A.101, 3747(2004)
%   Watts and Strogatz (1998) Nature 393:440-442

%Inputs:
%   W    the weighted or unweighted connectivity matrix
%   met   a string, to represent the method to be used for computing
%         the clustering coefficient
%         possible strings: 'O' (Onnela), 'Z' (Zhang), 'B' (Barrat), 'bin'
%         (binary)
%         default if none is chosen is Onnela


%code originally written by Mika Rubinov, UNSW, 2007-2010
%modified/written by Eric Bridgeford

if met == 'O'
    K=sum(W~=0,2);
    W = double(W);
    W2 = W/max(max(W));
    cyc3=diag(W2.^(1/3)^3);
    K(cyc3==0)=inf;             %if no 3-cycles exist, make C=0 (via K=inf)
    C=cyc3./(K.*(K-1));
end
if met == 'bin'
    G = double(W>0);
    n=length(G);
    C=zeros(n,1);
    for u=1:n
        V=find(G(u,:));
        k=length(V);
        if k>=2;                %degree must be at least 2
            S=G(V,V);
            C(u)=sum(S(:))/(k^2-k);
        end
    end
end

if met == 'Z'
    K=sum(W~=0,2);
    W = double(W);
    W2 = W/max(max((W)));
    cyc3=diag((W2)^3);
    denom = zeros(length(W),1);
    for i = 1:length(W)
        denom(i) = (sum(W2(i,:))^2-sum(W2(i,:).^2));
    end
    C = cyc3./denom;
end

if met == 'B'
    A = double(W>0);
    C = zeros(length(W),1);
    for i = 1:length(W)
        sum1 = 0;
        for j = 1:length(W)
            for k = 1:length(W)
                sum1 = ((W(i,j)+W(i,k))/2)*A(i,j)*A(j,k)*A(i,k)+sum1;
            end
        end
        C(i) = 1/(sum(W(i,:))*(sum(A(i,:))-1))*sum1;
    end
end


end


function A_rand=randomize_matrix(A);

%This code creates a random undirected network from the connectivity
%distribution of an undirected adjacency matrix, ie, the intital matrix
%must be symmetric.

% INPUTS:
%   A: an undirected adjacency matrix (symmetric) with no self connections

% OUTPUTS:
%   A_rand: a comparable random network with same number of nodes and
%       connectivity distribution

% written by Sarah F. Muldoon

num_nodes=length(A);
A_rand=zeros(num_nodes);
mask=triu(ones(num_nodes),1);
grab_indices=find(mask > 0);

orig_edges=A(grab_indices);
num_edges=length(orig_edges);

rand_index=randperm(num_edges);
randomized_edges=orig_edges(rand_index);

edge=1;
for i=1:num_nodes-1
    for j=i+1:num_nodes
        A_rand(i,j)=randomized_edges(edge);
        A_rand(j,i)=randomized_edges(edge);
        edge=edge+1;
    end
end
end

        
function M = regular_matrix_generator(G,r)
%generates a regular matrix, with weights obtained form the 
%original adjacency matrix representation of the network

% note that all inputs should be symmeterized prior to forming a regular
% matrix, since otherwise half of the connnections will be trashed. This
% can be accomplished with the built in symm_matrix function, however,
% the function is not performed here so that users can use their own 
% symmeterization procedure.

%Inputs:
%   G    the adjacency matrix for the given network; must be symmmeterized
%   r    the approximate radius of the regular network 

%Outputs:
%   M    the regular matrix for the given network, where all 
%        weights are sorted such that the inner radius has the
%        highest weights randomly distributed across the nodes, 
%        and so on

%written by Eric W. Bridgeford 

n = length(G);
G = triu(G);
%reshape the matrix G into an array, B
B = reshape(G,[length(G)^2,1]);
%sorts the array in descending order
B = sort(B,'descend');
%computes the number of connections and adds zeros if 
%numel(G) < 2*n*r
num_els =ceil(numel(G)/(2*n));
num_zeros = 2*n*num_els - numel(G);
%adds zeros to the remaineder of the list, so length(B) = 2*n*r
B = cat(1,B,zeros(num_zeros,1));
%reshapes B into a matrix, where the values descend top to
%bottom, as well as left to right. The greatest value in each 
%column is less than the smallest value of the column to its left.
B = reshape(B,[n],[]);

M = zeros(length(G));

%distributes the connections into a regular network, M, where
%the innermost radius represents the highest values
for i = 1:length(G)
    for z = 1:r
        a = randi([1,n]);

        %random integer chosen to take a value from B
        while (B(a,z) == 0 && z ~= r) || (B(a,z) == 0 && z == r && ~isempty(find(B(:,r),1)))
            a = randi([1,n]);
        end
        %finds the two nodes a distance of z from the origin node
        %and places the entries from the matrix B
        y_coor_1 = mod(i+z-1,length(G))+1;
        %y_coor_2 = mod(i-z-1,length(G))+1;
        M(i,y_coor_1) = B(a,z);
        M(y_coor_1,i) = B(a,z);
        %removes the weights from the matrix B so they cannot be
        %reused
        B(a,z) = 0;      
        
    end
end
end


function [W] = symm_matrix(A, bin_key)

% a function to symmetrize an input matrix. The procedure
% by which items are symmetrized such that:
%   in the binary case:
%       if a(i,j) || a(j,i) == 1 for i,j in A
%           w(i,j) && w(j,i) == 1
%   in  the weighted case:
%       if (a(i,j) || a(j,i) > 0 for i,j in A
%           w(i,j) && w(j,i) == (a(i,j) + a(j,i) )/ 2

% Inputs:
%   A:          The binary or weighted input matrix
%   bin_key:    the key to indicate whether weighted or binary analysis
%               will take place
%               1 indicates binarized, 0 indicates weighted

% Outputs
%   W:          The symmeterized matrix

% if binary analysis is specified, let binary symmeterization take place 

% written by Eric W. Bridgeford

W = zeros(length(A));

if bin_key == 1
    A = A > 0; % verify that the input matrix is binary
    for i = 1:length(A)
        for j = i:length(A)
            if A(i,j) || A(j,i)
                W(i,j) = 1;
                W(j,i) = 1;
            end
        end
    end
else
    for i = 1:length(A)
        for j = i:length(A)
            if A(i,j) || A(j,i)
                val = (A(i,j) + A(j,i)) / 2;
                W(i,j) = val;
                W(j,i) = val;
            end
        end
    end
end
end
    





