Home‎ > ‎

All help headers

adjacency_plot_und
 ADJACENCY_PLOT_UND     Quick visualization tool
 
    [X,Y,Z] = ADJACENCY_PLOT(AIJ,COOR) takes adjacency matrix AIJ and node
    spatial coordinates COOR and generates three vectors that can be used
    for quickly plotting the edges in AIJ. If no coordinates are specified,
    then each node is assigned a position on a circle. COOR can, in 
    general, be 2D or 3D.
 
    Example:
 
    >> load AIJ;                                % load your adjacency matrix
    >> load COOR;                               % load 3D coordinates for each node
    >> [x,y,z] = adjacency_plot_und(AIJ,COOR);  % call function
    >> plot3(x,y,z);                            % plots network as a single line object
 
    If COOR were 2D, the PLOT3 command changes to a PLOT command.
 
    NOTE: This function is similar to MATLAB's GPLOT command.
 
    Richard Betzel, Indiana University, 2013

agreement
 AGREEMENT      Agreement matrix from clusters
 
    D = AGREEMENT(CI) takes as input a set of vertex partitions CI of
    dimensions [vertex x partition]. Each column in CI contains the
    assignments of each vertex to a class/community/module. This function
    aggregates the partitions in CI into a square [vertex x vertex]
    agreement matrix D, whose elements indicate the number of times any two
    vertices were assigned to the same class.
 
    In the case that the number of nodes and partitions in CI is large
    (greater than ~1000 nodes or greater than ~1000 partitions), the script
    can be made faster by computing D in pieces. The optional input BUFFSZ
    determines the size of each piece. Trial and error has found that
    BUFFSZ ~ 150 works well.
 
    Inputs,     CI,     set of (possibly) degenerate partitions
                BUFFSZ, optional second argument to set buffer size
 
    Outputs:    D,      agreement matrix
 
    Richard Betzel, Indiana University, 2012

