From 45a46d6d338566937ab776951ab30c0a03a9d3a4 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Wed, 7 Dec 2022 12:45:11 -0500
Subject: [PATCH 01/11] create models file

---
 bct/algorithms/models.py | 7 +++++++
 1 file changed, 7 insertions(+)
 create mode 100644 bct/algorithms/models.py

diff --git a/bct/algorithms/models.py b/bct/algorithms/models.py
new file mode 100644
index 0000000..bcb64f7
--- /dev/null
+++ b/bct/algorithms/models.py
@@ -0,0 +1,7 @@
+from __future__ import division, print_function
+import numpy as np
+
+from bct.utils import BCTParamError, get_rng
+from ..due import due, BibTex
+
+

From 5e9e49faea44ef5f14b2fb16705f4a1939c5331f Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Wed, 7 Dec 2022 13:24:22 -0500
Subject: [PATCH 02/11] create function stub for mleme_constraint_model

---
 bct/algorithms/generative.py |  6 ++++-
 bct/algorithms/models.py     | 50 ++++++++++++++++++++++++++++++++++++
 2 files changed, 55 insertions(+), 1 deletion(-)

diff --git a/bct/algorithms/generative.py b/bct/algorithms/generative.py
index 3b40c1b..ad92ebc 100644
--- a/bct/algorithms/generative.py
+++ b/bct/algorithms/generative.py
@@ -598,7 +598,8 @@ def kstats(x, y):
 
     return np.max(K, axis=1)
 
-def generate_fc(sc, beta, ed=None, pred_var=(), model='linear'):
+def generate_fc(sc, beta, ed=None, pred_var=(), model='linear', seed=None):
+
     '''
     Uses a vector beta of regression coefficients from the model
     FC = pred_var*beta to predict functional connectivity.
@@ -641,6 +642,9 @@ def generate_fc(sc, beta, ed=None, pred_var=(), model='linear'):
             'quadratic'. Depends on some matlab calls and 
             #TODO figure out what implementations exist
             The default value is 'linear'
+    seed : hashable, optional
+        If None (default), use the np.random's global random state to generate random numbers.
+        Otherwise, use a new np.random.RandomState instance seeded with the given value.
 
     In the matlab function there is an additional parameter for inserting a
     real functional connectivity matrix. If supplied, the function calculates
diff --git a/bct/algorithms/models.py b/bct/algorithms/models.py
index bcb64f7..7827ac8 100644
--- a/bct/algorithms/models.py
+++ b/bct/algorithms/models.py
@@ -4,4 +4,54 @@
 from bct.utils import BCTParamError, get_rng
 from ..due import due, BibTex
 
+def mleme_constraint_model(nr_samples, W, ci=None, lo=None, li=None, lm=None, 
+                           seed=None):
+    '''
+    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.
 
+    Parameters
+    ----------
+    W : np.ndarray
+        NxN square directed weighted connectivity matrix. All inputs must be
+        nonnegative integers. Real valued weights could be converted to
+        integers through rescaling and rounding.
+    ci : np.ndarray
+        Nx1 module affiliation vector. Can be None if there are no module
+        constraints. Must contain nonnegative integers. The default value
+        is None.
+    lo : np.ndarray
+        Nx1 out strength constraing logical vector. This vector specifies
+        out strength constraints for each node. Alternately, it can be
+        True to constrain all out-strengths or None for no constraints.
+        The default value is None.
+    li : np.ndarray
+        Nx1 in strength constraing logical vector. This vector specifies
+        in strength constraints for each node. Alternately, it can be
+        True to constrain all in-strengths or None for no constraints.
+        The default value is None.
+    lm : np.ndarray
+        Mx1 module-weight constraint logical matrix where M is the number of
+        modules. Specifies module-weight constraints for all pairs of modules.
+        Can be True, 'all', or 2, to constrain all inter-module and 
+        intra-module weights, 'intra' or 1 to constrain all intra-module
+        weights only, or None for no constraints. The default value is None.
+    seed : hashable, optional
+        If None (default), use the np.random's global random state to generate random numbers.
+        Otherwise, use a new np.random.RandomState instance seeded with the given value.
+
+    Returns
+    -------
+    W0 : np.ndarray
+        NxNxnr_samples an ensemble of sampled networks with constraints
+    E0 : np.ndarray
+        expected weights matrix
+    P0 : np.ndarray
+        probability matrix
+    delt0 : float
+        algorithm convergence error
+    '''
+    rng = get_rng(seed)
+    raise NotImplementedError()

From 779e350b1ee920c6874420725bbb61f32132deac Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Wed, 7 Dec 2022 13:52:23 -0500
Subject: [PATCH 03/11] implement diffusion_efficiency()

---
 bct/algorithms/distance.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/bct/algorithms/distance.py b/bct/algorithms/distance.py
index 32b4651..1456b9f 100644
--- a/bct/algorithms/distance.py
+++ b/bct/algorithms/distance.py
@@ -1091,3 +1091,31 @@ def mean_first_passage_time(adjacency):
     mfpt = (np.repeat(np.atleast_2d(np.diag(Z)), n, 0) - Z) / W
 
     return mfpt
+
+def 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.
+
+    Parameters
+    ----------
+    adj : np.ndarray
+        weighted/unweighted, directed/undirected adjacency matrix
+
+    Returns
+    -------
+    gediff : float
+        mean global diffusion efficiency
+    ediff : np.ndarray
+        pairwise NxN diffusion efficiency matrix
+    '''
+    n = len(adj)
+    mfpt = mean_first_passage_time(adj)
+    ediff = 1 / mpft
+    np.fill_diagonal(ediff, 0)
+    gediff = np.sum(ediff) / (n ** 2 - n)
+    return gediff, ediff

From 79c58260a160aae889628166b9648668b38bb555 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Wed, 7 Dec 2022 14:40:10 -0500
Subject: [PATCH 04/11] add test case for diffusion_efficiency, added part of
 code for clique_communities not currently working or formatted at all

---
 bct/algorithms/core.py     | 70 ++++++++++++++++++++++++++++++++++++++
 bct/algorithms/distance.py |  3 +-
 test/distance_test.py      |  7 ++++
 3 files changed, 79 insertions(+), 1 deletion(-)

diff --git a/bct/algorithms/core.py b/bct/algorithms/core.py
index 2bba5a0..8814936 100644
--- a/bct/algorithms/core.py
+++ b/bct/algorithms/core.py
@@ -653,3 +653,73 @@ def score_wu(CIJ, s):
 
     sn = np.sum(str > 0)
     return CIJscore, sn
