gh-40284: PNC k shortest simple path (for directed graphs)
    
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->

Implement PNC k shortest simple path for directed graphs.
I implemented for only directed graphs for simplicity first.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->

#40248
#40217
    
URL: #40284
Reported by: Yuta Inoue
Reviewer(s): David Coudert

Author: Release Manager

src/doc/en/reference/references/index.rst

@@ -108,6 +108,10 @@ REFERENCES:
          *Coxeter submodular functions and deformations of Coxeter permutahedra*,
          Advances in Mathematics, Volume 365, 13 May 2020.
 
+.. [ACN2023] Ali Al Zoobi, David Coudert, Nicolas Nisse
+             Finding the k Shortest Simple Paths: Time and Space trade-offs
+             ACM Journal of Experimental Algorithmics, 2023, 28, pp.23 :doi:`10.1145/3626567`.
+
 .. [ALL2002] P. Auger, G. Labelle and P. Leroux, *Combinatorial
              addition formulas and applications*, Advances in Applied
              Mathematics 28 (2002) 302-342.

src/sage/graphs/path_enumeration.pyx

@@ -13,6 +13,7 @@ This module is meant for all functions related to path enumeration in graphs.
     :func:`all_paths` | Return the list of all paths between a pair of vertices.
     :func:`yen_k_shortest_simple_paths` | Return an iterator over the simple paths between a pair of vertices in increasing order of weights.
     :func:`feng_k_shortest_simple_paths` | Return an iterator over the simple paths between a pair of vertices in increasing order of weights.
+    :func:`pnc_k_shortest_simple_paths` | Return an iterator over the simple paths between a pair of vertices in increasing order of weights.
     :func:`all_paths_iterator` | Return an iterator over the paths of ``self``.
     :func:`all_simple_paths` | Return a list of all the simple paths of ``self`` starting with one of the given vertices.
     :func:`shortest_simple_paths` | Return an iterator over the simple paths between a pair of vertices.