agreement_weighted
  AGREEMENT_WEIGHTED     Weights agreement matrix
 
    D = AGREEMENT_WEIGHTED(CI,WTS) is identical to AGREEMENT, with the 
    exception that each partitions contribution is weighted according to 
    the corresponding scalar value stored in the vector WTS. As an example,
    suppose CI contained partitions obtained using some heuristic for 
    maximizing modularity. A possible choice for WTS might be the Q metric
    (Newman's modularity score). Such a choice would add more weight to 
    higher modularity partitions.
 
    NOTE: Unlike AGREEMENT, this script does not have the input argument
    BUFFSZ.
 
    Inputs:     CI,     set of partitions
                WTS,    relative weight of importance of each paritition
 
    Outputs:    D,      weighted agreement matrix
 
    Richard Betzel, Indiana University, 2013

align_matrices
 ALIGN_MATRICES         Alignment of two matrices
 
    [Mreordered,Mindices,cost] = align_matrices(M1,M2,dfun,flag)
 
    This function aligns two matrices relative to one another by reordering
    the nodes in M2.  The function uses a version of simulated annealing.
 
    Inputs:     M1             = first connection matrix (square)
                M2             = second connection matrix (square)
                dfun           = distance metric to use for matching:
                                 'absdff' = absolute difference
                                 'sqrdff' = squared difference
                                 'cosang' = cosine of vector angle
 
                Mreordered     = reordered connection matrix M2
                Mindices       = reordered indices
                cost           = distance between M1 and Mreordered
 
    Connection matrices can be weighted or binary, directed or undirected.
    They must have the same number of nodes.  M1 can be entered in any
    node ordering.
 
    Note that in general, the outcome will depend on the initial condition
    (the setting of the random number seed).  Also, there is no good way to 
    determine optimal annealing parameters in advance - these parameters 
    will need to be adjusted "by hand" (particularly H, Texp, T0, and Hbrk).  
    For large and/or dense matrices, it is highly recommended to perform 
    exploratory runs varying the settings of 'H' and 'Texp' and then select 
    the best values.
 
    Based on extensive testing, it appears that T0 and Hbrk can remain
    unchanged in most cases.  Texp may be varied from 1-1/H to 1-10/H, for
    example.  H is the most important parameter - set to larger values as
    the problem size increases.  Good solutions can be obtained for
    matrices up to about 100 nodes.  It is advisable to run this function
    multiple times and select the solution(s) with the lowest 'cost'.
 
    If the two matrices are related it may be very helpful to pre-align them
    by reordering along their largest eigenvectors:
        [v,~] = eig(M1); v1 = abs(v(:,end)); [a1,b1] = sort(v1);
        [v,~] = eig(M2); v2 = abs(v(:,end)); [a2,b2] = sort(v2);
        [a,b,c] = overlapMAT2(M1(b1,b1),M2(b2,b2),'dfun',1);
 
    Setting 'Texp' to zero cancels annealing and uses a greedy algorithm
    instead.
 
    Yusuke Adachi, University of Tokyo, 2010
    Olaf Sporns, Indiana University, 2010

assortativity_bin
  ASSORTATIVITY_BIN      Assortativity coefficient
 
    r = assortativity(CIJ,flag);
 
    The assortativity coefficient is a correlation coefficient between the
    degrees of all nodes on two opposite ends of a link. A positive
    assortativity coefficient indicates that nodes tend to link to other
    nodes with the same or similar degree.
 
    Inputs:     CIJ,    binary directed/undirected connection matrix
                flag,   0, undirected graph: degree/degree correlation
                        1, directed graph: out-degree/in-degree correlation
                        2, directed graph: in-degree/out-degree correlation
                        3, directed graph: out-degree/out-degree correlation
                        4, directed graph: in-degree/in-degree correlation
 
    Outputs:    r,      assortativity coefficient
 
    Notes: The function accepts weighted networks, but all connection
    weights are ignored. The main diagonal should be empty. For flag 1
    the function computes the directed assortativity described in Rubinov
    and Sporns (2010) NeuroImage.
 
    Reference:  Newman (2002) Phys Rev Lett 89:208701
                Foster et al. (2010) PNAS 107:10815�10820
 
    Olaf Sporns, Indiana University, 2007/2008
    Vassilis Tsiaras, University of Crete, 2009
    Murray Shanahan, Imperial College London, 2012
    Mika Rubinov, University of Cambridge, 2012

assortativity_wei
  ASSORTATIVITY_WEI      Assortativity coefficient
 
    r = assortativity_wei(CIJ,flag);
 
    The assortativity coefficient is a correlation coefficient between the
    strengths (weighted degrees) of all nodes on two opposite ends of a link.
    A positive assortativity coefficient indicates that nodes tend to link to
    other nodes with the same or similar strength.
 
    Inputs:     CIJ,    weighted directed/undirected connection matrix
                flag,   0, undirected graph: strength/strength correlation
                        1, directed graph: out-strength/in-strength correlation
                        2, directed graph: in-strength/out-strength correlation
                        3, directed graph: out-strength/out-strength correlation
                        4, directed graph: in-strength/in-strength correlation
 
    Outputs:    r,      assortativity coefficient
 
    Notes: The main diagonal should be empty. For flag 1 the function computes 
    the directed assortativity described in Rubinov and Sporns (2010) NeuroImage.
 
    Reference:  Newman (2002) Phys Rev Lett 89:208701
                Foster et al. (2010) PNAS 107:10815-10820
 
    Olaf Sporns, Indiana University, 2007/2008
    Vassilis Tsiaras, University of Crete, 2009
    Murray Shanahan, Imperial College London, 2012
    Mika Rubinov, University of Cambridge, 2012

backbone_wu
 BACKBONE_WU        Backbone
 
    [CIJtree,CIJclus] = backbone_wu(CIJ,avgdeg)
 
    The network backbone contains the dominant connections in the network
    and may be used to aid network visualization. This function computes
    the backbone of a given weighted and undirected connection matrix CIJ, 
    using a minimum-spanning-tree based algorithm.
 
    input:      CIJ,    connection/adjacency matrix (weighted, undirected)
             avgdeg,    desired average degree of backbone
    output: 
            CIJtree,    connection matrix of the minimum spanning tree of CIJ
            CIJclus,    connection matrix of the minimum spanning tree plus
                        strongest connections up to an average degree 'avgdeg'
 
    NOTE: nodes with zero strength are discarded.
    NOTE: CIJclus will have a total average degree exactly equal to 
          (or very close to) 'avgdeg'.
    NOTE: 'avgdeg' backfill is handled slightly differently than in Hagmann
          et al 2008.
 
    Reference: Hidalgo et al. (2007) Science 317, 482.
               Hagmann et al. (2008) PLoS Biol
 
    Olaf Sporns, Indiana University, 2007/2008/2010/2012

betweenness_bin
 BETWEENNESS_BIN    Node betweenness centrality
 
    BC = betweenness_bin(A);
 
    Node betweenness centrality is the fraction of all shortest paths in 
    the network that contain a given node. Nodes with high values of 
    betweenness centrality participate in a large number of shortest paths.
 
    Input:      A,      binary (directed/undirected) connection matrix.
 
    Output:     BC,     node betweenness centrality vector.
 
    Note: Betweenness centrality may be normalised to the range [0,1] as
    BC/[(N-1)(N-2)], where N is the number of nodes in the network.
 
    Reference: Kintali (2008) arXiv:0809.1906v2 [cs.DS]
               (generalization to directed and disconnected graphs)
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2012

betweenness_wei
 BETWEENNESS_WEI    Node betweenness centrality
 
    BC = betweenness_wei(L);
 
    Node betweenness centrality is the fraction of all shortest paths in 
    the network that contain a given node. Nodes with high values of 
    betweenness centrality participate in a large number of shortest paths.
 
    Input:      L,      Directed/undirected connection-length matrix.
 
    Output:     BC,     node betweenness centrality vector.
 
    Notes:
        The input matrix must be a connection-length matrix, typically
    obtained via a mapping from weight to length. For instance, in a
    weighted correlation network higher correlations are more naturally
    interpreted as shorter distances and the input matrix should
    consequently be some inverse of the connectivity matrix. 
        Betweenness centrality may be normalised to the range [0,1] as
    BC/[(N-1)(N-2)], where N is the number of nodes in the network.
 
    Reference: Brandes (2001) J Math Sociol 25:163-177.
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2012

breadth
 BREADTH        Auxiliary function for breadthdist.m
 
    [distance,branch] = breadth(CIJ,source);
 
    Implementation of breadth-first search.
 
    Input:      CIJ,        binary (directed/undirected) connection matrix
                source,     source vertex
 
    Outputs:    distance,   distance between 'source' and i'th vertex
                            (0 for source vertex)
                branch,     vertex that precedes i in the breadth-first search tree
                            (-1 for source vertex)
         
    Notes: Breadth-first search tree does not contain all paths (or all 
    shortest paths), but allows the determination of at least one path with
    minimum distance. The entire graph is explored, starting from source 
    vertex 'source'.
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

breadthdist
 BREADTHDIST      Reachability and distance matrices
 
    [R,D] = breadthdist(CIJ);
 
    The binary reachability matrix describes reachability between all pairs
    of nodes. An entry (u,v)=1 means that there exists a path from node u
    to node v; alternatively (u,v)=0.
 
    The distance matrix contains lengths of shortest paths between all
    pairs of nodes. An entry (u,v) represents the length of shortest path 
    from node u to  node v. The average shortest path length is the 
    characteristic path length of the network.
 
    Input:      CIJ,     binary (directed/undirected) connection matrix
 
    Outputs:    R,       reachability matrix
                D,       distance matrix
 
    Note: slower but less memory intensive than "reachdist.m".
 
    Algorithm: Breadth-first search.
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

charpath
 CHARPATH       Characteristic path length, global efficiency and related statistics
 
    lambda                                  = charpath(D);
    lambda                                  = charpath(D);
    [lambda,efficiency]                     = charpath(D);
    [lambda,efficiency,ecc,radius,diameter] = charpath(D,diagonal_dist,infinite_dist);
 
    The network characteristic path length is the average shortest path
    length between all pairs of nodes in the network. The global efficiency
    is the average inverse shortest path length in the network. The nodal
    eccentricity is the maximal path length between a node and any other
    node in the network. The radius is the minimal eccentricity, and the
    diameter is the maximal eccentricity.
 
    Input:      D,              distance matrix
                diagonal_dist   optional argument
                                include distances on the main diagonal
                                    (default: diagonal_dist=0)
                infinite_dist   optional argument
                                include infinite distances in calculation
                                    (default: infinite_dist=1)
 
    Outputs:    lambda,         network characteristic path length
                efficiency,     network global efficiency
                ecc,            nodal eccentricity
                radius,         network radius
                diameter,       network diameter
 
    Notes:
        The input distance matrix may be obtained with any of the distance
    functions, e.g. distance_bin, distance_wei.
        Characteristic path length is defined here as the mean shortest
    path length between all pairs of nodes, for consistency with common
    usage. Note that characteristic path length is also defined as the
    median of the mean shortest path length from each node to all other
    nodes.
        Infinitely long paths (i.e. paths between disconnected nodes) are
    included in computations by default. This behavior may be modified with
    via the infinite_dist argument.
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008
    Mika Rubinov, U Cambridge, 2010/2015

clique_communities
  CLIQUE_COMMUNITIES     Overlapping community structure via clique percolation
 
    M = clique_communities(A, cq_thr)
 
    The optimal community structure is a subdivision of the network into
    groups of nodes which have a high number of within-group connections
    and a low number of between group connections.
 
    This algorithm uncovers overlapping community structure in binary
    undirected networks via the clique percolation method.
 
    Inputs:
        A,          Binary undirected connection matrix.
 
       	cq_thr,     Clique size threshold (integer). Larger clique size
                    thresholds potentially result in larger communities.
 
    Output:     
        M,          Overlapping community-affiliation matrix
                    Binary matrix of size CxN [communities x nodes]
 
    Algorithms:
        Bron–Kerbosch algorithm for detection of maximal cliques.
        Dulmage-Mendelsohn decomposition for detection of components
                    (implemented in get_components.m)
 
 
    Note: This algorithm can be slow and memory intensive in large
    matrices. The algorithm requires the function get_components.m
 
    Reference: Palla et al. (2005) Nature 435, 814-818.
 
    Mika Rubinov, Janelia HHMI, 2017

clustering_coef_bd
 CLUSTERING_COEF_BD     Clustering coefficient
 
    C = clustering_coef_bd(A);
 
    The clustering coefficient is the fraction of triangles around a node
    (equiv. the fraction of node's neighbors that are neighbors of each other).
 
    Input:      A,      binary directed connection matrix
 
    Output:     C,      clustering coefficient vector
 
    Reference: Fagiolo (2007) Phys Rev E 76:026107.
 
 
    Mika Rubinov, UNSW, 2007-2010

clustering_coef_bu
 CLUSTERING_COEF_BU     Clustering coefficient
 
    C = clustering_coef_bu(A);
 
    The clustering coefficient is the fraction of triangles around a node
    (equiv. the fraction of node's neighbors that are neighbors of each other).
 
    Input:      A,      binary undirected connection matrix
 
    Output:     C,      clustering coefficient vector
 
    Reference: Watts and Strogatz (1998) Nature 393:440-442.
 
 
    Mika Rubinov, UNSW, 2007-2010

clustering_coef_wd
 CLUSTERING_COEF_WD     Clustering coefficient
 
    C = clustering_coef_wd(W);
 
    The weighted clustering coefficient is the average "intensity"
    (geometric mean) of all triangles associated with each node.
 
    Input:      W,      weighted directed connection matrix
                        (all weights must be between 0 and 1)
 
    Output:     C,      clustering coefficient vector
 
    Reference: Fagiolo (2007) Phys Rev E 76:026107.
 
    Note:   All weights must be between 0 and 1.
            This may be achieved using the weight_conversion.m function,
            W_nrm = weight_conversion(W, 'normalize');
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

clustering_coef_wu
 CLUSTERING_COEF_WU     Clustering coefficient
 
    C = clustering_coef_wu(W);
 
    The weighted clustering coefficient is the average "intensity"
    (geometric mean) of all triangles associated with each node.
 
    Input:      W,      weighted undirected connection matrix
                        (all weights must be between 0 and 1)
 
    Output:     C,      clustering coefficient vector
 
    Note:   All weights must be between 0 and 1.
            This may be achieved using the weight_conversion.m function,
            W_nrm = weight_conversion(W, 'normalize');
 
    Reference: Onnela et al. (2005) Phys Rev E 71:065103
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

clustering_coef_wu_sign
 CLUSTERING_COEF_WU_SIGN     Multiple generalizations of the clustering coefficient 
 
    [C_pos,C_neg,Ctot_pos,Ctot_neg] = clustering_coef_wu_sign(W,coef_type);
 
    The weighted clustering coefficient is the average weight or intensity
    of all triangles associated with each node.
 
    Inputs:
        W,          
            Weighted undirected connection matrix
 
        corr_type,
            Desired type of clustering coefficient.
            Options:  
            1,  (default) Onnela et al. formula, used in original
                clustering_coef_wu.m. Computed separately for positive &
                negative weights.
            2,  Zhang & Horvath formula, similar to Onnela formula except
                denominator of Onnela formula relies on binarizing the
                network whereas this denominator is based on weight value,
                which reduces the sensitivity of this measure to the
                weights directly connected to the node of interest.
                Computed separately for positive & negative weights.
            3,  Constantini & Perugini's generalization of the Zhang &
                Horvath formula. This formula takes both positive &
                negative weights into account simultaneously, & is
                particularly sensitive to non-redundancy in path
                information based on sign (i.e., when two weights are
                positive & one negative, or all three are negative, both of
                which indicate that the weight of the third path is not
                redundant information). Produces only one value.
 
 
    Outputs: 
        C_pos/C_neg,
            Clustering coefficient vector for positive/negative weights.
            For the third option, only one vector is outputted (as C_pos). 
        Ctot_pos/Ctot_neg,
            Mean clustering coefficient for positive and negative weights.
 
    References: 
        Onnela et al. (2005) Phys Rev E 71:065103
        Zhang & Horvath (2005) Stat Appl Genet Mol Biol 41:1544-6115
        Costantini & Perugini (2014) PLOS ONE 9:e88669
 
 
    Contributor: Jeff Spielberg, Boston University, 2014-2015
                 (script based on clustering_coef_wu.m)

community_louvain
 COMMUNITY_LOUVAIN     Optimal community structure
 
    M     = community_louvain(W);
    [M,Q] = community_louvain(W,gamma);
    [M,Q] = community_louvain(W,gamma,M0);
    [M,Q] = community_louvain(W,gamma,M0,'potts');
    [M,Q] = community_louvain(W,gamma,M0,'negative_asym');
    [M,Q] = community_louvain(W,[],[],B);
 
    The optimal community structure is a subdivision of the network into
    nonoverlapping groups of nodes which maximizes the number of within-
 	group edges, and minimizes the number of between-group edges.
 
    This function is a fast and accurate multi-iterative generalization of
    the Louvain community detection algorithm. This function subsumes and
    improves upon,
 		modularity_louvain_und.m, modularity_finetune_und.m,
 		modularity_louvain_dir.m, modularity_finetune_dir.m,
        modularity_louvain_und_sign.m
 	and additionally allows to optimize other objective functions (includes
 	built-in Potts-model Hamiltonian, allows for custom objective-function
 	matrices).
 
    Inputs:
        W,
            directed/undirected weighted/binary connection matrix with
            positive and possibly negative weights.
        gamma,
            resolution parameter (optional)
                gamma>1,        detects smaller modules
                0<=gamma<1,     detects larger modules
                gamma=1,        classic modularity (default)
        M0,
            initial community affiliation vector (optional)
        B,
            objective-function type or custom objective matrix (optional)
            'modularity',       modularity (default)
            'potts',            Potts-model Hamiltonian (for binary networks)
            'negative_sym',     symmetric treatment of negative weights
            'negative_asym',    asymmetric treatment of negative weights
            B,                  custom objective-function matrix
 
            Note: see Rubinov and Sporns (2011) for a discussion of
            symmetric vs. asymmetric treatment of negative weights.
 
    Outputs:
        M,
            community affiliation vector
        Q,
            optimized community-structure statistic (modularity by default)
 
    Example:
        % Iterative community finetuning.
        % W is the input connection matrix.
        n  = size(W,1);             % number of nodes
        M  = 1:n;                   % initial community affiliations
        Q0 = -1; Q1 = 0;            % initialize modularity values
        while Q1-Q0>1e-5;           % while modularity increases
            Q0 = Q1;                % perform community detection
            [M, Q1] = community_louvain(W, [], M);
        end
 
    References:
        Blondel et al. (2008)  J. Stat. Mech. P10008.
        Reichardt and Bornholdt (2006) Phys. Rev. E 74, 016110.
        Ronhovde and Nussinov (2008) Phys. Rev. E 80, 016109
        Sun et al. (2008) Europhysics Lett 86, 28004.
        Rubinov and Sporns (2011) Neuroimage 56:2068-79.
 
    Mika Rubinov, U Cambridge 2015-2016

consensus_und
  CONSENSUS_UND      Consensus clustering
 
    CIU = CONSENSUS(D,TAU,REPS) seeks a consensus partition of the 
    agreement matrix D. The algorithm used here is almost identical to the
    one introduced in Lancichinetti & Fortunato (2012): The agreement
    matrix D is thresholded at a level TAU to remove an weak elements. The
    resulting matrix is then partitions REPS number of times using the
    Louvain algorithm (in principle, any clustering algorithm that can
    handle weighted matrixes is a suitable alternative to the Louvain
    algorithm and can be substituted in its place). This clustering
    produces a set of partitions from which a new agreement is built. If
    the partitions have not converged to a single representative partition,
    the above process repeats itself, starting with the newly built
    agreement matrix.
 
    NOTE: In this implementation, the elements of the agreement matrix must
    be converted into probabilities.
 
    NOTE: This implementation is slightly different from the original
    algorithm proposed by Lanchichinetti & Fortunato. In its original
    version, if the thresholding produces singleton communities, those
    nodes are reconnected to the network. Here, we leave any singleton
    communities disconnected.
 
    Inputs:     D,      agreement matrix with entries between 0 and 1
                        denoting the probability of finding node i in the
                        same cluster as node j
                TAU,    threshold which controls the resolution of the
                        reclustering
                REPS,   number of times that the clustering algorithm is
                        reapplied
 
    Outputs:    CIU,    consensus partition
 
    References: Lancichinetti & Fortunato (2012). Consensus clustering in
    complex networks. Scientific Reports 2, Article number: 336.
 
    Richard Betzel, Indiana University, 2012
 
    modified on 3/2014 to include "unique_partitions"

core_periphery_dir
 CORE_PERIPHERY_DIR     Core/periphery structure and core-ness statistic
 
    C     = core_periphery_dir(W)
    [C,q] = core_periphery_dir(W,gamm,C0)
 
    The optimal core/periphery subdivision is a partition of the network
    into two non-overlapping groups of nodes, a core group and a periphery
    group, in a way that maximizes the number/weight of within core-group
    edges, and minimizes the number/weight of within periphery-group edges.
 
    The core-ness is a statistic which quantifies the goodness of the
    optimal core/periphery subdivision.
 
    Input:      W       directed (weighted or binary) connection matrix.
                gamma,  core-ness resolution parameter (optional)
                        gamma>1     detects small core/large periphery
                        0<=gamma<1  detects large core/small periphery
                        default is gamma=1
 
    Outputs:    C,      binary vector of optimal core structure
                        C = 1 represents nodes in the core
                        C = 0 represents nodes in the periphery
                q,      maximized core-ness statistic
 
    Algorithm: A version of Kernighan-Lin algorithm for graph partitioning
    used in community detection (Newman, 2006) applied to optimize a
    core-structure objective described in Borgatti and Everett (2000).
 
    Reference: Borgatti and Everett (2000) Soc Networks 21:375–395.
               Newman (2006) Phys Rev E 74:036104, PNAS 23:8577-8582.
               Rubinov, Ypma et al. (2015) PNAS 112:10032-7
 
    2015, Mika Rubinov, U Cambridge

cycprob
  CYCPROB       Cycle probability
 
    [fcyc,pcyc] = cycprob(Pq);
 
    Cycles are paths which begin and end at the same node. Cycle 
    probability for path length d, is the fraction of all paths of length 
    d-1 that may be extended to form cycles of length d.
 
    Input:      Pq,     3D matrix, with Pq(i,j,q) = number of paths from 
                        'i' to 'j' of length 'q' (produced by 'findpaths')
 
    Outputs:    fcyc,   fraction of all paths that are cycles for each path
                        length 'q'. 
                pcyc,   probability that a non-cyclic path of length 'q-1' 
                        can be extended to form a cycle of length 'q', for 
                        each path length 'q', 
 
 
  Olaf Sporns, Indiana University, 2002/2007/2008

degrees_dir
 DEGREES_DIR        Indegree and outdegree
 
    [id,od,deg] = degrees_dir(CIJ);
 
    Node degree is the number of links connected to the node. The indegree 
    is the number of inward links and the outdegree is the number of 
    outward links.
 
    Input:      CIJ,    directed (binary/weighted) connection matrix
 
    Output:     id,     node indegree
                od,     node outdegree
                deg,    node degree (indegree + outdegree)
 
    Notes:  Inputs are assumed to be on the columns of the CIJ matrix.
            Weight information is discarded.
 
 
    Olaf Sporns, Indiana University, 2002/2006/2008

degrees_und
 DEGREES_UND        Degree
 
    deg = degrees_und(CIJ);
 
    Node degree is the number of links connected to the node.
 
    Input:      CIJ,    undirected (binary/weighted) connection matrix
 
    Output:     deg,    node degree
 
    Note: Weight information is discarded.
 
 
    Olaf Sporns, Indiana University, 2002/2006/2008

density_dir
  DENSITY_DIR        Density
 
    kden = density_dir(CIJ);
    [kden,N,K] = density_dir(CIJ);
 
    Density is the fraction of present connections to possible connections.
 
    Input:      CIJ,    directed (weighted/binary) connection matrix
 
    Output:     kden,   density
                N,      number of vertices
                K,      number of edges
 
    Notes:  Assumes CIJ is directed and has no self-connections.
            Weight information is discarded.
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

density_und
  DENSITY_UND        Density
 
    kden = density_und(CIJ);
    [kden,N,K] = density_und(CIJ);
 
    Density is the fraction of present connections to possible connections.
 
    Input:      CIJ,    undirected (weighted/binary) connection matrix
 
    Output:     kden,   density
                N,      number of vertices
                K,      number of edges
 
    Notes:  Assumes CIJ is undirected and has no self-connections.
            Weight information is discarded.
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

diffusion_efficiency
  DIFFUSION_EFFICIENCY      Global mean and pair-wise diffusion efficiency
 
    [GEdiff,Ediff] = diffusion_efficiency(adj);
 
    The diffusion efficiency between nodes i and j is the inverse of the
    mean first passage time from i to j, that is the expected number of
    steps it takes a random walker starting at node i to arrive for the
    first time at node j. Note that the mean first passage time is not a
    symmetric measure -- mfpt(i,j) may be different from mfpt(j,i) -- and
    the pair-wise diffusion efficiency matrix is hence also not symmetric.
 
 
    Input:
        adj,    Weighted/Unweighted, directed/undirected adjacency matrix
 
 
    Outputs:
        GEdiff, Mean Global diffusion efficiency (scalar)
        Ediff,  Pair-wise diffusion efficiency (matrix)
 
 
    References: Goñi J, et al (2013) PLoS ONE
 
    Joaquin Goñi and Andrea Avena-Koenigsberger, IU Bloomington, 2012

distance_bin
 DISTANCE_BIN       Distance matrix
 
    D = distance_bin(A);
 
    The distance matrix contains lengths of shortest paths between all
    pairs of nodes. An entry (u,v) represents the length of shortest path 
    from node u to node v. The average shortest path length is the 
    characteristic path length of the network.
 
    Input:      A,      binary directed/undirected connection matrix
 
    Output:     D,      distance matrix
 
    Notes: 
        Lengths between disconnected nodes are set to Inf.
        Lengths on the main diagonal are set to 0.
 
    Algorithm: Algebraic shortest paths.
 
 
    Mika Rubinov, U Cambridge
    Jonathan Clayden, UCL
    2007-2013

distance_wei
  DISTANCE_WEI       Distance matrix (Dijkstra's algorithm)
 
    D = distance_wei(L);
    [D,B] = distance_wei(L);
 
    The distance matrix contains lengths of shortest paths between all
    pairs of nodes. An entry (u,v) represents the length of shortest path 
    from node u to node v. The average shortest path length is the 
    characteristic path length of the network.
 
    Input:      L,      Directed/undirected connection-length matrix.
    *** NB: The length matrix L isn't the weights matrix W (see below) ***
 
    Output:     D,      distance (shortest weighted path) matrix
                B,      number of edges in shortest weighted path matrix
 
    Notes:
        The input matrix must be a connection-length matrix, typically
    obtained via a mapping from weight to length. For instance, in a
    weighted correlation network higher correlations are more naturally
    interpreted as shorter distances and the input matrix should
    consequently be some inverse of the connectivity matrix. 
        The number of edges in shortest weighted paths may in general 
    exceed the number of edges in shortest binary paths (i.e. shortest
    paths computed on the binarized connectivity matrix), because shortest 
    weighted paths have the minimal weighted distance, but not necessarily 
    the minimal number of edges.
        Lengths between disconnected nodes are set to Inf.
        Lengths on the main diagonal are set to 0.
 
    Algorithm: Dijkstra's algorithm.
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2012.
    Rick Betzel and Andrea Avena, IU, 2012

distance_wei_floyd
  DISTANCE_WEI_FLOYD        Distance matrix (Floyd-Warshall algorithm)
 
    [SPL,hops,Pmat] = distance_wei_floyd(D,transform)
 
    Computes the topological length of the shortest possible path
    connecting every pair of nodes in the network.
 
    Inputs:
 
        D,
            Weighted/unweighted directed/undirected 
            connection *weight* OR *length* matrix.
 
        transform,
            If the input matrix is a connection *weight* matrix, specify a
            transform that map input connection weights to connection
            lengths. Two transforms are available.
                'log' -> l_ij = -log(w_ij)
                'inv' -> l_ij =    1/w_ij
 
            If the input matrix is a connection *length* matrix, do not
            specify a transform (or specify an empty transform argument).
 
 
    Outputs:
 
        SPL,
            Unweighted/Weighted shortest path-length matrix.
            If W is directed matrix, then SPL is not symmetric.
 
        hops,
            Number of edges in the shortest path matrix. If W is
            unweighted, SPL and hops are identical.
 
        Pmat,
            Elements {i,j} of this matrix indicate the next node in the
            shortest path between i and j. This matrix is used as an input
            argument for function 'retrieve_shortest_path.m', which returns
            as output the sequence of nodes comprising the shortest path
            between a given pair of nodes.
 
 
    Notes:
 
        There may be more than one shortest path between any pair of nodes
        in the network. Non-unique shortest paths are termed shortest path
        degeneracies, and are most likely to occur in unweighted networks.
        When the shortest-path is degenerate, The elements of matrix Pmat
        correspond to the first shortest path discovered by the algorithm.
 
        The input matrix may be either a connection weight matrix, or a
        connection length matrix. The connection length matrix is typically
        obtained with a mapping from weight to length, such that higher
        weights are mapped to shorter lengths (see above).
 
 
    Algorithm:  Floyd–Warshall Algorithm
 
 
    Andrea Avena-Koenigsberger, IU, 2012

diversity_coef_sign
 DIVERSITY_COEF_SIGN     Shannon-entropy based diversity coefficient
 
    [Hpos Hneg] = diversity_coef_sign(W,Ci);
 
    The Shannon-entropy based diversity coefficient measures the diversity
    of intermodular connections of individual nodes and ranges from 0 to 1.
 
    Inputs:     W,      undirected connection matrix with positive and
                        negative weights
 
                Ci,     community affiliation vector
 
    Output:     Hpos,   diversity coefficient based on positive connections
                Hneg,   diversity coefficient based on negative connections
 
    References: Shannon CE (1948) Bell Syst Tech J 27, 379-423.
                Rubinov and Sporns (2011) NeuroImage.
 
 
    2011-2012, Mika Rubinov, U Cambridge

edge_betweenness_bin
 EDGE_BETWEENNESS_BIN    Edge betweenness centrality
 
    EBC = edge_betweenness_bin(A);
    [EBC BC] = edge_betweenness_bin(A);
 
    Edge betweenness centrality is the fraction of all shortest paths in 
    the network that contain a given edge. Edges with high values of 
    betweenness centrality participate in a large number of shortest paths.
 
    Input:      A,      binary (directed/undirected) connection matrix.
 
    Output:     EBC,    edge betweenness centrality matrix.
                BC,     node betweenness centrality vector.
 
    Note: Betweenness centrality may be normalised to the range [0,1] as
    BC/[(N-1)(N-2)], where N is the number of nodes in the network.
 
    Reference: Brandes (2001) J Math Sociol 25:163-177.
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2012

edge_betweenness_wei
 EDGE_BETWEENNESS_WEI    Edge betweenness centrality
 
    EBC = edge_betweenness_wei(L);
    [EBC BC] = edge_betweenness_wei(L);
 
    Edge betweenness centrality is the fraction of all shortest paths in 
    the network that contain a given edge. Edges with high values of 
    betweenness centrality participate in a large number of shortest paths.
 
    Input:      L,      Directed/undirected connection-length matrix.
 
    Output:     EBC,    edge betweenness centrality matrix.
                BC,     nodal betweenness centrality vector.
 
    Notes:
        The input matrix must be a connection-length matrix, typically
    obtained via a mapping from weight to length. For instance, in a
    weighted correlation network higher correlations are more naturally
    interpreted as shorter distances and the input matrix should
    consequently be some inverse of the connectivity matrix. 
        Betweenness centrality may be normalised to the range [0,1] as
    BC/[(N-1)(N-2)], where N is the number of nodes in the network.
 
    Reference: Brandes (2001) J Math Sociol 25:163-177.
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2012

edge_nei_overlap_bd
  EDGE_NEI_OVERLAP_BD     Overlap amongst neighbors of two adjacent nodes
 
    [EC,ec,degij] = edge_nei_bd(CIJ);
 
    This function determines the neighbors of two nodes that are linked by 
    an edge, and then computes their overlap.  Connection matrix must be
    binary and directed.  Entries of 'EC' that are 'inf' indicate that no
    edge is present.  Entries of 'EC' that are 0 denote "local bridges",
    i.e. edges that link completely non-overlapping neighborhoods.  Low
    values of EC indicate edges that are "weak ties".
 
    If CIJ is weighted, the weights are ignored. Neighbors of a node can be
    linked by incoming, outgoing, or reciprocal connections.
 
    Inputs:     CIJ,      directed (binary/weighted) connection matrix
   
    Outputs:    EC,     edge neighborhood overlap matrix
                ec,     edge neighborhood overlap per edge, in vector format
                degij,  degrees of node pairs connected by each edge
 
    Reference:
 
        Easley and Kleinberg (2010) Networks, Crowds, and Markets. 
            Cambridge University Press, Chapter 3
 
    Olaf Sporns, Indiana University, 2012

edge_nei_overlap_bu
  EDGE_NEI_OVERLAP_BU        Overlap amongst neighbors of two adjacent nodes
 
    [EC,ec,degij] = edge_nei_bu(CIJ);
 
    This function determines the neighbors of two nodes that are linked by 
    an edge, and then computes their overlap.  Connection matrix must be
    binary and directed.  Entries of 'EC' that are 'inf' indicate that no
    edge is present.  Entries of 'EC' that are 0 denote "local bridges", i.e.
    edges that link completely non-overlapping neighborhoods.  Low values
    of EC indicate edges that are "weak ties".
 
    If CIJ is weighted, the weights are ignored.
 
    Inputs:     CIJ,    undirected (binary/weighted) connection matrix
   
    Outputs:    EC,     edge neighborhood overlap matrix
                ec,     edge neighborhood overlap per edge, in vector format
                degij,  degrees of node pairs connected by each edge
 
    Reference: Easley and Kleinberg (2010) Networks, Crowds, and Markets. 
               Cambridge University Press, Chapter 3.
 
    Olaf Sporns, Indiana University, 2012

efficiency_bin
 EFFICIENCY_BIN     Global efficiency, local efficiency.
 
    Eglob = efficiency_bin(A);
    Eloc = efficiency_bin(A,1);
 
    The global efficiency is the average of inverse shortest path length,
    and is inversely related to the characteristic path length.
 
    The local efficiency is the global efficiency computed on the
    neighborhood of the node, and is related to the clustering coefficient.
 
    Inputs:     A,              binary undirected or directed connection matrix
                local,          optional argument
                                    local=0 computes global efficiency (default)
                                    local=1 computes local efficiency
 
    Output:     Eglob,          global efficiency (scalar)
                Eloc,           local efficiency (vector)
 
 
    Algorithm: algebraic path count
 
    Reference: Latora and Marchiori (2001) Phys Rev Lett 87:198701.
               Fagiolo (2007) Phys Rev E 76:026107.
               Rubinov M, Sporns O (2010) NeuroImage 52:1059-69
 
 
    Mika Rubinov, U Cambridge
    Jonathan Clayden, UCL
    2008-2013

efficiency_wei
 EFFICIENCY_WEI     Global efficiency, local efficiency.
 
    Eglob = efficiency_wei(W);
    Eloc = efficiency_wei(W,2);
 
    The global efficiency is the average of inverse shortest path length,
    and is inversely related to the characteristic path length.
 
    The local efficiency is the global efficiency computed on the
    neighborhood of the node, and is related to the clustering coefficient.
 
    Inputs:     W,
                    weighted undirected or directed connection matrix
 
                local,
                    optional argument
                    local=0  computes the global efficiency (default).
                   	local=1  computes the original version of the local
                                efficiency.
                	local=2  computes the modified version of the local
                                efficiency (recommended, see below). 
 
    Output:     Eglob,
                    global efficiency (scalar)
                Eloc,
                    local efficiency (vector)
 
    Notes:
        The  efficiency is computed using an auxiliary connection-length
    matrix L, defined as L_ij = 1/W_ij for all nonzero L_ij; This has an
    intuitive interpretation, as higher connection weights intuitively
    correspond to shorter lengths.
        The weighted local efficiency broadly parallels the weighted
    clustering coefficient of Onnela et al. (2005) and distinguishes the
    influence of different paths based on connection weights of the
    corresponding neighbors to the node in question. In other words, a path
    between two neighbors with strong connections to the node in question
    contributes more to the local efficiency than a path between two weakly
    connected neighbors. Note that the original weighted variant of the
    local efficiency (described in Rubinov and Sporns, 2010) is not a
    true generalization of the binary variant, while the modified variant
    (described in Wang et al., 2016) is a true generalization.
        For ease of interpretation of the local efficiency it may be
    advantageous to rescale all weights to lie between 0 and 1.
 
    Algorithm:  Dijkstra's algorithm
 
    References: Latora and Marchiori (2001) Phys Rev Lett 87:198701.
                Onnela et al. (2005) Phys Rev E 71:065103
                Fagiolo (2007) Phys Rev E 76:026107.
                Rubinov M, Sporns O (2010) NeuroImage 52:1059-69
                Wang Y et al. (2016) Neural Comput 21:1-19.
 
    Mika Rubinov, U Cambridge/Janelia HHMI, 2011-2017

eigenvector_centrality_und
 EIGENVECTOR_CENTRALITY_UND      Spectral measure of centrality
 
    v = eigenvector_centrality_und(CIJ)
 
    Eigenector centrality is a self-referential measure of centrality:
    nodes have high eigenvector centrality if they connect to other nodes
    that have high eigenvector centrality. The eigenvector centrality of
    node i is equivalent to the ith element in the eigenvector 
    corresponding to the largest eigenvalue of the adjacency matrix.
 
    Inputs:     CIJ,        binary/weighted undirected adjacency matrix.
 
    Outputs:      v,        eigenvector associated with the largest
                            eigenvalue of the adjacency matrix CIJ.
 
    Reference: Newman, MEJ (2002). The mathematics of networks.
 
    Contributors:
    Xi-Nian Zuo, Chinese Academy of Sciences, 2010
    Rick Betzel, Indiana University, 2012
    Mika Rubinov, University of Cambridge, 2015

erange
 ERANGE     Shortcuts
 
    [Erange,eta,Eshort,fs] = erange(CIJ);
 
    Shorcuts are central edges which significantly reduce the
    characteristic path length in the network.
 
    Input:      CIJ,        binary directed connection matrix
 
    Outputs:    Erange,     range for each edge, i.e. the length of the
                            shortest path from i to j for edge c(i,j) AFTER
                            the edge has been removed from the graph.
                eta         average range for entire graph.
                Eshort      entries are ones for shortcut edges.
                fs          fraction of shortcuts in the graph.
 
    Follows the treatment of 'shortcuts' by Duncan Watts
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

evaluate_generative_model
  EVALUATE_GENERATIVE_MODEL     Generation and evaluation of synthetic networks
 
    [B,E,K] = EVALUATE_GENERATIVE_MODEL(A,Atgt,D,m,modeltype,modelvar,params) 
 
    Generates synthetic networks and evaluates their energy function (see
    below) using the models described in the study by Betzel et al (2016)
    in Neuroimage.
 
    Inputs:
            A,          binary network of seed connections
            Atgt,       binary network against which synthetic networks are
                        compared
            D,          Euclidean distance/fiber length matrix
            m,          number of connections that should be present in
                        final synthetic network
            modeltype,  specifies the generative rule (see below)
            modelvar,   specifies whether the generative rules are based on
                        power-law or exponential relationship
                        ({'powerlaw'}|{'exponential})
            params,     either a vector (in the case of the geometric
                        model) or a matrix (for all other models) of
                        parameters at which the model should be evaluated.
 
    Outputs:
            B,          m x number of networks matrix of connections
            E,          energy for each synthetic network
            K,          Kolmogorov-Smirnov statistics for each synthetic
                        network.
 
    Full list of model types:
    (each model type realizes a different generative rule)
 
        1.  'sptl'          spatial model
        2.  'neighbors'     number of common neighbors
        3.  'matching'      matching index
        4.  'clu-avg'       average clustering coeff.
        5.  'clu-min'       minimum clustering coeff.
        6.  'clu-max'       maximum clustering coeff.
        7.  'clu-diff'      difference in clustering coeff.
        8.  'clu-prod'      product of clustering coeff.
        9.  'deg-avg'       average degree
        10. 'deg-min'       minimum degree
        11. 'deg-max'       maximum degree
        12. 'deg-diff'      difference in degree
        13. 'deg-prod'      product of degree
 
    Note: Energy is calculated in exactly the same way as in Betzel et
    al (2016). There are four components to the energy are KS statistics
    comparing degree, clustering coefficient, betweenness centrality, and 
    edge length distributions. Energy is calculated as the maximum across
    all four statistics.
 
    Reference: Betzel et al (2016) Neuroimage 124:1054-64.
 
    Richard Betzel, Indiana University/University of Pennsylvania, 2015

find_motif34
 FIND_MOTIF34       Motif legend
 
    Motif_matrices = find_motif34(Motif_id,Motif_class);
    Motif_id = find_motif34(Motif_matrix);
 
    This function returns all motif isomorphs for a given motif id and 
    class (3 or 4). The function also returns the motif id for a given
    motif matrix
 
    1. Input:       Motif_id,           e.g. 1 to 13, if class is 3
                    Motif_class,        number of nodes, 3 or 4.
 
       Output:      Motif_matrices,     all isomorphs for the given motif
 
    2. Input:       Motif_matrix        e.g. [0 1 0; 0 0 1; 1 0 0]
 
       Output       Motif_id            e.g. 1 to 13, if class is 3
 
 
 Mika Rubinov, UNSW, 2007-2008

findpaths
 FINDPATHS      Network paths
 
    [Pq,tpath,plq,qstop,allpths,util] = findpaths(CIJ,sources,qmax,savepths);
 
    Paths are sequences of linked nodes, that never visit a single node
    more than once. This function finds all paths that start at a set of 
    source nodes, up to a specified length. Warning: very memory-intensive.
 
    Inputs:     CIJ,        binary (directed/undirected) connection matrix
                qmax,       maximal path length
                sources,    source units from which paths are grown
                savepths,   set to 1 if all paths are to be collected in
                            'allpths'
 
    Outputs:    Pq,         3D matrix, with P(i,j,q) = number of paths from
                            'i' to 'j' of length 'q'.
                tpath,      total number of paths found (lengths 1 to 'qmax')
                plq,        path length distribution as a function of 'q'
                qstop,      path length at which 'findpaths' is stopped
                allpths,    a matrix containing all paths up to 'qmax'
                util,       node use index
 
    Note that Pq(:,:,N) can only carry entries on the diagonal, as all
    "legal" paths of length N-1 must terminate.  Cycles of length N are
    possible, with all vertices visited exactly once (except for source and
    target). 'qmax = N' can wreak havoc (due to memory problems).
 
    Note: Weights are discarded.
    Note: I am certain that this algorithm is rather inefficient -
    suggestions for improvements are welcome.
 
    Olaf Sporns, Indiana University, 2002/2007/2008/2010

findwalks
 FINDWALKS      Network walks
 
    [Wq,twalk,wlq] = findwalks(CIJ);
 
    Walks are sequences of linked nodes, that may visit a single node more
    than once. This function finds the number of walks of a given length, 
    between any two nodes.
 
    Input:      CIJ         binary (directed/undirected) connection matrix
 
    Outputs:    Wq          3D matrix, Wq(i,j,q) is the number of walks
                            from 'i' to 'j' of length 'q'.
                twalk       total number of walks found
                wlq         walk length distribution as function of 'q'
 
    Notes: Wq grows very quickly for larger N,K,q. Weights are discarded.
 
    Algorithm: algebraic path count
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

flow_coef_bd
 FLOW_COEF_BD        Node-wise flow coefficients
 
    [hc,HC,total_flo] = flow_coef_bd(CIJ)
 
    Computes the flow coefficient for each node and averaged over the
    network, as described in Honey et al. (2007) PNAS. The flow coefficient
    is similar to betweenness centrality, but works on a local
    neighborhood. It is mathematically related to the clustering
    coefficient  (cc) at each node as, fc+cc <= 1.
 
    input:      CIJ,	connection/adjacency matrix (binary, directed)
    output:     fc,     flow coefficient for each node
                FC,     average flow coefficient over the network
         total_flo,     number of paths that "flow" across the central node
 
    Reference:  Honey et al. (2007) Proc Natl Acad Sci U S A
 
    Olaf Sporns, Indiana University, 2007/2010/2012

gateway_coef_sign
 GATEWAY_COEF_SIGN     Gateway coefficient
 
    [Gpos,Gneg] = gateway_coef_sign(W,Ci,centtype);
 
    Gateway coefficient is a variant of participation coefficient. Similar 
    to participation coefficient, gateway coefficient measures the 
    diversity of intermodular connections of individual nodes, but this is
    weighted by how critical these connections are to intermodular
    connectivity (e.g., if a node is the only connection between it's
    module and another module, it will have a higher gateway coefficient). 
 
    Inputs:     W,        undirected connection matrix with positive and
                          negative weights
 
                Ci,       community affiliation vector
 
                centtype, centrality measure to use
                          1 = Node Strength
                          2 = Betweenness Centrality
 
    Output:     Gpos,     gateway coefficient for positive weights
 
                Gneg,     gateway coefficient for negative weights
 
    Reference: Vargas ER, Wahl LM. Eur Phys J B (2014) 87:1-10.
 
    Jeff Spielberg, Boston University

generate_fc
    GENERATE_FC     Generation of synthetic functional connectivity matrices
 
    [FCpre,pred_data,Fcorr] = generate_fc(SC,beta,ED,{'SPLwei_log','SIwei_log'},FC)
    [FCpre,pred_data] = generate_fc(SC,beta,[],{'SPLwei_log','SIwei_log'})
 
    Uses a vector beta of regression coefficients from the model
    FC = pred_var*beta to predict FC. pred_var are structural-based network
    measures derived from the structural connectivity network.
 
    Inputs:
 
        SC,
            Weighted/unweighted undirected NxN Structural Connectivity matrix.
 
        beta,
            Regression coefficients (vector). These may be obtained as an
            output parameter from function predict_fc.m
 
        ED,
            Euclidean distance matrix or upper triangular vector of the
            matrix (optional)
 
        pred_var,
            Set of M predictors. These can be given as an KxM array where
            K = ((N*(N-1))/2) and M is the number of predictors.
            Alternatively, pred_var can be a cell with the names of network
            measures to be used as predictors. Accepted network measure
            names are:
                SPLbin        - Shortest-path length (binary)
                SPLwei_inv    - Shortest-path length computed with an inv transform
              	SPLwei_log    - Shortest-path length computed with a log transform
              	SPLdist       - Shortest-path length computed with no transform
               	SIbin         - Search Information of binary shortest-paths
              	SIwei_inv     - Search Information of shortest-paths computed with an inv transform
            	SIwei_log     - Search Information of shortest-paths computed with a log transform
               	SIdist        - Search Information of shortest-paths computed with no transform
               	T             - Path Transitivity
               	deltaMFPT     - Column-wise z-scored mean first passage time
                neighOverlap  - Neighborhood Overlap
               	MI            - Matching Index
 
            Predictors must be specified in the order that matches the
            given beta values.
 
       	model,
            Specifies the order of the regression model used within
            matlab's function regstats.m. 'model' can be any option
            accepted by matlab's regstats.m function (e.g.'linear',
            'interaction' 'quadratic', etc). If no model is specified,
            'linear' is the default.
 
     	FC,
            Functional connections. FC can be a NxN symmetric matrix or a
            ((N*(N-1))/2) x 1 vector containing the upper triangular
            elements of the square FC matrix (excluding diagonal elements).
            This argument is optional and only used to compute the
            correlation between the predicted FC and empirical FC.
 
 
    Outputs:
 
        FCpre,
            Predicted NxN Functional Connectivity matrix
 
       	pred_data,
            KxM array of predictors.
 
        FCcorr,
            Pearson Correlation between FCpred and FC
 
 
    Reference: Goñi et al. (2014) PNAS 111: 833–838
 
 
    Andrea Avena-Koenigsberger, Joaquin Goñi and Olaf Sporns; IU Bloomington, 2016

generative_model
  GENERATIVE_MODEL          Run generative model code
 
    B = GENERATIVE_MODEL(A,D,m,modeltype,modelvar,params)
 
    Generates synthetic networks using the models described in the study by
    Betzel et al (2016) in Neuroimage.
 
    Inputs:
            A,          binary network of seed connections
            D,          Euclidean distance/fiber length matrix
            m,          number of connections that should be present in
                        final synthetic network
            modeltype,  specifies the generative rule (see below)
            modelvar,   specifies whether the generative rules are based on
                        power-law or exponential relationship
                        ({'powerlaw'}|{'exponential})
            params,     either a vector (in the case of the geometric
                        model) or a matrix (for all other models) of
                        parameters at which the model should be evaluated.
            epsilon,    the baseline probability of forming a particular
                        connection (should be a very small number
                        {default = 1e-5}).
 
    Output:
            B,          m x number of networks matrix of connections
 
 
    Full list of model types:
    (each model type realizes a different generative rule)
 
        1.  'sptl'          spatial model
        2.  'neighbors'     number of common neighbors
        3.  'matching'      matching index
        4.  'clu-avg'       average clustering coeff.
        5.  'clu-min'       minimum clustering coeff.
        6.  'clu-max'       maximum clustering coeff.
        7.  'clu-diff'      difference in clustering coeff.
        8.  'clu-prod'      product of clustering coeff.
        9.  'deg-avg'       average degree
        10. 'deg-min'       minimum degree
        11. 'deg-max'       maximum degree
        12. 'deg-diff'      difference in degree
        13. 'deg-prod'      product of degree
 
 
    Example usage:
 
        load demo_generative_models_data
 
        % get number of bi-directional connections
        m = nnz(A)/2;
  
        % get cardinality of network
        n = length(A);
  
        % set model type
        modeltype = 'neighbors';
  
        % set whether the model is based on powerlaw or exponentials
        modelvar = [{'powerlaw'},{'powerlaw'}];
  
        % choose some model parameters
        params = [-2,0.2; -5,1.2; -1,1.5];
        nparams = size(params,1);
  
        % generate synthetic networks
        B = generative_model(Aseed,D,m,modeltype,modelvar,params);
 
        % store them in adjacency matrix format
        Asynth = zeros(n,n,nparams);
        for i = 1:nparams; 
            a = zeros(n); a(B(:,i)) = 1; a = a + a'; 
            Asynth(:,:,i) = a; 
        end
 
    Reference: Betzel et al (2016) Neuroimage 124:1054-64.
 
    Richard Betzel, Indiana University/University of Pennsylvania, 2015

get_components
  GET_COMPONENTS        Connected components
 
    [comps,comp_sizes] = get_components(adj);
 
    Returns the components of an undirected graph specified by the binary and 
    undirected adjacency matrix adj. Components and their constitutent nodes are 
    assigned the same index and stored in the vector, comps. The vector, comp_sizes,
    contains the number of nodes beloning to each component.
 
    Inputs:         adj,    binary and undirected adjacency matrix
 
    Outputs:      comps,    vector of component assignments for each node
             comp_sizes,    vector of component sizes
 
    Note: disconnected nodes will appear as components of size 1
 
    J Goni, University of Navarra and Indiana University, 2009/2011

grid_communities
  GRID_COMMUNITIES       Outline communities along diagonal
 
    [X Y INDSORT] = GRID_COMMUNITIES(C) takes a vector of community
    assignments C and returns three output arguments for visualizing the
    communities. The third is INDSORT, which is an ordering of the vertices
    so that nodes with the same community assignment are next to one
    another. The first two arguments are vectors that, when overlaid on the
    adjacency matrix using the PLOT function, highlight the communities.
 
    Example:
 
    >> load AIJ;                                % load adjacency matrix
    >> [C,Q] = modularity_louvain_und(AIJ);     % get community assignments
    >> [X,Y,INDSORT] = fcn_grid_communities(C); % call function
    >> imagesc(AIJ(INDSORT,INDSORT));           % plot ordered adjacency matrix
    >> hold on;                                 % hold on to overlay community visualization
    >> plot(X,Y,'r','linewidth',2);             % plot community boundaries
 
    Inputs:     C,       community assignments
 
    Outputs:    X,       x coor
                Y,       y coor
                INDSORT, indices
 
    Richard Betzel, Indiana University, 2012
 

gtom
 GTOM       Generalized topological overlap measure
 
    gt = gtom(adj,numSteps);
 
    The m-th step generalized topological overlap measure (GTOM) quantifies
    the extent to which a pair of nodes have similar m-th step neighbors.
    Mth-step neighbors are nodes that are reachable by a path of at most
    length m.
 
    This function computes the the M x M generalized topological overlap
    measure (GTOM) matrix for number of steps, numSteps. 
 
    Inputs:       adj,    adjacency matrix (binary,undirected)
             numSteps,    number of steps
 
    Outputs:       gt,    GTOM matrix
 
    NOTE: When numSteps is equal to 1, GTOM is identical to the topological
    overlap measure (TOM) from reference [2]. In that case the 'gt' matrix
    records, for each pair of nodes, the fraction of neighbors the two
    nodes share in common, where "neighbors" are one step removed. As
    'numSteps' is increased, neighbors that are furter out are considered.
    Elements of 'gt' are bounded between 0 and 1.  The 'gt' matrix can be
    converted from a similarity to a distance matrix by taking 1-gt.
 
    References: [1] Yip & Horvath (2007) BMC Bioinformatics 2007, 8:22
                [2] Ravasz et al (2002) Science 297 (5586), 1551.
 
    J Goni, University of Navarra and Indiana University, 2009/2011

jdegree
 JDEGREE        Joint degree distribution
 
    [J,J_od,J_id,J_bl] = jdegree(CIJ);
 
    This function returns a matrix in which the value of each element (u,v)
    corresponds to the number of nodes that have u outgoing connections 
    and v incoming connections.
 
    Input:      CIJ,    directed (weighted/binary) connection matrix
 
    Outputs:    J,      joint degree distribution matrix (shifted by one)
                J_od,   number of vertices with od>id.
                J_id,   number of vertices with id>od.
                J_bl,   number of vertices with id=od.
 
    Note: Weights are discarded.
 
 
    Olaf Sporns, Indiana University, 2002/2006/2008

kcore_bd
 KCORE_BD       K-core
 
    [CIJkcore,kn,peelorder,peellevel] = kcore_bd(CIJ,k);
 
    The k-core is the largest subnetwork comprising nodes of degree at
    least k. This function computes the k-core for a given binary directed
    connection matrix by recursively peeling off nodes with degree lower
    than k, until no such nodes remain.
 
    input:          CIJ,        connection/adjacency matrix (binary, directed)
                      k,        level of k-core
 
    output:    CIJkcore,        connection matrix of the k-core.  This matrix
                                only contains nodes of degree at least k.
                     kn,        size of k-core
                     peelorder, indices in the order in which they were
                                peeled away during k-core decomposition
                     peellevel, corresponding level - nodes at the same
                                level have been peeled away at the same time
 
    'peelorder' and 'peellevel' are similar the the k-core sub-shells
    described in Modha and Singh (2010).
 
    References: e.g. Hagmann et al. (2008) PLoS Biology
 
    Olaf Sporns, Indiana University, 2007/2008/2010/2012

kcore_bu
 KCORE_BU       K-core
 
    [CIJkcore,kn,peelorder,peellevel] = kcore_bu(CIJ,k);
 
    The k-core is the largest subnetwork comprising nodes of degree at
    least k. This function computes the k-core for a given binary
    undirected connection matrix by recursively peeling off nodes with
    degree lower than k, until no such nodes remain.
 
    input:          CIJ,        connection/adjacency matrix (binary, undirected)
                      k,        level of k-core
 
    output:    CIJkcore,        connection matrix of the k-core.  This matrix
                                only contains nodes of degree at least k.
                     kn,        size of k-core
                     peelorder, indices in the order in which they were
                                peeled away during k-core decomposition
                     peellevel, corresponding level - nodes at the same
                                level were peeled away at the same time
 
    'peelorder' and 'peellevel' are similar the the k-core sub-shells
    described in Modha and Singh (2010).
 
    References: e.g. Hagmann et al. (2008) PLoS Biology
 
    Olaf Sporns, Indiana University, 2007/2008/2010/2012

kcoreness_centrality_bd
 KCORENESS_CENTRALITY_BD       K-coreness centrality
 
    [coreness,kn] = kcoreness_centrality_bd(CIJ)
 
    The k-core is the largest subgraph comprising nodes of degree at least
    k. The coreness of a node is k if the node belongs to the k-core but
    not to the (k+1)-core. This function computes k-coreness of all nodes
    for a given binary directed connection matrix.
 
    input:          CIJ,        connection/adjacency matrix (binary, directed)
 
    output:    coreness,        node coreness.
                     kn,        size of k-core
 
    References: e.g. Hagmann et al. (2008) PLoS Biology
 
    Olaf Sporns, Indiana University, 2007/2008/2010/2012

kcoreness_centrality_bu
 KCORENESS_CENTRALITY_BU       K-coreness centrality
 
    [coreness,kn] = kcoreness_centrality_bu(CIJ)
 
    The k-core is the largest subgraph comprising nodes of degree at least
    k. The coreness of a node is k if the node belongs to the k-core but
    not to the (k+1)-core. This function computes the coreness of all nodes
    for a given binary undirected connection matrix.
 
    input:          CIJ,        connection/adjacency matrix (binary, undirected)
 
    output:    coreness,        node coreness.
                     kn,        size of k-core
 
    References: e.g. Hagmann et al. (2008) PLoS Biology
 
    Olaf Sporns, Indiana University, 2007/2008/2010/2012

latmio_dir
 LATMIO_DIR     Lattice with preserved in/out degree distribution
 
    [Rlatt,Rrp,ind_rp,eff] = latmio_dir(R,ITER,D);
 
    This function "latticizes" a directed network, while preserving the in-
    and out-degree distributions. In weighted networks, the function
    preserves the out-strength but not the in-strength distributions.
 
    Input:      R,      directed (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
                D,      distance-to-diagonal matrix
 
    Output:     Rlatt,  latticized network in original node ordering
                Rrp,    latticized network in node ordering used for
                        latticization
                ind_rp, node ordering used for latticization
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
                Sporns and Zwi (2004) Neuroinformatics 2:145
 
    Mika Rubinov, UNSW, 2007-2010
    Olaf Sporns, IU, 2012

latmio_dir_connected
 LATMIO_DIR_CONNECTED     Lattice with preserved in/out degree distribution
 
    [Rlatt,Rrp,ind_rp,eff] = latmio_dir_connected(R,ITER,D);
 
    This function "latticizes" a directed network, while preserving the in-
    and out-degree distributions. In weighted networks, the function
    preserves the out-strength but not the in-strength distributions. The 
    function also ensures that the randomized network maintains
    connectedness, the ability for every node to reach every other node in
    the network. The input network for this function must be connected.
 
    Input:      R,      directed (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
                D,      distance-to-diagonal matrix
 
    Output:     Rlatt,  latticized network in original node ordering
                Rrp,    latticized network in node ordering used for
                        latticization
                ind_rp, node ordering used for latticization
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
                Sporns and Zwi (2004) Neuroinformatics 2:145
 
    Mika Rubinov, UNSW, 2007-2010
    Olaf Sporns, Indiana University, 2012

latmio_und
 LATMIO_UND     Lattice with preserved degree distribution
 
    [Rlatt,Rrp,ind_rp,eff] = latmio_und(R,ITER,D);
 
    This function "latticizes" an undirected network, while preserving the 
    degree distribution. The function does not preserve the strength 
    distribution in weighted networks.
 
    Input:      R,      undirected (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
                D,      distance-to-diagonal matrix
 
    Output:     Rlatt,  latticized network in original node ordering
                Rrp,    latticized network in node ordering used for
                        latticization
                ind_rp, node ordering used for latticization
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
                Sporns and Zwi (2004) Neuroinformatics 2:145
 
    2007-2012
    Mika Rubinov, UNSW
    Jonathan Power, WUSTL
    Olaf Sporns, IU

latmio_und_connected
 LATMIO_UND_CONNECTED     Lattice with preserved degree distribution
 
    [Rlatt,Rrp,ind_rp,eff] = latmio_und_connected(R,ITER,D);
 
    This function "latticizes" an undirected network, while preserving the 
    degree distribution. The function does not preserve the strength 
    distribution in weighted networks. The function also ensures that the 
    randomized network maintains connectedness, the ability for every node 
    to reach every other node in the network. The input network for this 
    function must be connected.
 
    Input:      R,      undirected (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
                D,      distance-to-diagonal matrix
 
    Output:     Rlatt,  latticized network in original node ordering
                Rrp,    latticized network in node ordering used for
                        latticization
                ind_rp, node ordering used for latticization
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
                Sporns and Zwi (2004) Neuroinformatics 2:145
 
    2007-2012
    Mika Rubinov, UNSW
    Jonathan Power, WUSTL
    Olaf Sporns, IU

link_communities
 LINK_COMMUNITIES     Optimal overlapping community structure
 
    M = link_communities(W)
    M = link_communities(W,'complete');
 
    The optimal community structure is a subdivision of the network into
    groups of nodes which have a high number of within-group connections
    and a low number of between group connections.
 
    This algorithm uncovers overlapping community structure via
    hierarchical clustering of network links. This algorith is generalized
    for weighted/directed/fully-connected networks.
 
    Input:      W,                  directed (weighted or binary) connection matrix.
                type_clustering,    type of hierarchical clustering (optional)
                                        'single'        single-linkage (default)
                                        'complete'      complete-linkage
 
    Output:     M,                  nodal community-affiliation matrix
                                    binary matrix of size CxN [communities x nodes]
 
    NB: The algorithm can be slow and memory intensive.
 
    Reference: Ahn, Bagrow and Lehmann (2010) Nature 466, 761–764.
 
    Mika Rubinov, U Cambridge, 2014-2015

local_assortativity_wu_sign
 LOCAL_ASSORTATIVITY_WU_SIGN     Local Assortativity
 
    [loc_assort_pos,loc_assort_neg] = local_assortativity_wu_sign(W);
 
    Local Assortativity measures the extent to which nodes are connected to
    nodes of similar strength (vs. higher or lower strength). Adapted from
    Thedchanamoorthy et al. (2014)'s formula to allow weighted/signed 
    networks (node degree replaced with node strength). Note, output values 
    sum to total assortativity. 
 
    Inputs:     W,        undirected connection matrix with positive and
                          negative weights
 
    Output:     loc_assort_pos, local assortativity from positive weights
 
                loc_assort_neg, local assortativity from negative weights
 
    Reference: Thedchanamoorthy G, Piraveenan M, Kasthuriratna D, 
               Senanayake U. Proc Comp Sci (2014) 29:2449-2461.
 
 
    Jeff Spielberg, Boston University

make_motif34lib
 MAKE_MOTIF34LIB        Auxiliary motif library function
 
    make_motif34lib;
 
    This function generates the motif34lib.mat library required for all
    other motif computations.
 
 
    Mika Rubinov, UNSW, 2007-2010

makeevenCIJ
 MAKEEVENCIJ        Synthetic modular small-world network
 
    CIJ = makeevenCIJ(N,K,sz_cl);
 
    This function generates a random, directed network with a specified 
    number of fully connected modules linked together by evenly distributed
    remaining random connections.
 
    Inputs:     N,      number of vertices (must be power of 2)
                K,      number of edges
                sz_cl,  size of clusters (power of 2)
 
    Outputs:    CIJ,    connection matrix
 
    Notes:  N must be a power of 2.
            A warning is generated if all modules contain more edges than K.
            Cluster size is 2^sz_cl;
 
 
    Olaf Sporns, Indiana University, 2005/2007

makefractalCIJ
 MAKEFRACTALCIJ     Synthetic hierarchical modular network
 
    [CIJ,K] = makefractalCIJ(mx_lvl,E,sz_cl);
 
    This function generates a directed network with a hierarchical modular
    organization. All modules are fully connected and connection density 
    decays as 1/(E^n), with n = index of hierarchical level.
 
    Inputs:     mx_lvl,     number of hierarchical levels, N = 2^mx_lvl
                E,          connection density fall-off per level
                sz_cl,      size of clusters (power of 2)
 
    Outputs:    CIJ,        connection matrix
                K,          number of connections present in the output CIJ
 
 
  Olaf Sporns, Indiana University, 2005/2007

makelatticeCIJ
 MAKELATTICECIJ     Synthetic lattice network
 
    CIJ = makelatticeCIJ(N,K);
 
    This function generates a directed lattice network without toroidal 
    boundary counditions (i.e. no ring-like "wrapping around").
 
    Inputs:     N,      number of vertices
                K,      number of edges
 
    Outputs:    CIJ,    connection matrix
 
    Note: The lattice is made by placing connections as close as possible 
    to the main diagonal, without wrapping around. No connections are made 
    on the main diagonal. In/Outdegree is kept approx. constant at K/N.
 
 
    Olaf Sporns, Indiana University, 2005/2007

makerandCIJ_dir
 MAKERANDCIJ_DIR        Synthetic directed random network
 
    CIJ = makerandCIJ_dir(N,K);
 
    This function generates a directed random network
 
    Inputs:     N,      number of vertices
                K,      number of edges
 
    Output:     CIJ,    directed random connection matrix
 
    Note: no connections are placed on the main diagonal.
 
 
  Olaf Sporns, Indiana University, 2007/2008

makerandCIJ_und
 MAKERANDCIJ_UND        Synthetic directed random network
 
    CIJ = makerandCIJ_und(N,K);
 
    This function generates an undirected random network
 
    Inputs:     N,      number of vertices
                K,      number of edges
 
    Output:     CIJ,    undirected random connection matrix
 
    Note: no connections are placed on the main diagonal.
 
 
  Olaf Sporns, Indiana University, 2007/2008

makerandCIJdegreesfixed
 MAKERANDCIJDEGREESFIXED        Synthetic directed random network
 
    CIJ = makerandCIJdegreesfixed(N,K);
 
    This function generates a directed random network with a specified 
    in-degree and out-degree sequence. The function returns a flag, 
    denoting whether the algorithm succeeded or failed.
 
    Inputs:     in,     indegree vector
                out,    outdegree vector
 
    Output:     CIJ,    binary directed connectivity matrix
                flag,   flag=1 if the algorithm succeeded; flag=0 otherwise
 
 
    Notes:  Necessary conditions include:
                length(in) = length(out) = n
                sum(in) = sum(out) = k
                in(i), out(i) < n-1
                in(i) + out(j) < n+2
                in(i) + out(i) < n
 
            No connections are placed on the main diagonal
 
 
  Aviad Rubinstein, Indiana University 2005/2007

makeringlatticeCIJ
 MAKERINGLATTICECIJ     Synthetic lattice network
 
    CIJ = makeringlatticeCIJ(N,K);
 
    This function generates a directed lattice network with toroidal 
    boundary counditions (i.e. with ring-like "wrapping around").
 
    Inputs:     N,      number of vertices
                K,      number of edges
 
    Outputs:    CIJ,    connection matrix
 
    Note: The lattice is made by placing connections as close as possible 
    to the main diagonal, with wrapping around. No connections are made 
    on the main diagonal. In/Outdegree is kept approx. constant at K/N.
 
 
    Olaf Sporns, Indiana University, 2005/2007

maketoeplitzCIJ
 MAKETOEPLITZCIJ    A synthetic directed network with Gaussian drop-off of
                    connectivity with distance
 
    CIJ = maketoeprandCIJ(N,K,s)
 
    This function generates a directed network with a Gaussian drop-off in
    edge density with increasing distance from the main diagonal. There are
    toroidal boundary counditions (i.e. no ring-like "wrapping around").
 
    Inputs:     N,      number of vertices
                K,      number of edges
                s,      standard deviation of toeplitz
 
    Output:     CIJ,    connection matrix
 
    Note: no connections are placed on the main diagonal.
 
 
  Olaf Sporns, Indiana University, 2005/2007

matching_ind
 MATCHING_IND       Matching index
 
    [Min,Mout,Mall] = matching_ind(CIJ);
 
    For any two nodes u and v, the matching index computes the amount of
    overlap in the connection patterns of u and v. Self-connections and
    u-v connections are ignored. The matching index is a symmetric 
    quantity, similar to a correlation or a dot product.
 
    Input:      CIJ,    connection/adjacency matrix
 
    Output:     Min,    matching index for incoming connections
                Mout,   matching index for outgoing connections
                Mall,   matching index for all connections
 
    Notes:
        Does not use self- or cross connections for comparison.
        Does not use connections that are not present in BOTH u and v.
        All output matrices are calculated for upper triangular only.
 
 
  Olaf Sporns, Indiana University, 2002/2007/2008

matching_ind_und
 MATCHING_IND_UND       matching index
 
    M0 = MATCHING_IND_UND(CIJ) computes matching index for undirected
    graph specified by adjacency matrix CIJ. Matching index is a measure of
    similarity between two nodes' connectivity profiles (excluding their
    mutual connection, should it exist).
 
    Inputs:     CIJ,    undirected adjacency matrix
 
    Outputs:    M0,     matching index matrix.
 
    Richard Betzel, Indiana University, 2013
 

mean_first_passage_time
  MEAN_FIRST_PASSAGE_TIME           Mean first passage time
 
    MFPT = mean_first_passage_time(adj)
 
    The first passage time (MFPT) from i to j is the expected number of
    steps it takes a random walker starting at node i to arrive for the
    first time at node j. The mean first passage time is not a
    symmetric measure: mfpt(i,j) may be different from mfpt(j,i).  
  
    Input:
        adj,    Weighted/Unweighted, directed/undirected adjacency matrix 
 
    Output:
        MFPT, 	Pairwise mean first passage time matrix.   
 
 
    References: Goñi J, et al (2013) PLoS ONE
                
    Joaquin Goñi, IU Bloomington, 2012

mleme_constraint_model
 MLEME_CONSTRAINT_MODEL     Unbiased sampling of networks with soft constraints
 
    W0 = mleme_constraint_model(samp, W);
    W0 = mleme_constraint_model(samp, W, M);
    W0 = mleme_constraint_model(samp, W, M, Lo, Li, Lm);
    [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W, M, Lo, Li, Lm, opts);
 
    This function returns an ensemble of unbiasedly sampled networks with
    weighted node-strength and module-weight constraints. These constraints
    are soft in that they are satisfied on average for the full network
    ensemble but not, in general, for each individual network.
 
    Inputs (for a network with n nodes, m modules and c constraints):
 
        samp,   Number of networks to sample.
 
        W,      (length n) square directed and weighted connectivity
                matrix. All weights must be nonnegative integers. Note that
                real-valued weights may be converted to integers with
                arbitrary precision through rescaling and rounding, e.g.
                W_int = round(10^precision * W_real).
 
        M,      (length n) module affiliation vector. This vector is often
                obtained as the output of a community detection algorithm.
                The vector must contain nonnegative integers, with zeros
                specifying nodes which are not part of any community. This
                input may be left empty if there are no module constraints.
 
        Lo,     (length n) out-strength constraint logical vector. This
                vector specifies out-strength constraints for each node.
                Alternatively, it is possible to specify 1 to constrain all
                out-strengths or 0 for no constraints. Empty or no input
                results in default behavour (no constraints).
 
        Lo,     (length n) in-strength constraint logical vector. This
                vector specifies in-strength constraints for each node.
                Alternatively, it is possible to specify 1 to constrain all
                in-strengths or 0 for no constraints. Empty or no input
                results in default behavour (no constraints).
 
        Lm,     (length m) module-weight constraint logical matrix. This
                matrix specifies module-weight constraints for all pairs of
                modules. Alternatively, it is possible to specify
                2 to constrain all inter-module and intra-module weights,
                1 to constrain all intra-module weights, or 0  for no
                constraints. Empty or no input results in default behavour
                (no constraints).
 
        opts,   optional argument: pass optimization and display options with optimset.
                Default: optimset('MaxFunEvals', 1e6*c, 'MaxIter', 1e6, 'Display', 'iter');
 
 
    Outputs:
        W0,     an ensemble of sampled networks with constraints.
 
        E0,     expected weights matrix.
 
        P0,     probability matrix.
 
        Delt0,  algorithm convergence error.
 
 
    Algorithm:
                Maximum-likelihood estimation of network probability
                distribution by numerical solution of systems of nonlinear
                equations, and sampling of individual networks directly
                from this distribution.
 
 
    Notes:
                Empirical connection weights are
                not preserved. Constraint errors are guaranteed to vanish
                in the limit of the full network ensemble.
 
 
    Examples:
                % get community structure of a weighted network W
                M = community_louvain(W, 2);
 
                % specify node and module constraints
                n = length(W);              % number of nodes
                m = max(M);                 % number of modules
                Lo = true(n, 1);            % out-strength constraints
                Li = true(n, 1);            % in-strength constraints
                Lm = eye(m);                % module-weight constraints
 
                % sample networks with the above constraints
                [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W, M, Lo, Li, Lm);
 
                % equivalent formulation
                [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W, M, 1, 1, 1);
 
                % alternative: sample networks with average weight constraints only
                [W0, E0, P0, Delt0] = mleme_constraint_model(samp, W);
 
 
    References: Squartini and Garlaschelli (2011) New J Phys 13:083001
                Rubinov (2016) Nat Commun 7:13812
 
 
    2016, Mika Rubinov, Janelia HHMI

modularity_dir
 MODULARITY_DIR     Optimal community structure and modularity
 
    Ci = modularity_dir(W);
    [Ci Q] = modularity_dir(W);
 
    The optimal community structure is a subdivision of the network into
    nonoverlapping groups of nodes in a way that maximizes the number of
    within-group edges, and minimizes the number of between-group edges. 
    The modularity is a statistic that quantifies the degree to which the
    network may be subdivided into such clearly delineated groups. 
 
    Inputs:
        W,
            directed weighted/binary connection matrix
        gamma,
            resolution parameter (optional)
                gamma>1,        detects smaller modules
                0<=gamma<1,     detects larger modules
                gamma=1,        classic modularity (default)
 
    Outputs:    
        Ci,     optimal community structure
        Q,      maximized modularity
 
    Note:
        This algorithm is essentially deterministic. The only potential
        source of stochasticity occurs at the iterative finetuning step, in
        the presence of non-unique optimal swaps. However, the present
        implementation always makes the first available optimal swap and
        is therefore deterministic.
 
    References:
        Leicht and Newman (2008) Phys Rev Lett 100:118703.
        Reichardt and Bornholdt (2006) Phys Rev E 74:016110.
 
    2008-2016
    Mika Rubinov, UNSW
    Jonathan Power, WUSTL
    Dani Bassett, UCSB
    Xindi Wang, Beijing Normal University
    Roan LaPlante, Martinos Center for Biomedical Imaging

modularity_und
 MODULARITY_UND     Optimal community structure and modularity
 
    Ci = modularity_und(W);
    [Ci Q] = modularity_und(W,gamma);
 
    The optimal community structure is a subdivision of the network into
    nonoverlapping groups of nodes in a way that maximizes the number of
    within-group edges, and minimizes the number of between-group edges.
    The modularity is a statistic that quantifies the degree to which the
    network may be subdivided into such clearly delineated groups.
 
    Inputs:
        W,
            undirected weighted/binary connection matrix
        gamma,
            resolution parameter (optional)
                gamma>1,        detects smaller modules
                0<=gamma<1,     detects larger modules
                gamma=1,        classic modularity (default)
 
    Outputs:    
        Ci,     optimal community structure
        Q,      maximized modularity
 
    Note:
        This algorithm is essentially deterministic. The only potential
        source of stochasticity occurs at the iterative finetuning step, in
        the presence of non-unique optimal swaps. However, the present
        implementation always makes the first available optimal swap and
        is therefore deterministic.
 
    References: 
        Newman (2006) -- Phys Rev E 74:036104, PNAS 23:8577-8582.
        Reichardt and Bornholdt (2006) Phys Rev E 74:016110.
 
    2008-2016
    Mika Rubinov, UNSW
    Jonathan Power, WUSTL
    Dani Bassett, UCSB
    Xindi Wang, Beijing Normal University
    Roan LaPlante, Martinos Center for Biomedical Imaging

module_degree_zscore
 MODULE_DEGREE_ZSCORE       Within-module degree z-score
 
    Z=module_degree_zscore(W,Ci,flag);
 
    The within-module degree z-score is a within-module version of degree
    centrality.
 
    Inputs:     W,      binary/weighted, directed/undirected connection matrix
                Ci,     community affiliation vector
                flag,   0, undirected graph (default)
                        1, directed graph: out-degree
                        2, directed graph: in-degree
                        3, directed graph: out-degree and in-degree
 
    Output:     Z,      within-module degree z-score.
 
    Reference: Guimera R, Amaral L. Nature (2005) 433:895-900.
 
 
    Mika Rubinov, UNSW, 2008-2010

motif3funct_bin
 MOTIF3FUNCT_BIN       Frequency of functional class-3 motifs
 
    [f,F] = motif3funct_bin(A);
 
    *Structural motifs* are patterns of local connectivity in complex
    networks. In contrast, *functional motifs* are all possible subsets of
    patterns of local connectivity embedded within structural motifs. Such
    patterns are particularly diverse in directed networks. The motif
    frequency of occurrence around an individual node is known as the motif
    fingerprint of that node. The total motif frequency of occurrence in
    the whole network is correspondingly known as the motif fingerprint of
    the network.
 
    Input:      A,      binary directed connection matrix
 
    Output:     F,      node motif frequency fingerprint
                f,      network motif frequency fingerprint
 
    Notes: 
        1. The function find_motif34.m outputs the motif legend.
        2. There is a source of possible confusion in motif terminology.
           Motifs ("structural" and "functional") are most frequently
           considered only in the context of anatomical brain networks
           (Sporns and Kötter, 2004). On the other hand, motifs are not
           commonly studied in undirected networks, due to the paucity of
           local undirected connectivity patterns.
 
    References: Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif3funct_wei
 MOTIF3FUNCT_WEI       Intensity and coherence of functional class-3 motifs
 
    [I,Q,F] = motif3funct_wei(W);
 
    *Structural motifs* are patterns of local connectivity in complex
    networks. In contrast, *functional motifs* are all possible subsets of
    patterns of local connectivity embedded within structural motifs. Such
    patterns are particularly diverse in directed networks. The motif
    frequency of occurrence around an individual node is known as the motif
    fingerprint of that node. The motif intensity and coherence are
    weighted generalizations of the motif frequency. The motif
    intensity is equivalent to the geometric mean of weights of links
    comprising each motif. The motif coherence is equivalent to the ratio
    of geometric and arithmetic means of weights of links comprising each
    motif.  
 
    Input:      W,      weighted directed connection matrix
                        (all weights must be between 0 and 1)
 
    Output:     I,      node motif intensity fingerprint
                Q,      node motif coherence fingerprint
                F,      node motif frequency fingerprint
 
    Notes: 
        1. The function find_motif34.m outputs the motif legend.
        2. Average intensity and coherence are given by I./F and Q./F
        3. All weights must be between 0 and 1. This may be achieved using
           the weight_conversion.m function, as follows: 
           W_nrm = weight_conversion(W, 'normalize');
    	4. There is a source of possible confusion in motif terminology.
           Motifs ("structural" and "functional") are most frequently
           considered only in the context of anatomical brain networks
           (Sporns and Kötter, 2004). On the other hand, motifs are not
           commonly studied in undirected networks, due to the paucity of
           local undirected connectivity patterns.
 
    References: Onnela et al. (2005), Phys Rev E 71:065103
                Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif3struct_bin
 MOTIF3STRUCT_BIN       Frequency of structural class-3 motifs
 
    [f,F] = motif3struct_bin(A);
 
    Structural motifs are patterns of local connectivity in complex
    networks. Such patterns are particularly diverse in directed networks.
    The motif frequency of occurrence around an individual node is known as
    the motif fingerprint of that node. The total motif frequency of
    occurrence in the whole network is correspondingly known as the
    motif fingerprint of the network. 
 
    Input:      A,      binary directed connection matrix
 
    Output:     F,      node motif frequency fingerprint
                f,      network motif frequency fingerprint
 
    Note: The function find_motif34.m outputs the motif legend.
 
    References: Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif3struct_wei
 MOTIF3STRUCT_WEI       Intensity and coherence of structural class-3 motifs
 
    [I,Q,F] = motif3struct_wei(W);
 
    Structural motifs are patterns of local connectivity in complex
    networks. Such patterns are particularly diverse in directed networks.
    The motif frequency of occurrence around an individual node is known as
    the motif fingerprint of that node. The motif intensity and coherence
    are weighted generalizations of the motif frequency. The motif
    intensity is equivalent to the geometric mean of weights of links
    comprising each motif. The motif coherence is equivalent to the ratio
    of geometric and arithmetic means of weights of links comprising each
    motif.
 
    Input:      W,      weighted directed connection matrix
                        (all weights must be between 0 and 1)
 
    Output:     I,      node motif intensity fingerprint
                Q,      node motif coherence fingerprint
                F,      node motif frequency fingerprint
 
    Notes: 
        1. The function find_motif34.m outputs the motif legend.
        2. Average intensity and coherence are given by I./F and Q./F
        3. All weights must be between 0 and 1. This may be achieved using
           the weight_conversion.m function, as follows:
           W_nrm = weight_conversion(W, 'normalize');
 
    References: Onnela et al. (2005), Phys Rev E 71:065103
                Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369%
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif4funct_bin
 MOTIF4FUNCT_BIN       Frequency of functional class-4 motifs
 
    [f,F] = motif4funct_bin(A);
 
    *Structural motifs* are patterns of local connectivity in complex
    networks. In contrast, *functional motifs* are all possible subsets of
    patterns of local connectivity embedded within structural motifs. Such
    patterns are particularly diverse in directed networks. The motif
    frequency of occurrence around an individual node is known as the motif
    fingerprint of that node. The total motif frequency of occurrence in
    the whole network is correspondingly known as the motif fingerprint of
    the network.
 
    Input:      A,      binary directed connection matrix
 
    Output:     F,      node motif frequency fingerprint
                f,      network motif frequency fingerprint
 
    Notes: 
        1. The function find_motif34.m outputs the motif legend.
        2. There is a source of possible confusion in motif terminology.
           Motifs ("structural" and "functional") are most frequently
           considered only in the context of anatomical brain networks
           (Sporns and Kötter, 2004). On the other hand, motifs are not
           commonly studied in undirected networks, due to the paucity of
           local undirected connectivity patterns.
 
    References: Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif4funct_wei
 MOTIF4FUNCT_WEI       Intensity and coherence of functional class-4 motifs
 
    [I,Q,F] = motif4funct_wei(W);
 
    *Structural motifs* are patterns of local connectivity in complex
    networks. In contrast, *functional motifs* are all possible subsets of
    patterns of local connectivity embedded within structural motifs. Such
    patterns are particularly diverse in directed networks. The motif
    frequency of occurrence around an individual node is known as the motif
    fingerprint of that node. The motif intensity and coherence are
    weighted generalizations of the motif frequency. The motif
    intensity is equivalent to the geometric mean of weights of links
    comprising each motif. The motif coherence is equivalent to the ratio
    of geometric and arithmetic means of weights of links comprising each
    motif.
 
    Input:      W,      weighted directed connection matrix
                        (all weights must be between 0 and 1)
 
    Output:     I,      node motif intensity fingerprint
                Q,      node motif coherence fingerprint
                F,      node motif frequency fingerprint
 
    Notes: 
        1. The function find_motif34.m outputs the motif legend.
        2. Average intensity and coherence are given by I./F and Q./F
        3. All weights must be between 0 and 1. This may be achieved using
           the weight_conversion.m function, as follows: 
           W_nrm = weight_conversion(W, 'normalize');
        4. There is a source of possible confusion in motif terminology.
           Motifs ("structural" and "functional") are most frequently
           considered only in the context of anatomical brain networks
           (Sporns and Kötter, 2004). On the other hand, motifs are not
           commonly studied in undirected networks, due to the paucity of
           local undirected connectivity patterns.
 
    References: Onnela et al. (2005), Phys Rev E 71:065103
                Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369%
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif4struct_bin
 MOTIF4STRUCT_BIN       Frequency of structural class-4 motifs
 
    [f,F] = motif4struct_bin(A);
 
    Structural motifs are patterns of local connectivity in complex
    networks. Such patterns are particularly diverse in directed networks.
    The motif frequency of occurrence around an individual node is known as
    the motif fingerprint of that node. The total motif frequency of
    occurrence in the whole network is correspondingly known as the
    motif fingerprint of the network. 
 
    Input:      A,      binary directed connection matrix
 
    Output:     F,      node motif frequency fingerprint
                f,      network motif frequency fingerprint
 
    Note: The function find_motif34.m outputs the motif legend.
 
    References: Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

motif4struct_wei
 MOTIF4STRUCT_WEI       Intensity and coherence of structural class-4 motifs
 
    [I,Q,F] = motif4struct_wei(W);
 
    Structural motifs are patterns of local connectivity in complex
    networks. Such patterns are particularly diverse in directed networks.
    The motif frequency of occurrence around an individual node is known as
    the motif fingerprint of that node. The motif intensity and coherence
    are weighted generalizations of the motif frequency. The motif
    intensity is equivalent to the geometric mean of weights of links
    comprising each motif. The motif coherence is equivalent to the ratio
    of geometric and arithmetic means of weights of links comprising each
    motif.
 
    Input:      W,      weighted directed connection matrix
                        (all weights must be between 0 and 1)
 
    Output:     I,      node motif intensity fingerprint
                Q,      node motif coherence fingerprint
                F,      node motif frequency fingerprint
 
    Notes: 
        1. The function find_motif34.m outputs the motif legend.
        2. Average intensity and coherence are given by I./F and Q./F
        3. All weights must be between 0 and 1. This may be achieved using
           the weight_conversion.m function, as follows: 
           W_nrm = weight_conversion(W, 'normalize');
 
    References: Onnela et al. (2005), Phys Rev E 71:065103
                Milo et al. (2002) Science 298:824-827
                Sporns O, Kötter R (2004) PLoS Biol 2: e369%
 
 
    Mika Rubinov, UNSW/U Cambridge, 2007-2015

null_model_dir_sign
 NULL_MODEL_DIR_SIGN     Directed random graphs with preserved weight,
                         degree and strength distributions
 
    W0 = null_model_dir_sign(W);
    W0 = null_model_dir_sign(W,bin_swaps);
    W0 = null_model_dir_sign(W,bin_swaps,wei_freq);
    [W0 R] = null_model_dir_sign(W,bin_swaps,wei_freq);
 
    This function randomizes an directed network with positive and
    negative weights, while preserving the degree and strength
    distributions. This function calls randmio_dir_signed.m
 
    Inputs: W,          Directed weighted connection matrix
            bin_swaps,  Average number of swaps of each edge in binary randomization.
                            bin_swap=5 is the default (each edge rewired 5 times)
                            bin_swap=0 implies no binary randomization 
            wei_freq,   Frequency of weight sorting in weighted randomization
                            wei_freq must be in the range of: 0 < wei_freq <= 1
                            wei_freq=1 implies that weights are sorted at each step
                                (default in older [<2011] versions of MATLAB)
                            wei_freq=0.1 implies that weights are sorted at each 10th step
                                (faster, default in newer versions of MATLAB)
 
    Output:     W0,     Randomized weighted connection matrix
                R,      Correlation coefficients between strength sequences
                            of input and output connection matrices
 
    Notes:
        The value of bin_swaps is ignored when binary topology is fully
    connected (e.g. when the network has no negative weights).
        Randomization may be better (and execution time will be slower) for
    higher values of bin_swaps and wei_freq. Higher values of bin_swaps may
    enable a more random binary organization, and higher values of wei_freq
    may enable a more accurate conservation of strength sequences.
        R are the correlation coefficients between positive and negative
    in-strength and out-strength sequences of input and output connection
    matrices and are used to evaluate the accuracy with which strengths
    were preserved. Note that correlation coefficients may be a rough
    measure of strength-sequence accuracy and one could implement more
    formal tests (such as the Kolmogorov-Smirnov test) if desired.
 
    Example usage:
 
    %Create random directed weights matrix
 
    W=randn(100);
 
    %Compute one instance of null model (slow execution time):
    %bin_swaps=5,   rewire each binary edge 5 times on average
    %wei_freq=1,    sort all edges at every step
 
    tic; [W0_slow R_slow]=null_model_dir_sign(W,5,1); R_slow, toc
 
    R_slow =
        0.9795    0.9724    0.9772    0.9773
    Elapsed time is 3.485388 seconds.
 
    %Compute another instance of of null model (fast execution time):
    %bin_swaps=5,   rewire each binary edge 5 times on average
    %wei_freq=0.1,  sort all edges at every 10th step (10=1/0.1)
 
    tic; [W0_fast R_fast]=null_model_dir_sign(W,5,0.1); R_fast, toc
 
    R_fast =
        0.9655    0.9652    0.9717    0.9804
    Elapsed time is 0.763831 seconds.
 
 
    Reference: Rubinov and Sporns (2011) Neuroimage 56:2068-79
 
 
    2011-2015, Mika Rubinov, U Cambridge

null_model_und_sign
 NULL_MODEL_UND_SIGN     Random graphs with preserved weight, degree and
                         strength distributions
 
    W0 = null_model_und_sign(W);
    W0 = null_model_und_sign(W,bin_swaps);
    W0 = null_model_und_sign(W,bin_swaps,wei_freq);
    [W0 R] = null_model_und_sign(W,bin_swaps,wei_freq);
 
    This function randomizes an undirected network with positive and
    negative weights, while preserving the degree and strength
    distributions. This function calls randmio_und_signed.m
 
    Inputs: W,          Undirected weighted connection matrix
            bin_swaps,  Average number of swaps of each edge in binary randomization.
                            bin_swap=5 is the default (each edge rewired 5 times)
                            bin_swap=0 implies no binary randomization
            wei_freq,   Frequency of weight sorting in weighted randomization
                            wei_freq must be in the range of: 0 < wei_freq <= 1
                            wei_freq=1 implies that weights are resorted at each step
                                    (default in older [<2011] versions of MATLAB)
                            wei_freq=0.1 implies that weights are sorted at each 10th step
                                    (faster, default in newer versions of Matlab)
 
    Output:     W0,     Randomized weighted connection matrix
                R,      Correlation coefficient between strength sequences
                            of input and output connection matrices
 
    Notes:
        The value of bin_swaps is ignored when binary topology is fully
    connected (e.g. when the network has no negative weights).
        Randomization may be better (and execution time will be slower) for
    higher values of bin_swaps and wei_freq. Higher values of bin_swaps may
    enable a more random binary organization, and higher values of wei_freq
    may enable a more accurate conservation of strength sequences.
        R are the correlation coefficients between positive and negative
    strength sequences of input and output connection matrices and are
    used to evaluate the accuracy with which strengths were preserved. Note
    that correlation coefficients may be a rough measure of
    strength-sequence accuracy and one could implement more formal tests
    (such as the Kolmogorov-Smirnov test) if desired.
 
    Example usage:
 
    %Create random weights matrix
 
    W=tril(randn(100),-1); W=W+W.';
 
    %Compute one instance of null model (slow execution time):
    %bin_swaps=5,   rewire each binary edge 5 times on average
    %wei_freq=1,    sort all edges at every step
 
    tic; [W0_slow R_slow]=null_model_und_sign(W,5,1); R_slow, toc
 
    R_slow =
        0.9720    0.9746
    Elapsed time is 2.112715 seconds.
 
    %Compute another instance of of null model (fast execution time):
    %bin_swaps=5,   rewire each binary edge 5 times on average
    %wei_freq=0.1,  sort all edges at every 10th step (10=1/0.1)
 
    tic; [W0_fast R_fast]=null_model_und_sign(W,5,0.1); R_fast, toc
 
    R_fast =
        0.9623    0.9789
    Elapsed time is 0.584797 seconds.
 
 
    Reference: Rubinov and Sporns (2011) Neuroimage 56:2068-79
 
 
    2011-2015, Mika Rubinov, U Cambridge

pagerank_centrality
 PAGERANK_CENTRALITY       PageRank centrality
 
    r = pagerank_centrality(A, d, falff)
 
    The PageRank centrality is a variant of eigenvector centrality. This
    function computes the PageRank centrality of each vertex in a graph.
 
    Formally, PageRank is defined as the stationary distribution achieved
    by instantiating a Markov chain on a graph. The PageRank centrality of
    a given vertex, then, is proportional to the number of steps (or amount
    of time) spent at that vertex as a result of such a process. 
 
    The PageRank index gets modified by the addition of a damping factor,
    d. In terms of a Markov chain, the damping factor specifies the
    fraction of the time that a random walker will transition to one of its
    current state's neighbors. The remaining fraction of the time the
    walker is restarted at a random vertex. A common value for the damping
    factor is d = 0.85.
 
    Inputs:     A,      adjacency matrix
                d,      damping factor
            falff,      initial page rank probability (non-negative)
 
    Outputs:    r,      vectors of page rankings
 
    Note: The algorithm will work well for smaller matrices (number of
    nodes around 1000 or less) 
 
    References:
 
    [1]. GeneRank: Using search engine technology for the analysis of
    microarray experiments, by Julie L. Morrison, Rainer Breitling, Desmond
    J. Higham and David R. Gilbert, BMC Bioinformatics, 6:233, 2005.
  	[2]. Boldi P, Santini M, Vigna S (2009) PageRank: Functional
  	dependencies. ACM Trans Inf Syst 27, 1-23.
 
    Xi-Nian Zuo, Institute of Psychology, Chinese Academy of Sciences, 2011
    Rick Betzel, Indiana University, 2012

participation_coef
 PARTICIPATION_COEF     Participation coefficient
 
    P = participation_coef(W,Ci);
 
    Participation coefficient is a measure of diversity of intermodular
    connections of individual nodes.
 
    Inputs:     W,      binary/weighted, directed/undirected connection matrix
                Ci,     community affiliation vector
                flag,   0, undirected graph (default)
                        1, directed graph: out-degree
                        2, directed graph: in-degree
 
    Output:     P,      participation coefficient
 
    Reference: Guimera R, Amaral L. Nature (2005) 433:895-900.
 
 
    2008-2015
    Mika Rubinov, UNSW/U Cambridge
    Alex Fornito, University of Melbourne

participation_coef_sign
 PARTICIPATION_COEF_SIGN     Participation coefficient
 
    [Ppos Pneg] = participation_coef_sign(W,Ci);
 
    Participation coefficient is a measure of diversity of intermodular
    connections of individual nodes.
 
    Inputs:     W,      undirected connection matrix with positive and
                        negative weights
 
                Ci,     community affiliation vector
 
    Output:     Ppos,   participation coefficient from positive weights
 
                Pneg,   participation coefficient from negative weights
 
    Reference: Guimera R, Amaral L. Nature (2005) 433:895-900.
 
 
    2011, Mika Rubinov, UNSW

partition_distance
 PARTITION_DISTANCE     Distance or similarity between community partitions
 
    This function quantifies information-theoretic distance (normalized
    variation of information) or similarity (normalized mutual information)
    between community partitions.
 
    VIn        = partition_distance(Cx);
    VIn        = partition_distance(Cx, Cy);
    [VIn, MIn] = partition_distance(Cx, Cy);
 
    Inputs:
        Cx,
            Community partition vector or matrix of n rows and p columns,
            n is the number of network nodes, and p is the number of input
            community partitions (in the case of vector input p=1).
 
        Cy (optional argument),
            Community partition vector or matrix of n rows and q columns. n
            is the number of nodes (must be equal to the number of nodes in
            Cq) and q is the number of input community partitions (may be
            different to the number of nodes in Cq). This argument may be
            omitted, in which case, the partition distance is computed
            between all pairwise partitions of Cx.
 
    Outputs:
        VIn,
            Normalized variation of information ([p, q] matrix)
 
        MIn,
            Normalized mutual information ([p, q] matrix)
 
    Notes:
        Mathematical definitions.
 
            VIn = [H(X) + H(Y) - 2MI(X, Y)]/log(n)
            MIn = 2MI(X, Y) / [H(X) + H(Y)]
 
            where H is the entropy and MI is the mutual information
 
 
    Reference: Meila M (2007) J Multivar Anal 98, 873-895.
 
 
    2011-2017, Mika Rubinov, UNSW, Janelia HHMI

path_transitivity
  PATH_TRANSITIVITY             Transitivity based on shortest paths
 
    T = path_transitivity(W,transform)
 
    This function computes the density of local detours (triangles) that
    are available along the shortest-paths between all pairs of nodes.
 
    Inputs:
 
        W,
            unweighted/weighted undirected connection *weight* OR *length*
            matrix.
 
 
        transform,
            If the input matrix is a connection *weight* matrix, specify a
            transform that map input connection weights to connection
            lengths. Two transforms are available.
                'log' -> l_ij = -log(w_ij)
                'inv' -> l_ij =    1/w_ij
 
            If the input matrix is a connection *length* matrix, do not
            specify a transform (or specify an empty transform argument).
 
 
    Output:
 
        T,
            matrix of pairwise path transitivity.
 
 
    Olaf Sporns, Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2014
 
    References: Goñi et al (2014) PNAS doi: 10.1073/pnas.131552911
 

predict_fc
  PREDICT_FC    Prediction of functional connectivity from structural connectivity
 
    [FCpre,FCcorr,beta,pred_data,R] = predict_fc(SC,FC,ED,pred_var,model)
    [FCpre,FCcorr,beta] = predict_fc(SC,FC,[],{'SPLwei_log','SIwei_log'},'quadratic')
 
    This function computes regression coefficients to predict FC from
    structural-based measures that are used as predictor variables.
 
    Inputs:
 
        SC,
            Weighted/unweighted undirected NxN Structural Connectivity matrix.
 
        FC,
            Functional connections. FC can be a NxN symmetric matrix or an
            ((N*(N-1))/2) x 1 vector containing the upper triangular
            elements of the square FC matrix (excluding diagonal elements).
 
        ED,
            Euclidean distance matrix or upper triangular vector of the
            matrix (optional)
 
        pred_var,
            Set of M predictors. These can be given as an KxM array where
            K = ((N*(N-1))/2) and M is the number of predictors.
            Alternatively, pred_var can be a cell with the names of network
            measures to be used as predictors. Accepted network measure
            names are:
                SPLbin        - Shortest-path length (binary)
                SPLwei_inv    - Shortest-path length computed with an inv transform
              	SPLwei_log    - Shortest-path length computed with a log transform
               	SPLdist       - Shortest-path length computed with no transform
               	SIbin         - Search Information of binary shortest-paths
               	SIwei_inv     - Search Information of shortest-paths computed with an inv transform
               	SIwei_log     - Search Information of shortest-paths computed with a log transform
               	SIdist        - Search Information of shortest-paths computed with no transform
               	T             - Path Transitivity
               	deltaMFPT     - Column-wise z-scored mean first passage time
                neighOverlap  - Neighborhood Overlap
               	MI            - Matching Index
 
        	If no predictors are specified, the defaults are {'SPLwei_log', 'SIwei_log'}.
 
    	model,
        	Specifies the order of the regression model used within
        	matlab's function regstats.m. 'model' can be any option
        	accepted by matlab's regstats.m function (e.g. 'linear',
        	'interaction', 'quadratic', etc.) If no model is specified,
        	'linear' is the default.
 
    Output:
 
        FCpre,
            Predicted NxN Functional Connectivity matrix
 
       	FCcorr,
            Pearson Correlation between PCpred and FC
 
       	beta,
            Regression Coefficients
 
       	pred_data,
            KxM array of predictors.
 
       	R,
            Output from regstats.m (e.g. 'beta', 'yhat', 'rsquare',
            'adjrsquare', 'tstat', 'r', 'mse', 'standres').
 
 
    References: Goñi et al (2014) PNAS,  833–838, doi: 10.1073/pnas.1315529111
 
 
    Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2014

randmio_dir
 RANDMIO_DIR    Random graph with preserved in/out degree distribution
 
    R = randmio_dir(W, ITER);
    [R eff] = randmio_dir(W, ITER);
 
    This function randomizes a directed network, while preserving the in-
    and out-degree distributions. In weighted networks, the function
    preserves the out-strength but not the in-strength distributions.
 
    Input:      W,      directed (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
 
    Output:     R,      randomized network
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
 
 
    2007-2012
    Mika Rubinov, UNSW
    Olaf Sporns, IU

randmio_dir_connected
 RANDMIO_DIR_CONNECTED    Random graph with preserved in/out degree distribution
 
    R = randmio_dir_connected(W, ITER);
    [R eff] = randmio_dir_connected(W, ITER);
 
    This function randomizes a directed network, while preserving the in-
    and out-degree distributions. In weighted networks, the function
    preserves the out-strength but not the in-strength distributions. The
    function also ensures that the randomized network maintains
    connectedness, the ability for every node to reach every other node in
    the network. The input network for this function must be connected.
 
    Input:      W,      directed (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
 
    Output:     R,      randomized network
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
 
 
    2007-2012
    Mika Rubinov, UNSW
    Olaf Sporns, IU

randmio_dir_signed
  RANDMIO_DIR_SIGNED        Random directed graph with preserved signed
                            in/out degree distribution
 
    R       = randmio_dir(W, ITER);
    [R,eff] = randmio_dir(W, ITER);
 
    This function randomizes a directed weighted network with positively and
    negatively signed connections, while preserving the positively and
    negatively signed in- and out-degree distributions. In weighted
    networks, the function preserves the out-strength but not the
    in-strength distributions.
 
    Input:      W,      directed (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
 
    Output:     R,      randomized network
                eff,    number of actual rewirings carried out
 
    Reference:  Maslov and Sneppen (2002) Science 296:910
 
 
    2011-2015
    Dani Bassett, UCSB
    Olaf Sporns,  Indiana U
    Mika Rubinov, U Cambridge

randmio_und
 RANDMIO_UND     Random graph with preserved degree distribution
 
    R = randmio_und(W,ITER);
    [R eff]=randmio_und(W, ITER);
 
    This function randomizes an undirected network, while preserving the
    degree distribution. The function does not preserve the strength
    distribution in weighted networks.
 
    Input:      W,      undirected (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
 
    Output:     R,      randomized network
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
 
 
    2007-2012
    Mika Rubinov, UNSW
    Jonathan Power, WUSTL
    Olaf Sporns, IU

randmio_und_connected
 RANDMIO_UND_CONNECTED     Random graph with preserved degree distribution
 
    R = randmio_und_connected(W,ITER);
    [R eff] = randmio_und_connected(W, ITER);
 
    This function randomizes an undirected network, while preserving the
    degree distribution. The function does not preserve the strength
    distribution in weighted networks. The function also ensures that the
    randomized network maintains connectedness, the ability for every node
    to reach every other node in the network. The input network for this
    function must be connected.
 
    Input:      W,      undirected (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
 
    Output:     R,      randomized network
                eff,    number of actual rewirings carried out
 
    References: Maslov and Sneppen (2002) Science 296:910
 
 
    2007-2012
    Mika Rubinov, UNSW
    Jonathan Power, WUSTL
    Olaf Sporns, IU

randmio_und_signed
  RANDMIO_UND_SIGNED	Random graph with preserved signed degree distribution
 
    R       = randmio_und_signed(W,ITER);
    [R,eff] = randmio_und_signed(W,ITER);
 
    This function randomizes an undirected network with positively and
    negatively signed connections, while preserving the positively and
    negatively signed degree distribution. The function does not preserve
    the strength distribution in weighted networks.
 
    Input:      W,      undirected (binary/weighted) connection matrix
                ITER,   rewiring parameter
                        (each edge is rewired approximately ITER times)
 
    Output:     R,      randomized network
                eff,    number of actual rewirings carried out
 
    Reference:  Maslov and Sneppen (2002) Science 296:910
 
 
    2011-2015
    Dani Bassett, UCSB
    Olaf Sporns,  Indiana U
    Mika Rubinov, U Cambridge

randomize_graph_partial_und
  RANDOMIZE_GRAPH_PARTIAL_UND    Swap edges with preserved degree sequence
 
    A = RANDOMIZE_GRAPH_PARTIAL_UND(A,B,MAXSWAP) takes adjacency matrices A 
    and B and attempts to randomize matrix A by performing MAXSWAP 
    rewirings. The rewirings will avoid any spots where matrix B is 
    nonzero.
 
    Inputs:       A,      undirected adjacency matrix
                  B,      edges to avoid
            MAXSWAP,      number of rewirings
 
    Outputs:      A,      randomized matrix
 
    Richard Betzel, Indiana University, 2013
 
    Notes:
    1. Based on the script randmio_und.m.
    2. Graph may become disconnected as a result of rewiring. Always
       important to check.
    3. A can be weighted, though the weighted degree sequence will not be
       preserved.
 

randomizer_bin_und
 RANDOMIZER_BIN_UND     Random graph with preserved in/out degree distribution
 
    R = randomizer_bin_und(A,alpha);
 
    This function randomizes a binary undirected network, while preserving 
    the degree distribution. The function directly searches for rewirable 
    edge pairs (rather than trying to rewire edge pairs at random), and 
    hence avoids long loops and works especially well in dense matrices.
 
    Inputs:     A,          binary undirected connection matrix
                alpha,      fraction of edges to rewire
 
    Outputs:    R,          randomized network
 
    References: Maslov and Sneppen (2002) Science 296:910
 
 
    Jonathan Power, WUSTL. 3/1/10.

reachdist
 REACHDIST      Reachability and distance matrices
 
    [R,D] = reachdist(CIJ);
 
    The binary reachability matrix describes reachability between all pairs
    of nodes. An entry (u,v)=1 means that there exists a path from node u
    to node v; alternatively (u,v)=0.
 
    The distance matrix contains lengths of shortest paths between all
    pairs of nodes. An entry (u,v) represents the length of shortest path 
    from node u to  node v. The average shortest path length is the 
    characteristic path length of the network.
 
    Input:      CIJ,     binary (directed/undirected) connection matrix
 
    Outputs:    R,       reachability matrix
                D,       distance matrix
 
    Note: faster but more memory intensive than "breadthdist.m".
 
    Algorithm: algebraic path count.
 
 
    Olaf Sporns, Indiana University, 2002/2007/2008

rentian_scaling_2d
  RENTIAN_SCALING_2D    Rentian scaling for networks embedded in two dimensions.
 
  [N,E] = rentian_scaling_2d(A,XY,n,tol)
 
  Physical Rentian scaling (or more simply Rentian scaling) is a property
  of systems that are cost-efficiently embedded into physical space. It is
  what is called a "topo-physical" property because it combines information
  regarding the topological organization of the graph with information
  about the physical placement of connections. Rentian scaling is present
  in very large scale integrated circuits, the C. elegans neuronal network,
  and morphometric and diffusion-based graphs of human anatomical networks.
  Rentian scaling is determined by partitioning the system into cubes,
  counting the number of nodes inside of each cube (N), and the number of
  edges traversing the boundary of each cube (E). If the system displays
  Rentian scaling, these two variables N and E will scale with one another
  in loglog space. The Rent's exponent is given by the slope of log10(E)
  vs. log10(N), and can be reported alone or can be compared to the
  theoretical minimum Rent's exponent to determine how cost efficiently the
  network has been embedded into physical space. Note: if a system displays
  Rentian scaling, it does not automatically mean that the system is
  cost-efficiently embedded (although it does suggest that). Validation
  occurs when comparing to the theoretical minimum Rent's exponent for that
  system.
 
  INPUTS:
 
  A:              MxM adjacency matrix.
                  Must be unweighted, binary, and symmetric.
  XY:             Matrix of node placement coordinates.
                  Must be in the form of an Mx2 matrix [x y], where M is the
                  number of nodes and x and y are column vectors of node
                  coordinates.
  n:              Number of partitions to compute. Each partition is a data
                  point. You want a large enough number to adequately
                  estimate the Rent's exponent.
  tol:            This should be a small value (for example 1e-6).
                  In order to mitigate the effects of boundary conditions due
                  to the finite size of the network, we only allow partitions
                  that are contained within the boundary of the network. This
                  is achieved by first computing the volume of the convex
                  hull of the node coordinates (V). We then ensure that the
                  volume of the convex hull computed on the original node
                  coordinates plus the coordinates of the randomly generated
                  partition (Vnew) is within a given tolerance of the
                  original (i.e. check abs(V - Vnew) < tol). Thus tol, should
                  be a small value in order to make sure the partitions are
                  contained largely within the boundary of the network, and
                  thus the number of nodes and edges within the box are not
                  skewed by finite size effects.
 
  OUTPUTS:
 
  N:              nx1 vector of the number of nodes in each of the n partitions.
  E:              nx1 vector of the number of edges crossing the boundary of
                  each partition.
 
  Subsequent Analysis:
 
      Rentian scaling plots are created by: figure; loglog(E,N,'*');
 
      To determine the Rent's exponent, p, we need to determine
      the slope of E vs. N in loglog space, which is the Rent's
      exponent. There are many ways of doing this with more or less
      statistical rigor. Robustfit in MATLAB is one such option:
 
          [b,stats] = robustfit(log10(N),log10(E))
 
      Then the Rent's exponent is b(1,2) and the standard error of the
      estimation is given by stats.se(1,2).
 
  Note: n=5000 was used in Bassett et al. 2010 in PLoS CB.
 
  Reference:
  Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg,
  Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient
  physical embedding of topologically complex information processing
  networks in brains and computer circuits. PLoS Comput Biol, 2010,
  6(4):e1000748.
 
  Modification History:
 
      2010:     Original (Dani Bassett)
      Dec 2016: Updated code so that both partition centers and partition
                sizes are chosen at random. Also added in a constraint on
                partition placement that prevents boxes from being located
                outside the edges of the network. This helps prevent skewed
                results due to boundary effects arising from the finite size
                of the network. (Lia Papadopoulos)

rentian_scaling_3d
  RENTIAN_SCALING_3D 	Rentian scaling for networks embedded in three dimensions.
 
  [N,E] = rentian_scaling_3d(A,XYZ,n,tol)
 
  Physical Rentian scaling (or more simply Rentian scaling) is a property
  of systems that are cost-efficiently embedded into physical space. It is
  what is called a "topo-physical" property because it combines information
  regarding the topological organization of the graph with information
  about the physical placement of connections. Rentian scaling is present
  in very large scale integrated circuits, the C. elegans neuronal network,
  and morphometric and diffusion-based graphs of human anatomical networks.
  Rentian scaling is determined by partitioning the system into cubes,
  counting the number of nodes inside of each cube (N), and the number of
  edges traversing the boundary of each cube (E). If the system displays
  Rentian scaling, these two variables N and E will scale with one another
  in loglog space. The Rent's exponent is given by the slope of log10(E)
  vs. log10(N), and can be reported alone or can be compared to the
  theoretical minimum Rent's exponent to determine how cost efficiently the
  network has been embedded into physical space. Note: if a system displays
  Rentian scaling, it does not automatically mean that the system is
  cost-efficiently embedded (although it does suggest that). Validation
  occurs when comparing to the theoretical minimum Rent's exponent for that
  system.
 
  INPUTS:
 
  A:              MxM adjacency matrix.
                  Must be unweighted, binary, and symmetric.
  XYZ:            Matrix of node placement coordinates.
                  Must be in the form of an Mx3 matrix [x y z], where M is
                  the number of nodes and x, y, z are column vectors of node
                  coordinates.
  n:              Number of partitions to compute. Each partition is a data
                  point. You want a large enough number to adequately
                  estimate the Rent's exponent.
  tol:            This should be a small value (for example 1e-6).
                  In order to mitigate the effects of boundary conditions due
                  to the finite size of the network, we only allow partitions
                  that are contained within the boundary of the network. This
                  is achieved by first computing the volume of the convex
                  hull of the node coordinates (V). We then ensure that the
                  volume of the convex hull computed on the original node
                  coordinates plus the coordinates of the randomly generated
                  partition (Vnew) is within a given tolerance of the
                  original (i.e. check abs(V - Vnew) < tol). Thus tol, should
                  be a small value in order to make sure the partitions are
                  contained largely within the boundary of the network, and
                  thus the number of nodes and edges within the box are not
                  skewed by finite size effects.
 
  OUTPUTS:
 
  N:              nx1 vector of the number of nodes in each of the n partitions.
  E:              nx1 vector of the number of edges crossing the boundary of
                  each partition.
 
  Subsequent Analysis:
 
      Rentian scaling plots are created by: figure; loglog(E,N,'*');
 
      To determine the Rent's exponent, p, we need to determine
      the slope of E vs. N in loglog space, which is the Rent's
      exponent. There are many ways of doing this with more or less
      statistical rigor. Robustfit in MATLAB is one such option:
 
          [b,stats] = robustfit(log10(N),log10(E))
 
      Then the Rent's exponent is b(1,2) and the standard error of the
      estimation is given by stats.se(1,2).
 
  Note: n=5000 was used in Bassett et al. 2010 in PLoS CB.
 
  Reference:
  Danielle S. Bassett, Daniel L. Greenfield, Andreas Meyer-Lindenberg,
  Daniel R. Weinberger, Simon W. Moore, Edward T. Bullmore. Efficient
  physical embedding of topologically complex information processing
  networks in brains and computer circuits. PLoS Comput Biol, 2010,
  6(4):e1000748.
 
  Modification History:
 
      2010:     Original (Dani Bassett)
      Dec 2016: Updated code so that both partition centers and partition
                sizes are chosen at random. Also added in a constraint on
                partition placement that prevents boxes from being located
                outside the edges of the network. This helps prevent skewed
                results due to boundary effects arising from the finite size
                of the network. (Lia Papadopoulos)
 

reorderMAT
 REORDERMAT         Reorder matrix for visualization
 
    [MATreordered,MATindices,MATcost] = reorderMAT(MAT,H,cost);
 
    This function reorders the connectivity matrix in order to place more
    edges closer to the diagonal. This often helps in displaying community
    structure, clusters, etc.
 
    Inputs:     MAT,            connection matrix
                H,              number of reordering attempts
                cost,           'line' or 'circ', for shape of lattice
                                (linear or ring lattice)
 
                MATreordered    reordered connection matrix
                MATindices      reordered indices
                MATcost         cost of reordered matrix
    
 
    Olaf Sporns, Indiana University

reorder_matrix
  REORDER_MATRIX         Matrix reordering for visualization
 
    [Mreordered,Mindices,cost] = reorder_matrix(M1,cost,flag)
 
    This function rearranges the nodes in matrix M1 such that the matrix
    elements are squeezed along the main diagonal.  The function uses a
    version of simulated annealing. 
 
    Inputs:     M1             = connection matrix (weighted or binary, 
                                 directed or undirected)
                cost           = 'line' or 'circ', for shape of lattice
                                 cost (linear or ring lattice)
 
                Mreordered     = reordered connection matrix
                Mindices       = reordered indices
                cost           = distance between M1 and Mreordered
 
    Note that in general, the outcome will depend on the initial condition
    (the setting of the random number seed).  Also, there is no good way to 
    determine optimal annealing parameters in advance - these paramters 
    will need to be adjusted "by hand" (particularly H, Texp, and T0).  
    For large and/or dense matrices, it is highly recommended to perform 
    exploratory runs varying the settings of 'H' and 'Texp' and then select 
    the best values.
 
    Based on extensive testing, it appears that T0 and Hbrk can remain
    unchanged in most cases.  Texp may be varied from 1-1/H to 1-10/H, for
    example.  H is the most important parameter - set to larger values as
    the problem size increases.  It is advisable to run this function
    multiple times and select the solution(s) with the lowest 'cost'.
 
    Setting 'Texp' to zero cancels annealing and uses a greedy algorithm
    instead.
 
    Yusuke Adachi, University of Tokyo 2010
    Olaf Sporns, Indiana University 2010

reorder_mod
 REORDER_MOD         Reorder connectivity matrix by modular structure
 
    On = reorder_mod(W,M);
    [On Wr] = reorder_mod(W,M);
 
    This function reorders the connectivity matrix by modular structure and
    may consequently be useful in visualization of modular structure.
 
    Inputs:
        W,      connectivity matrix (binary/weighted undirected/directed)
        M,      module affiliation vector
 
    Outputs:
        On,     new node order
        Wr,     reordered connectivity matrix
 
 
    Used in: Rubinov and Sporns (2011) NeuroImage; Zingg et al. (2014) Cell.
 
 
    2011, Mika Rubinov, UNSW/U Cambridge

resource_efficiency_bin
  RESOURCE_EFFICIENCY_BIN       Resource efficiency and shortest-path probability
 
    [Eres,prob_SPL] = resource_efficiency_bin(adj,lambda,SPL,M)
 
    The resource efficiency between nodes i and j is inversly proportional
    to the amount of resources (i.e. number of particles or messages)
    required to ensure with probability 0 < lambda < 1 that at least one of
    them will arrive at node j in exactly SPL steps, where SPL is the
    length of the shortest-path between i and j.
 
    The shortest-path probability between nodes i and j is the probability
    that a single random walker starting at node i will arrive at node j by
    following (one of) the shortest path(s).
 
    Inputs:
 
        adj,
            Unweighted, undirected adjacency matrix
 
        lambda,
            Probability (0 < lambda < 1)
           	set to NAN if computation of Eres is not desired
 
        SPL,
            Shortest-path length matrix (optional)
 
      	M,
            Transition probability matrix (optional)
 
 
    Outputs:
 
        Eres,
            Resource efficiency matrix.
 
        prob_SPL,
            Shortest-path probability matrix
 
    Notes:
 
        Global measures for both Eres and prob_SPL are well defined and can
        be computed as the average across the off-diagonal elements:
            GEres = mean(Eres(~eye(N)>0));
            Gprob_SPL = mean(rob_SPL(~eye(N)>0));
 
 
    Reference: Goñi J, et al (2013) PLoS ONE
 
 
    Joaquin Goñi, IU Bloomington, 2012

retrieve_shortest_path
  RETRIEVE_SHORTEST_PATH        Retrieval of shortest path
 
    This function finds the sequence of nodes that comprise the shortest
    path between a given source and target node.
 
    Inputs:
        s,
            Source node: i.e. node where the shortest path begins. 
     	t,
            Target node: i.e. node where the shortest path ends.
        hops,
            Number of edges in the path. This matrix may be obtained as the
            second output argument of the function "distance_wei_floyd.m".
        Pmat,
            Pmat is a matrix whose elements {k,t} indicate the next node in
            the shortest path between k and t. This matrix may be obtained
            as the third output of the function "distance_wei_floyd.m"
 
    Output:
        path,
            Nodes comprising the shortest path between nodes s and t.
 
 
    Andrea Avena-Koenigsberger and Joaquin Goñi, IU, 2012

rich_club_bd
 RICH_CLUB_BD        Rich club coefficients (binary directed graph)
 
    R = rich_club_bd(CIJ)
    [R,Nk,Ek] = rich_club_bd(CIJ,klevel)
 
    The rich club coefficient, R, at level k is the fraction of edges that
    connect nodes of degree k or higher out of the maximum number of edges
    that such nodes might share.
 
    Input:      CIJ,        connection matrix, binary and directed
             klevel,        optional input argument. klevel sets the
                               maximum level at which the rich club
                               coefficient will be calculated. If klevel is
                               not included the the maximum level will be
                               set to the maximum degree of CIJ.
 
    Output:       R,        vector of rich-club coefficients for levels
                               1 to klevel.
                 Nk,        number of nodes with degree>k
                 Ek,        number of edges remaining in subgraph with
                               degree>k
 
    Reference: Colizza et al. (2006) Nat. Phys. 2:110.
 
    Martijn van den Heuvel, University Medical Center Utrecht, 2011

rich_club_bu
 RICH_CLUB_BU        Rich club coefficients (binary undirected graph)
 
    R = rich_club_bu(CIJ)
    [R,Nk,Ek] = rich_club_bu(CIJ,klevel)
 
    The rich club coefficient, R, at level k is the fraction of edges that
    connect nodes of degree k or higher out of the maximum number of edges
    that such nodes might share.
 
    Input:      CIJ,        connection matrix, binary and undirected
             klevel,        optional input argument. klevel sets the
                               maximum level at which the rich club
                               coefficient will be calculated. If klevel is
                               not included the the maximum level will be
                               set to the maximum degree of CIJ.
 
    Output:       R,        vector of rich-club coefficients for levels
                               1 to klevel.
                 Nk,        number of nodes with degree>k
                 Ek,        number of edges remaining in subgraph with
                               degree>k
 
    Reference: Colizza et al. (2006) Nat. Phys. 2:110.
 
    Martijn van den Heuvel, University Medical Center Utrecht, 2011

rich_club_wd
 RICH_CLUB_WD	Rich club coefficients curve (weighted directed graph)
 
    Rw = rich_club_wd(CIJ,varargin)
 
    The weighted rich club coefficient, Rw, at level k is the fraction of
    edge weights that connect nodes of degree k or higher out of the
    maximum edge weights that such nodes might share.
 
    Inputs:
        CIJ:        weighted directed connection matrix
 
        k-level:    (optional) max level of RC(k).
                    (by default k-level quals the maximal degree of CIJ)
 
    Output:
        Rw:         rich-club curve
 
 
    References:
        T Opsahl et al. Phys Rev Lett, 2008, 101(16)
        M van den Heuvel, O Sporns, J Neurosci 2011 31(44)
 
    Martijn van den Heuvel, University Medical Center Utrecht, 2011

rich_club_wu
 RICH_CLUB_WU 	Rich club coefficients curve (weighted undirected graph)
 
    Rw = rich_club_wu(CIJ,varargin) % rich club curve for weighted graph
 
    The weighted rich club coefficient, Rw, at level k is the fraction of
    edge weights that connect nodes of degree k or higher out of the
    maximum edge weights that such nodes might share.
 
    Inputs:
        CIJ:        weighted directed connection matrix
 
        k-level:    (optional) max level of RC(k).
                    (by default k-level quals the maximal degree of CIJ)
 
    Output:
        Rw:         rich-club curve
 
 
    References:
        T Opsahl et al. Phys Rev Lett, 2008, 101(16)
        M van den Heuvel, O Sporns, J Neurosci 2011 31(44)
 
    Martijn van den Heuvel, University Medical Center Utrecht, 2011

rout_efficiency
  ROUT_EFFICIENCY       Mean, pair-wise and local routing efficiency
 
    [GErout,Erout,Eloc] = rout_efficiency(D,transform);
 
    The routing efficiency is the average of inverse shortest path length.
 
    The local routing efficiency of a node u is the routing efficiency
    computed on the subgraph formed by the neighborhood of node u
    (excluding node u).
 
 
    Inputs:
 
        D,
            Weighted/unweighted directed/undirected
            connection *weight* OR *length* matrix.
 
        transform,
            If the input matrix is a connection *weight* matrix, specify a
            transform that map input connection weights to connection
            lengths. Two transforms are available.
                'log' -> l_ij = -log(w_ij)
                'inv' -> l_ij =    1/w_ij
 
            If the input matrix is a connection *length* matrix, do not
            specify a transform (or specify an empty transform argument).
 
 
    Outputs:
 
        GErout,
            Mean global routing efficiency (scalar).
 
    	Erout,
            Pair-wise routing efficiency (matrix).
 
     	Eloc,
            Local efficiency (vector)
 
 
    Note:
 
        The input matrix may be either a connection weight matrix, or a
        connection length matrix. The connection length matrix is typically
        obtained with a mapping from weight to length, such that higher
        weights are mapped to shorter lengths (see above).
 
 
    Algorithm:  Floyd–Warshall Algorithm
 
 
    References:
        Latora and Marchiori (2001) Phys Rev Lett
        Goñi et al (2013) PLoS ONE
        Avena-Koenigsberger et al (2016) Brain Structure and Function
 
 
    Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2012
 

score_wu
 SCORE_WU       S-score
 
    [CIJscore,sn] = score_wu(CIJ,s);
 
    The s-core is the largest subnetwork comprising nodes of strength at
    least s. This function computes the s-core for a given weighted
    undirected connection matrix. Computation is analogous to the more
    widely used k-core, but is based on node strengths instead of node
    degrees. 
 
    input:          CIJ,	connection/adjacency matrix (weighted, undirected)
                      s,    level of s-core. Note: s can take on any fractional value
 
    output:    CIJscore,    connection matrix of the s-core.  This matrix 
                            contains only nodes with a strength of at least s.
                     sn,    size of s-score
 
    Olaf Sporns, Indiana University, 2007/2008/2010/2012

search_information
  SEARCH_INFORMATION                    Search information
 
    SI = search_information(adj,transform,has_memory)
 
    Computes the amount of information (measured in bits) that a random
    walker needs to follow the shortest path between a given pair of nodes.
 
    Inputs:
 
        adj,
            Weighted/unweighted directed/undirected
            connection *weight* OR *length* matrix.
 
        transform,
            If the input matrix is a connection *weight* matrix, specify a
            transform that map input connection weights to connection
            lengths. Two transforms are available.
                'log' -> l_ij = -log(w_ij)
                'inv' -> l_ij =    1/w_ij
 
            If the input matrix is a connection *length* matrix, do not
            specify a transform (or specify an empty transform argument).
 
       	has_memory,
            This flag defines whether or not the random walker "remembers"
            its previous step, which has the effect of reducing the amount
            of information needed to find the next state. If this flag is
            not set, the walker has no memory by default.
 
 
    Outputs:
 
        SI,
            pair-wise search information (matrix). Note that SI(i,j) may be
            different from SI(j,i), hense, SI is not a symmetric matrix
            even when adj is symmetric.
 
 
    References: Rosvall et al. (2005) Phys Rev Lett 94, 028701
                Goñi et al (2014) PNAS doi: 10.1073/pnas.131552911
 
 
    Andrea Avena-Koenigsberger and Joaquin Goñi, IU Bloomington, 2014

strengths_dir
 STRENGTHS_DIR      In-strength and out-strength
 
    [is,os,str] = strengths_dir(CIJ);
 
    Node strength is the sum of weights of links connected to the node. The
    instrength is the sum of inward link weights and the outstrength is the
    sum of outward link weights.
 
    Input:      CIJ,    directed weighted connection matrix
 
    Output:     is,     node instrength
                os,     node outstrength
                str,    node strength (instrength + outstrength)
 
    Notes:  Inputs are assumed to be on the columns of the CIJ matrix.
 
 
    Olaf Sporns, Indiana University, 2002/2006/2008

strengths_und
 STRENGTHS_UND        Strength
 
    str = strengths_und(CIJ);
 
    Node strength is the sum of weights of links connected to the node.
 
    Input:      CIJ,    undirected weighted connection matrix
 
    Output:     str,    node strength
 
 
    Olaf Sporns, Indiana University, 2002/2006/2008

strengths_und_sign
 STRENGTHS_UND_SIGN        Strength and weight
 
    [Spos Sneg] = strengths_und_sign(W);
    [Spos Sneg vpos vneg] = strengths_und_sign(W);
 
    Node strength is the sum of weights of links connected to the node.
 
    Inputs:     W,              undirected connection matrix with positive
                                and negative weights
 
    Output:     Spos/Sneg,      nodal strength of positive/negative weights
                vpos/vneg,      total positive/negative weight
 
 
    2011, Mika Rubinov, UNSW

subgraph_centrality
  SUBGRAPH_CENTRALITY     Subgraph centrality of a network
 
    Cs = subgraph_centrality(CIJ)
 
    The subgraph centrality of a node is a weighted sum of closed walks of
    different lengths in the network starting and ending at the node. This
    function returns a vector of subgraph centralities for each node of the
    network.
 
    Inputs:     CIJ,        adjacency matrix (binary)
 
    Outputs:     Cs,        subgraph centrality
 
    Reference: Estrada and Rodriguez-Velasquez (2005) Phys Rev E 71, 056103
               Estrada and Higham (2010) SIAM Rev 52, 696.
 
    Xi-Nian Zuo, Chinese Academy of Sciences, 2010
    Rick Betzel, Indiana University, 2012

threshold_absolute
  THRESHOLD_ABSOLUTE    Absolute thresholding
  
    W_thr = threshold_absolute(W, thr);
 
    This function thresholds the connectivity matrix by absolute weight
    magnitude. All weights below the given threshold, and all weights
    on the main diagonal (self-self connections) are set to 0.
 
    Inputs: W           weighted or binary connectivity matrix
            thr         weight treshold
 
    Output: W_thr       thresholded connectivity matrix
 
 
    Mika Rubinov, UNSW, 2009-2010

threshold_proportional
 THRESHOLD_PROPORTIONAL     Proportional thresholding
 
    W_thr = threshold_proportional(W, p);
 
    This function "thresholds" the connectivity matrix by preserving a
    proportion p (0<p<1) of the strongest weights. All other weights, and
    all weights on the main diagonal (self-self connections) are set to 0.
 
    Inputs: W,      weighted or binary connectivity matrix
            p,      proportion of weights to preserve
                        range:  p=1 (all weights preserved) to
                                p=0 (no weights preserved)
 
    Output: W_thr,  thresholded connectivity matrix
 
 
    Mika Rubinov, U Cambridge,
    Roan LaPlante, Martinos Center, MGH
    Zitong Zhang, Penn Engineering

transitivity_bd
 TRANSITIVITY_BD    Transitivity
 
    T = transitivity_bd(A);
 
    Transitivity is the ratio of 'triangles to triplets' in the network.
    (A classical version of the clustering coefficient).
 
    Input:      A       binary directed connection matrix
 
    Output:     T       transitivity scalar
 
    Reference:  Rubinov M, Sporns O (2010) NeuroImage 52:1059-69
                based on Fagiolo (2007) Phys Rev E 76:026107.
 
 
    Contributors:
    Mika Rubinov, UNSW/University of Cambridge
    Christoph Schmidt, Friedrich Schiller University Jena
    Andrew Zalesky, University of Melbourne
    2007-2015

transitivity_bu
 TRANSITIVITY_BU    Transitivity
 
    T = transitivity_bu(A);
 
    Transitivity is the ratio of 'triangles to triplets' in the network.
    (A classical version of the clustering coefficient).
 
    Input:      A       binary undirected connection matrix
 
    Output:     T       transitivity scalar
 
    Reference: e.g. Humphries et al. (2008) Plos ONE 3: e0002051
 
 
    Alexandros Goulas, Maastricht University, 2010

transitivity_wd
 TRANSITIVITY_WD    Transitivity
 
    T = transitivity_wd(W);
 
    Transitivity is the ratio of 'triangles to triplets' in the network.
    (A classical version of the clustering coefficient).
 
    Input:      W       weighted directed connection matrix
 
    Output:     T       transitivity scalar
 
    Note:       All weights must be between 0 and 1.
                This may be achieved using the weight_conversion.m function,
                W_nrm = weight_conversion(W, 'normalize');
 
    Reference:  Rubinov M, Sporns O (2010) NeuroImage 52:1059-69
                based on Fagiolo (2007) Phys Rev E 76:026107.
 
 
    Contributors:
    Mika Rubinov, UNSW/University of Cambridge
    Christoph Schmidt, Friedrich Schiller University Jena
    Andrew Zalesky, University of Melbourne
    2007-2015

transitivity_wu
 TRANSITIVITY_WU    Transitivity
 
    T = transitivity_wu(W);
 
    Transitivity is the ratio of 'triangles to triplets' in the network.
    (A classical version of the clustering coefficient).
 
    Input:      W       weighted undirected connection matrix
 
    Output:     T       transitivity scalar
 
    Note:      All weights must be between 0 and 1.
               This may be achieved using the weight_conversion.m function,
               W_nrm = weight_conversion(W, 'normalize');
 
    Reference: Rubinov M, Sporns O (2010) NeuroImage 52:1059-69
               based on Onnela et al. (2005) Phys Rev E 71:065103
 
 
    Mika Rubinov, UNSW/U Cambridge, 2010-2015

weight_conversion
  WEIGHT_CONVERSION    Conversion of weights in input matrix
 
    W_bin = weight_conversion(W, 'binarize');
    W_nrm = weight_conversion(W, 'normalize');
    L = weight_conversion(W, 'lengths');
    W_fix = weight_conversion(W, 'autofix');
 
    This function may either binarize an input weighted connection matrix,
    normalize an input weighted connection matrix, convert an input
    weighted connection matrix to a weighted connection-length matrix, or
    fix common connection problems in binary or weighted connection matrices.
 
        Binarization converts all present connection weights to 1.
 
        Normalization rescales all weight magnitudes to the range [0,1] and
    should be done prior to computing some weighted measures, such as the
    weighted clustering coefficient.
 
        Conversion of connection weights to connection lengths is needed
    prior to computation of weighted distance-based measures, such as
    distance and betweenness centrality. In a weighted connection network,
    higher weights are naturally interpreted as shorter lengths. The
    connection-lengths matrix here is defined as the inverse of the
    connection-weights matrix.
 
        Autofix removes all Inf and NaN values, remove all self connections
    (sets all weights on the main diagonal to 0), ensures that symmetric matrices
    are exactly symmetric (by correcting for round-off error), and ensures that
    binary matrices are exactly binary (by correcting for round-off error).
 
    Inputs: W           binary or weighted connectivity matrix
            wcm         weight-conversion command - possible values:
                            'binarize'      binarize weights
                            'normalize'     normalize weights

writetoPAJ
 WRITETOPAJ         Write to Pajek
 
    writetoPAJ(CIJ, fname, arcs);
 
    This function writes a Pajek .net file from a MATLAB matrix
 
    Inputs:     CIJ,        adjacency matrix
                fname,      filename minus .net extension
                arcs,       1 for directed network
                            0 for an undirected network
 
    Chris Honey, Indiana University, 2007