+
+def 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.
+
+    Parameters
+    ----------
+    A : np.ndaray
+        binary undirected connection matrix
+    cq_thr : int
+        Clique size threshold. Larger clique size thresholds potentially result
+        in larger communities
+
+    Returns
+    -------
+    M : np.ndarray
+        MxN Overlapping community affiliation matrix, each node can participate
+        in arbitrarily many of the M communities
+    
+    Note: This algorithm can be slow and memory intensive in large matrices.
+    Requires get_components
+    '''
+    if not np.all(A == A.T):
+        raise BCTParamError('Input must be undirected')
+    elif not a.ndim == 2:
+        raise BCTParamError('Input must be 2 dimensional NxN matrix')
+    elif a.shape[0] != a.shape[1]:
+        raise BCTParamError('Input must be square')
+
+    def maximal_cliques(A, n):
+        #Bron-Kerbosch algorithm
+        q = 0
+        R = np.zeros((n, 1))        #current
+        P = np.ones((n, 1))         #prospective
+        X = np.zeros((n, 1))        #processed
+
+        MQ = np.zeros((n, 1000*n))
+
+        def bk(R, P, X):
+            if not (np.sum(P) or np.sum(X)):
+                MQ[:, q] = R
+                q += 1
+            else:
+                U_p = np.where()
+
+    n = len(a)
+    A = A.copy()
+    np.fill_diagonal(A, 0)
+    cq = maximal_cliques(A, n)
+    #remove subthreshold cliques
+    cq = mq[np.where(np.sum(mq, axis=1) >= cq_thr), :]
+    #compute clique overlap
+    ov = np.matmul(cq, cq.T)
+    #keep percolating cliques
+    ov_thr = (ov >= cq_thr - 1)
+
+    cq_components = get_components(ov_thr)
+    #get number of components
+    nr_comopnents = np.max(cq_components)
+    M = np.zeros((nr_components, n))
+    for i in range(nr_components):
+        M[i, np.where(cq[cq_components==i, :])] = 1
+
+    return M
+    
diff --git a/bct/algorithms/distance.py b/bct/algorithms/distance.py
index 1456b9f..f1423ae 100644
--- a/bct/algorithms/distance.py
+++ b/bct/algorithms/distance.py
@@ -1114,8 +1114,9 @@ def diffusion_efficiency(adj):
         pairwise NxN diffusion efficiency matrix
     '''
     n = len(adj)
+    adj = adj.copy()
     mfpt = mean_first_passage_time(adj)
-    ediff = 1 / mpft
+    ediff = 1 / mfpt
     np.fill_diagonal(ediff, 0)
     gediff = np.sum(ediff) / (n ** 2 - n)
     return gediff, ediff
diff --git a/test/distance_test.py b/test/distance_test.py
index 033d416..6ffb155 100644
--- a/test/distance_test.py
+++ b/test/distance_test.py
@@ -54,3 +54,10 @@ def test_charpath():
     assert np.any(np.isinf(d))
     assert not np.isnan(radius)
     assert not np.isnan(diameter)
+
+def test_diffusion_efficiency():
+    x = load_sample(thres=.23)
+    gde, ed = bct.diffusion_efficiency(x)
+    print(gde, np.sum(ed)) 
+    assert np.allclose(gde, .0069472)
+    assert np.allclose(np.sum(ed), 131.34, atol=.01)

From f22f27263ecd75b9693ce46145ded9b083781e63 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Wed, 7 Dec 2022 17:22:20 -0500
Subject: [PATCH 05/11] correct code for clique_communities() and add test case

---
 bct/algorithms/core.py | 63 +++++++++++++++++++++++++++++-------------
 test/core_test.py      | 11 ++++++++
 2 files changed, 55 insertions(+), 19 deletions(-)

diff --git a/bct/algorithms/core.py b/bct/algorithms/core.py
index 8814936..8408c04 100644
--- a/bct/algorithms/core.py
+++ b/bct/algorithms/core.py
@@ -2,8 +2,10 @@
 import numpy as np
 
 from ..utils.miscellaneous_utilities import get_rng, BCTParamError
+from ..utils.other import binarize
 from .degree import degrees_dir, degrees_und, strengths_dir, strengths_und
 from .degree import strengths_und_sign
+from .clustering import get_components
 
 from ..due import due, BibTeX
 from ..citations import NEWMAN2002, FOSTER2010, HAGMANN2008, COLIZZA2006, OPSAHL2008, HEUVEL2011