@@ -1383,6 +1384,264 @@ def feng_k_shortest_simple_paths(self, source, target, weight_function=None,
             reduced_cost[e[0], e[1]] = temp_dict[e[0], e[1]]
 
 
+def pnc_k_shortest_simple_paths(self, source, target, weight_function=None,
+                                by_weight=False, check_weight=True,
+                                report_edges=False,
+                                labels=False, report_weight=False):
+    r"""
+    Return an iterator over the simple paths between a pair of vertices in
+    increasing order of weights.
+
+    Works only for directed graphs.
+
+    In case of weighted graphs, negative weights are not allowed.
+
+    If ``source`` is the same vertex as ``target``, then ``[[source]]`` is
+    returned -- a list containing the 1-vertex, 0-edge path ``source``.
+
+    The loops and the multiedges if present in the given graph are ignored and
+    only minimum of the edge labels is kept in case of multiedges.
+
+    INPUT:
+
+    - ``source`` -- a vertex of the graph, where to start
+
+    - ``target`` -- a vertex of the graph, where to end
+
+    - ``weight_function`` -- function (default: ``None``); a function that
+      takes as input an edge ``(u, v, l)`` and outputs its weight. If not
+      ``None``, ``by_weight`` is automatically set to ``True``. If ``None``
+      and ``by_weight`` is ``True``, we use the edge label ``l`` as a
+      weight.
+
+    - ``by_weight`` -- boolean (default: ``False``); if ``True``, the edges
+      in the graph are weighted, otherwise all edges have weight 1
+
+    - ``check_weight`` -- boolean (default: ``True``); whether to check that
+      the ``weight_function`` outputs a number for each edge
+
+    - ``report_edges`` -- boolean (default: ``False``); whether to report
+      paths as list of vertices (default) or list of edges, if ``False``
+      then ``labels`` parameter is ignored
+
+    - ``labels`` -- boolean (default: ``False``); if ``False``, each edge
+      is simply a pair ``(u, v)`` of vertices. Otherwise a list of edges
+      along with its edge labels are used to represent the path.
+
+    - ``report_weight`` -- boolean (default: ``False``); if ``False``, just
+      a path is returned. Otherwise a tuple of path length and path is
+      returned.
+
+    ALGORITHM:
+
+    This algorithm is based on the ``feng_k_shortest_simple_paths`` algorithm
+    in [Feng2014]_, but postpones the shortest path tree computation when non-simple
+    deviations occur. See Postponed Node Classification algorithm in [ACN2023]_
+    for the algorithm description.
+
+    EXAMPLES::
+
+        sage: from sage.graphs.path_enumeration import pnc_k_shortest_simple_paths
+        sage: g = DiGraph([(1, 2, 20), (1, 3, 10), (1, 4, 30), (2, 5, 20), (3, 5, 10), (4, 5, 30)])
+        sage: list(pnc_k_shortest_simple_paths(g, 1, 5, by_weight=True, report_weight=True))
+        [(20.0, [1, 3, 5]), (40.0, [1, 2, 5]), (60.0, [1, 4, 5])]
+        sage: list(pnc_k_shortest_simple_paths(g, 1, 5, report_weight=True))
+        [(2.0, [1, 2, 5]), (2.0, [1, 4, 5]), (2.0, [1, 3, 5])]
+
+    TESTS::
+
+        sage: from sage.graphs.path_enumeration import pnc_k_shortest_simple_paths
+        sage: g = DiGraph([(0, 1, 9), (0, 3, 1), (0, 4, 2), (1, 6, 4),
+        ....:              (1, 7, 1), (2, 0, 5), (2, 1, 4), (2, 7, 1),
+        ....:              (3, 1, 7), (3, 2, 4), (3, 4, 2), (4, 0, 8),
+        ....:              (4, 1, 10), (4, 3, 3), (4, 7, 10), (5, 2, 5),
+        ....:              (5, 4, 9), (6, 2, 9)], weighted=True)
+        sage: list(pnc_k_shortest_simple_paths(g, 5, 1, by_weight=True, report_weight=True,
+        ....:                                  labels=True, report_edges=True))
+        [(9.0, [(5, 2, 5), (2, 1, 4)]),
+         (18.0, [(5, 2, 5), (2, 0, 5), (0, 3, 1), (3, 1, 7)]),
+         (19.0, [(5, 2, 5), (2, 0, 5), (0, 1, 9)]),
+         (19.0, [(5, 4, 9), (4, 1, 10)]),
+         (19.0, [(5, 4, 9), (4, 3, 3), (3, 1, 7)]),
+         (20.0, [(5, 4, 9), (4, 3, 3), (3, 2, 4), (2, 1, 4)]),
+         (22.0, [(5, 2, 5), (2, 0, 5), (0, 4, 2), (4, 1, 10)]),
+         (22.0, [(5, 2, 5), (2, 0, 5), (0, 4, 2), (4, 3, 3), (3, 1, 7)]),
+         (23.0, [(5, 2, 5), (2, 0, 5), (0, 3, 1), (3, 4, 2), (4, 1, 10)]),
+         (25.0, [(5, 4, 9), (4, 0, 8), (0, 3, 1), (3, 1, 7)]),
+         (26.0, [(5, 4, 9), (4, 0, 8), (0, 1, 9)]),
+         (26.0, [(5, 4, 9), (4, 0, 8), (0, 3, 1), (3, 2, 4), (2, 1, 4)]),
+         (30.0, [(5, 4, 9), (4, 3, 3), (3, 2, 4), (2, 0, 5), (0, 1, 9)])]
+        sage: g = DiGraph(graphs.Grid2dGraph(2, 6).relabel(inplace=False))
+        sage: for u, v in g.edge_iterator(labels=False):
+        ....:     g.set_edge_label(u, v, 1)
+        sage: [w for  w, P in pnc_k_shortest_simple_paths(g, 5, 1, by_weight=True, report_weight=True)]
+        [4.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 8.0, 8.0,
+         8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 10.0, 10.0, 10.0, 10.0]
+
+    Same tests as ``yen_k_shortest_simple_paths``::
+
+        sage: g = DiGraph([(1, 2, 1), (2, 3, 1), (3, 4, 1), (4, 5, 1),
+        ....:              (1, 7, 1), (7, 8, 1), (8, 5, 1), (1, 6, 1),
+        ....:              (6, 9, 1), (9, 5, 1), (4, 2, 1), (9, 3, 1),
+        ....:              (9, 10, 1), (10, 5, 1), (9, 11, 1), (11, 10, 1)])
+        sage: [w for w, P in pnc_k_shortest_simple_paths(g, 1, 5, by_weight=True, report_weight=True)]
+        [3.0, 3.0, 4.0, 4.0, 5.0, 5.0]
+
+    More tests::
+
+        sage: D = graphs.Grid2dGraph(5, 5).relabel(inplace=False).to_directed()
+        sage: A = [w for w, P in pnc_k_shortest_simple_paths(D, 0, 24, report_weight=True)]
+        sage: assert len(A) == 8512
+        sage: for i in range(len(A) - 1):
+        ....:     assert A[i] <= A[i + 1]
+    """
+    if not self.is_directed():
+        raise ValueError("this algorithm works only for directed graphs")
+
+    if source not in self:
+        raise ValueError("vertex '{}' is not in the graph".format(source))
+    if target not in self:
+        raise ValueError("vertex '{}' is not in the graph".format(target))
+    if source == target:
+        P = [] if report_edges else [source]
+        yield (0, P) if report_weight else P
+        return
+
+    if self.has_loops() or self.allows_multiple_edges():
+        G = self.to_simple(to_undirected=False, keep_label='min', immutable=False)
+    else:
+        G = self.copy(immutable=False)
+
+    G.delete_edges(G.incoming_edges(source, labels=False))
+    G.delete_edges(G.outgoing_edges(target, labels=False))
+
+    by_weight, weight_function = G._get_weight_function(by_weight=by_weight,
+                                                        weight_function=weight_function,
+                                                        check_weight=check_weight)
+
+    def reverse_weight_function(e):
+        return weight_function((e[1], e[0], e[2]))
+
+    cdef dict edge_labels
+    edge_labels = {(e[0], e[1]): e for e in G.edge_iterator()}
+
+    cdef dict edge_wt
+    edge_wt = {(e[0], e[1]): weight_function(e) for e in G.edge_iterator()}
+
+    # The first shortest path tree T_0
+    from sage.graphs.base.boost_graph import shortest_paths
+    cdef dict dist
+    cdef dict successor
+    reverse_graph = G.reverse()
+    dist, successor = shortest_paths(reverse_graph, target, weight_function=reverse_weight_function,
+                                     algorithm='Dijkstra_Boost')
+    cdef set unnecessary_vertices = set(G) - set(dist)  # no path to target
+    if source in unnecessary_vertices:  # no path from source to target
+        return
+
+    # sidetrack cost
+    cdef dict sidetrack_cost = {(e[0], e[1]): weight_function(e) + dist[e[1]] - dist[e[0]]
+                                for e in G.edge_iterator()
+                                if e[0] in dist and e[1] in dist}
+
+    def sidetrack_length(path):
+        return sum(sidetrack_cost[e] for e in zip(path, path[1:]))
+
+    # v-t path in the first shortest path tree T_0
+    def tree_path(v):
+        path = [v]
+        while v != target:
+            v = successor[v]
+            path.append(v)
+        return path
+
+    # shortest path
+    shortest_path = tree_path(source)
+    cdef double shortest_path_length = dist[source]
+
+    # idx of paths
+    cdef dict idx_to_path = {0: shortest_path}
+    cdef int idx = 1
+
+    # candidate_paths collects (cost, path_idx, dev_idx, is_simple)
+    # + cost is sidetrack cost from the first shortest path tree T_0
+    #   (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.bidirectional_dijkstra_special
+
+    candidate_paths.push(((0, True), (0, 0)))
+    while candidate_paths.size():
+        (negative_cost, is_simple), (path_idx, dev_idx) = candidate_paths.top()
+        cost = -negative_cost
+        candidate_paths.pop()
+
+        path = idx_to_path[path_idx]
+        del idx_to_path[path_idx]
+
+        # ancestor_idx_dict[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]``.
+        ancestor_idx_dict = {v: i for i, v in enumerate(path)}
+
+        def ancestor_idx_func(v, t, len_path):
+            if v not in successor:
+                # target vertex is not reachable from v
+                return -1
+            if v in ancestor_idx_dict:
+                if ancestor_idx_dict[v] <= t or ancestor_idx_dict[v] == len_path - 1:
+                    return ancestor_idx_dict[v]
+            ancestor_idx_dict[v] = ancestor_idx_func(successor[v], t, len_path)
+            return ancestor_idx_dict[v]
+
+        if is_simple:
+            # output
+            if report_edges and labels:
+                P = [edge_labels[e] for e in zip(path, path[1:])]
+            elif report_edges:
+                P = list(zip(path, path[1:]))
+            else:
+                P = path
+            if report_weight:
+                yield (shortest_path_length + cost, P)
+            else:
+                yield P
+
+            # GET DEVIATION PATHS
+            original_cost = cost
+            for deviation_i in range(len(path) - 1, dev_idx - 1, -1):
+                for e in G.outgoing_edge_iterator(path[deviation_i]):
+                    if e[1] in path[:deviation_i + 2]:  # 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
+                    new_path = path[:deviation_i + 1] + tree_path(e[1])
+                    new_path_idx = idx
+                    idx_to_path[new_path_idx] = new_path
+                    idx += 1
+                    new_cost = original_cost + sidetrack_cost[(e[0], e[1])]
+                    new_is_simple = ancestor_idx > deviation_i
+                    candidate_paths.push(((-new_cost, new_is_simple), (new_path_idx, deviation_i + 1)))
+                if deviation_i == dev_idx:
+                    continue
+                original_cost -= sidetrack_cost[(path[deviation_i - 1], path[deviation_i])]
+        else:
+            # get a path to target in G \ path[:dev_idx]
+            deviation = shortest_path_func(path[dev_idx], target,
+                                           exclude_vertices=unnecessary_vertices.union(path[:dev_idx]),
+                                           reduced_weight=sidetrack_cost)
+            if not deviation:
+                continue  # no path to target in G \ path[:dev_idx]
+            new_path = path[:dev_idx] + deviation
+            new_path_idx = idx
+            idx_to_path[new_path_idx] = new_path
+            idx += 1
+            new_cost = sidetrack_length(new_path)
+            candidate_paths.push(((-new_cost, True), (new_path_idx, dev_idx)))
+
+
 def _all_paths_iterator(self, vertex, ending_vertices=None,
                         simple=False, max_length=None, trivial=False,
                         use_multiedges=False, report_edges=False,