diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index 42d3ae2d886..747bb35033e 100644 --- a/src/doc/en/reference/references/index.rst +++ b/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. diff --git a/src/sage/graphs/path_enumeration.pyx b/src/sage/graphs/path_enumeration.pyx index 52643122e56..a592d0c9b54 100644 --- a/src/sage/graphs/path_enumeration.pyx +++ b/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,