@@ -682,44 +684,67 @@ def clique_communities(A, cq_thr):
     '''
     if not np.all(A == A.T):
         raise BCTParamError('Input must be undirected')
-    elif not a.ndim == 2:
+    elif not A.ndim == 2:
         raise BCTParamError('Input must be 2 dimensional NxN matrix')
-    elif a.shape[0] != a.shape[1]:
+    elif A.shape[0] != A.shape[1]:
         raise BCTParamError('Input must be square')
 
-    def maximal_cliques(A, n):
+    n = len(A)
+    A = binarize(A, copy=True)
+    np.fill_diagonal(A, 0)
+
+    def maximal_cliques(A):
         #Bron-Kerbosch algorithm
-        q = 0
         R = np.zeros((n, 1))        #current
         P = np.ones((n, 1))         #prospective
         X = np.zeros((n, 1))        #processed
 
-        MQ = np.zeros((n, 1000*n))
+        MQ = []
 
-        def bk(R, P, X):
+        def bk(R, P, X, nrec):
+            #print(nrec)
+            nrec += 1
             if not (np.sum(P) or np.sum(X)):
-                MQ[:, q] = R
-                q += 1
+                MQ.append(R)
             else:
-                U_p = np.where()
+                U_p, = np.where(np.any(np.hstack((P, X)), axis=1))
+                ix = np.argmax(np.dot(A[:, U_p].T, P))
+                u_p = U_p[ix]
 
-    n = len(a)
-    A = A.copy()
-    np.fill_diagonal(A, 0)
-    cq = maximal_cliques(A, n)
+                Aup = np.reshape(np.logical_not(A[:,u_p]), (n, 1))
+                U, = np.where(np.all(np.hstack((P.astype(bool), Aup)), axis=1))
+                for u in U:
+                    Nu = np.reshape(A[:, u], (n, 1))
+                    P[u] = 0
+                    Rnew = R.copy()
+                    Rnew[u] = 1
+
+                    Pnew = np.reshape(np.all(np.hstack((P, Nu)), axis=1), (n, 1))
+                    Xnew = np.reshape(np.all(np.hstack((X, Nu)), axis=1), (n, 1))
+
+                    bk(Rnew, Pnew, Xnew, nrec)
+                    X[u] = 1
+
+        bk(R, P, X, 0)
+        return np.squeeze(MQ)
+
+    cq = maximal_cliques(A)
+    print(np.shape(cq))
     #remove subthreshold cliques
-    cq = mq[np.where(np.sum(mq, axis=1) >= cq_thr), :]
+    ix, = np.where(np.sum(cq, axis=1) >= cq_thr)
+    cq = cq[ix, :]
     #compute clique overlap
-    ov = np.matmul(cq, cq.T)
+    ov = np.dot(cq, cq.T)
+    print(ov)
     #keep percolating cliques
     ov_thr = (ov >= cq_thr - 1)
+    print(ov_thr.shape)
 
-    cq_components = get_components(ov_thr)
+    cq_components, _ = get_components(ov_thr)
     #get number of components
-    nr_comopnents = np.max(cq_components)
+    nr_components = np.max(cq_components)
     M = np.zeros((nr_components, n))
     for i in range(nr_components):
-        M[i, np.where(cq[cq_components==i, :])] = 1
+        M[i, np.where(np.any(cq[cq_components==i+1, :], axis=0))] = 1
 
     return M
-    
diff --git a/test/core_test.py b/test/core_test.py
index 8151785..b5a0930 100644
--- a/test/core_test.py
+++ b/test/core_test.py
@@ -17,3 +17,14 @@ def test_core_periphery_dir():
     assert np.sum(c) == 57
     assert np.sum(np.cumsum(c)) == 4170
     assert np.allclose(q, .3086, atol=.0001)
+
+def test_clique_communities():
+    x = load_sample(thres=.23)
+
+    print(np.sum(bct.binarize(x)))
+
+    cis = bct.clique_communities(x, 9)
+    print(cis.shape, np.max(np.sum(cis, axis=0)))
+    print(np.sum(cis, axis=1))
+    assert np.sum(cis) == 199
+    assert 4 == 8

From a8ade7e5be3b806ecf3e4f478554d37f58cca010 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Sun, 18 Dec 2022 23:54:20 -0500
Subject: [PATCH 06/11] add path_transitivity() and test case

---
 bct/algorithms/clustering.py | 65 +++++++++++++++++++++++++++++++++++-
 bct/algorithms/distance.py   | 39 +++++++++++-----------
 bct/utils/other.py           | 27 +++++++++++++++
 test/clustering_test.py      | 10 ++++++
 test/distance_test.py        |  8 +++++
 5 files changed, 128 insertions(+), 21 deletions(-)

diff --git a/bct/algorithms/clustering.py b/bct/algorithms/clustering.py
index 484cc23..9b06b03 100644
--- a/bct/algorithms/clustering.py
+++ b/bct/algorithms/clustering.py
@@ -2,7 +2,7 @@
 import numpy as np
 from .modularity import modularity_louvain_und_sign
 from bct.utils import cuberoot, BCTParamError, dummyvar, binarize, get_rng
-from .distance import breadthdist
+from .distance import breadthdist, distance_wei_floyd, retrieve_shortest_path
 from ..due import due, BibTeX
 from ..citations import (
     WATTS1998, ONNELA2005, FAGIOLO2007, LANCICHINETTI2012,
@@ -719,3 +719,66 @@ def transitivity_wu(W):
     ws = cuberoot(W)
     cyc3 = np.diag(np.dot(ws, np.dot(ws, ws)))
     return np.sum(cyc3, axis=0) / np.sum(K * (K - 1), axis=0)
+
+def path_transitivity(W, transform=None):
+    '''
+    This function computes the density of local detours (triangles) that
+    are available along the shortest-paths between all pairs of nodes.
+
+    Parameters
+    ----------
+    W : NxN np.ndarray
+        weighted or unweighted undirected connection weight or length matrix
+    transform : None or enum
+        if the input is a connection length matrix, no transform is needed
+        if the input is a connection weight matrix, use 'log' for
+        log transform l_ij = -log(w_ij)
+        or 'inv' for inversion l_ij = 1/w_ij
+        The default value is None
+
+    Returns
+    -------
+    T : NxN
+        matrix of pairwise path transitivity
+    '''
+    n = len(W)
+    m = np.zeros((n, n))
+    T = np.zeros((n, n))
+
+    for i in range(n-1):
+        for j in range(i+1, n):
+            x = 0
+            y = 0
+            z = 0
+            
+            for k in range(n):
+                if W[i, k] != 0 and W[j, k] != 0 and k not in (i, j):
+                    x += W[i, k] + W[j, k]
+                if k != j:
+                    y += W[i, k]
+                if k != i:
+                    z += W[j, k]
+
+            m[i,j] = x/(y+z)
+
+    m = m + m.T
+
+    _, hops, pmat = distance_wei_floyd(W, transform=transform)
+
+    for i in range(n-1):
+        for j in range(i+1, n):
+
+            x = 0
+            path = retrieve_shortest_path(i, j, hops, pmat)
+            k = len(path)
+            print(path)
+            print(k)
+
+            for t in range(k-1):
+                for l in range(t+1, k):
+                    x += m[path[t], path[l]]
+
+            T[i, j] = 2 * x / (k * (k - 1))
+
+    T = T + T.T
+    return T
diff --git a/bct/algorithms/distance.py b/bct/algorithms/distance.py
index f1423ae..12d097e 100644
--- a/bct/algorithms/distance.py
+++ b/bct/algorithms/distance.py
@@ -1,6 +1,6 @@
 from __future__ import division, print_function
 import numpy as np
-from bct.utils import cuberoot, binarize, invert, BCTParamError
+from bct.utils import cuberoot, binarize, invert, BCTParamError, logtransform
 from ..due import due, BibTeX
 from ..citations import LATORA2001, ONNELA2005, FAGIOLO2007, RUBINOV2010
 
@@ -381,44 +381,43 @@ def distance_wei_floyd(adjacency, transform=None):
     .. [4] https://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm
     """
 
+    #it is important not to do these transformations safely, to allow infinity
     if transform is not None:
-        if transform == 'log':
-            if np.logical_or(adjacency > 1, adjacency < 0).any():
-                raise ValueError("Connection strengths must be in the " +
-                                 "interval [0,1) to use the transform " +
-                                 "-log(w_ij).")
-            SPL = -np.log(adjacency)
-        elif transform == 'inv':
-            SPL = 1. / adjacency
-        else:
-            raise ValueError("Unexpected transform type. Only 'log' and " +
-                             "'inv' are accepted")
+        with np.errstate(divide='ignore'):
+            if transform == 'log':
+                #SPL = logtransform(adjacency)
+                SPL = -np.log(adjacency)
+            elif transform == 'inv':
+                #SPL = invert(adjacency)
+                SPL = 1 / adjacency
+            else:
+                raise ValueError("Unexpected transform type. Only 'log' and " +
+                                 "'inv' are accepted")
     else:
         SPL = adjacency.copy().astype('float')
         SPL[SPL == 0] = np.inf
 
     n = adjacency.shape[1]
 
-    flag_find_paths = True
     hops = np.array(adjacency != 0).astype('float')
     Pmat = np.repeat(np.atleast_2d(np.arange(0, n)), n, 0)
 
+    #print(SPL)
+
     for k in range(n):
         i2k_k2j = np.repeat(SPL[:, [k]], n, 1) + np.repeat(SPL[[k], :], n, 0)
 
-        if flag_find_paths:
-            path = SPL > i2k_k2j
-            i, j = np.where(path)
-            hops[path] = hops[i, k] + hops[k, j]
-            Pmat[path] = Pmat[i, k]
+        path = SPL > i2k_k2j
+        i, j = np.where(path)
+        hops[path] = hops[i, k] + hops[k, j]
+        Pmat[path] = Pmat[i, k]
 
         SPL = np.min(np.stack([SPL, i2k_k2j], 2), 2)
 
     I = np.eye(n) > 0
     SPL[I] = 0
 
-    if flag_find_paths:
-        hops[I], Pmat[I] = 0, 0
+    hops[I], Pmat[I] = 0, 0
 
     return SPL, hops, Pmat
 
diff --git a/bct/utils/other.py b/bct/utils/other.py
index 7542bda..44704cf 100644
--- a/bct/utils/other.py
+++ b/bct/utils/other.py
@@ -241,6 +241,33 @@ def invert(W, copy=True):
     W[E] = 1. / W[E]
     return W
 
+def logtransform(W, copy=True):
+    '''
+    Makes a log transformation of the weights of an input matrix, such that each
+    value W[i,j] will be -log(W[i,l])
+
+    If copy is not set, this function will *modify W in place.*
+
+    Parameters
+    ----------
+    W : np.ndarray
+        weighted connectivity matrix
+    copy : bool
+        if True, returns a copy of the matrix. Otherwise, modifies the matrix
+        in place. Default value=True
+
+    Returns
+    -------
+    W : np.ndarray
+        Log transformed connectivity matrix
+    '''
+    if copy:
+        W = W.copy()
+    if np.logical_or(W > 1, W <= 0).any():
+        raise ValueError("Connection strengths must be in the interval (0,1] "
+                         "to use the transform -log(w_ij)")
+    W = -np.log(W)
+    return W
 
 def autofix(W, copy=True):
     '''
