faster bidirectional dijkstra using pairing heap

src/sage/data_structures/pairing_heap.pxd

@@ -30,6 +30,7 @@ cdef extern from "./pairing_heap.h" namespace "pairing_heap":
         void decrease(TypeOfItem, TypeOfValue) except +
         bint contains(TypeOfItem)
         TypeOfValue value(TypeOfItem) except +
+        size_t size()
 
     cdef cppclass PairingHeapNodePy:
         PyObject * value           # value associated with the item

src/sage/graphs/base/c_graph.pyx

@@ -52,6 +52,7 @@ from libcpp.pair cimport pair
 from sage.rings.integer_ring import ZZ
 from cysignals.memory cimport check_allocarray, sig_free
 from sage.data_structures.bitset cimport FrozenBitset
+from sage.data_structures.pairing_heap cimport PairingHeap
 
 
 cdef extern from "Python.h":
@@ -3899,7 +3900,7 @@ cdef class CGraphBackend(GenericGraphBackend):
 
             sage: G = Graph(graphs.PetersenGraph())
             sage: for (u, v) in G.edges(sort=True, labels=None):
-            ....:    G.set_edge_label(u, v, 1)
+            ....:     G.set_edge_label(u, v, 1)
             sage: G.shortest_path(0, 1, by_weight=True)
             [0, 1]
             sage: G.shortest_path_length(0, 1, by_weight=True)
@@ -3928,7 +3929,7 @@ cdef class CGraphBackend(GenericGraphBackend):
 
             sage: G = DiGraph({0: [1, 2], 1: [4], 2: [3, 4], 4: [5], 5: [6]}, multiedges=True)
             sage: for u, v in list(G.edges(labels=None, sort=False)):
