Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
dag_constrained_shortest_path.h
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14#ifndef OR_TOOLS_GRAPH_DAG_CONSTRAINED_SHORTEST_PATH_H_
15#define OR_TOOLS_GRAPH_DAG_CONSTRAINED_SHORTEST_PATH_H_
16
17#include <cmath>
18#include <limits>
19#include <optional>
20#include <vector>
21
22#include "absl/algorithm/container.h"
23#include "absl/base/log_severity.h"
24#include "absl/log/check.h"
25#include "absl/strings/str_format.h"
26#include "absl/types/span.h"
29#include "ortools/graph/graph.h"
30
31namespace operations_research {
32
33// This library provides APIs to compute the constrained shortest path (CSP) on
34// a given directed acyclic graph (DAG) with resources on each arc. A CSP is a
35// shortest path on a DAG which does not exceed a set of maximum resources
36// consumption. The algorithm is exponential and has no guarantee to finish. It
37// is based on bi-drectionnal search. First is a forward pass from the source to
38// nodes “somewhere in the middle” to generate forward labels, just as the
39// onedirectional labeling algorithm we discussed; then a symmetric backward
40// pass from the destination generates backward labels; and finally at each node
41// with both forward and backward labels, it joins any pair of labels to form a
42// feasible complete path. Intuitively, the number of labels grows exponentially
43// with the number of arcs in the path. The overall number of labels are then
44// expected to be smaller with shorter paths. For DAG with a topological
45// ordering, we can pick any node (usually right in the middle) as a *midpoint*
46// to stop each pass at. Then labels can be joined at only one half of the nodes
47// by considering all edges between each half.
48//
49// In the DAG, multiple arcs between the same pair of nodes is allowed. However,
50// self-loop arcs are not allowed.
51//
52// Note that we use the length formalism here, but the arc lengths can represent
53// any numeric physical quantity. A shortest path will just be a path minimizing
54// this quantity where the length/resources of a path is the sum of the
55// length/resources of its arcs. An arc length can be negative, or +inf
56// (indicating that it should not be used). An arc length cannot be -inf or nan.
57//
58// Resources on each arc must be non-negative and cannot be +inf or nan.
59
60// -----------------------------------------------------------------------------
61// Basic API.
62// -----------------------------------------------------------------------------
63
64// `tail` and `head` should both be in [0, num_nodes)
65// If the length is +inf, then the arc is not used.
67 int from = 0;
68 int to = 0;
69 double length = 0.0;
70 std::vector<double> resources;
71};
73// Returns {+inf, {}, {}} if there is no path of finite length from the source
74// to the destination. Dies if `arcs_with_length_and_resources` has a cycle.
76 int num_nodes,
77 absl::Span<const ArcWithLengthAndResources> arcs_with_length_and_resources,
78 int source, int destination, const std::vector<double>& max_resources);
79
80// -----------------------------------------------------------------------------
81// Advanced API.
82// -----------------------------------------------------------------------------
83template <class GraphType>
85 double length = 0.0;
86 // The returned arc indices points into the `arcs_with_length` passed to the
87 // function below.
88 std::vector<typename GraphType::ArcIndex> arc_path;
89 std::vector<typename GraphType::NodeIndex>
90 node_path; // includes the source node.
91};
93// A wrapper that holds the memory needed to run many constrained shortest path
94// computations efficiently on the given DAG (on which resources do not change).
95// `GraphType` can use one of the interfaces defined in `util/graph/graph.h`.
96template <class GraphType>
98 public:
99 using NodeIndex = typename GraphType::NodeIndex;
100 using ArcIndex = typename GraphType::ArcIndex;
102 // IMPORTANT: All arguments must outlive the class.
103 //
104 // The vectors of `arc_lengths` and `arc_resources[i]` (for all resource i)
105 // *must* be of size `graph.num_arcs()` and indexed the same way as in
106 // `graph`. The vector `arc_resources` and `max_resources` *must* be of same
107 // size.
108 //
109 // You *must* provide a topological order. You can use
110 // `util::graph::FastTopologicalSort(graph)` to compute one if you don't
111 // already have one. An invalid topological order results in an upper bound
112 // for all shortest path computations. For maximum performance, you can
113 // further reindex the nodes under the topological order so that the memory
114 // access pattern is generally forward instead of random. For example, if the
115 // topological order for a graph with 4 nodes is [2,1,0,3], you can re-label
116 // the nodes 2, 1, and 0 to 0, 1, and 2 (and updates arcs accordingly).
117 //
118 // Validity of arcs and topological order are DCHECKed.
119 //
120 // If the number of labels in memory exceeds `max_num_created_labels / 2` at
121 // any point in each pass of the algorithm, new labels are not generated
122 // anymore and it returns the best path found so far, most particularly the
123 // empty path if none were found.
124 //
125 // IMPORTANT: You cannot modify anything except `arc_lengths` between calls to
126 // the `RunConstrainedShortestPathOnDag()` function.
128 const GraphType* graph, const std::vector<double>* arc_lengths,
129 const std::vector<std::vector<double>>* arc_resources,
130 absl::Span<const NodeIndex> topological_order,
131 absl::Span<const NodeIndex> sources,
132 absl::Span<const NodeIndex> destinations,
133 const std::vector<double>* max_resources,
134 int max_num_created_labels = 1e9);
135
136 // Returns {+inf, {}, {}} if there is no constrained path of finite length
137 // wihtin resources constraints from one node in `sources` to one node in
138 // `destinations`.
140
141 // For benchmarking and informational purposes, returns the number of labels
142 // generated in the call of `RunConstrainedShortestPathOnDag()`.
143 int label_count() const {
144 return lengths_from_sources_[FORWARD].size() +
145 lengths_from_sources_[BACKWARD].size();
146 }
148 private:
149 enum Direction {
150 FORWARD = 0,
151 BACKWARD = 1,
152 };
153
154 inline static Direction Reverse(Direction d) {
155 return d == FORWARD ? BACKWARD : FORWARD;
156 }
157
158 // A LabelPair includes the `length` of a path that can be constructed by
159 // merging the paths from two *linkable* labels corresponding to
160 // `label_index`.
161 struct LabelPair {
162 double length = 0.0;
163 int label_index[2];
164 };
165
166 void RunHalfConstrainedShortestPathOnDag(
167 const GraphType& reverse_graph, absl::Span<const double> arc_lengths,
168 absl::Span<const std::vector<double>> arc_resources,
169 absl::Span<const std::vector<double>> min_arc_resources,
170 absl::Span<const double> max_resources, int max_num_created_labels,
171 std::vector<double>& lengths_from_sources,
172 std::vector<std::vector<double>>& resources_from_sources,
173 std::vector<ArcIndex>& incoming_arc_indices_from_sources,
174 std::vector<int>& incoming_label_indices_from_sources,
175 std::vector<int>& first_label, std::vector<int>& num_labels);
176
177 // Returns the arc index linking two nodes from each pass forming the best
178 // path. Returns -1 if no better path than the one found from
179 // `best_label_pair` is found.
180 ArcIndex MergeHalfRuns(
181 const GraphType& graph, absl::Span<const double> arc_lengths,
182 absl::Span<const std::vector<double>> arc_resources,
183 absl::Span<const double> max_resources,
184 const std::vector<NodeIndex> sub_node_indices[2],
185 const std::vector<double> lengths_from_sources[2],
186 const std::vector<std::vector<double>> resources_from_sources[2],
187 const std::vector<int> first_label[2],
188 const std::vector<int> num_labels[2], LabelPair& best_label_pair);
189
190 // Returns the path as list of arc indices that starts from a node in
191 // `sources` (if `direction` iS FORWARD) or `destinations` (if `direction` is
192 // BACKWARD) and ends in node represented by `best_label_index`.
193 std::vector<ArcIndex> ArcPathTo(
194 int best_label_index,
195 absl::Span<const ArcIndex> incoming_arc_indices_from_sources,
196 absl::Span<const int> incoming_label_indices_from_sources) const;
197
198 // Returns the list of all the nodes implied by a given `arc_path`.
199 std::vector<NodeIndex> NodePathImpliedBy(absl::Span<const ArcIndex> arc_path,
200 const GraphType& graph) const;
201
202 static constexpr double kTolerance = 1e-6;
203
204 const GraphType* const graph_;
205 const std::vector<double>* const arc_lengths_;
206 const std::vector<std::vector<double>>* const arc_resources_;
207 const std::vector<double>* const max_resources_;
208 absl::Span<const NodeIndex> sources_;
209 absl::Span<const NodeIndex> destinations_;
210 const int num_resources_;
211
212 // Set to a node if and only if this node is in both `sources_` and
213 // `destinations_`.
214 std::optional<NodeIndex> source_is_destination_ = std::nullopt;
215 // Data about *reachable* sub-graphs split in two for bidirectional search.
216 // Reachable nodes are nodes that can be reached given the resources
217 // constraints, i.e., for each resource, the sum of the minimum resource to
218 // get to a node from a node in `sources` and to get from a node to a node in
219 // `destinations` should be less than the maximum resource. Reachable arcs are
220 // arcs linking reachable nodes.
221 //
222 // `sub_reverse_graph_[dir]` is the reachable sub-graph split in *half* with
223 // an additional linked to sources (resp. destinations) for the forward (resp.
224 // backward) direction. For the forward (resp. backward) direction, nodes are
225 // indexed using the original (resp. reverse) topological order.
226 GraphType sub_reverse_graph_[2];
227 std::vector<std::vector<double>> sub_arc_resources_[2];
228 // `sub_full_arc_indices_[dir]` has size `sub_reverse_graph_[dir].num_arcs()`
229 // such that `sub_full_arc_indices_[dir][sub_arc] = arc` where `sub_arc` is
230 // the arc in the reachable sub-graph for direction `dir` (i.e.
231 // `sub_reverse_graph[dir]`) and `arc` is the arc in the original graph (i.e.
232 // `graph`).
233 std::vector<NodeIndex> sub_full_arc_indices_[2];
234 // `sub_node_indices_[dir]` has size `graph->num_nodes()` such that
235 // `sub_node_indices[dir][node] = sub_node` where `node` is the node in the
236 // original graph (i.e. `graph`) and `sub_node` is the node in the reachable
237 // sub-graph for direction `dir` (i.e. `sub_reverse_graph[dir]`) and -1 if
238 // `node` is not present in reachable sub-graph.
239 std::vector<NodeIndex> sub_node_indices_[2];
240 // `sub_is_source_[dir][sub_dir]` has size
241 // `sub_reverse_graph_[dir].num_nodes()` such that
242 // `sub_is_source_[dir][sub_dir][sub_node]` is true if `sub_node` is a node in
243 // the reachable sub-graph for direction `dir` (i.e. `sub_reverse_graph[dir]`)
244 // which is a source (resp. destination) is `sub_dir` is FORWARD (resp.
245 // BACKWARD).
246 std::vector<bool> sub_is_source_[2][2];
247 // `sub_min_arc_resources_[dir]` has size `max_resources->size()` and
248 // `sub_min_arc_resources_[dir][r]`, `sub_reverse_graph_[dir].num_nodes()`
249 // such that `sub_min_arc_resources_[dir][r][sub_node]` is the minimum of
250 // resource r needed to get to a destination (resp. come from a source) if
251 // `dir` is FORWARD (resp. BACKWARD).
252 std::vector<std::vector<double>> sub_min_arc_resources_[2];
253 // Maximum number of labels created for each sub-graph.
254 int max_num_created_labels_[2];
255
256 // Data about the last call of the RunConstrainedShortestPathOnDag()
257 // function. A path is only added to the following vectors if and only if
258 // it is feasible with respect to all resources.
259 // A Label includes the cumulative length, resources and the previous arc used
260 // in the path to get to this node.
261 // Instead of having a single vector of `Label` objects (cl/590819865), we
262 // split them into 3 vectors of more fundamental types as this improves
263 // push_back operations and memory release.
264 std::vector<double> lengths_from_sources_[2];
265 std::vector<std::vector<double>> resources_from_sources_[2];
266 std::vector<ArcIndex> incoming_arc_indices_from_sources_[2];
267 std::vector<int> incoming_label_indices_from_sources_[2];
268 std::vector<int> node_first_label_[2];
269 std::vector<int> node_num_labels_[2];
270};
271
272namespace internal {
273template <typename T>
274std::vector<T> GetInversePermutation(const std::vector<T>& permutation);
275} // namespace internal
276
277// -----------------------------------------------------------------------------
278// Implementation.
279// -----------------------------------------------------------------------------
280
281template <class GraphType>
284 const GraphType* graph, const std::vector<double>* arc_lengths,
285 const std::vector<std::vector<double>>* arc_resources,
286 absl::Span<const NodeIndex> topological_order,
287 absl::Span<const NodeIndex> sources,
288 absl::Span<const NodeIndex> destinations,
289 const std::vector<double>* max_resources, int max_num_created_labels)
290 : graph_(graph),
291 arc_lengths_(arc_lengths),
292 arc_resources_(arc_resources),
293 max_resources_(max_resources),
294 sources_(sources),
295 destinations_(destinations),
296 num_resources_(max_resources->size()) {
297 CHECK(graph_ != nullptr);
298 CHECK(arc_lengths_ != nullptr);
299 CHECK(arc_resources_ != nullptr);
300 CHECK(!sources_.empty());
301 CHECK(!destinations_.empty());
302 CHECK(max_resources_ != nullptr);
303 CHECK(!max_resources_->empty())
304 << "max_resources cannot be empty. Use "
305 "ortools/graph/dag_shortest_path.h instead";
306 if (DEBUG_MODE) {
307 CHECK_EQ(arc_lengths->size(), graph->num_arcs());
308 CHECK_EQ(arc_resources->size(), max_resources->size());
309 for (absl::Span<const double> arcs_resource : *arc_resources) {
310 CHECK_EQ(arcs_resource.size(), graph->num_arcs());
311 for (const double arc_resource : arcs_resource) {
312 CHECK(arc_resource >= 0 &&
313 arc_resource != std::numeric_limits<double>::infinity() &&
314 !std::isnan(arc_resource))
315 << absl::StrFormat("resource cannot be negative nor +inf nor NaN");
316 }
317 }
318 for (const double arc_length : *arc_lengths) {
319 CHECK(arc_length != -std::numeric_limits<double>::infinity() &&
320 !std::isnan(arc_length))
321 << absl::StrFormat("length cannot be -inf nor NaN");
322 }
323 CHECK_OK(TopologicalOrderIsValid(*graph, topological_order))
324 << "Invalid topological order";
325 for (const double max_resource : *max_resources) {
326 CHECK(max_resource >= 0 &&
327 max_resource != std::numeric_limits<double>::infinity() &&
328 !std::isnan(max_resource))
329 << absl::StrFormat(
330 "max_resource cannot be negative not +inf nor NaN");
331 }
332 }
333 std::vector<bool> is_source(graph->num_nodes(), false);
334 for (const NodeIndex source : sources) {
335 is_source[source] = true;
336 }
337 for (const NodeIndex destination : destinations) {
338 if (is_source[destination]) {
339 source_is_destination_ = destination;
340 return;
341 }
342 }
343
344 // Full graphs.
345 const GraphType* full_graph[2];
346 const std::vector<std::vector<double>>* full_arc_resources[2];
347 absl::Span<const NodeIndex> full_topological_order[2];
348 absl::Span<const NodeIndex> full_sources[2];
349 // Forward.
350 const int num_nodes = graph->num_nodes();
351 const int num_arcs = graph->num_arcs();
352 full_graph[FORWARD] = graph;
353 full_arc_resources[FORWARD] = arc_resources;
354 full_topological_order[FORWARD] = topological_order;
355 full_sources[FORWARD] = sources;
356 // Backward.
357 GraphType full_backward_graph(num_nodes, num_arcs);
358 for (ArcIndex arc_index = 0; arc_index < num_arcs; ++arc_index) {
359 full_backward_graph.AddArc(graph->Head(arc_index), graph->Tail(arc_index));
360 }
361 std::vector<ArcIndex> full_permutation;
362 full_backward_graph.Build(&full_permutation);
363 const std::vector<ArcIndex> full_inverse_arc_indices =
364 internal::GetInversePermutation(full_permutation);
365 std::vector<std::vector<double>> backward_arc_resources(num_resources_);
366 for (int r = 0; r < num_resources_; ++r) {
367 backward_arc_resources[r] = (*arc_resources)[r];
368 util::Permute(full_permutation, &backward_arc_resources[r]);
369 }
370 std::vector<NodeIndex> full_backward_topological_order;
371 full_backward_topological_order.reserve(num_nodes);
372 for (int i = num_nodes - 1; i >= 0; --i) {
373 full_backward_topological_order.push_back(topological_order[i]);
374 }
375 full_graph[BACKWARD] = &full_backward_graph;
376 full_arc_resources[BACKWARD] = &backward_arc_resources;
377 full_topological_order[BACKWARD] = full_backward_topological_order;
378 full_sources[BACKWARD] = destinations;
379
380 // Get the minimum resources sources -> node and node -> destination for each
381 // node.
382 std::vector<std::vector<double>> full_min_arc_resources[2];
383 for (const Direction dir : {FORWARD, BACKWARD}) {
384 full_min_arc_resources[dir].reserve(num_resources_);
385 std::vector<double> full_arc_resource = full_arc_resources[dir]->front();
386 ShortestPathsOnDagWrapper<GraphType> shortest_paths_on_dag(
387 full_graph[dir], &full_arc_resource, full_topological_order[dir]);
388 for (int r = 0; r < num_resources_; ++r) {
389 full_arc_resource = (*(full_arc_resources[dir]))[r];
390 shortest_paths_on_dag.RunShortestPathOnDag(full_sources[dir]);
391 full_min_arc_resources[dir].push_back(shortest_paths_on_dag.LengthTo());
392 }
393 }
394
395 // Get reachable subgraph.
396 std::vector<bool> is_reachable(num_nodes, true);
397 std::vector<NodeIndex> sub_topological_order;
398 sub_topological_order.reserve(num_nodes);
399 for (const NodeIndex node_index : topological_order) {
400 for (int r = 0; r < num_resources_; ++r) {
401 if (full_min_arc_resources[FORWARD][r][node_index] +
402 full_min_arc_resources[BACKWARD][r][node_index] >
403 (*max_resources)[r]) {
404 is_reachable[node_index] = false;
405 break;
406 }
407 }
408 if (is_reachable[node_index]) {
409 sub_topological_order.push_back(node_index);
410 }
411 }
412 const int reachable_node_count = sub_topological_order.size();
413
414 // We split the number of labels evenly between each search (+1 for the
415 // additional source node).
416 max_num_created_labels_[BACKWARD] = max_num_created_labels / 2 + 1;
417 max_num_created_labels_[FORWARD] =
418 max_num_created_labels - max_num_created_labels / 2 + 1;
419
420 // Split sub-graphs and related information.
421 // The split is based on the number of paths. This is used as a simple proxy
422 // for the number of labels.
423 int mid_index = 0;
424 {
425 // We use double to avoid overflow. Note that this is an heuristic, so we
426 // don't care too much if we are not precise enough.
427 std::vector<double> path_count[2];
428 for (const Direction dir : {FORWARD, BACKWARD}) {
429 const GraphType& reverse_full_graph = *(full_graph[Reverse(dir)]);
430 path_count[dir].resize(num_nodes);
431 for (const NodeIndex source : full_sources[dir]) {
432 ++path_count[dir][source];
433 }
434 for (const NodeIndex to : full_topological_order[dir]) {
435 if (!is_reachable[to]) continue;
436 for (const ArcIndex arc : reverse_full_graph.OutgoingArcs(to)) {
437 const NodeIndex from = reverse_full_graph.Head(arc);
438 if (!is_reachable[from]) continue;
439 path_count[dir][to] += path_count[dir][from];
440 }
441 }
442 }
443 for (const NodeIndex node_index : sub_topological_order) {
444 if (path_count[FORWARD][node_index] > path_count[BACKWARD][node_index]) {
445 break;
446 }
447 ++mid_index;
448 }
449 if (mid_index == reachable_node_count) {
450 mid_index = reachable_node_count / 2;
451 }
452 }
453
454 for (const Direction dir : {FORWARD, BACKWARD}) {
455 absl::Span<const NodeIndex> const sub_nodes =
456 dir == FORWARD
457 ? absl::MakeSpan(sub_topological_order).subspan(0, mid_index)
458 : absl::MakeSpan(sub_topological_order)
459 .subspan(mid_index, reachable_node_count - mid_index);
460 sub_node_indices_[dir].assign(num_nodes, -1);
461 sub_min_arc_resources_[dir].resize(num_resources_);
462 for (int r = 0; r < num_resources_; ++r) {
463 sub_min_arc_resources_[dir][r].resize(sub_nodes.size());
464 }
465 for (NodeIndex i = 0; i < sub_nodes.size(); ++i) {
466 const NodeIndex sub_node_index =
467 dir == FORWARD ? i : sub_nodes.size() - 1 - i;
468 sub_node_indices_[dir][sub_nodes[i]] = sub_node_index;
469 for (int r = 0; r < num_resources_; ++r) {
470 sub_min_arc_resources_[dir][r][sub_node_index] =
471 full_min_arc_resources[Reverse(dir)][r][sub_nodes[i]];
472 }
473 }
474 // IMPORTANT: The sub-graph has an additional node linked to sources (resp.
475 // destinations) for the forward (resp. backward) direction. This additional
476 // node is indexed with the last index. All added arcs are given to have an
477 // arc index in the original graph of -1.
478 const int sub_arcs_count = num_arcs + full_sources[dir].size();
479 sub_reverse_graph_[dir] = GraphType(sub_nodes.size() + 1, sub_arcs_count);
480 sub_arc_resources_[dir].resize(num_resources_);
481 for (int r = 0; r < num_resources_; ++r) {
482 sub_arc_resources_[dir][r].reserve(sub_arcs_count);
483 }
484 sub_full_arc_indices_[dir].reserve(sub_arcs_count);
485 const GraphType& reverse_full_graph = *(full_graph[Reverse(dir)]);
486 for (ArcIndex arc_index = 0; arc_index < num_arcs; ++arc_index) {
487 const NodeIndex from =
488 sub_node_indices_[dir][reverse_full_graph.Tail(arc_index)];
489 const NodeIndex to =
490 sub_node_indices_[dir][reverse_full_graph.Head(arc_index)];
491 if (from == -1 || to == -1) {
492 continue;
493 }
494 sub_reverse_graph_[dir].AddArc(from, to);
495 ArcIndex sub_full_arc_index;
496 if (dir == FORWARD && !full_inverse_arc_indices.empty()) {
497 sub_full_arc_index = full_inverse_arc_indices[arc_index];
498 } else {
499 sub_full_arc_index = arc_index;
500 }
501 for (int r = 0; r < num_resources_; ++r) {
502 sub_arc_resources_[dir][r].push_back(
503 (*arc_resources_)[r][sub_full_arc_index]);
504 }
505 sub_full_arc_indices_[dir].push_back(sub_full_arc_index);
506 }
507 for (const NodeIndex source : full_sources[dir]) {
508 const NodeIndex sub_source = sub_node_indices_[dir][source];
509 if (sub_source == -1) {
510 continue;
511 }
512 sub_reverse_graph_[dir].AddArc(sub_source, sub_nodes.size());
513 for (int r = 0; r < num_resources_; ++r) {
514 sub_arc_resources_[dir][r].push_back(0.0);
515 }
516 sub_full_arc_indices_[dir].push_back(-1);
517 }
518 std::vector<ArcIndex> sub_permutation;
519 sub_reverse_graph_[dir].Build(&sub_permutation);
520 for (int r = 0; r < num_resources_; ++r) {
521 util::Permute(sub_permutation, &sub_arc_resources_[dir][r]);
522 }
523 util::Permute(sub_permutation, &sub_full_arc_indices_[dir]);
524 }
525
526 // Memory allocation is done here and only once in order to avoid
527 // reallocation at each call of `RunConstrainedShortestPathOnDag()` for
528 // better performance.
529 for (const Direction dir : {FORWARD, BACKWARD}) {
530 resources_from_sources_[dir].resize(num_resources_);
531 node_first_label_[dir].resize(sub_reverse_graph_[dir].size());
532 node_num_labels_[dir].resize(sub_reverse_graph_[dir].size());
533 }
534}
535
536template <class GraphType>
538 GraphType>::RunConstrainedShortestPathOnDag() {
539 if (source_is_destination_.has_value()) {
540 return {
541 .length = 0, .arc_path = {}, .node_path = {*source_is_destination_}};
542 }
543 // Assign lengths on sub-relevant graphs.
544 std::vector<double> sub_arc_lengths[2];
545 for (const Direction dir : {FORWARD, BACKWARD}) {
546 sub_arc_lengths[dir].reserve(sub_reverse_graph_[dir].num_arcs());
547 for (ArcIndex sub_arc_index = 0;
548 sub_arc_index < sub_reverse_graph_[dir].num_arcs(); ++sub_arc_index) {
549 const ArcIndex arc_index = sub_full_arc_indices_[dir][sub_arc_index];
550 if (arc_index == -1) {
551 sub_arc_lengths[dir].push_back(0.0);
552 continue;
553 }
554 sub_arc_lengths[dir].push_back((*arc_lengths_)[arc_index]);
555 }
556 }
557
558 {
559 ThreadPool search_threads(2);
560 search_threads.StartWorkers();
561 for (const Direction dir : {FORWARD, BACKWARD}) {
562 search_threads.Schedule([this, dir, &sub_arc_lengths]() {
563 RunHalfConstrainedShortestPathOnDag(
564 /*reverse_graph=*/sub_reverse_graph_[dir],
565 /*arc_lengths=*/sub_arc_lengths[dir],
566 /*arc_resources=*/sub_arc_resources_[dir],
567 /*min_arc_resources=*/sub_min_arc_resources_[dir],
568 /*max_resources=*/*max_resources_,
569 /*max_num_created_labels=*/max_num_created_labels_[dir],
570 /*lengths_from_sources=*/lengths_from_sources_[dir],
571 /*resources_from_sources=*/resources_from_sources_[dir],
572 /*incoming_arc_indices_from_sources=*/
573 incoming_arc_indices_from_sources_[dir],
574 /*incoming_label_indices_from_sources=*/
575 incoming_label_indices_from_sources_[dir],
576 /*first_label=*/node_first_label_[dir],
577 /*num_labels=*/node_num_labels_[dir]);
578 });
579 }
580 }
581
582 // Check destinations within relevant half sub-graphs.
583 LabelPair best_label_pair = {
584 .length = std::numeric_limits<double>::infinity(),
585 .label_index = {-1, -1}};
586 for (const Direction dir : {FORWARD, BACKWARD}) {
587 absl::Span<const NodeIndex> destinations =
588 dir == FORWARD ? destinations_ : sources_;
589 for (const NodeIndex dst : destinations) {
590 const NodeIndex sub_dst = sub_node_indices_[dir][dst];
591 if (sub_dst == -1) {
592 continue;
593 }
594 const int num_labels_dst = node_num_labels_[dir][sub_dst];
595 if (num_labels_dst == 0) {
596 continue;
597 }
598 const int first_label_dst = node_first_label_[dir][sub_dst];
599 for (int label_index = first_label_dst;
600 label_index < first_label_dst + num_labels_dst; ++label_index) {
601 const double length_dst = lengths_from_sources_[dir][label_index];
602 if (length_dst < best_label_pair.length) {
603 best_label_pair.length = length_dst;
604 best_label_pair.label_index[dir] = label_index;
605 }
606 }
607 }
608 }
609
610 const ArcIndex merging_arc_index = MergeHalfRuns(
611 /*graph=*/*graph_, /*arc_lengths=*/*arc_lengths_,
612 /*arc_resources=*/*arc_resources_,
613 /*max_resources=*/*max_resources_,
614 /*sub_node_indices=*/sub_node_indices_,
615 /*lengths_from_sources=*/lengths_from_sources_,
616 /*resources_from_sources=*/resources_from_sources_,
617 /*first_label=*/node_first_label_,
618 /*num_labels=*/node_num_labels_, /*best_label_pair=*/best_label_pair);
619
620 std::vector<ArcIndex> arc_path;
621 for (const Direction dir : {FORWARD, BACKWARD}) {
622 for (const ArcIndex sub_arc_index : ArcPathTo(
623 /*best_label_index=*/best_label_pair.label_index[dir],
624 /*incoming_arc_indices_from_sources=*/
625 incoming_arc_indices_from_sources_[dir],
626 /*incoming_label_indices_from_sources=*/
627 incoming_label_indices_from_sources_[dir])) {
628 const ArcIndex arc_index = sub_full_arc_indices_[dir][sub_arc_index];
629 if (arc_index == -1) {
630 break;
631 }
632 arc_path.push_back(arc_index);
633 }
634 if (dir == FORWARD && merging_arc_index != -1) {
635 absl::c_reverse(arc_path);
636 arc_path.push_back(merging_arc_index);
637 }
638 }
639
640 // Clear all labels from the next run.
641 for (const Direction dir : {FORWARD, BACKWARD}) {
642 lengths_from_sources_[dir].clear();
643 for (int r = 0; r < num_resources_; ++r) {
644 resources_from_sources_[dir][r].clear();
645 }
646 incoming_arc_indices_from_sources_[dir].clear();
647 incoming_label_indices_from_sources_[dir].clear();
648 }
649 return {.length = best_label_pair.length,
650 .arc_path = arc_path,
651 .node_path = NodePathImpliedBy(arc_path, *graph_)};
652}
653
654template <class GraphType>
657 const GraphType& reverse_graph, absl::Span<const double> arc_lengths,
658 absl::Span<const std::vector<double>> arc_resources,
659 absl::Span<const std::vector<double>> min_arc_resources,
660 absl::Span<const double> max_resources,
661 const int max_num_created_labels,
662 std::vector<double>& lengths_from_sources,
663 std::vector<std::vector<double>>& resources_from_sources,
664 std::vector<ArcIndex>& incoming_arc_indices_from_sources,
665 std::vector<int>& incoming_label_indices_from_sources,
666 std::vector<int>& first_label, std::vector<int>& num_labels) {
667 // Initialize source node.
668 const NodeIndex source_node = reverse_graph.num_nodes() - 1;
669 first_label[source_node] = 0;
670 num_labels[source_node] = 1;
671 lengths_from_sources.push_back(0);
672 for (int r = 0; r < num_resources_; ++r) {
673 resources_from_sources[r].push_back(0);
674 }
675 incoming_arc_indices_from_sources.push_back(-1);
676 incoming_label_indices_from_sources.push_back(-1);
677
678 std::vector<double> lengths_to;
679 std::vector<std::vector<double>> resources_to(num_resources_);
680 std::vector<ArcIndex> incoming_arc_indices_to;
681 std::vector<int> incoming_label_indices_to;
682 std::vector<int> label_indices_to;
683 std::vector<double> resources(num_resources_);
684 for (NodeIndex to = 0; to < source_node; ++to) {
685 lengths_to.clear();
686 for (int r = 0; r < num_resources_; ++r) {
687 resources_to[r].clear();
688 }
689 incoming_arc_indices_to.clear();
690 incoming_label_indices_to.clear();
691 for (const ArcIndex reverse_arc_index : reverse_graph.OutgoingArcs(to)) {
692 const NodeIndex from = reverse_graph.Head(reverse_arc_index);
693 const double arc_length = arc_lengths[reverse_arc_index];
694 DCHECK(arc_length != -std::numeric_limits<double>::infinity());
695 if (arc_length == std::numeric_limits<double>::infinity()) {
696 continue;
697 }
698 for (int label_index = first_label[from];
699 label_index < first_label[from] + num_labels[from]; ++label_index) {
700 bool path_is_feasible = true;
701 for (int r = 0; r < num_resources_; ++r) {
702 DCHECK_GE(arc_resources[r][reverse_arc_index], 0.0);
703 resources[r] = resources_from_sources[r][label_index] +
704 arc_resources[r][reverse_arc_index];
705 if (resources[r] + min_arc_resources[r][to] > max_resources[r]) {
706 path_is_feasible = false;
707 break;
708 }
709 }
710 if (!path_is_feasible) {
711 continue;
712 }
713 lengths_to.push_back(lengths_from_sources[label_index] + arc_length);
714 for (int r = 0; r < num_resources_; ++r) {
715 resources_to[r].push_back(resources[r]);
716 }
717 incoming_arc_indices_to.push_back(reverse_arc_index);
718 incoming_label_indices_to.push_back(label_index);
719 }
720 }
721 // Sort labels lexicographically with lengths then resources.
722 label_indices_to.clear();
723 label_indices_to.reserve(lengths_to.size());
724 for (int i = 0; i < lengths_to.size(); ++i) {
725 label_indices_to.push_back(i);
726 }
727 absl::c_sort(label_indices_to, [&](const int i, const int j) {
728 if (lengths_to[i] < lengths_to[j]) return true;
729 if (lengths_to[i] > lengths_to[j]) return false;
730 for (int r = 0; r < num_resources_; ++r) {
731 if (resources_to[r][i] < resources_to[r][j]) return true;
732 if (resources_to[r][i] > resources_to[r][j]) return false;
733 }
734 return i < j;
735 });
736
737 first_label[to] = lengths_from_sources.size();
738 int& num_labels_to = num_labels[to];
739 // Reset the number of labels to zero otherwise it holds the previous run
740 // result.
741 num_labels_to = 0;
742 for (int i = 0; i < label_indices_to.size(); ++i) {
743 // Check if label "i" on node `to` is dominated by any other label.
744 const int label_i_index = label_indices_to[i];
745 bool label_i_is_dominated = false;
746 for (int j = 0; j < i - 1; ++j) {
747 const int label_j_index = label_indices_to[j];
748 if (lengths_to[label_i_index] <= lengths_to[label_j_index]) continue;
749 bool label_j_dominates_label_i = true;
750 for (int r = 0; r < num_resources_; ++r) {
751 if (resources_to[r][label_i_index] <=
752 resources_to[r][label_j_index]) {
753 label_j_dominates_label_i = false;
754 break;
755 }
756 }
757 if (label_j_dominates_label_i) {
758 label_i_is_dominated = true;
759 break;
760 }
761 }
762 if (label_i_is_dominated) continue;
763 lengths_from_sources.push_back(lengths_to[label_i_index]);
764 for (int r = 0; r < num_resources_; ++r) {
765 resources_from_sources[r].push_back(resources_to[r][label_i_index]);
766 }
767 incoming_arc_indices_from_sources.push_back(
768 incoming_arc_indices_to[label_i_index]);
769 incoming_label_indices_from_sources.push_back(
770 incoming_label_indices_to[label_i_index]);
771 ++num_labels_to;
772 if (lengths_from_sources.size() >= max_num_created_labels) {
773 return;
774 }
775 }
776 }
777}
778
779template <class GraphType>
780typename GraphType::ArcIndex
781ConstrainedShortestPathsOnDagWrapper<GraphType>::MergeHalfRuns(
782 const GraphType& graph, absl::Span<const double> arc_lengths,
783 absl::Span<const std::vector<double>> arc_resources,
784 absl::Span<const double> max_resources,
785 const std::vector<NodeIndex> sub_node_indices[2],
786 const std::vector<double> lengths_from_sources[2],
787 const std::vector<std::vector<double>> resources_from_sources[2],
788 const std::vector<int> first_label[2], const std::vector<int> num_labels[2],
789 LabelPair& best_label_pair) {
790 const std::vector<NodeIndex>& forward_sub_node_indices =
791 sub_node_indices[FORWARD];
792 absl::Span<const double> forward_lengths = lengths_from_sources[FORWARD];
793 const std::vector<std::vector<double>>& forward_resources =
794 resources_from_sources[FORWARD];
795 absl::Span<const int> forward_first_label = first_label[FORWARD];
796 absl::Span<const int> forward_num_labels = num_labels[FORWARD];
797 const std::vector<NodeIndex>& backward_sub_node_indices =
798 sub_node_indices[BACKWARD];
799 absl::Span<const double> backward_lengths = lengths_from_sources[BACKWARD];
800 const std::vector<std::vector<double>>& backward_resources =
801 resources_from_sources[BACKWARD];
802 absl::Span<const int> backward_first_label = first_label[BACKWARD];
803 absl::Span<const int> backward_num_labels = num_labels[BACKWARD];
804 ArcIndex merging_arc_index = -1;
805 for (ArcIndex arc_index = 0; arc_index < graph.num_arcs(); ++arc_index) {
806 const NodeIndex sub_from = forward_sub_node_indices[graph.Tail(arc_index)];
807 if (sub_from == -1) {
808 continue;
809 }
810 const NodeIndex sub_to = backward_sub_node_indices[graph.Head(arc_index)];
811 if (sub_to == -1) {
812 continue;
813 }
814 const int num_labels_from = forward_num_labels[sub_from];
815 if (num_labels_from == 0) {
816 continue;
817 }
818 const int num_labels_to = backward_num_labels[sub_to];
819 if (num_labels_to == 0) {
820 continue;
821 }
822 const double arc_length = arc_lengths[arc_index];
823 DCHECK(arc_length != -std::numeric_limits<double>::infinity());
824 if (arc_length == std::numeric_limits<double>::infinity()) {
825 continue;
826 }
827 const int first_label_from = forward_first_label[sub_from];
828 const int first_label_to = backward_first_label[sub_to];
829 for (int label_to_index = first_label_to;
830 label_to_index < first_label_to + num_labels_to; ++label_to_index) {
831 const double length_to = backward_lengths[label_to_index];
832 for (int label_from_index = first_label_from;
833 label_from_index < first_label_from + num_labels_from;
834 ++label_from_index) {
835 const double length_from = forward_lengths[label_from_index];
836 if (length_from + arc_length + length_to >= best_label_pair.length) {
837 continue;
838 }
839 bool path_is_feasible = true;
840 for (int r = 0; r < num_resources_; ++r) {
841 DCHECK_GE(arc_resources[r][arc_index], 0.0);
842 if (forward_resources[r][label_from_index] +
843 arc_resources[r][arc_index] +
844 backward_resources[r][label_to_index] >
845 max_resources[r]) {
846 path_is_feasible = false;
847 break;
848 }
849 }
850 if (!path_is_feasible) {
851 continue;
852 }
853 best_label_pair.length = length_from + arc_length + length_to;
854 best_label_pair.label_index[FORWARD] = label_from_index;
855 best_label_pair.label_index[BACKWARD] = label_to_index;
856 merging_arc_index = arc_index;
857 }
858 }
859 }
860 return merging_arc_index;
861}
862
863template <class GraphType>
864std::vector<typename GraphType::ArcIndex>
865ConstrainedShortestPathsOnDagWrapper<GraphType>::ArcPathTo(
866 const int best_label_index,
867 absl::Span<const ArcIndex> incoming_arc_indices_from_sources,
868 absl::Span<const int> incoming_label_indices_from_sources) const {
869 int current_label_index = best_label_index;
870 std::vector<ArcIndex> arc_path;
871 for (int i = 0; i < graph_->num_nodes(); ++i) {
872 if (current_label_index == -1) {
873 break;
874 }
875 arc_path.push_back(incoming_arc_indices_from_sources[current_label_index]);
876 current_label_index =
877 incoming_label_indices_from_sources[current_label_index];
878 }
879 return arc_path;
880}
881
882template <typename GraphType>
883std::vector<typename GraphType::NodeIndex>
884ConstrainedShortestPathsOnDagWrapper<GraphType>::NodePathImpliedBy(
885 absl::Span<const ArcIndex> arc_path, const GraphType& graph) const {
886 if (arc_path.empty()) {
887 return {};
888 }
889 std::vector<NodeIndex> node_path;
890 node_path.reserve(arc_path.size() + 1);
891 for (const ArcIndex arc_index : arc_path) {
892 node_path.push_back(graph.Tail(arc_index));
893 }
894 node_path.push_back(graph.Head(arc_path.back()));
895 return node_path;
896}
897
898namespace internal {
899
900template <typename T>
901std::vector<T> GetInversePermutation(const std::vector<T>& permutation) {
902 std::vector<T> inverse_permutation(permutation.size());
903 if (!permutation.empty()) {
904 for (T i = 0; i < permutation.size(); ++i) {
905 inverse_permutation[permutation[i]] = i;
906 }
908 return inverse_permutation;
909}
910
911} // namespace internal
912
913} // namespace operations_research
914
915#endif // OR_TOOLS_GRAPH_DAG_CONSTRAINED_SHORTEST_PATH_H_
ConstrainedShortestPathsOnDagWrapper(const GraphType *graph, const std::vector< double > *arc_lengths, const std::vector< std::vector< double > > *arc_resources, absl::Span< const NodeIndex > topological_order, absl::Span< const NodeIndex > sources, absl::Span< const NodeIndex > destinations, const std::vector< double > *max_resources, int max_num_created_labels=1e9)
double LengthTo(NodeIndex node) const
Returns the length of the shortest path from node's source to node.
void RunShortestPathOnDag(absl::Span< const NodeIndex > sources)
End of the interface. Below is the implementation.
std::vector< T > GetInversePermutation(const std::vector< T > &permutation)
In SWIG mode, we don't want anything besides these top-level includes.
absl::Status TopologicalOrderIsValid(const GraphType &graph, absl::Span< const typename GraphType::NodeIndex > topological_order)
BlossomGraph::NodeIndex NodeIndex
const bool DEBUG_MODE
Definition radix_sort.h:266
std::vector< typename GraphType::NodeIndex > NodePathImpliedBy(absl::Span< const typename GraphType::ArcIndex > arc_path, const GraphType &graph)
PathWithLength ConstrainedShortestPathsOnDag(const int num_nodes, absl::Span< const ArcWithLengthAndResources > arcs_with_length_and_resources, int source, int destination, const std::vector< double > &max_resources)
void Permute(const IntVector &permutation, Array *array_to_permute)
Definition graph.h:1105
trees with all degrees equal to
std::vector< typename GraphType::ArcIndex > arc_path
std::vector< typename GraphType::NodeIndex > node_path