Improve PNC algorithm

src/sage/graphs/path_enumeration.pyx

@@ -37,7 +37,10 @@ from sage.misc.misc_c import prod
 from libcpp.queue cimport priority_queue
 from libcpp.pair cimport pair
 from libcpp.vector cimport vector
+
+from sage.data_structures.pairing_heap cimport PairingHeap
 from sage.rings.integer_ring import ZZ
+
 import copy
 
 
@@ -1546,10 +1549,6 @@ def pnc_k_shortest_simple_paths(self, source, target, weight_function=None,
     reverse_graph = G.reverse()
     dist, successor = shortest_paths(reverse_graph, id_target, weight_function=reverse_weight_function,
                                      algorithm='Dijkstra_Boost')
-    cdef dict predecessor = {v: [] for v in dist}
-    for v, u in successor.items():
-        if u is not None:
-            predecessor[u].append(v)
     cdef set unnecessary_vertices = set(G) - set(dist)  # no path to target
     if id_source in unnecessary_vertices:  # no path from source to target
         return
@@ -1580,40 +1579,65 @@ def pnc_k_shortest_simple_paths(self, source, target, weight_function=None,
     #   (i.e. real length = cost + shortest_path_length in T_0)
     cdef priority_queue[pair[pair[double, bint], pair[int, int]]] candidate_paths
 
-    # shortest path function for weighted/unweighted graph using reduced weights
-    shortest_path_func = G._backend.shortest_path_to_set
-
-    # See explanation of ancestor_idx_vec below
+    # ancestor_idx_vec[v] := the first vertex of ``path[:t+1]`` or ``id_target`` reachable by
+    #                    edges of first shortest path tree from v.
     cdef vector[int] ancestor_idx_vec = [-1 for _ in range(len(G))]
 
-    def ancestor_idx_func(v, t, len_path):
-        if v not in successor:
-            # target vertex is not reachable from v
-            return -1
+    def ancestor_idx_func(v, t, target_idx):
         if ancestor_idx_vec[v] != -1:
-            if ancestor_idx_vec[v] <= t or ancestor_idx_vec[v] == len_path - 1:
+            if ancestor_idx_vec[v] <= t or ancestor_idx_vec[v] == target_idx:
                 return ancestor_idx_vec[v]
-        ancestor_idx_vec[v] = ancestor_idx_func(successor[v], t, len_path)
+        ancestor_idx_vec[v] = ancestor_idx_func(successor[v], t, target_idx)
         return ancestor_idx_vec[v]
 
-    def green_vertices(path):
-        # Find all vertices that are reachable to target in T_0
-        # without crossing vertices of path
-        green = set()
-        in_path = [False for _ in range(len(G))]
-        for v in path:
-            in_path[v] = True
-
-        def dfs(v):
-            if in_path[v]:
-                return
-            green.add(v)
-            for u in predecessor[v]:
-                dfs(u)
-            return
-        dfs(id_target)
-        return green
+    # used inside shortest_path_func
+    cdef PairingHeap[int, double] pq = PairingHeap[int, double]()
+    cdef dict dist_in_func = {}
+    cdef dict pred = {}
+
+    def shortest_path_func(dev, exclude_vertices):
+        t = len(exclude_vertices)
+        ancestor_idx_vec[id_target] = t + 1
+        # clear
+        while not pq.empty():
+            pq.pop()
+        dist_in_func.clear()
+        pred.clear()
+
+        pq.push(dev, 0)
+        dist_in_func[dev] = 0
+
+        while not pq.empty():
+            v, d = pq.top()
+            pq.pop()
+
+            if ancestor_idx_func(v, t, t + 1) == t + 1:  # green
+                path = []
+                while v in pred:
+                    path.append(v)
+                    v = pred[v]
+                path.append(dev)
+                path.reverse()
+                return (d, path)
+
+            if d > dist_in_func.get(v, float('inf')):
+                continue  # already found a better path
+
+            for _, u, _ in G.outgoing_edge_iterator(v):
+                if u in exclude_vertices:
+                    continue
+                new_dist = d + sidetrack_cost[(v, u)]
+                if new_dist < dist_in_func.get(u, float('inf')):
+                    dist_in_func[u] = new_dist
+                    pred[u] = v
+                    if pq.contains(u):
+                        if pq.value(u) > new_dist:
+                            pq.decrease(u, new_dist)
+                    else:
+                        pq.push(u, new_dist)
+        return
 
+    cdef int i, deviation_i
     candidate_paths.push(((0, True), (0, 0)))
     while candidate_paths.size():
         (negative_cost, is_simple), (path_idx, dev_idx) = candidate_paths.top()
@@ -1623,6 +1647,11 @@ def pnc_k_shortest_simple_paths(self, source, target, weight_function=None,
         path = idx_to_path[path_idx]
         del idx_to_path[path_idx]
 
+        for i in range(ancestor_idx_vec.size()):
+            ancestor_idx_vec[i] = -1
+        for i, v in enumerate(path):
+            ancestor_idx_vec[v] = i
+
         if is_simple:
             # output
             if report_edges and labels:
@@ -1636,24 +1665,14 @@ def pnc_k_shortest_simple_paths(self, source, target, weight_function=None,
             else:
                 yield P
 
-            # ancestor_idx_vec[v] := the first vertex of ``path[:t+1]`` or ``path[-1]`` reachable by
-            #                    edges of first shortest path tree from v when enumerating deviating edges
-            #                    from ``path[t]``.
-            for i in range(ancestor_idx_vec.size()):
-                ancestor_idx_vec[i] = -1
-            for i, v in enumerate(path):
-                ancestor_idx_vec[v] = i
-
             # GET DEVIATION PATHS
             original_cost = cost
             former_part = set(path)
             for deviation_i in range(len(path) - 2, dev_idx - 1, -1):
                 for e in G.outgoing_edge_iterator(path[deviation_i]):
                     if e[1] in former_part:  # e[1] is red or e in path
                         continue
-                    ancestor_idx = ancestor_idx_func(e[1], deviation_i, len(path))
-                    if ancestor_idx == -1:
-                        continue
+                    ancestor_idx = ancestor_idx_func(e[1], deviation_i, len(path) - 1)
                     new_is_simple = ancestor_idx > deviation_i
                     # no need to compute tree_path if new_is_simple is False
                     new_path = path[:deviation_i + 1] + (tree_path(e[1]) if new_is_simple else [e[1]])
@@ -1667,18 +1686,10 @@ def pnc_k_shortest_simple_paths(self, source, target, weight_function=None,
                 original_cost -= sidetrack_cost[(path[deviation_i - 1], path[deviation_i])]
                 former_part.remove(path[deviation_i + 1])
         else:
-            green = green_vertices(path[:dev_idx + 1])
-            # get a path to target in G \ path[:dev_idx] to one of green vertices
-            try:
-                deviation_weight, deviation = shortest_path_func(path[dev_idx], green,
-                                                                 report_weight=True,
-                                                                 exclude_vertices=set(path[:dev_idx]),
-                                                                 edge_weight=sidetrack_cost)
-            except ValueError as e:
-                if str(e) == "no path found from source to targets.":
-                    continue  # no path to target in G \ path[:dev_idx]
-                else:
-                    raise
+            deviations = shortest_path_func(path[dev_idx], set(path[:dev_idx]))
+            if not deviations:
+                continue  # no path to target in G \ path[:dev_idx]
+            deviation_weight, deviation = deviations
             new_path = path[:dev_idx] + deviation[:-1] + tree_path(deviation[-1])
             new_path_idx = idx
             idx_to_path[new_path_idx] = new_path