-            ....:    G.set_edge_label(u, v, 1)
+            ....:     G.set_edge_label(u, v, 1)
             sage: G.distance(0, 5, by_weight=true)
             3
         """
@@ -3954,119 +3955,123 @@ cdef class CGraphBackend(GenericGraphBackend):
         cdef dict pred_x = {}
         cdef dict pred_y = {}
         cdef dict pred_current
-        cdef dict pred_other
 
         # Stores the distances from x and y
         cdef dict dist_x = {}
         cdef dict dist_y = {}
         cdef dict dist_current
         cdef dict dist_other
 
-        # Lists of vertices who are left to be explored. They are represented
-        # as pairs of pair and pair: ((distance, side), (predecessor, name)).
-        # 1 indicates x's side, -1 indicates y's, the distance being
-        # defined relatively.
-        cdef priority_queue[pair[pair[double, int], pair[int, int]]] pq
-        pq.push(((0, 1), (x_int, x_int)))
-        pq.push(((0, -1), (y_int, y_int)))
-        cdef list neighbors
+        # We use 2 min-heap data structures (pairing heaps), one for the
+        # exploration from x and the other for the reverse exploration to y.
+        # Each heap associates to a vertex a pair (distance, pred).
+        cdef PairingHeap[int, pair[double, int]] px = PairingHeap[int, pair[double, int]]()
+        cdef PairingHeap[int, pair[double, int]] py = PairingHeap[int, pair[double, int]]()
+        cdef PairingHeap[int, pair[double, int]] * ptmp
+        px.push(x_int, (0, x_int))
+        py.push(y_int, (0, y_int))
 
-        cdef list shortest_path = []
+        cdef list neighbors
 
         # Meeting_vertex is a vertex discovered through x and through y
         # which defines the shortest path found
         # (of length shortest_path_length).
         cdef int meeting_vertex = -1
+        cdef double shortest_path_length
+        cdef double f_tmp
 
         if weight_function is None:
             def weight_function(e):
                 return 1 if e[2] is None else e[2]
 
         # As long as the current side (x or y) is not totally explored ...
-        while not pq.empty():
-            (distance, side), (pred, v) = pq.top()
-            # priority_queue by default is max heap
-            # negative value of distance is stored in priority_queue to get
-            # minimum distance
-            distance = -distance
-            pq.pop()
+        while not (px.empty() and py.empty()):
+            if (px.empty() or
+                    (not py.empty() and px.top_value().first > py.top_value().first)):
+                side = -1
+                ptmp = &py
+            else:  # px is not empty
+                side = 1
+                ptmp = &px
+            v, (distance, pred) = ptmp.top()
             if meeting_vertex != -1 and distance > shortest_path_length:
                 break
+            ptmp.pop()
 
             if side == 1:
                 dist_current, dist_other = dist_x, dist_y
-                pred_current, pred_other = pred_x, pred_y
+                pred_current = pred_x
+                neighbors = self.cg().out_neighbors(v)
             else:
                 dist_current, dist_other = dist_y, dist_x
-                pred_current, pred_other = pred_y, pred_x
-
-            if v not in dist_current:
-                if not distance_flag:
-                    pred_current[v] = pred
-                dist_current[v] = distance
-
-                if v in dist_other:
-                    f_tmp = distance + dist_other[v]
-                    if meeting_vertex == -1 or f_tmp < shortest_path_length:
-                        meeting_vertex = v
-                        shortest_path_length = f_tmp
-
-                if side == 1:
-                    neighbors = self.cg().out_neighbors(v)
-                else:
-                    neighbors = self.cg().in_neighbors(v)
-                for w in neighbors:
-                    # If the neighbor is new, adds its non-found neighbors to
-                    # the queue.
-                    if w not in dist_current:
-                        v_obj = self.vertex_label(v)
-                        w_obj = self.vertex_label(w)
-                        if side == -1:
-                            v_obj, w_obj = w_obj, v_obj
-                        if self._multiple_edges:
-                            edge_label = min(weight_function((v_obj, w_obj, l)) for l in self.get_edge_label(v_obj, w_obj))
-                        else:
-                            edge_label = weight_function((v_obj, w_obj, self.get_edge_label(v_obj, w_obj)))
-                        if edge_label < 0:
-                            raise ValueError("the graph contains an edge with negative weight")
-                        # priority_queue is by default max_heap
-                        # negative value of distance + edge_label is stored in
-                        # priority_queue to get minimum distance
-                        pq.push(((-(distance + edge_label), side), (v, w)))
+                pred_current = pred_y
+                neighbors = self.cg().in_neighbors(v)
+
+            dist_current[v] = distance
+            if not distance_flag:
+                pred_current[v] = pred
+
+            if v in dist_other:
+                f_tmp = distance + dist_other[v]
+                if meeting_vertex == -1 or f_tmp < shortest_path_length:
+                    meeting_vertex = v
+                    shortest_path_length = f_tmp
+
+            for w in neighbors:
+                # If w has not yet been extracted from the heap, we check if we
+                # can improve its path
+                if w not in dist_current:
+                    v_obj = self.vertex_label(v)
+                    w_obj = self.vertex_label(w)
+                    if side == -1:
+                        v_obj, w_obj = w_obj, v_obj
+                    if self._multiple_edges:
+                        edge_label = min(weight_function((v_obj, w_obj, l)) for l in self.get_edge_label(v_obj, w_obj))
+                    else:
+                        edge_label = weight_function((v_obj, w_obj, self.get_edge_label(v_obj, w_obj)))
+                    if edge_label < 0:
+                        raise ValueError("the graph contains an edge with negative weight")
+                    f_tmp = distance + edge_label
+                    if ptmp.contains(w):
+                        if ptmp.value(w).first > f_tmp:
+                            ptmp.decrease(w, (f_tmp, v))
+                    else:
+                        ptmp.push(w, (f_tmp, v))
 
         # No meeting point has been found
         if meeting_vertex == -1:
             if distance_flag:
                 from sage.rings.infinity import Infinity
                 return Infinity
             return []
-        else:
-            # build the shortest path and returns it.
-            if distance_flag:
-                if shortest_path_length in ZZ:
-                    return int(shortest_path_length)
-                else:
-                    return shortest_path_length
-            w = meeting_vertex
-
-            while w != x_int:
-                shortest_path.append(self.vertex_label(w))
-                w = pred_x[w]
 
-            shortest_path.append(x)
-            shortest_path.reverse()
+        if distance_flag:
+            if shortest_path_length in ZZ:
+                return int(shortest_path_length)
+            return shortest_path_length
 
-            if meeting_vertex == y_int:
-                return shortest_path
+        # build the shortest path and returns it.
+        cdef list shortest_path = []
+        w = meeting_vertex
+        while w != x_int:
+            shortest_path.append(self.vertex_label(w))
+            w = pred_x[w]
 
-            w = pred_y[meeting_vertex]
-            while w != y_int:
-                shortest_path.append(self.vertex_label(w))
-                w = pred_y[w]
-            shortest_path.append(y)
+        shortest_path.append(x)
+        shortest_path.reverse()
 
+        if meeting_vertex == y_int:
             return shortest_path
 
+        w = pred_y[meeting_vertex]
+        while w != y_int:
+            shortest_path.append(self.vertex_label(w))
+            w = pred_y[w]
+
+        shortest_path.append(y)
+
+        return shortest_path
+
     def shortest_path_all_vertices(self, v, cutoff=None,
                                    distance_flag=False):
         r"""

src/sage/graphs/generic_graph.py

@@ -17350,7 +17350,7 @@ def shortest_path(self, u, v, by_weight=False, algorithm=None,
             sage: D.shortest_path(4, 8, algorithm='Dijkstra_Bid_NetworkX')              # needs networkx
             [4, 3, 2, 1, 8]
             sage: D.shortest_path(4, 9, algorithm='Dijkstra_Bid')
-            [4, 3, 19, 0, 10, 9]
+            [4, 3, 2, 1, 8, 9]
             sage: D.shortest_path(5, 5)
             [5]
             sage: D.delete_edges(D.edges_incident(13))