diff --git a/test/clustering_test.py b/test/clustering_test.py
index 211c692..26dda21 100644
--- a/test/clustering_test.py
+++ b/test/clustering_test.py
@@ -146,3 +146,13 @@ def test_agreement():
         output_2_buffered = bct.agreement(input_2, buffsz=buffsz)
         print(output_2_buffered)
         assert (output_2_buffered == correct_2).all()
+
+def test_path_transitivity():
+    x = load_sample(thres=.31)
+    xn = bct.normalize(x)
+    Ti = bct.path_transitivity(x, transform='inv')
+    Tl = bct.path_transitivity(xn, transform='log')
+
+    print(np.sum(Ti), np.sum(Tl))
+    assert np.allclose(np.sum(Ti), 8340.8, atol=.01)
+    assert np.allclose(np.sum(Tl), 8621.1, atol=.01)
diff --git a/test/distance_test.py b/test/distance_test.py
index 6ffb155..1fb965e 100644
--- a/test/distance_test.py
+++ b/test/distance_test.py
@@ -61,3 +61,11 @@ def test_diffusion_efficiency():
     print(gde, np.sum(ed)) 
     assert np.allclose(gde, .0069472)
     assert np.allclose(np.sum(ed), 131.34, atol=.01)
+
+def test_distance_floyd():
+    x = load_sample(thres=.31)
+    print('guomish')
+    spli, hopsi, pmati = bct.distance_wei_floyd(x, transform='inv')
+    print(np.sum(spli))
+    assert np.allclose(np.sum(spli), 11536.1, atol=.01)
+    

From d2c03a157b2513e6a7641b2805da3b27431ff0a4 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Tue, 20 Dec 2022 00:16:24 -0500
Subject: [PATCH 07/11] add resource_efficiency, move some other functions to
 new file efficiency.py and create test cases

---
 bct/algorithms/__init__.py   |   1 +
 bct/algorithms/distance.py   | 241 +----------------------
 bct/algorithms/efficiency.py | 361 +++++++++++++++++++++++++++++++++++
 test/distance_test.py        |   7 -
 4 files changed, 363 insertions(+), 247 deletions(-)
 create mode 100644 bct/algorithms/efficiency.py

diff --git a/bct/algorithms/__init__.py b/bct/algorithms/__init__.py
index dee6cb5..9879bfa 100644
--- a/bct/algorithms/__init__.py
+++ b/bct/algorithms/__init__.py
@@ -1,5 +1,6 @@
 from .centrality import *
 from .clustering import *
+from .efficiency import *
 from .generative import *
 from .modularity import *
 from .core import *
diff --git a/bct/algorithms/distance.py b/bct/algorithms/distance.py
index 12d097e..781b7dc 100644
--- a/bct/algorithms/distance.py
+++ b/bct/algorithms/distance.py
@@ -1,6 +1,6 @@
 from __future__ import division, print_function
 import numpy as np
-from bct.utils import cuberoot, binarize, invert, BCTParamError, logtransform
+from bct.utils import binarize, BCTParamError, invert, logtransform
 from ..due import due, BibTeX
 from ..citations import LATORA2001, ONNELA2005, FAGIOLO2007, RUBINOV2010
 
@@ -467,217 +467,6 @@ def retrieve_shortest_path(s, t, hops, Pmat):
     return path
 
 
