New Feature - Resistance Distance (#3385)

* First pass of resistance distance logic

* Made inverting the weight optional.
Removed zero_approx parameter
Added name to contributors

* csr/csc slice and tranformation more effecient and simple for single row/col deletion

* * Split permulation count into seperate funciton
* added tests
* Faster submatrix approach

* Made tests more readable

* Added whitespace

* Clean up styling, pep8

* Fixed nx test class import

* Simplied test imports

* Fixed prior merge error

* Fixed prior merge issue in docstring

* Couple of final pep8 changes

* Added testing for non-scipy builds

* minor tweaks of if logic and comments

5 files changed, 270 insertions, 10 deletions

diff --git a/.gitignore b/.gitignore
index 549df698..8a63b3b1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -47,3 +47,6 @@ networkx.egg-info/
 
 # PyCharm project file
 .idea
+
+# VS Code settings
+.vscode
\ No newline at end of file
diff --git a/CONTRIBUTORS.rst b/CONTRIBUTORS.rst
index 995a8a11..f8631917 100644
--- a/CONTRIBUTORS.rst
+++ b/CONTRIBUTORS.rst
@@ -123,6 +123,7 @@ is partially historical, and now, mostly arbitrary.
 - Mike Trenfield
 - Jon Crall, Github: `Erotemic <https://github.com/Erotemic>`_
 - Issa Moradnejad, `Github <https://github.com/Moradnejad>`_, `LinkedIn <https://linkedin.com/in/moradnejad/>`_
+- Brian Kiefer, Github: `bkief <https://github.com/bkief>`_
 - Julien Klaus
 - Weisheng Si, Github: `ws4u <https://github.com/ws4u>`_
 
diff --git a/doc/reference/algorithms/distance_measures.rst b/doc/reference/algorithms/distance_measures.rst
index ee3f1fca..1a6a3c71 100644
--- a/doc/reference/algorithms/distance_measures.rst
+++ b/doc/reference/algorithms/distance_measures.rst
@@ -13,5 +13,6 @@ Distance Measures
    extrema_bounding
    periphery
    radius
+   resistance_distance
 
 
diff --git a/networkx/algorithms/distance_measures.py b/networkx/algorithms/distance_measures.py
index a59c209b..e2b7a960 100644
--- a/networkx/algorithms/distance_measures.py
+++ b/networkx/algorithms/distance_measures.py
@@ -9,11 +9,15 @@
 #
 # Authors: Aric Hagberg (hagberg@lanl.gov)
 #          Dan Schult (dschult@colgate.edu)
+#          Brian Kiefer (bkiefer@asu.edu)
 """Graph diameter, radius, eccentricity and other properties."""
-import networkx
+
+import networkx as nx
+from networkx.utils import not_implemented_for
 
 __all__ = ['extrema_bounding', 'eccentricity', 'diameter',
-           'radius', 'periphery', 'center', 'barycenter']
+           'radius', 'periphery', 'center', 'barycenter',
+           'resistance_distance']
 
 
 def extrema_bounding(G, compute="diameter"):
@@ -90,10 +94,10 @@ def extrema_bounding(G, compute="diameter"):
         high = not high
 
         # get distances from/to current node and derive eccentricity
-        dist = dict(networkx.single_source_shortest_path_length(G, current))
+        dist = dict(nx.single_source_shortest_path_length(G, current))
         if len(dist) != N:
             msg = ('Cannot compute metric because graph is not connected.')
-            raise networkx.NetworkXError(msg)
+            raise nx.NetworkXError(msg)
         current_ecc = max(dist.values())
 
         # print status update
@@ -224,14 +228,14 @@ def eccentricity(G, v=None, sp=None):
     e = {}
     for n in G.nbunch_iter(v):
         if sp is None:
-            length = networkx.single_source_shortest_path_length(G, n)
+            length = nx.single_source_shortest_path_length(G, n)
             L = len(length)
         else:
             try:
                 length = sp[n]
                 L = len(length)
             except TypeError:
