Merge pull request #185 from jGaboardi/ring_bug

Debug of data structure/algo that led to ring/graph errors

spaghetti/network.py

@@ -46,11 +46,11 @@ class Network:
         nodes. Default is ``True``.
     
     w_components : bool
-        Set to ``True`` to record connected components from a
+        Set to ``False`` to not record connected components from a
         `libpysal.weights.weights.W 
         <https://libpysal.readthedocs.io/en/latest/generated/
         libpysal.weights.W.html#libpysal.weights.W>`_
-        object. Default is False.
+        object. Default is True.
         
         
     weightings : {dict, bool}
@@ -102,6 +102,10 @@ class Network:
         Keys are the graph edge ids (tuple). Values are the graph edge
         length (``float``).
     
+    non_articulation_points : list
+        All vertices with degree 2 that are not in an isolated
+        island ring (loop) component
+    
     w_network : `libpysal.weights.weights.W <https://libpysal.readthedocs.io/en/latest/generated/libpysal.weights.W.html#libpysal.weights.W>`_
         Weights object created from the network arcs
     
@@ -161,7 +165,7 @@ class Network:
     """
 
     def __init__(self, in_data=None, vertex_sig=11, unique_arcs=True,
-                 extractgraph=True, w_components=False, weightings=False):
+                 extractgraph=True, w_components=True, weightings=False):
 
         if in_data is not None:
             self.in_data = in_data
@@ -347,15 +351,13 @@ def extractgraph(self):
         self.edges = []
         self.edge_lengths = {}
 
-        # Find all nodes with cardinality 2. (non-articulation points)
-        arc_vertices = []
-        for k, v in self.adjacencylist.items():
-            # len(v) == 1 #cul-de-sac
-            # len(v) == 2 #bridge segment
-            # len(v) > 2 #intersection
-            if len(v) == 2:
-                arc_vertices.append(k)
-
+        # Find all vertices with degree 2 that are not in an isolated
+        # island ring (loop) component. These are non-articulation
+        # points on the graph representation.
+        non_articulation_points = self._yield_napts()
+        # retain non_articulation_points as an attribute
+        self.non_articulation_points = list(non_articulation_points)
+
         # Start with a copy of the spatial representation and
         # iteratively remove edges deemed to be segments.
         self.edges = copy.deepcopy(self.arcs)
@@ -367,21 +369,23 @@ def extractgraph(self):
 
         # build up bridges "rooted" on the initial
         # non-articulation points
-        bridges = []
-        for s in arc_vertices:
+        bridge_roots = []
+        for s in non_articulation_points:
             bridge = [s]
-            neighbors = self._yieldneighbor(s, arc_vertices, bridge)
+            neighbors = self._yieldneighbor(s,
+                                            non_articulation_points,
+                                            bridge)
             while neighbors:
                 cnode = neighbors.pop()
-                arc_vertices.remove(cnode)
+                non_articulation_points.remove(cnode)
                 bridge.append(cnode)
                 newneighbors = self._yieldneighbor(cnode,
-                                                   arc_vertices,
+                                                   non_articulation_points,
                                                    bridge)
                 neighbors += newneighbors
-            bridges.append(bridge)
+            bridge_roots.append(bridge)
 
-        for bridge in bridges:
+        for bridge in bridge_roots:
             if len(bridge) == 1:
                 n = self.adjacencylist[bridge[0]]
                 new_edge = tuple(sorted([n[0], n[1]]))
@@ -426,9 +430,104 @@ def extractgraph(self):
                     self.edge_lengths.pop(r, None)
                     self.edges_to_arcs[r] = new_edge
                 self.edge_lengths[new_edge] = cumulative_length
-
+
+            # add the updated graph edge
             self.edges.append(new_edge)
-        self.edges = sorted(self.edges)
+
+        # converted the graph edges into a sorted set to prune out
+        # duplicate graph edges created during simplification
+        self.edges = sorted(set(self.edges))
+
+
+    def _yield_napts(self):
+        """Find all nodes with degree 2 that are not in an isolated
+        island ring (loop) component. These are non-articulation
+        points on the graph representation.
+        
+        Returns
+        -------
+        
+        napts : list
+            non-articulation points on a graph representation
+        
+        """
+
+        # non-articulation points
+        napts = set()
+
+        # network vertices remaining to evaluate
+        unvisted = set(self.vertices.values())
+
+        while unvisted:
+
+            # iterate over each component
+            for component_id, ring in self.network_component_is_ring.items():
+
+                # evaluate for non-articulation points
+                napts, unvisted = self._evaluate_napts(napts, unvisted,
+                                                       component_id, ring)
+
+        napts = list(napts)
+
+        return napts
+
+
+    def _evaluate_napts(self, napts, unvisited, component_id, ring):
+        """Evaluate one connected component in a network for
+        non-articulation points (napts) and return an updated set of
+        napts and unvisted vertices.
+        
+        Parameters
+        ----------
+        
+        napts : set
+            Non-articulation points (napts) in the network. The
+            'napts' here do not include those within an isolated
+            loop island.
+        
+        unvisited : set
+            Vertices left to evaluate in the network.
+        
+        component_id : int
+            ID for the network connected component for the
+            current iteration of the algorithm.
+        
+        ring : bool
+            Network component is isolated island loop ``True`` or
+            not ``False``.
+        
+        Returns
+        -------
+        
+        napts : set
+            Updated 'napts' object.
+        
+        unvisited : set
+            Updated 'napts' object.
+        
+        """
+
+        # iterate over each `edge` of the `component`
+        for component in self.network_component2arc[component_id]:
+
+            # each `component` has two vertices
+            for vertex in component:
+
+                # if `component` is not an isolated island
+                # and `vertex` has exactly 2 neighbors,
+                # add `vertex` to `napts`
+                if not ring:
+                    if len(self.adjacencylist[vertex]) == 2:
+                        napts.add(vertex)
+
+                # remove `vertex` from `unvisited` if
+                # it is still in the set
+                try:
+                    unvisited.remove(vertex)
+                except KeyError:
+                    pass
+
+        return napts, unvisited
 
 
     def _yieldneighbor(self, vtx, arc_vertices, bridge):

spaghetti/tests/test_network.py

@@ -45,7 +45,8 @@ def test_network_data_read(self):
     def test_network_from_geopandas(self):
         # network instantiated from geodataframe
         gdf = geopandas.read_file(self.path_to_shp)
-        self.ntw_from_gdf = spgh.Network(in_data=gdf)
+        self.ntw_from_gdf = spgh.Network(in_data=gdf,
+                                         w_components=True)
 
         # gdf test against known
         self.assertEqual(len(self.ntw_from_gdf.arcs),
@@ -64,7 +65,7 @@ def test_contiguity_weights(self):
         observed_network_histo = self.ntw_from_shp.w_network.histogram
         self.assertEqual(known_network_histo, observed_network_histo)
 
-        known_graph_histo = [(2, 2), (3, 2), (4, 45), (5, 82), (6, 48)]
+        known_graph_histo = [(2, 2), (3, 2), (4, 47), (5, 80), (6, 48)]
         observed_graph_histo = self.ntw_from_shp.w_graph.histogram
         self.assertEqual(observed_graph_histo, known_graph_histo)
 
@@ -73,7 +74,7 @@ def test_components(self):
         observed_network_arc = self.ntw_from_shp.network_component2arc[0][-1]
         self.assertEqual(observed_network_arc, known_network_arc)
 
-        known_graph_edge = (206, 207)
+        known_graph_edge = (207, 208)
         observed_graph_edge = self.ntw_from_shp.graph_component2edge[0][-1]
         self.assertEqual(observed_graph_edge, known_graph_edge)
 

spaghetti/tests/test_network_api.py

@@ -44,7 +44,8 @@ def test_network_data_read(self):
     def test_network_from_geopandas(self):
         # network instantiated from geodataframe
         gdf = geopandas.read_file(self.path_to_shp)
-        self.ntw_from_gdf = spgh.Network(in_data=gdf)
+        self.ntw_from_gdf = spgh.Network(in_data=gdf,
+                                         w_components=True)
 
         # gdf test against known
         self.assertEqual(len(self.ntw_from_gdf.arcs),
@@ -63,7 +64,7 @@ def test_contiguity_weights(self):
         observed_network_histo = self.ntw_from_shp.w_network.histogram
         self.assertEqual(known_network_histo, observed_network_histo)
 
-        known_graph_histo = [(2, 2), (3, 2), (4, 45), (5, 82), (6, 48)]
+        known_graph_histo = [(2, 2), (3, 2), (4, 47), (5, 80), (6, 48)]
         observed_graph_histo = self.ntw_from_shp.w_graph.histogram
         self.assertEqual(observed_graph_histo, known_graph_histo)
 
@@ -72,7 +73,7 @@ def test_components(self):
         observed_network_arc = self.ntw_from_shp.network_component2arc[0][-1]
         self.assertEqual(observed_network_arc, known_network_arc)
 
-        known_graph_edge = (206, 207)
+        known_graph_edge = (207, 208)
         observed_graph_edge = self.ntw_from_shp.graph_component2edge[0][-1]
         self.assertEqual(observed_graph_edge, known_graph_edge)