-@due.dcite(BibTeX(LATORA2001), description="Unweighted global efficiency")
-@due.dcite(BibTeX(ONNELA2005), description="Unweighted global efficiency")
-@due.dcite(BibTeX(FAGIOLO2007), description="Unweighted global efficiency")
-@due.dcite(BibTeX(RUBINOV2010), description="Unweighted global efficiency")
-def efficiency_bin(G, local=False):
-    '''
-    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.
-
-    Parameters
-    ----------
-    A : NxN np.ndarray
-        binary undirected connection matrix
-    local : bool
-        If True, computes local efficiency instead of global efficiency.
-        Default value = False.
-
-    Returns
-    -------
-    Eglob : float
-        global efficiency, only if local=False
-    Eloc : Nx1 np.ndarray
-        local efficiency, only if local=True
-    '''
-    def distance_inv(g):
-        D = np.eye(len(g))
-        n = 1
-        nPATH = g.copy()
-        L = (nPATH != 0)
-
-        while np.any(L):
-            D += n * L
-            n += 1
-            nPATH = np.dot(nPATH, g)
-            L = (nPATH != 0) * (D == 0)
-        D[np.logical_not(D)] = np.inf
-        D = 1 / D
-        np.fill_diagonal(D, 0)
-        return D
-
-    G = binarize(G)
-    n = len(G)  # number of nodes
-    if local:
-        E = np.zeros((n,))  # local efficiency
-
-        for u in range(n):
-            # V,=np.where(G[u,:])			#neighbors
-            # k=len(V)					#degree
-            # if k>=2:					#degree must be at least 2
-            #	e=distance_inv(G[V].T[V])
-            #	E[u]=np.sum(e)/(k*k-k)	#local efficiency computation
-
-            # find pairs of neighbors
-            V, = np.where(np.logical_or(G[u, :], G[u, :].T))
-            # inverse distance matrix
-            e = distance_inv(G[np.ix_(V, V)])
-            # symmetrized inverse distance matrix
-            se = e + e.T
-
-            # symmetrized adjacency vector
-            sa = G[u, V] + G[V, u].T
-            numer = np.sum(np.outer(sa.T, sa) * se) / 2
-            if numer != 0:
-                denom = np.sum(sa)**2 - np.sum(sa * sa)
-                E[u] = numer / denom  # local efficiency
-
-    else:
-        e = distance_inv(G)
-        E = np.sum(e) / (n * n - n)  # global efficiency
-    return E
-
-
-@due.dcite(BibTeX(LATORA2001), description="Weighted global efficiency")
-@due.dcite(BibTeX(ONNELA2005), description="Weighted global efficiency")
-@due.dcite(BibTeX(FAGIOLO2007), description="Weighted global efficiency")
-@due.dcite(BibTeX(RUBINOV2010), description="Weighted global efficiency")
-def efficiency_wei(Gw, local=False):
-    '''
-    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.
-
-    Parameters
-    ----------
-    W : NxN np.ndarray
-        undirected weighted connection matrix
-        (all weights in W must be between 0 and 1)
-    local = bool or enum
-        If True or 'local', computes local efficiency instead of global efficiency.
-        If False or 'global', uses the global efficiency
-        If 'original', will use the original algorithm provided by (Rubinov
-        & Sporns 2010). This version is not recommended. The local efficiency
-        calculation was improved in (Wang et al. 2016) as a true generalization
-        of the binary variant.
-        
-    Returns
-    -------
-    Eglob : float
-        global efficiency, only if local in (False, 'global')
-    Eloc : Nx1 np.ndarray
-        local efficiency, only if local in (True, 'local', 'original')
-
-    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 this weighted variant of the local
-    efficiency is hence not a strict generalization of the binary variant.
-
-    Algorithm:  Dijkstra's algorithm
-    '''
-    if local not in (True, False, 'local', 'global', 'original'):
-        raise BCTParamError("local param must be any of True, False, "
-            "'local', 'global', or 'original'")
-
-    def distance_inv_wei(G):
-        n = len(G)
-        D = np.zeros((n, n))  # distance matrix
-        D[np.logical_not(np.eye(n))] = np.inf
-
-        for u in range(n):
-            # distance permanence (true is temporary)
-            S = np.ones((n,), dtype=bool)
-            G1 = G.copy()
-            V = [u]
-            while True:
-                S[V] = 0  # distance u->V is now permanent
-                G1[:, V] = 0  # no in-edges as already shortest
-                for v in V:
-                    W, = np.where(G1[v, :])  # neighbors of smallest nodes
-                    td = np.array(
-                        [D[u, W].flatten(), (D[u, v] + G1[v, W]).flatten()])
-                    D[u, W] = np.min(td, axis=0)
-
-                if D[u, S].size == 0:  # all nodes reached
-                    break
-                minD = np.min(D[u, S])
-                if np.isinf(minD):  # some nodes cannot be reached
-                    break
-                V, = np.where(D[u, :] == minD)
-
-        np.fill_diagonal(D, 1)
-        D = 1 / D
-        np.fill_diagonal(D, 0)
-        return D
-
-    n = len(Gw)
-    Gl = invert(Gw, copy=True)  # connection length matrix
-    A = np.array((Gw != 0), dtype=int)
-    #local efficiency algorithm described by Rubinov and Sporns 2010, not recommended
-    if local == 'original':
-        E = np.zeros((n,))
-        for u in range(n):
-            # V,=np.where(Gw[u,:])		#neighbors
-            # k=len(V)					#degree
-            # if k>=2:					#degree must be at least 2
-            #	e=(distance_inv_wei(Gl[V].T[V])*np.outer(Gw[V,u],Gw[u,V]))**1/3
-            #	E[u]=np.sum(e)/(k*k-k)
-
-            # find pairs of neighbors
-            V, = np.where(np.logical_or(Gw[u, :], Gw[:, u].T))
-            # symmetrized vector of weights
-            sw = cuberoot(Gw[u, V]) + cuberoot(Gw[V, u].T)
-            # inverse distance matrix
-            e = distance_inv_wei(Gl[np.ix_(V, V)])
-            # symmetrized inverse distance matrix
-            se = cuberoot(e) + cuberoot(e.T)
-
-            numer = np.sum(np.outer(sw.T, sw) * se) / 2
-            if numer != 0:
-                # symmetrized adjacency vector
-                sa = A[u, V] + A[V, u].T
-                denom = np.sum(sa)**2 - np.sum(sa * sa)
-                # print numer,denom
-                E[u] = numer / denom  # local efficiency
-
-    #local efficiency algorithm described by Wang et al 2016, recommended
-    elif local in (True, 'local'):
-        E = np.zeros((n,))
-        for u in range(n):
-            V, = np.where(np.logical_or(Gw[u, :], Gw[:, u].T))
-            sw = cuberoot(Gw[u, V]) + cuberoot(Gw[V, u].T)
-            e = distance_inv_wei(cuberoot(Gl)[np.ix_(V, V)])
-            se = e+e.T
-         
-            numer = np.sum(np.outer(sw.T, sw) * se) / 2
-            if numer != 0:
-                # symmetrized adjacency vector
-                sa = A[u, V] + A[V, u].T
-                denom = np.sum(sa)**2 - np.sum(sa * sa)
-                # print numer,denom
-                E[u] = numer / denom  # local efficiency
-
-    elif local in (False, 'global'):
-        e = distance_inv_wei(Gl)
-        E = np.sum(e) / (n * n - n)
-    return E
-
 
 def findpaths(CIJ, qmax, sources, savepths=False):
     '''
@@ -1091,31 +880,3 @@ def mean_first_passage_time(adjacency):
 
     return mfpt
 
-def 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.
-
-    Parameters
-    ----------
-    adj : np.ndarray
-        weighted/unweighted, directed/undirected adjacency matrix
-
-    Returns
-    -------
-    gediff : float
-        mean global diffusion efficiency
-    ediff : np.ndarray
-        pairwise NxN diffusion efficiency matrix
-    '''
-    n = len(adj)
-    adj = adj.copy()
-    mfpt = mean_first_passage_time(adj)
-    ediff = 1 / mfpt
-    np.fill_diagonal(ediff, 0)
-    gediff = np.sum(ediff) / (n ** 2 - n)
-    return gediff, ediff
diff --git a/bct/algorithms/efficiency.py b/bct/algorithms/efficiency.py
new file mode 100644
index 0000000..2a7cdf1
--- /dev/null
+++ b/bct/algorithms/efficiency.py
@@ -0,0 +1,361 @@
+from __future__ import division, print_function
+import numpy as np
+from bct.utils import cuberoot, BCTParamError, binarize, invert
+from .distance import distance_wei_floyd, mean_first_passage_time
+from ..due import due, BibTeX
+from ..citations import LATORA2001, ONNELA2005, FAGIOLO2007, RUBINOV2010
+
+@due.dcite(BibTeX(LATORA2001), description="Unweighted global efficiency")
+@due.dcite(BibTeX(ONNELA2005), description="Unweighted global efficiency")
+@due.dcite(BibTeX(FAGIOLO2007), description="Unweighted global efficiency")
+@due.dcite(BibTeX(RUBINOV2010), description="Unweighted global efficiency")
+def efficiency_bin(G, local=False):
+    '''
+    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.
+
+    Parameters
+    ----------
+    A : NxN np.ndarray
+        binary undirected connection matrix
+    local : bool
+        If True, computes local efficiency instead of global efficiency.
+        Default value = False.
+
+    Returns
+    -------
+    Eglob : float
+        global efficiency, only if local=False
+    Eloc : Nx1 np.ndarray
+        local efficiency, only if local=True
+    '''
+    def distance_inv(g):
+        D = np.eye(len(g))
+        n = 1
+        nPATH = g.copy()
+        L = (nPATH != 0)
+
+        while np.any(L):
+            D += n * L
+            n += 1
+            nPATH = np.dot(nPATH, g)
+            L = (nPATH != 0) * (D == 0)
+        D[np.logical_not(D)] = np.inf
+        D = 1 / D
+        np.fill_diagonal(D, 0)
+        return D
+
+    G = binarize(G)
+    n = len(G)  # number of nodes
+    if local:
+        E = np.zeros((n,))  # local efficiency
+
+        for u in range(n):
+            # V,=np.where(G[u,:])			#neighbors
+            # k=len(V)					#degree
+            # if k>=2:					#degree must be at least 2
+            #	e=distance_inv(G[V].T[V])
+            #	E[u]=np.sum(e)/(k*k-k)	#local efficiency computation
+
+            # find pairs of neighbors
+            V, = np.where(np.logical_or(G[u, :], G[u, :].T))
+            # inverse distance matrix
+            e = distance_inv(G[np.ix_(V, V)])
+            # symmetrized inverse distance matrix
+            se = e + e.T
+
+            # symmetrized adjacency vector
+            sa = G[u, V] + G[V, u].T
+            numer = np.sum(np.outer(sa.T, sa) * se) / 2
+            if numer != 0:
+                denom = np.sum(sa)**2 - np.sum(sa * sa)
+                E[u] = numer / denom  # local efficiency
+
+    else:
+        e = distance_inv(G)
+        E = np.sum(e) / (n * n - n)  # global efficiency
+    return E
+
+
+@due.dcite(BibTeX(LATORA2001), description="Weighted global efficiency")
+@due.dcite(BibTeX(ONNELA2005), description="Weighted global efficiency")
+@due.dcite(BibTeX(FAGIOLO2007), description="Weighted global efficiency")
+@due.dcite(BibTeX(RUBINOV2010), description="Weighted global efficiency")
+def efficiency_wei(Gw, local=False):
+    '''
+    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.
+
+    Parameters
+    ----------
+    W : NxN np.ndarray
+        undirected weighted connection matrix
+        (all weights in W must be between 0 and 1)
+    local = bool or enum
+        If True or 'local', computes local efficiency instead of global efficiency.
+        If False or 'global', uses the global efficiency
+        If 'original', will use the original algorithm provided by (Rubinov
+        & Sporns 2010). This version is not recommended. The local efficiency
+        calculation was improved in (Wang et al. 2016) as a true generalization
+        of the binary variant.
+        
+    Returns
+    -------
+    Eglob : float
+        global efficiency, only if local in (False, 'global')
+    Eloc : Nx1 np.ndarray
+        local efficiency, only if local in (True, 'local', 'original')
+
+    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 this weighted variant of the local
+    efficiency is hence not a strict generalization of the binary variant.
+
+    Algorithm:  Dijkstra's algorithm
+    '''
+    if local not in (True, False, 'local', 'global', 'original'):
+        raise BCTParamError("local param must be any of True, False, "
+            "'local', 'global', or 'original'")
+
+    def distance_inv_wei(G):
+        n = len(G)
+        D = np.zeros((n, n))  # distance matrix
+        D[np.logical_not(np.eye(n))] = np.inf
+
+        for u in range(n):
+            # distance permanence (true is temporary)
+            S = np.ones((n,), dtype=bool)
+            G1 = G.copy()
+            V = [u]
+            while True:
+                S[V] = 0  # distance u->V is now permanent
+                G1[:, V] = 0  # no in-edges as already shortest
+                for v in V:
+                    W, = np.where(G1[v, :])  # neighbors of smallest nodes
+                    td = np.array(
+                        [D[u, W].flatten(), (D[u, v] + G1[v, W]).flatten()])
+                    D[u, W] = np.min(td, axis=0)
+
+                if D[u, S].size == 0:  # all nodes reached
+                    break
+                minD = np.min(D[u, S])
+                if np.isinf(minD):  # some nodes cannot be reached
+                    break
+                V, = np.where(D[u, :] == minD)
+
+        np.fill_diagonal(D, 1)
+        D = 1 / D
+        np.fill_diagonal(D, 0)
+        return D
+
+    n = len(Gw)
+    Gl = invert(Gw, copy=True)  # connection length matrix
+    A = np.array((Gw != 0), dtype=int)
+    #local efficiency algorithm described by Rubinov and Sporns 2010, not recommended
+    if local == 'original':
+        E = np.zeros((n,))
+        for u in range(n):
+            # V,=np.where(Gw[u,:])		#neighbors
+            # k=len(V)					#degree
+            # if k>=2:					#degree must be at least 2
+            #	e=(distance_inv_wei(Gl[V].T[V])*np.outer(Gw[V,u],Gw[u,V]))**1/3
+            #	E[u]=np.sum(e)/(k*k-k)
+
+            # find pairs of neighbors
+            V, = np.where(np.logical_or(Gw[u, :], Gw[:, u].T))
+            # symmetrized vector of weights
+            sw = cuberoot(Gw[u, V]) + cuberoot(Gw[V, u].T)
+            # inverse distance matrix
+            e = distance_inv_wei(Gl[np.ix_(V, V)])
+            # symmetrized inverse distance matrix
+            se = cuberoot(e) + cuberoot(e.T)
+
+            numer = np.sum(np.outer(sw.T, sw) * se) / 2
+            if numer != 0:
+                # symmetrized adjacency vector
+                sa = A[u, V] + A[V, u].T
+                denom = np.sum(sa)**2 - np.sum(sa * sa)
+                # print numer,denom
+                E[u] = numer / denom  # local efficiency
+
+    #local efficiency algorithm described by Wang et al 2016, recommended
+    elif local in (True, 'local'):
+        E = np.zeros((n,))
+        for u in range(n):
+            V, = np.where(np.logical_or(Gw[u, :], Gw[:, u].T))
+            sw = cuberoot(Gw[u, V]) + cuberoot(Gw[V, u].T)
+            e = distance_inv_wei(cuberoot(Gl)[np.ix_(V, V)])
+            se = e+e.T
+         
+            numer = np.sum(np.outer(sw.T, sw) * se) / 2
+            if numer != 0:
+                # symmetrized adjacency vector
+                sa = A[u, V] + A[V, u].T
+                denom = np.sum(sa)**2 - np.sum(sa * sa)
+                # print numer,denom
+                E[u] = numer / denom  # local efficiency
+
+    elif local in (False, 'global'):
+        e = distance_inv_wei(Gl)
+        E = np.sum(e) / (n * n - n)
+    return E
+
+
+def 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.
+
+    Parameters
+    ----------
+    adj : np.ndarray
+        weighted/unweighted, directed/undirected adjacency matrix
+
+    Returns
+    -------
+    gediff : float
+        mean global diffusion efficiency
+    ediff : np.ndarray
+        pairwise NxN diffusion efficiency matrix
+    '''
+    n = len(adj)
+    adj = adj.copy()
+    mfpt = mean_first_passage_time(adj)
+    ediff = 1 / mfpt
+    np.fill_diagonal(ediff, 0)
+    gediff = np.sum(ediff) / (n ** 2 - n)
+    return gediff, ediff
+
+
+def resource_efficiency_bin(adj, lamb, spl=None, m=None):
+    '''
+    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).
+
+    Parameters 
+    ----------
+    adj : NxN np.ndarray
+        Unweighted, undirected adjacency matrix
+    lamb : float
+        Probability parameter of finding a path, set to np.nan if computation
+        of resource efficiency matrix is not desired.
+        Thanks Guido for making calling a variable lambda a syntax error in 
+        python while still not providing inline syntax for a fully first class
+        function.
+    spl : NxN np.ndarray
+        Shortest path length matrix (optional)
+    m : NxN np.ndarray
+        Transition probability matrix (optional)
+
+    Returns
+    -------
+    Eres : NxN np.ndarray
+        resource efficiency matrix. If lambda was provided as NaN, then Eres
+        will be np.nan
+    prob_spl : NxN np.ndarray
+        probabilistic shortest path matrix
+
+    Note that global measures for the resource efficiency and probabilistic
+    shortest paths exist and are well defined, they are the mean values of the
+    off diagonal elements
+    GEres = mean(Eres(~eye(N) > 0))
+    Gprob_SPL = mean(prob_SPL(~eye(N) > 0))
+    '''
+    n = len(adj)
+
+    if not np.isnan(lamb):
+        if lamb <= 0 or lamb >= 1:
+            raise BCTParamError("Lambda must be a nonzero probability")
+
+        z = np.zeros((n, n))
+
+    def prob_first_particle_arrival(m, l, hvec, lamb):
+        prob = np.zeros((n, n))
+        
+        if not np.isnan(lamb):
+            resources = np.zeros((n, n))    
+        else:
+            resources = None
+
+        # for each destination node h
+        for h in hvec:
+            B_h = m.copy()
+            B_h[h, :] = 0
+            B_h[h, h] = 1  
+            # h becomes absorbant state
+
+            B_h_L = B_h ** l
+
+            term = 1 - B_h_L[:, h]
+
+            prob[:, h] = 1-term
+
+            if not np.isnan(lamb):
+                resources[:, h] = np.repeat( np.log(1 - lamb), n, 0) / np.log(term)
+
+        return prob, resources
+
+    if spl is None:
+        spl, _, _ = distance_wei_floyd(adj)
+
+    if m is None:
+        m = np.linalg.solve( np.diag( np.sum(adj, axis=1) ), adj)
+
+    lvalues = np.unique(spl)
+    lvalues = lvalues[np.where(lvalues != 0)]
+
+    #a priori zero probability of going through SPL among nodes
+    prob_spl = np.zeros((n, n))
+    
+    for splvalue in lvalues:
+        hcols = np.where(spl == splvalue)
+        hvec = np.unique(hcols)
+        entries = spl == splvalue
+
+        prob_aux, z_aux = prob_first_particle_arrival(m, splvalue, hvec, lamb)
+
+        prob_aux[entries == 0] = 0
+        prob_spl += prob_aux
+
+        if not np.isnan(lamb):
+            z_aux[entries == 0] = 0
+            z += z_aux
+
+    np.fill_diagonal(prob_spl, 0)
+
+    if not np.isnan(lamb):
+        z[prob_spl == 1] = 1
+        with np.errstate(divide='ignore'):
+            Eres = 1 / z
+        np.fill_diagonal(Eres, 0)
+    else:
+        Eres = np.nan
+
+    return Eres, prob_spl
diff --git a/test/distance_test.py b/test/distance_test.py
index 1fb965e..f73d070 100644
--- a/test/distance_test.py
+++ b/test/distance_test.py
@@ -55,13 +55,6 @@ def test_charpath():
     assert not np.isnan(radius)
     assert not np.isnan(diameter)
 
-def test_diffusion_efficiency():
-    x = load_sample(thres=.23)
-    gde, ed = bct.diffusion_efficiency(x)
-    print(gde, np.sum(ed)) 
-    assert np.allclose(gde, .0069472)
-    assert np.allclose(np.sum(ed), 131.34, atol=.01)
-
 def test_distance_floyd():
     x = load_sample(thres=.31)
     print('guomish')

From 619a5bab173576ca60b745844499a67e50252742 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Tue, 20 Dec 2022 11:55:20 -0500
Subject: [PATCH 08/11] add rout_efficiency

---
 bct/algorithms/efficiency.py | 57 ++++++++++++++++++++++++++++++++++--
 1 file changed, 55 insertions(+), 2 deletions(-)

diff --git a/bct/algorithms/efficiency.py b/bct/algorithms/efficiency.py
index 2a7cdf1..a7c6919 100644
--- a/bct/algorithms/efficiency.py
+++ b/bct/algorithms/efficiency.py
@@ -241,7 +241,8 @@ def diffusion_efficiency(adj):
     n = len(adj)
     adj = adj.copy()
     mfpt = mean_first_passage_time(adj)
-    ediff = 1 / mfpt
+    with np.errstate(divide='ignore'):
+        ediff = 1 / mfpt
     np.fill_diagonal(ediff, 0)
     gediff = np.sum(ediff) / (n ** 2 - n)
     return gediff, ediff
@@ -318,7 +319,9 @@ def prob_first_particle_arrival(m, l, hvec, lamb):
             prob[:, h] = 1-term
 
             if not np.isnan(lamb):
-                resources[:, h] = np.repeat( np.log(1 - lamb), n, 0) / np.log(term)
+                with np.errstate(divide='ignore'):
+                    resources[:, h] = (np.repeat( np.log(1 - lamb), n, 0) / 
+                                       np.log(term))
 
         return prob, resources
 
@@ -359,3 +362,53 @@ def prob_first_particle_arrival(m, l, hvec, lamb):
         Eres = np.nan
 
     return Eres, prob_spl
+
+def rout_efficiency(D, transform=None):
+    '''
+    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).
+
+    Parameters
+    ----------
+    D : NxN np.ndarray
+        Weighted/unweighted directed/undirected connection weight or length
+        matrix
+    transform : str or None, optional
+        If `adjacency` is a connection weight array, specify a transform to map
+        input connection weights to connection lengths. Options include ['log',
+        'inv'], where 'log' is `-np.log(adjacency)` and 'inv' is `1/adjacency`.
+        Default: None
+
+    Returns
+    -------
+    GErout : float
+        Mean global routing efficiency
+    Erout : NxN np.ndarray
+        Pairwise routing efficiency matrix
+    Eloc : Nx1 np.ndarray
+        Local efficiency vector
+    '''
+    n = len(D)
+    Erout, _, _ = distance_wei_floyd(D, transform=transform)
+    with np.errstate(divide='ignore'):
+        Erout = 1 / Erout
+    np.fill_diagonal(Erout, 0)
+    GErout = (np.sum(Erout[np.where(np.logical_not(np.isnan(Erout)))]) / 
+              (n ** 2 - n))
+
+    Eloc = np.zeros((n,))
+    for u in range(n):
+        Gu, = np.where(np.logical_or(D[u, :], D[:, u].T))
+        nGu = len(Gu)
+        e, _, _ = distance_wei_floyd(D[Gu, :][:, Gu], transform=transform)
+        with np.errstate(divide='ignore'):
+            e = 1 / e
+        np.fill_diagonal(e, 0)
+        Eloc[u] = np.sum(e) / nGu
+
+    return GErout, Erout, Eloc
+
+    

From 832b1422cdf3a54347912a31d3b1d959c88bffa8 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Tue, 20 Dec 2022 11:56:13 -0500
Subject: [PATCH 09/11] add efficiency test cases

---
 test/efficiency_test.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 test/efficiency_test.py

diff --git a/test/efficiency_test.py b/test/efficiency_test.py
new file mode 100644
index 0000000..c303979
--- /dev/null
+++ b/test/efficiency_test.py
@@ -0,0 +1,28 @@
+from .load_samples import load_sample, load_directed_sample
+import numpy as np
+import bct
+
+def test_diffusion_efficiency():
+    x = load_sample(thres=.23)
+    gde, ed = bct.diffusion_efficiency(x)
+    print(gde, np.sum(ed)) 
+    assert np.allclose(gde, .0069472)
+    assert np.allclose(np.sum(ed), 131.34, atol=.01)
+
+def test_resource_efficiency():
+    x = load_sample(thres=.39)
+    x = bct.binarize(x)
+
+    eres, prob = bct.resource_efficiency_bin(x, .35)
+
+    assert np.allclose(np.sum(eres), 323.5398, atol=.0001)
+    assert np.allclose(np.sum(prob), 138.0000, atol=.0001)
+
+def test_rout_efficiency():
+    x = load_directed_sample(thres=1)
+    GErout, Erout, Eloc = bct.rout_efficiency(x, 'inv')
+
+    assert np.allclose(np.sum(Erout), 9515.25, atol=.01)
+    assert np.allclose(GErout, 1.0655, atol=.0001)
+    assert np.allclose(np.sum(Eloc), 2906.574, atol=.001)
+    

From f42fbb83e5ff8058601935e61cf9a06e48213910 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Sat, 11 Feb 2023 15:47:09 -0500
Subject: [PATCH 10/11] add navigation_wu, test case not working yet because
 uses 138 parcellation, old parcellation, need real parcellation with average
 talairach centroids

---
 bct/algorithms/distance.py | 121 +++++++++++++++++++++++++++++++++++++
 test/core_test.py          |   1 -
 test/distance_test.py      |  27 ++++++++-
 3 files changed, 147 insertions(+), 2 deletions(-)

diff --git a/bct/algorithms/distance.py b/bct/algorithms/distance.py
index 781b7dc..a75b338 100644
--- a/bct/algorithms/distance.py
+++ b/bct/algorithms/distance.py
@@ -880,3 +880,124 @@ def mean_first_passage_time(adjacency):
 
     return mfpt
 
+
+def navigation_wu(L, D, max_hops=None):
+    '''
+    Navigation of connectivity length matrix L guided by nodal distance D
+   
+    % Navigation
+    [sr, PL_bin, PL_wei] = navigation_wu(L,D);
+    % Binary shortest path length
+    sp_PL_bin = distance_bin(L);
+    % Weighted shortest path length
+    sp_PL_wei = distance_wei_floyd(L);
+    % Binary efficiency ratio
+    er_bin = mean(mean(sp_PL_bin./PL_bin));
+    % Weighted efficiency ratio
+    er_wei = mean(mean(sp_PL_wei./PL_wei));
+   
+    Parameters
+    ----------
+    L : NxN np.ndarray
+        Weighted/unweighted directed/undirected NxN SC matrix of connection
+        *lengths*, L(i,j) is the strength-to-length remapping of the connection
+        weight between i and j. L(i,j) = 0 denotes the lack of a connection 
+        between i and j.
+   
+    D : NxN np.ndarray
+        Symmetric NxN nodal distance matrix (e.g., Euclidean distance between 
+        node centroids)
+
+    max_hops : int | None
+        Limits the maximum number of hops of navigation paths
+   
+    Returns
+    ------- 
+    sr : int
+        Success ratio scalar, proportion of node pairs successfully reached by
+        navigation
+    
+    PL_bin : NxN np.ndarray
+        NxN matrix of binary navigation path length (i.e., number of hops in 
+        navigation paths). Infinite values indicate failed navigation paths
+    
+    PL_wei : NxN np.ndarray
+        NxN matrix of weighted navigation path length (i.e., sum of connection
+        weights as defined by C along navigation path). Infinite values
+        indicate failed paths.
+   
+    PL_dis : NxN np.ndarray
+        NxN matrix of distance-based navigation path length (i.e., sum of
+        connection distances as defined by D along navigation paths. Infinite 
+        values indicate failed paths.
+   
+    paths - dict(tuple -> list)
+        array of nodes comprising navigation paths. The key (i,j) specifies
+        the path from i to j, and the value is a list of all nodes traveled
+        between i and j.
+    '''
+
+    n = len(L)
+    PL_bin = np.zeros((n, n))
+    PL_wei = np.zeros((n, n))
+    PL_dis = np.zeros((n, n))
+    paths = {}
+
+    for i in range(n):
+        for j in range(n):
+            if i == j:
+                continue
+            
+            curr_node = i
+            last_node = curr_node
+            target = j
+            curr_paths = [curr_node]
+
+            pl_bin = 0
+            pl_wei = 0
+            pl_dis = 0
+
+            while curr_node != target:
+                #print(curr_node, "WHEEF")
+                #print(np.where(L[curr_node, :] != 0))         
+                #print(np.shape(np.where(L[curr_node, :] != 0)))
+
+                neighbors, = np.where(L[curr_node, :] != 0)
+                if len(neighbors) == 0:
+                    pl_bin = np.inf
+                    pl_wei = np.inf
+                    pl_dis = np.inf
+                    break
+
+                min_ix = np.argmin(D[target, neighbors])
+                next_node = neighbors[min_ix]
+
+                if (next_node == last_node or 
+                    (max_hops is not None and pl_bin > max_hops)):
+
+                    pl_bin = np.inf
+                    pl_wei = np.inf
+                    pl_dis = np.inf
+                    break
+
+                curr_paths.append(next_node)
+                pl_bin += 1
+                pl_wei += L[curr_node, next_node]
+                pl_dis += D[curr_node, next_node]
+
+                last_node = curr_node
+                curr_node = next_node
+                
+            PL_bin[i, j] = pl_bin
+            PL_wei[i, j] = pl_wei
+            PL_dis[i, j] = pl_dis
+            paths[(i, j)] = curr_paths
+
+    np.fill_diagonal(PL_bin, np.inf)
+    np.fill_diagonal(PL_wei, np.inf)
+    np.fill_diagonal(PL_dis, np.inf)
+
+    inf_ixes, = np.where(PL_bin.flat == np.inf)
+    sr = 1 - (len(inf_ixes) - n)/(n**2 - n)
+
+    return sr, PL_bin, PL_wei, PL_dis, paths
diff --git a/test/core_test.py b/test/core_test.py
index b5a0930..29d5553 100644
--- a/test/core_test.py
+++ b/test/core_test.py
@@ -27,4 +27,3 @@ def test_clique_communities():
     print(cis.shape, np.max(np.sum(cis, axis=0)))
     print(np.sum(cis, axis=1))
     assert np.sum(cis) == 199
-    assert 4 == 8
diff --git a/test/distance_test.py b/test/distance_test.py
index f73d070..611b9fc 100644
--- a/test/distance_test.py
+++ b/test/distance_test.py
@@ -57,8 +57,33 @@ def test_charpath():
 
 def test_distance_floyd():
     x = load_sample(thres=.31)
-    print('guomish')
     spli, hopsi, pmati = bct.distance_wei_floyd(x, transform='inv')
     print(np.sum(spli))
     assert np.allclose(np.sum(spli), 11536.1, atol=.01)
     
+def test_navigation_wu():
+    x = load_sample(thres=.24)
+    x_len = bct.invert(x)
+
+    #randomly generate distances for testing purposes
+    n = len(x)
+    while True:
+        centroids = np.random.randint(512, size=(n, 3))
+        #make sure every centroid is unique
+        if len(np.unique(centroids, axis=0)) == n:
+            break
+
+    d = np.zeros((n, n))
+    for i in range(n):
+        for j in range(n):
+            d[i, j] = np.linalg.norm(centroids[i, :] - centroids[j, :])
+            
+        
+    sr, plbin, plwei, pldis, paths = bct.navigation_wu(x_len, d, max_hops=14)
+
+    sr2, plbin2, plwei2, pldis2, paths2 = bct.navigation_wu(x_len, d, max_hops=None)
+
+    #TODO create the centroids for an actual bit of sample data and converge the matlab algorithm
+    #this procedure of random centroid generation generates a random reachability which is usually around 45-60%
+    #but not guaranteed
+

From 0906c856fa52255fc4845c386674bf8c775f2e71 Mon Sep 17 00:00:00 2001
From: aestrivex <aestrivex@gmail.com>
Date: Wed, 26 Jul 2023 12:39:28 -0400
Subject: [PATCH 11/11] corrected version

---
 README.md | 2 +-
 setup.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index ebb0aa1..f9aa136 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# Brain Connectivity Toolbox for Python version 0.6.0
+# Brain Connectivity Toolbox for Python version 0.6.1
 
 Author: Roan LaPlante <rlaplant@nmr.mgh.harvard.edu>
 
diff --git a/setup.py b/setup.py
index 8677f93..be26be7 100644
--- a/setup.py
+++ b/setup.py
@@ -7,7 +7,7 @@ def read(fname):
 
 setuptools.setup(
     name="bctpy",
-    version="0.6.0",
+    version="0.6.1",
     maintainer="Roan LaPlante",
     maintainer_email="rlaplant@nmr.mgh.harvard.edu",
     description=("Brain Connectivity Toolbox for Python"),