-                raise networkx.NetworkXError('Format of "sp" is invalid.')
+                raise nx.NetworkXError('Format of "sp" is invalid.')
         if L != order:
             if G.is_directed():
                 msg = ('Found infinite path length because the digraph is not'
@@ -239,7 +243,7 @@ def eccentricity(G, v=None, sp=None):
             else:
                 msg = ('Found infinite path length because the graph is not'
                        ' connected')
-            raise networkx.NetworkXError(msg)
+            raise nx.NetworkXError(msg)
 
         e[n] = max(length.values())
 
@@ -416,7 +420,7 @@ def barycenter(G, weight=None, attr=None, sp=None):
     periphery
     """
     if sp is None:
-        sp = networkx.shortest_path_length(G, weight=weight)
+        sp = nx.shortest_path_length(G, weight=weight)
     else:
         sp = sp.items()
         if weight is not None:
@@ -424,7 +428,7 @@ def barycenter(G, weight=None, attr=None, sp=None):
     smallest, barycenter_vertices, n = float('inf'), [], len(G)
     for v, dists in sp:
         if len(dists) < n:
-            raise networkx.NetworkXNoPath(
+            raise nx.NetworkXNoPath(
                 ("Input graph %r is disconnected, so every induced subgraph "
                  "has infinite barycentricity.") % G)
         barycentricity = sum(dists.values())
@@ -436,3 +440,151 @@ def barycenter(G, weight=None, attr=None, sp=None):
         elif barycentricity == smallest:
             barycenter_vertices.append(v)
     return barycenter_vertices
+
+
+def _laplacian_submatrix(node, mat, node_list):
+    """Removes row/col from a sparse matrix and returns the submatrix
+    """
+    j = node_list.index(node)
+    n = list(range(len(node_list)))
+    n.pop(j)
+
+    if mat.shape[0] != mat.shape[1]:
+        raise nx.NetworkXError('Matrix must be square')
+    elif len(node_list) != mat.shape[0]:
+        msg = "Node list length does not match matrix dimentions"
+        raise nx.NetworkXError(msg)
+
+    mat = mat.tocsr()
+    mat = mat[n, :]
+
+    mat = mat.tocsc()
+    mat = mat[:, n]
+
+    node_list.pop(j)
+
+    return mat, node_list
+
+
+def _count_lu_permutations(perm_array):
+    """Counts the number of permutations in SuperLU perm_c or perm_r
+    """
+    perm_cnt = 0
+    arr = perm_array.tolist()
+    for i in range(len(arr)):
+        if i != arr[i]:
+            perm_cnt += 1
+            n = arr.index(i)
+            arr[n] = arr[i]
+            arr[i] = i
+
+    return perm_cnt
+
+
+@not_implemented_for('directed')
+def resistance_distance(G, nodeA, nodeB, weight=None, invert_weight=True):
+    """Returns the resistance distance between node A and node B on graph G.
+
+    The resistance distance between two nodes of a graph is akin to treating
+    the graph as a grid of resistorses with a resistance equal to the provided
+    weight.
+
+    If weight is not provided, then a weight of 1 is used for all edges.
+
+    Parameters
+    ----------
+    G : NetworkX graph
+       A graph
+
+    nodeA : node
+      A node within graph G.
+
+    nodeB : node
+      A node within graph G, exclusive of Node A.
+
+    weight : string or None, optional (default=None)
+       The edge data key used to compute the resistance distance.
+       If None, then each edge has weight 1.
+
+    invert_weight : boolean (default=True)
+        Proper calculation of resistance distance requires building the
+        Laplacian matrix with the reciprocal of the weight. Not required
+        if the weight is already inverted. Weight cannot be zero.
+
+    Returns
+    -------
+    rd : float
+       Value of effective resistance distance
+
+    Notes
+    -----
+    Overview discussion:
+    * https://en.wikipedia.org/wiki/Resistance_distance
+    * http://mathworld.wolfram.com/ResistanceDistance.html
+
+    Additional details:
+    Vaya Sapobi Samui Vos, “Methods for determining the effective resistance,” M.S.,
+    Mathematisch Instituut, Universiteit Leiden, Leiden, Netherlands, 2016
+    Available: `Link to thesis <https://www.universiteitleiden.nl/binaries/content/assets/science/mi/scripties/master/vos_vaya_master.pdf>`_
+    """
+    import numpy as np
+    import scipy.sparse
+
+    if not nx.is_connected(G):
+        msg = ('Graph G must be strongly connected.')
+        raise nx.NetworkXError(msg)
+    elif nodeA not in G:
+        msg = ('Node A is not in graph G.')
+        raise nx.NetworkXError(msg)
+    elif nodeB not in G:
+        msg = ('Node B is not in graph G.')
+        raise nx.NetworkXError(msg)
+    elif nodeA == nodeB:
+        msg = ('Node A and Node B cannot be the same.')
+        raise nx.NetworkXError(msg)
+
+    G = G.copy()
+    node_list = list(G)
+
+    if invert_weight and weight is not None:
+        if G.is_multigraph():
+            for (u, v, k, d) in G.edges(keys=True, data=True):
+                d[weight] = 1/d[weight]
+        else:
+           for (u, v, d) in G.edges(data=True):
+               d[weight] = 1/d[weight]
+    # Replace with collapsing topology or approximated zero?
+
+    # Using determinants to compute the effective resistance is more memory
+    # efficent than directly calculating the psuedo-inverse
+    L = nx.laplacian_matrix(G, node_list, weight=weight)
+
+    Lsub_a, node_list_a = _laplacian_submatrix(nodeA, L.copy(),
+                                               node_list[:])
+    Lsub_ab, node_list_ab = _laplacian_submatrix(nodeB, Lsub_a.copy(),
+                                                 node_list_a[:])
+
+    # Factorize Laplacian submatrixes and extract diagonals
+    # Order the diagonals to minimize the likelihood over overflows
+    # during computing the determinant
+    lu_a = scipy.sparse.linalg.splu(Lsub_a, options=dict(SymmetricMode=True))
+    LdiagA = lu_a.U.diagonal()
+    LdiagA_s = np.product(np.sign(LdiagA)) * np.product(lu_a.L.diagonal())
+    LdiagA_s *= (-1)**_count_lu_permutations(lu_a.perm_r)
+    LdiagA_s *= (-1)**_count_lu_permutations(lu_a.perm_c)
+    LdiagA = np.absolute(LdiagA)
+    LdiagA = np.sort(LdiagA)
+
+    lu_ab = scipy.sparse.linalg.splu(Lsub_ab, options=dict(SymmetricMode=True))
+    LdiagAB = lu_ab.U.diagonal()
+    LdiagAB_s = np.product(np.sign(LdiagAB)) * np.product(lu_ab.L.diagonal())
+    LdiagAB_s *= (-1)**_count_lu_permutations(lu_ab.perm_r)
+    LdiagAB_s *= (-1)**_count_lu_permutations(lu_ab.perm_c)
+    LdiagAB = np.absolute(LdiagAB)
+    LdiagAB = np.sort(LdiagAB)
+
+    # Calculate the ratio of determinant, rd = det(Lsub_ab)/det(Lsub_a)
+    Ldet = np.product(np.divide(np.append(LdiagAB, [1]), LdiagA))
+    rd = Ldet * LdiagAB_s / LdiagA_s
+
+    return rd
diff --git a/networkx/algorithms/tests/test_distance_measures.py b/networkx/algorithms/tests/test_distance_measures.py
index 66597704..945818c3 100644
--- a/networkx/algorithms/tests/test_distance_measures.py
+++ b/networkx/algorithms/tests/test_distance_measures.py
@@ -1,7 +1,9 @@
 #!/usr/bin/env python
 from random import Random
+from nose import SkipTest
 from nose.tools import assert_equal, assert_is_instance, \
-        assert_raises, raises, assert_less_equal, assert_false
+        assert_raises, raises, assert_less_equal, assert_false, \
+        assert_true
 
 import networkx as nx
 from networkx import convert_node_labels_to_integers as cnlti
@@ -91,6 +93,107 @@ class TestDistance:
         nx.eccentricity(DG)
 
 
+class TestResistanceDistance:
+    @classmethod
+    def setupClass(cls):
+        global np
+        global sp_sparse
+        try:
+            import numpy as np
+        except ImportError:
+            raise SkipTest('NumPy not available.')
+        try:
+            import scipy.sparse as sp_sparse
+        except ImportError:
+            raise SkipTest('SciPy Sparse not available.')
+
+    def setUp(self):
+        G = nx.Graph()
+        G.add_edge(1, 2, weight=2)
+        G.add_edge(2, 3, weight=4)
+        G.add_edge(3, 4, weight=1)
+        G.add_edge(1, 4, weight=3)
+        self.G = G
+
+    def test_laplacian_submatrix(self):
+        from networkx.algorithms.distance_measures import _laplacian_submatrix
+        M = sp_sparse.csr_matrix([[1, 2, 3],
+                                     [4, 5, 6],
+                                     [7, 8, 9]], dtype=np.float32)
+        N = sp_sparse.csr_matrix([[5, 6],
+                                     [8, 9]], dtype=np.float32)
+        Mn, Mn_nodelist = _laplacian_submatrix(1, M, [1, 2, 3])
+        assert_equal(Mn_nodelist, [2, 3])
+        assert_true(np.allclose(Mn.toarray(), N.toarray()))
+
+    @raises(nx.NetworkXError)
+    def test_laplacian_submatrix_square(self):
+        from networkx.algorithms.distance_measures import _laplacian_submatrix
+        M = sp_sparse.csr_matrix([[1, 2],
+                                     [4, 5],
+                                     [7, 8]], dtype=np.float32)
+        _laplacian_submatrix(1, M, [1, 2, 3])
+
+    @raises(nx.NetworkXError)
+    def test_laplacian_submatrix_matrix_node_dim(self):
+        from networkx.algorithms.distance_measures import _laplacian_submatrix
+        M = sp_sparse.csr_matrix([[1, 2, 3],
+                                     [4, 5, 6],
+                                     [7, 8, 9]], dtype=np.float32)
+        _laplacian_submatrix(1, M, [1, 2, 3, 4])
+
+    def test_resistance_distance(self):
+        rd = nx.resistance_distance(self.G, 1, 3, 'weight', True)
+        test_data = 1/(1/(2+4) + 1/(1+3))
+        assert_equal(round(rd, 5), round(test_data, 5))
+
+    def test_resistance_distance_noinv(self):
+        rd = nx.resistance_distance(self.G, 1, 3, 'weight', False)
+        test_data = 1/(1/(1/2+1/4) + 1/(1/1+1/3))
+        assert_equal(round(rd, 5), round(test_data, 5))
+
+    def test_resistance_distance_no_weight(self):
+        rd = nx.resistance_distance(self.G, 1, 3)
+        assert_equal(round(rd, 5), 1)
+
+    def test_resistance_distance_neg_weight(self):
+        self.G[2][3]['weight'] = -4
+        rd = nx.resistance_distance(self.G, 1, 3, 'weight', True)
+        test_data = 1/(1/(2+-4) + 1/(1+3))
+        assert_equal(round(rd, 5), round(test_data, 5))
+
+    def test_multigraph(self):
+        G = nx.MultiGraph()
+        G.add_edge(1, 2, weight=2)
+        G.add_edge(2, 3, weight=4)
+        G.add_edge(3, 4, weight=1)
+        G.add_edge(1, 4, weight=3)
+        rd = nx.resistance_distance(G, 1, 3, 'weight', True)
+        assert_true(np.isclose(rd, 1/(1/(2+4) + 1/(1+3))))
+
+    @raises(ZeroDivisionError)
+    def test_resistance_distance_div0(self):
+        self.G[1][2]['weight'] = 0
+        nx.resistance_distance(self.G, 1, 3, 'weight')
+
+    @raises(nx.NetworkXError)
+    def test_resistance_distance_not_connected(self):
+        self.G.add_node(5)
+        nx.resistance_distance(self.G, 1, 5)
+
+    @raises(nx.NetworkXError)
+    def test_resistance_distance_same_node(self):
+        nx.resistance_distance(self.G, 1, 1)
+
+    @raises(nx.NetworkXError)
+    def test_resistance_distance_nodeA_not_in_graph(self):
+        nx.resistance_distance(self.G, 9, 1)
+
+    @raises(nx.NetworkXError)
+    def test_resistance_distance_nodeB_not_in_graph(self):
+        nx.resistance_distance(self.G, 1, 9)
+
+
 class TestBarycenter(object):
     """Test :func:`networkx.algorithms.distance_measures.barycenter`."""
     def barycenter_as_subgraph(self, g, **kwargs):