Repository: torressa/cspy (language: python)

Issue:
Preprocessing fails to find nodes and may find incorrect nodes
**Describe the bug**
In [`preprocessing.py`](https://github.com/torressa/cspy/blob/master/src/python/preprocessing.py#L37-L51), the check for exceptions when where are no nodes violating the resource limit is tripped early by the "Source" or "Sink" node itself.  The paths for these are necessarily 1 long, and may come up in the list comprehension before other actually resource constraint violated paths.  The exception catches the index error, meaning we never add some of the violated nodes.

In large graphs, assuming a random order of keys in the dict, there is a smaller probability of this happening.  But in small graphs, the chance of missing a node is much higher.

A somewhat related issue is that the node selected for removal may be wrong.  By picking the next-to-last node in a path, we may be removing a node that is on some other path between the Source and the Sink.  I think the correct node to remove should be the last node in the path, the one we cannot reach.  There may be cases where we wouldn't even want to remove this node, I haven't examined all the possibilities.

**To Reproduce**
```python
from cspy import BiDirectional
import numpy as np
import networkx as nx
from networkx import DiGraph

# make a simple graph, with labeled resource usage
# B clearly superfluous
#            1       1
#   Source ---> A ----> Sink
#               ^ 
#            10 \---->B
H = DiGraph(n_res=1)
max_res = [5]
min_res = [1]
H.add_edge('Source', 'A', res_cost=np.array([1]), weight=1)
H.add_edge('A', 'B', res_cost=np.array([10]), weight=1)
H.add_edge('B', 'A', res_cost=np.array([10]), weight=1)
H.add_edge('A', 'Sink', res_cost=np.array([1]), weight=1)

# seems to works, but error is masked
bidirec = BiDirectional(H, max_res, min_res, direction='forward', seed=42, preprocess=False)
bidirec.run()
print(bidirec.path, bidirec.total_cost, bidirec.consumed_resources)
bidirec = BiDirectional(H, max_res, min_res, direction='forward', seed=42, preprocess=True)
bidirec.run()
print(bidirec.path, bidirec.total_cost, bidirec.consumed_resources)

length_s, path_s = nx.single_source_bellman_ford(H, 'Source', weight=getw)
length_t, path_t = nx.single_source_bellman_ford(H.reverse(copy=True), 'Sink', weight=getw)

nodes_source = {}
nodes_sink = {}

# B should be identified (or rather A because of B)
# But get IndexError
try:
    nodes_source.update({
	path_s[key][-2]: (val, r)
	for key, val in length_s.items()
	if val > max_res[r] or val < min_res[r]
    })
    nodes_sink.update({
	path_t[key][-2]: (val, r)
	for key, val in length_t.items()
	if val > max_res[r] or val < min_res[r]
    })
except IndexError:
    print(nodes_source, nodes_sink)

# Fix to check if key is already at Source or Sink
nodes_source.update({
	path_s[key][-2]: (val, r)                              # I think this index should be -1
	for key, val in length_s.items()
	if key != 'Source' and (val > max_res[r] or val < min_res[r])
    })
nodes_sink.update({
	path_t[key][-2]: (val, r)                              # I think this index should be -1
	for key, val in length_t.items()
	if key != 'Sink' and (val > max_res[r] or val < min_res[r])
    })
print(nodes_source, nodes_sink)

## But now node 'A' is slated to be removed, causing graph to be disconnected
H.remove('A')
bidirec = BiDirectional(H, max_res, min_res, direction='forward', seed=42, preprocess=False)
bidirec.run() # failure
```

**Expected behavior**
`preprocessing.py` should ignore "Source" or "Sink" itself when checking related length keys.  Further, the node to be removed should be the final node in the path, not the next to last, as that node may be required for other paths.  There may be a more theoretically robust preprocessing that identifies offending edges and removes those; any nodes that become disconnected could then be safely removed.

**Desktop (please complete the following information):**
 - OS: CentOS
 - `cspy` 1.0.0a0
 - `numpy` 1.19.2
 - `networkx` 2.5
 - Python 3.8.5

**Additional context**
This bug can happen in modest size (say 25) random graphs with reasonable probability.  I'd be happy to make a pull request if that's easiest, but I don't have a build system setup to run all the tests before submitting.



Discussion / hints:
Hey again @dvbuntu, thanks for bringing up another bug!

There are some preprocessing techniques but most of them use the source-to-all shortest paths which has to run twice (one for each direction) and then times the number of resources. So it can get really slow. 
Also, they are only effective on certain types of networks (e.g. doesn't work on connected graphs like the ones you get in VRP instances).

In a recent paper, link found in #67, has a much faster and elegant solution to this. They also have a crazy good implementation
Haven't gotten round to implementing it yet, but was definitely one of my priorities.

You are very welcome to open a PR with the fix for this :) the workflows will show if the test pass.
The build process is pretty straight forward if you set up docker, or alternatively install cmake on your machine.
Looking at the preprocessing more closely, I'm not sure that it's checking the `min_res` correctly.  It looks for minimum weight paths and checks if the shortest path (for a given resource) is less than the minimum constraint.  But if the shortest path is too short, there could  be some longer path which is valid.

Ideally, we'd check if the longest path is not long enough.  But we can't just make the weights negative and look for the shortest since that will likely introduce negative-weight cycles (for the original weights being positive).  But the [Longest Path Problem](https://en.wikipedia.org/wiki/Longest_path_problem) is NP-Hard, so we're not going to have a quick way to do this check.  We may want to just avoid checking the minimum constraint in the preprocessor.

Fix patch:
diff --git a/src/python/preprocessing.py b/src/python/preprocessing.py
index f22f595..ae688ed 100644
--- a/src/python/preprocessing.py
+++ b/src/python/preprocessing.py
@@ -12,12 +12,16 @@ def prune_graph(G, max_res, min_res):
 
     Note
     -----
-    path_s and path_t contains all partial paths e.g.
-        {node_i: [Source, ..., node_k, node_i]}.
-    or  {node_i: [Sink, ..., node_k, node_i]}
-    Hence, if node_i is found to violate a resource bound, it is
-    because of node_k, therefore, we add path[key][-2] = node_k to
-    the dictionary of nodes to remove.
+    If the shortest path from Source to node_i has cost greater
+    than the resource limit, then node_i cannot be on the shortest
+    resource-constrained path from Source to Sink.  So we remove
+    node_i.
+
+    We cannot, however, make determinations about the minimum
+    resource constraints.  We would need to find the longest path
+    (which is NP-Hard) and compare that.  Because while the shortest
+    path might not use enough of a resource, a longer path could exist
+    that does.
     """
 
     def _check_resource(r):
@@ -34,18 +38,19 @@ def __get_weight(i, j, attr_dict):
         length_t, path_t = single_source_bellman_ford(G.reverse(copy=True),
                                                       'Sink',
                                                       weight=__get_weight)
+
         try:
             # Collect nodes in paths that violate the resource bounds
             # see note above
             nodes_source.update({
-                path_s[key][-2]: (val, r)
+                path_s[key][-1]: (val, r)
                 for key, val in length_s.items()
-                if val > max_res[r] or val < min_res[r]
+                if (key != 'Source') and (val > max_res[r])
             })
             nodes_sink.update({
-                path_t[key][-2]: (val, r)
+                path_t[key][-1]: (val, r)
                 for key, val in length_t.items()
-                if val > max_res[r] or val < min_res[r]
+                if (key != 'Sink') and (val > max_res[r])
             })
         except IndexError:  # No nodes violate resource limits
             pass
@@ -55,11 +60,13 @@ def __get_weight(i, j, attr_dict):
     n_nodes = len(G.nodes())
     # Map function for each resource
     list(map(_check_resource, range(0, G.graph['n_res'])))
-    if nodes_source and nodes_sink:  # if there are nodes to remove
-        # Filter out source or sink
+    if nodes_source or nodes_sink:  # if there are nodes to remove
+        # Remove node if either direction is bad; it can't be on shortest path
         nodes_to_remove = [
-            node for node in G.nodes() if node in nodes_source and nodes_sink
+                node for node in G.nodes() if (node in nodes_source) or
+                (node in nodes_sink)
         ]
+        # Filter out source or sink
         if "Source" in nodes_to_remove or "Sink" in nodes_to_remove:
             if "Sink" in nodes_source:
                 unreachable_res = nodes_source["Sink"][1]


Test patch:
diff --git a/test/python/tests_preprocessing.py b/test/python/tests_preprocessing.py
index 314b151..16ff60b 100644
--- a/test/python/tests_preprocessing.py
+++ b/test/python/tests_preprocessing.py
@@ -14,13 +14,14 @@ class TestsPreprocessing(unittest.TestCase):
     """
 
     def setUp(self):
-        # Create digraph with negative resource costs with unreachable node 'B'
+        # Create digraph with negative resource costs with unreachable node 'E'
         self.G = DiGraph(directed=True, n_res=2)
         self.G.add_edge('Source', 'A', res_cost=array([1, 2]), weight=0)
         self.G.add_edge('A', 'C', res_cost=array([-1, 0.3]), weight=0)
         self.G.add_edge('A', 'B', res_cost=array([-1, 3]), weight=0)
-        # Unreachable node B
         self.G.add_edge('B', 'D', res_cost=array([-1, 2]), weight=0)
+        # Unreachable node E
+        self.G.add_edge('B', 'E', res_cost=array([10, 1]), weight=0)
         self.G.add_edge('C', 'D', res_cost=array([1, 0.1]), weight=0)
         self.G.add_edge('D', 'Sink', res_cost=array([1, 0.1]), weight=0)
 
@@ -32,12 +33,26 @@ def setUp(self):
         self.H.add_edge('Sink', 'Source', weight=-1)
         self.max_res, self.min_res = ["foo"], ["bar"]
 
+        # Create digraph to test issue 72
+        self.F = DiGraph(n_res=1)
+        self.F.add_edge('Source', 'A', res_cost=array([1]), weight=1)
+        self.F.add_edge('A', 'B', res_cost=array([10]), weight=1)
+        self.F.add_edge('A', 'Sink', res_cost=array([1]), weight=1)
+
+
     def testUnreachable(self):
         """
-        Tests if the unreachable node 'B' is removed.
+        Tests if the unreachable node 'E' is removed.
         """
         self.G = preprocess_graph(self.G, [5, 20], [0, 0], True)
-        self.assertTrue('B' not in self.G.nodes())
+        self.assertTrue('E' not in self.G.nodes())
+
+    def testFindBadNode(self):
+        """
+        Tests if expensive node 'B' is removed
+        """
+        self.F = preprocess_graph(self.F, [5], [1], True)
+        self.assertTrue('B' not in self.F.nodes())
 
     def testFormatting(self):
         """