Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
routing_cuts.cc
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
15
16#include <algorithm>
17#include <atomic>
18#include <cmath>
19#include <cstdint>
20#include <cstdlib>
21#include <functional>
22#include <limits>
23#include <memory>
24#include <numeric>
25#include <queue>
26#include <stack>
27#include <string>
28#include <tuple>
29#include <utility>
30#include <vector>
31
32#include "absl/algorithm/container.h"
33#include "absl/cleanup/cleanup.h"
34#include "absl/container/flat_hash_map.h"
35#include "absl/container/flat_hash_set.h"
36#include "absl/flags/flag.h"
37#include "absl/log/check.h"
38#include "absl/log/log.h"
39#include "absl/log/vlog_is_on.h"
40#include "absl/numeric/bits.h"
41#include "absl/numeric/int128.h"
42#include "absl/random/distributions.h"
43#include "absl/strings/str_cat.h"
44#include "absl/types/span.h"
50#include "ortools/graph/graph.h"
52#include "ortools/sat/clause.h"
53#include "ortools/sat/cp_model.pb.h"
55#include "ortools/sat/cuts.h"
56#include "ortools/sat/integer.h"
60#include "ortools/sat/model.h"
62#include "ortools/sat/routes_support_graph.pb.h"
64#include "ortools/sat/sat_parameters.pb.h"
66#include "ortools/sat/util.h"
68
69ABSL_FLAG(bool, cp_model_dump_routes_support_graphs, false,
70 "DEBUG ONLY. When set to true, SolveCpModel() dumps the arcs with "
71 "non-zero LP values of the routes constraints, at decision level 0, "
72 "which are used to subsequently generate cuts. The values are "
73 "written as a SupportGraphProto in text format to "
74 "'FLAGS_cp_model_dump_prefix'support_graph_{counter}.pb.txt.");
75
76namespace operations_research {
77namespace sat {
78
79namespace {
80
81// If we have many sets of (var >= values relations) and we know at least one
82// set is true, then we can derive valid global lower bounds on the variable
83// that appear in all such sets.
84//
85// This represents "one step" of computing such bounds by adding the sets one at
86// the time. When we process a set 'new_lbs':
87// - For the first set, then 'new_lbs' is globally valid.
88// - For a new set, we need to erase variables not appearing in new_lbs and
89// take the min for the ones appearing in both.
90//
91// TODO(user): The operation is symmetric, so we could swap both hash_map if
92// new_lbs is smaller than current_lbs as a small optimization.
93void ComputeMinLowerBoundOfSharedVariables(
94 const absl::flat_hash_map<IntegerVariable, IntegerValue>& new_lbs,
95 absl::flat_hash_map<IntegerVariable, IntegerValue>* current_lbs) {
96 if (new_lbs.empty()) current_lbs->clear();
97 for (auto it = current_lbs->begin(); it != current_lbs->end();) {
98 const IntegerVariable var = it->first;
99 auto other_it = new_lbs.find(var);
100 if (other_it == new_lbs.end()) {
101 // We have no bounds about var in the new_fact, we need to erase it.
102 //
103 // Note: it = current_lbs->erase(it) does not work for flat_hash_map, this
104 // is the recommended way.
105 current_lbs->erase(it++);
106 } else {
107 // Update.
108 it->second = std::min(it->second, other_it->second);
109 ++it;
110 }
111 }
112}
113
114} // namespace
115
117 int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
118 const std::vector<Literal>& literals, Model* model)
119 : tails_(tails),
120 heads_(heads),
121 literals_(literals),
122 binary_relation_repository_(
123 *model->GetOrCreate<BinaryRelationRepository>()),
124 trail_(*model->GetOrCreate<Trail>()),
125 integer_trail_(*model->GetOrCreate<IntegerTrail>()),
126 integer_encoder_(*model->GetOrCreate<IntegerEncoder>()),
127 shared_stats_(model->GetOrCreate<SharedStatistics>()),
128 in_subset_(num_nodes, false),
129 index_in_subset_(num_nodes, -1),
130 incoming_arc_indices_(num_nodes),
131 outgoing_arc_indices_(num_nodes),
132 reachable_(num_nodes, false),
133 next_reachable_(num_nodes, false),
134 node_var_lower_bounds_(num_nodes),
135 next_node_var_lower_bounds_(num_nodes),
136 bin_packing_helper_(
137 model->GetOrCreate<SatParameters>()->routing_cut_dp_effort()) {}
138
140 if (!VLOG_IS_ON(1)) return;
141 std::vector<std::pair<std::string, int64_t>> stats;
142 stats.push_back({"RoutingDp/num_full_dp_calls", num_full_dp_calls_});
143 stats.push_back({"RoutingDp/num_full_dp_skips", num_full_dp_skips_});
144 stats.push_back(
145 {"RoutingDp/num_full_dp_early_abort", num_full_dp_early_abort_});
146 stats.push_back(
147 {"RoutingDp/num_full_dp_work_abort", num_full_dp_work_abort_});
148 stats.push_back({"RoutingDp/num_full_dp_rc_skip", num_full_dp_rc_skip_});
149 stats.push_back(
150 {"RoutingDp/num_full_dp_unique_arc", num_full_dp_unique_arc_});
151 for (const auto& [name, count] : num_by_type_) {
152 stats.push_back({absl::StrCat("RoutingDp/num_bounds_", name), count});
153 }
154 shared_stats_->AddStats(stats);
155}
156
157// TODO(user): This is costly and we might do it more than once per subset.
158// fix.
159void MinOutgoingFlowHelper::PrecomputeDataForSubset(
160 absl::Span<const int> subset) {
161 cannot_be_first_.clear();
162 cannot_be_last_.clear();
163
164 const int num_nodes = in_subset_.size();
165 in_subset_.assign(num_nodes, false);
166 index_in_subset_.assign(num_nodes, -1);
167 for (int i = 0; i < subset.size(); ++i) {
168 const int n = subset[i];
169 in_subset_[n] = true;
170 index_in_subset_[n] = i;
171 }
172
173 has_incoming_arcs_from_outside_.assign(num_nodes, false);
174 has_outgoing_arcs_to_outside_.assign(num_nodes, false);
175 incoming_arcs_from_outside_.assign(num_nodes, {});
176 outgoing_arcs_to_outside_.assign(num_nodes, {});
177
178 for (auto& v : incoming_arc_indices_) v.clear();
179 for (auto& v : outgoing_arc_indices_) v.clear();
180 for (int i = 0; i < tails_.size(); ++i) {
181 const int tail = tails_[i];
182 const int head = heads_[i];
183
184 // we always ignore self-arcs here.
185 if (tail == head) continue;
186 const bool tail_in = in_subset_[tail];
187 const bool head_in = in_subset_[head];
188
189 if (tail_in && head_in) {
190 outgoing_arc_indices_[tail].push_back(i);
191 incoming_arc_indices_[head].push_back(i);
192 } else if (tail_in && !head_in) {
193 outgoing_arcs_to_outside_[tail].Add(i);
194 has_outgoing_arcs_to_outside_[tail] = true;
195 } else if (!tail_in && head_in) {
196 incoming_arcs_from_outside_[head].Add(i);
197 has_incoming_arcs_from_outside_[head] = true;
198 }
199 }
200}
201
203 absl::Span<const int> subset, const RouteRelationsHelper& helper,
204 BestBoundHelper* best_bound) {
205 PrecomputeDataForSubset(subset);
206 DCHECK_EQ(helper.num_nodes(), in_subset_.size());
207 DCHECK_EQ(helper.num_arcs(), tails_.size());
208
209 int min_outgoing_flow = 1;
210 std::string best_name;
211 for (int d = 0; d < helper.num_dimensions(); ++d) {
212 for (const bool negate_expressions : {false, true}) {
213 for (const bool use_incoming : {true, false}) {
214 std::string info;
215 const absl::Span<SpecialBinPackingHelper::ItemOrBin> objects =
216 RelaxIntoSpecialBinPackingProblem(subset, d, negate_expressions,
217 use_incoming, helper);
218 const int bound = bin_packing_helper_.ComputeMinNumberOfBins(
219 objects, objects_that_cannot_be_bin_and_reach_minimum_, info);
220 if (bound >= min_outgoing_flow) {
221 min_outgoing_flow = bound;
222 best_name = absl::StrCat((use_incoming ? "in" : "out"), info,
223 (negate_expressions ? "_neg" : ""), "_", d);
224 cannot_be_first_.clear();
225 cannot_be_last_.clear();
226 auto& vec = use_incoming ? cannot_be_first_ : cannot_be_last_;
227 for (const int i : objects_that_cannot_be_bin_and_reach_minimum_) {
228 vec.push_back(objects[i].node);
229 }
230 best_bound->Update(bound, "AutomaticDimension", cannot_be_first_,
231 cannot_be_last_);
232 }
233 }
234 }
235 }
236
237 if (min_outgoing_flow > 1) num_by_type_[best_name]++;
238 return min_outgoing_flow;
239}
240
241absl::Span<SpecialBinPackingHelper::ItemOrBin>
242MinOutgoingFlowHelper::RelaxIntoSpecialBinPackingProblem(
243 absl::Span<const int> subset, int dimension, bool negate_expressions,
244 bool use_incoming, const RouteRelationsHelper& helper) {
245 tmp_bin_packing_problem_.clear();
246
247 // Computes UB/LB. The derivation in the .h used a simple version, but we can
248 // be a bit tighter. We explain for the use_incoming case only:
249 //
250 // In this setting, recall that bins correspond to the first node of each
251 // route covering the subset. The max_upper_bound is a maximum over the last
252 // node of each route, and as such it does not need to consider nodes with no
253 // arc leaving the subset: such nodes cannot be last.
254 //
255 // Moreover, if the bin is not empty, then the bound for that bin does not
256 // need to consider the variable of the "bin" node itself. If the new bound is
257 // non-negative, then it is also valid for the empty bin. If it is negative we
258 // have an issue as using it will force the bin to not be empty. But we can
259 // always use std::max(0, new_bound) which is always valid and still tighter
260 // than the bound explained in the .h which is always non-negative.
261 //
262 // This way we can have a different bound for each bin, which is slightly
263 // stronger than using the same bound for all of them.
264 const int num_nodes = in_subset_.size();
265 std::vector<IntegerValue> min_lower_bound_of_others(num_nodes,
267 std::vector<IntegerValue> max_upper_bound_of_others(num_nodes,
269
270 // Lower bound on the node expression when it is first to serve the subset.
271 const auto lb_when_first = [&helper, this](int n, int dimension,
272 bool negate_expressions) {
273 AffineExpression expr = helper.GetNodeExpression(n, dimension);
274 if (negate_expressions) expr = expr.Negated();
275 IntegerValue lb = integer_trail_.LevelZeroLowerBound(expr);
276 if (incoming_arcs_from_outside_[n].IsUnique()) {
277 const int a = incoming_arcs_from_outside_[n].Get();
278 AffineExpression tail_expr =
279 helper.GetNodeExpression(tails_[a], dimension);
280 if (negate_expressions) tail_expr = tail_expr.Negated();
281 const IntegerValue offset =
282 helper.GetArcOffsetLowerBound(a, dimension, negate_expressions);
283 lb = std::max(lb, integer_trail_.LevelZeroLowerBound(tail_expr) + offset);
284 }
285 return lb;
286 };
287
288 // Upper bound on the node expression when it is last to serve the subset.
289 const auto ub_when_last = [&helper, this](int n, int dimension,
290 bool negate_expressions) {
291 AffineExpression expr = helper.GetNodeExpression(n, dimension);
292 if (negate_expressions) expr = expr.Negated();
293 IntegerValue ub = integer_trail_.LevelZeroUpperBound(expr);
294 if (outgoing_arcs_to_outside_[n].IsUnique()) {
295 const int a = outgoing_arcs_to_outside_[n].Get();
296 AffineExpression head_expr =
297 helper.GetNodeExpression(heads_[a], dimension);
298 if (negate_expressions) head_expr = head_expr.Negated();
299 const IntegerValue offset =
300 helper.GetArcOffsetLowerBound(a, dimension, negate_expressions);
301 ub = std::min(ub, integer_trail_.LevelZeroUpperBound(head_expr) - offset);
302 }
303 return ub;
304 };
305
306 // We do a forward and a backward pass in order to compute the min/max
307 // of all the other nodes for each node.
308 for (const bool forward_pass : {true, false}) {
309 IntegerValue local_min = kMaxIntegerValue;
310 IntegerValue local_max = kMinIntegerValue;
311 const int size = subset.size();
312 for (int i = 0; i < size; ++i) {
313 const int n = forward_pass ? subset[i] : subset[size - 1 - i];
314
315 // The local_min/max contains the min/max of all nodes strictly before 'n'
316 // in the forward pass (resp. striclty after). So after two passes,
317 // min/max_..._of_others[n] will contains the min/max of all nodes
318 // different from n.
319 min_lower_bound_of_others[n] =
320 std::min(min_lower_bound_of_others[n], local_min);
321 max_upper_bound_of_others[n] =
322 std::max(max_upper_bound_of_others[n], local_max);
323
324 // Update local_min/max.
325 if (has_incoming_arcs_from_outside_[n]) {
326 local_min = std::min(local_min,
327 lb_when_first(n, dimension, negate_expressions));
328 }
329 if (has_outgoing_arcs_to_outside_[n]) {
330 local_max =
331 std::max(local_max, ub_when_last(n, dimension, negate_expressions));
332 }
333 }
334 }
335
336 for (const int n : subset) {
337 const absl::Span<const int> arcs =
338 use_incoming ? incoming_arc_indices_[n] : outgoing_arc_indices_[n];
339 IntegerValue demand = kMaxIntegerValue;
340 for (const int a : arcs) {
341 demand = std::min(demand, helper.GetArcOffsetLowerBound(
342 a, dimension, negate_expressions));
343 }
344
345 SpecialBinPackingHelper::ItemOrBin obj;
346 obj.node = n;
347 obj.demand = demand;
348
349 // We compute the capacity like this to avoid overflow and always have
350 // a non-negative capacity (see above).
351 obj.capacity = 0;
352 if (use_incoming) {
353 const IntegerValue lb = lb_when_first(n, dimension, negate_expressions);
354 if (max_upper_bound_of_others[n] > lb) {
355 obj.capacity = max_upper_bound_of_others[n] - lb;
356 }
357 } else {
358 const IntegerValue ub = ub_when_last(n, dimension, negate_expressions);
359 if (ub > min_lower_bound_of_others[n]) {
360 obj.capacity = ub - min_lower_bound_of_others[n];
361 }
362 }
363
364 // Note that we don't explicitly deal with the corner case of a subset node
365 // with no arcs. This corresponds to INFEASIBLE problem and should be dealt
366 // with elsewhere. Being "less restrictive" will still result in a valid
367 // bound and that is enough here.
368 if ((use_incoming && !has_incoming_arcs_from_outside_[n]) ||
369 (!use_incoming && !has_outgoing_arcs_to_outside_[n])) {
371 obj.capacity = 0;
372 } else if (arcs.empty()) {
374 obj.demand = 0; // We don't want kMaxIntegerValue !
375 } else {
377 }
378
379 tmp_bin_packing_problem_.push_back(obj);
380 }
381
382 return absl::MakeSpan(tmp_bin_packing_problem_);
383}
384
386 absl::Span<ItemOrBin> objects,
387 std::vector<int>& objects_that_cannot_be_bin_and_reach_minimum,
388 std::string& info) {
389 objects_that_cannot_be_bin_and_reach_minimum.clear();
390 if (objects.empty()) return 0;
391
392 IntegerValue sum_of_demands(0);
393 int64_t gcd = 0;
394 int64_t max_capacity = 0;
395 int max_num_items = 0;
396 bool all_demands_are_non_negative = true;
397 for (const ItemOrBin& obj : objects) {
398 sum_of_demands = CapAddI(sum_of_demands, obj.demand);
399 if (obj.type != MUST_BE_BIN) {
400 ++max_num_items;
401 if (obj.demand < 0) {
402 all_demands_are_non_negative = false;
403 }
404 gcd = std::gcd(gcd, std::abs(obj.demand.value()));
405 }
406 if (obj.type != MUST_BE_ITEM) {
407 max_capacity = std::max(max_capacity, obj.capacity.value());
408 }
409 }
410
411 // TODO(user): we can probably handle a couple of extra case rather than just
412 // bailing out here and below.
413 if (AtMinOrMaxInt64I(sum_of_demands)) return 0;
414
415 // If the gcd of all the demands term is positive, we can divide everything.
416 if (gcd > 1) {
417 absl::StrAppend(&info, "_gcd");
418 for (ItemOrBin& obj : objects) {
419 obj.demand /= gcd;
420 obj.capacity = FloorRatio(obj.capacity, gcd);
421 }
422 max_capacity = MathUtil::FloorOfRatio(max_capacity, gcd);
423 }
424
425 const int result = ComputeMinNumberOfBinsInternal(
426 objects, objects_that_cannot_be_bin_and_reach_minimum);
427
428 // We only try DP if it can help.
429 //
430 // Tricky: For the GreedyPackingWorks() to make sense, we use the fact that
431 // ComputeMinNumberOfBinsInternal() moves the best bins first. In any case
432 // the code will still be correct if it is not the case as we just use this
433 // as an heuristic to do less work.
434 if (result > 1 && all_demands_are_non_negative &&
435 max_bounded_subset_sum_exact_.ComplexityEstimate(
436 max_num_items, max_capacity) < dp_effort_ &&
437 !GreedyPackingWorks(result, objects) &&
438 UseDpToTightenCapacities(objects)) {
439 const int new_result = ComputeMinNumberOfBinsInternal(
440 objects, objects_that_cannot_be_bin_and_reach_minimum);
441 CHECK_GE(new_result, result);
442 if (new_result > result) {
443 absl::StrAppend(&info, "_dp");
444 }
445 return new_result;
446 }
447
448 return result;
449}
450
451// TODO(user): We could be more tight here and also compute the max_reachable
452// under smaller capacity than the max_capacity.
454 absl::Span<ItemOrBin> objects) {
455 tmp_demands_.clear();
456 int64_t max_capacity = 0;
457 for (const ItemOrBin& obj : objects) {
458 if (obj.type != MUST_BE_BIN) {
459 tmp_demands_.push_back(obj.demand.value());
460 }
461 if (obj.type != MUST_BE_ITEM) {
462 max_capacity = std::max(max_capacity, obj.capacity.value());
463 }
464 }
465 const int64_t max_reachable =
466 max_bounded_subset_sum_exact_.MaxSubsetSum(tmp_demands_, max_capacity);
467 bool some_change = false;
468 if (max_reachable < max_capacity) {
469 for (ItemOrBin& obj : objects) {
470 if (obj.capacity > max_reachable) {
471 some_change = true;
472 obj.capacity = max_reachable;
473 }
474 }
475 }
476 return some_change;
477}
478
479int SpecialBinPackingHelper::ComputeMinNumberOfBinsInternal(
480 absl::Span<ItemOrBin> objects,
481 std::vector<int>& objects_that_cannot_be_bin_and_reach_minimum) {
482 objects_that_cannot_be_bin_and_reach_minimum.clear();
483
484 // For a given choice of bins (set B), a feasible problem must satisfy.
485 // sum_{i \notin B} demands_i <= sum_{i \in B} capacity_i.
486 //
487 // Using this we can compute a lower bound on the number of bins needed
488 // using a greedy algorithm that chooses B in order to make the above
489 // inequality as "loose" as possible.
490 //
491 // This puts 'a' before 'b' if we get more unused capacity by using 'a' as a
492 // bin and 'b' as an item rather than the other way around. If we call the
493 // other items demands D and capacities C, the two options are:
494 // - option1: a.demand + D <= b.capacity + C
495 // - option2: b.demand + D <= a.capacity + C
496 //
497 // The option2 above leads to more unused capacity:
498 // (a.capacity - b.demand) > (b.capacity - a.demand).
499 std::stable_sort(objects.begin(), objects.end(),
500 [](const ItemOrBin& a, const ItemOrBin& b) {
501 if (a.type != b.type) {
502 // We want in order:
503 // MUST_BE_BIN, ITEM_OR_BIN, MUST_BE_ITEM.
504 return a.type > b.type;
505 }
506 return a.capacity + a.demand > b.capacity + b.demand;
507 });
508
509 // Note that we already checked for overflow.
510 IntegerValue sum_of_demands(0);
511 for (const ItemOrBin& obj : objects) {
512 sum_of_demands += obj.demand;
513 }
514
515 // We start with no bins (sum_of_demands=everything, sum_of_capacity=0) and
516 // add the best bins one by one until we have sum_of_demands <=
517 // sum_of_capacity.
518 int num_bins = 0;
519 IntegerValue sum_of_capacity(0);
520 for (; num_bins < objects.size(); ++num_bins) {
521 // Use obj as a bin instead of as an item, if possible (i.e., unless
522 // obj.type is MUST_BE_ITEM).
523 const ItemOrBin& obj = objects[num_bins];
524 if (obj.type != MUST_BE_BIN && sum_of_demands <= sum_of_capacity) {
525 break;
526 }
527 if (obj.type == MUST_BE_ITEM) {
528 // Because of our order, we only have objects of type MUST_BE_ITEM left.
529 // Hence we can no longer change items to bins, and since the demands
530 // exceed the capacity, the problem is infeasible. We just return the
531 // (number of objects + 1) in this case (any bound is valid).
532 DCHECK_GT(sum_of_demands, sum_of_capacity);
533 return objects.size() + 1;
534 }
535 sum_of_capacity = CapAddI(sum_of_capacity, obj.capacity);
536 if (AtMinOrMaxInt64I(sum_of_capacity)) return num_bins;
537 sum_of_demands -= obj.demand;
538 }
539
540 // We can identify items that cannot be a bin in a solution reaching the
541 // minimum 'num_bins'.
542 //
543 // TODO(user): Somewhat related are the items that can be removed while still
544 // having the same required 'num_bins' to solve the problem.
545 if (num_bins > 0 && num_bins != objects.size()) {
546 const int worst_used_bin = num_bins - 1;
547 if (objects[worst_used_bin].type == MUST_BE_BIN) {
548 // All first object must be bins, we don't have a choice if we want to
549 // reach the lower bound. The others cannot be bins.
550 for (int i = num_bins; i < objects.size(); ++i) {
551 objects_that_cannot_be_bin_and_reach_minimum.push_back(i);
552 }
553 } else {
554 // We revert the worst_used_bin and tries to use the other instead.
555 CHECK_EQ(objects[worst_used_bin].type, ITEM_OR_BIN);
556 sum_of_demands += objects[worst_used_bin].demand;
557 sum_of_capacity -= objects[worst_used_bin].capacity;
558 const IntegerValue slack = sum_of_demands - sum_of_capacity;
559 CHECK_GE(slack, 0) << num_bins << " " << objects.size() << " "
560 << objects[worst_used_bin].type;
561
562 for (int i = num_bins; i < objects.size(); ++i) {
563 if (objects[i].type == MUST_BE_ITEM) continue;
564
565 // If we use this as a bin instead of "worst_used_bin" can we still
566 // fulfill the demands?
567 const IntegerValue new_demand = sum_of_demands - objects[i].demand;
568 const IntegerValue new_capacity = sum_of_capacity + objects[i].capacity;
569 if (new_demand > new_capacity) {
570 objects_that_cannot_be_bin_and_reach_minimum.push_back(i);
571 }
572 }
573 }
574 }
575
576 return num_bins;
577}
578
580 int num_bins, absl::Span<const ItemOrBin> objects) {
581 if (num_bins >= objects.size()) return true;
582 CHECK_GE(num_bins, 1);
583 tmp_capacities_.resize(num_bins);
584 for (int i = 0; i < num_bins; ++i) {
585 tmp_capacities_[i] = objects[i].capacity;
586 }
587
588 // We place the objects in order in the first bin that fits.
589 for (const ItemOrBin object : objects.subspan(num_bins)) {
590 bool placed = false;
591 for (int b = 0; b < num_bins; ++b) {
592 if (object.demand <= tmp_capacities_[b]) {
593 placed = true;
594 tmp_capacities_[b] -= object.demand;
595 break;
596 }
597 }
598 if (!placed) return false;
599 }
600 return true;
601}
602
604 absl::Span<const int> subset) {
605 DCHECK_GE(subset.size(), 1);
606 DCHECK(absl::c_all_of(next_reachable_, [](bool b) { return !b; }));
607 DCHECK(absl::c_all_of(node_var_lower_bounds_,
608 [](const auto& m) { return m.empty(); }));
609 DCHECK(absl::c_all_of(next_node_var_lower_bounds_,
610 [](const auto& m) { return m.empty(); }));
611 PrecomputeDataForSubset(subset);
612
613 // Conservatively assume that each subset node is reachable from outside.
614 // TODO(user): use has_incoming_arcs_from_outside_[] to be more precise.
615 reachable_ = in_subset_;
616
617 // Maximum number of nodes of a feasible path inside the subset.
618 int longest_path_length = 1;
619 absl::flat_hash_map<IntegerVariable, IntegerValue> tmp_lbs;
620 while (longest_path_length < subset.size()) {
621 bool at_least_one_next_node_reachable = false;
622 for (const int head : subset) {
623 bool is_reachable = false;
624 for (const int incoming_arc_index : incoming_arc_indices_[head]) {
625 const int tail = tails_[incoming_arc_index];
626 const Literal lit = literals_[incoming_arc_index];
627 if (!reachable_[tail]) continue;
628
629 // If this arc cannot be taken skip.
630 tmp_lbs.clear();
631 if (!binary_relation_repository_.PropagateLocalBounds(
632 integer_trail_, lit, node_var_lower_bounds_[tail], &tmp_lbs)) {
633 continue;
634 }
635
636 if (!is_reachable) {
637 // This is the first arc that reach this node.
638 is_reachable = true;
639 next_node_var_lower_bounds_[head] = tmp_lbs;
640 } else {
641 // Combine information with previous one.
642 ComputeMinLowerBoundOfSharedVariables(
643 tmp_lbs, &next_node_var_lower_bounds_[head]);
644 }
645 }
646
647 next_reachable_[head] = is_reachable;
648 if (next_reachable_[head]) at_least_one_next_node_reachable = true;
649 }
650 if (!at_least_one_next_node_reachable) {
651 break;
652 }
653 std::swap(reachable_, next_reachable_);
654 std::swap(node_var_lower_bounds_, next_node_var_lower_bounds_);
655 for (const int n : subset) {
656 next_node_var_lower_bounds_[n].clear();
657 }
658 ++longest_path_length;
659 }
660
661 // The maximum number of distinct paths of length `longest_path_length`.
662 int max_longest_paths = 0;
663 for (const int n : subset) {
664 if (reachable_[n]) ++max_longest_paths;
665
666 // Reset the temporary data structures for the next call.
667 next_reachable_[n] = false;
668 node_var_lower_bounds_[n].clear();
669 next_node_var_lower_bounds_[n].clear();
670 }
671 return GetMinOutgoingFlow(subset.size(), longest_path_length,
672 max_longest_paths);
673}
674
675int MinOutgoingFlowHelper::GetMinOutgoingFlow(int subset_size,
676 int longest_path_length,
677 int max_longest_paths) {
678 if (max_longest_paths * longest_path_length < subset_size) {
679 // If `longest_path_length` is 1, there should be one such path per node,
680 // and thus `max_longest_paths` should be the subset size. This branch can
681 // thus not be taken in this case.
682 DCHECK_GT(longest_path_length, 1);
683 // TODO(user): Consider using information about the number of shorter
684 // paths to derive an even better bound.
685 return max_longest_paths +
687 subset_size - max_longest_paths * longest_path_length,
688 longest_path_length - 1);
689 }
690 return MathUtil::CeilOfRatio(subset_size, longest_path_length);
691}
692
693namespace {
694struct Path {
695 uint32_t node_set; // Bit i is set iif node subset[i] is in the path.
696 int last_node; // The last node in the path.
697
698 bool operator==(const Path& p) const {
699 return node_set == p.node_set && last_node == p.last_node;
700 }
701
702 template <typename H>
703 friend H AbslHashValue(H h, const Path& p) {
704 return H::combine(std::move(h), p.node_set, p.last_node);
705 }
706};
707
708} // namespace
709
711 absl::Span<const int> subset) {
712 DCHECK_GE(subset.size(), 1);
713 DCHECK_LE(subset.size(), 32);
714 PrecomputeDataForSubset(subset);
715
716 std::vector<int> longest_path_length_by_end_node(subset.size(), 1);
717
718 absl::flat_hash_map<IntegerVariable, IntegerValue> tmp_lbs;
719 absl::flat_hash_map<Path, absl::flat_hash_map<IntegerVariable, IntegerValue>>
720 path_var_bounds;
721 std::vector<Path> paths;
722 std::vector<Path> next_paths;
723 for (int i = 0; i < subset.size(); ++i) {
724 // TODO(user): use has_incoming_arcs_from_outside_[] to skip some nodes.
725 paths.push_back(
726 {.node_set = static_cast<uint32_t>(1 << i), .last_node = subset[i]});
727 path_var_bounds[paths.back()] = {}; // LevelZero bounds.
728 }
729 int longest_path_length = 1;
730 for (int path_length = 1; path_length <= subset.size(); ++path_length) {
731 for (const Path& path : paths) {
732 // We remove it from the hash_map since this entry should no longer be
733 // used as we create path of increasing length.
734 DCHECK(path_var_bounds.contains(path));
735 const absl::flat_hash_map<IntegerVariable, IntegerValue> path_bounds =
736 std::move(path_var_bounds.extract(path).mapped());
737
738 // We have found a feasible path, update the path length statistics...
739 longest_path_length = path_length;
740 longest_path_length_by_end_node[index_in_subset_[path.last_node]] =
741 path_length;
742
743 // ... and try to extend the path with each arc going out of its last
744 // node.
745 for (const int outgoing_arc_index :
746 outgoing_arc_indices_[path.last_node]) {
747 const int head = heads_[outgoing_arc_index];
748 const int head_index_in_subset = index_in_subset_[head];
749 DCHECK_NE(head_index_in_subset, -1);
750 if (path.node_set & (1 << head_index_in_subset)) continue;
751 const Path new_path = {
752 .node_set = path.node_set | (1 << head_index_in_subset),
753 .last_node = head};
754
755 // If this arc cannot be taken skip.
756 tmp_lbs.clear();
757 if (!binary_relation_repository_.PropagateLocalBounds(
758 integer_trail_, literals_[outgoing_arc_index], path_bounds,
759 &tmp_lbs)) {
760 continue;
761 }
762
763 const auto [it, inserted] = path_var_bounds.insert({new_path, tmp_lbs});
764 if (inserted) {
765 // We have a feasible path to a new state.
766 next_paths.push_back(new_path);
767 } else {
768 // We found another way to reach this state, only keep common best
769 // bounds.
770 ComputeMinLowerBoundOfSharedVariables(tmp_lbs, &it->second);
771 }
772 }
773 }
774 std::swap(paths, next_paths);
775 next_paths.clear();
776 }
777
778 int max_longest_paths = 0;
779 for (int i = 0; i < subset.size(); ++i) {
780 if (longest_path_length_by_end_node[i] == longest_path_length) {
781 ++max_longest_paths;
782 }
783 }
784
785 return GetMinOutgoingFlow(subset.size(), longest_path_length,
786 max_longest_paths);
787}
788
790 int k, absl::Span<const int> subset, RouteRelationsHelper* helper,
791 LinearConstraintManager* manager, int special_node,
792 bool use_forward_direction) {
793 cannot_be_first_.clear();
794 cannot_be_last_.clear();
795 if (k >= subset.size()) return true;
796 if (subset.size() > 31) return true;
797
798 // If we have a "special node" from which we must start, make sure it is
799 // always first. This simplifies the code a bit, and we don't care about the
800 // order of nodes in subset.
801 if (special_node >= 0) {
802 tmp_subset_.assign(subset.begin(), subset.end());
803 bool seen = false;
804 for (int i = 0; i < tmp_subset_.size(); ++i) {
805 if (tmp_subset_[i] == special_node) {
806 seen = true;
807 std::swap(tmp_subset_[0], tmp_subset_[i]);
808 break;
809 }
810 }
811 CHECK(seen);
812
813 // Make the span point to the new order.
814 subset = tmp_subset_;
815 }
816
817 PrecomputeDataForSubset(subset);
818 ++num_full_dp_calls_;
819
820 // We need a special graph here to handle both options.
821 // TODO(user): Maybe we should abstract that in a better way.
823 adjacency.reserve(subset.size());
824 for (int i = 0; i < subset.size(); ++i) {
825 adjacency.Add({});
826 const int node = subset[i];
827 if (helper != nullptr) {
828 CHECK(use_forward_direction);
829 for (int j = 0; j < subset.size(); ++j) {
830 if (i == j) continue;
831 if (helper->PathExists(node, subset[j])) {
832 adjacency.AppendToLastVector({subset[j], Literal(kNoLiteralIndex)});
833 }
834 }
835 } else {
836 if (use_forward_direction) {
837 for (const int arc : outgoing_arc_indices_[node]) {
838 adjacency.AppendToLastVector({heads_[arc], literals_[arc]});
839 }
840 } else {
841 for (const int arc : incoming_arc_indices_[node]) {
842 adjacency.AppendToLastVector({tails_[arc], literals_[arc]});
843 }
844 }
845 }
846 }
847
848 struct State {
849 // Bit i is set iif node subset[i] is in one of the current routes.
850 uint32_t node_set;
851
852 // The last nodes of each of the k routes. If the hamming weight is less
853 // that k, then at least one route is still empty.
854 uint32_t last_nodes_set;
855
856 // Valid lower bounds for this state.
857 //
858 // Note that unlike the other algorithm here, we keep the collective bounds
859 // of all the nodes so far, so this is likely in
860 // O(longest_route * num_dimensions) which can take quite a lot of space.
861 //
862 // By "dimensions", we mean the number of variables appearing in binary
863 // relation controlled by an arc literal. See for instance
864 // RouteRelationsHelper that also uses a similar definition.
865 //
866 // Hopefully the DFS order limit the number of entry to O(n^2 * k), so still
867 // somewhat reasonable for small values.
868 absl::flat_hash_map<IntegerVariable, IntegerValue> lbs;
869
870 // The sum of the reduced costs of all the literals selected to form the
871 // route(s) encoded by this state. If this becomes too large we know that
872 // this state can never be used to get to a better solution than the current
873 // best one and is thus "infeasible" for the problem of finding a better
874 // solution.
875 absl::int128 sum_of_reduced_costs = 0;
876 };
877
878 const int size = subset.size();
879 const uint32_t final_mask = (1 << size) - 1;
880
881 const auto& can_be_last_set = use_forward_direction
882 ? has_outgoing_arcs_to_outside_
883 : has_incoming_arcs_from_outside_;
884 const auto& can_be_first_set = use_forward_direction
885 ? has_incoming_arcs_from_outside_
886 : has_outgoing_arcs_to_outside_;
887
888 uint32_t can_be_last_mask = 0;
889 for (int i = 0; i < subset.size(); ++i) {
890 if (can_be_last_set[subset[i]]) {
891 can_be_last_mask |= (1 << i);
892 }
893 }
894
895 // This is also correlated to the work done, and we abort if we starts to
896 // do too much work on one instance.
897 int64_t allocated_memory_estimate = 0;
898
899 // We just do a DFS from the initial state.
900 std::vector<State> states;
901 states.push_back(State());
902
903 // If a state has an unique arc to enter/leave it, we can update it further.
904 const auto update_using_unique_arc = [manager, this](UniqueArc unique_arc,
905 State& state) {
906 if (!unique_arc.IsUnique()) return true;
907 if (manager == nullptr) return true;
908
909 ++num_full_dp_unique_arc_;
910 const Literal unique_lit = literals_[unique_arc.Get()];
911 state.sum_of_reduced_costs += manager->GetLiteralReducedCost(unique_lit);
912 if (state.sum_of_reduced_costs > manager->ReducedCostsGap()) {
913 ++num_full_dp_rc_skip_;
914 return false;
915 }
916
917 absl::flat_hash_map<IntegerVariable, IntegerValue> copy = state.lbs;
918 return binary_relation_repository_.PropagateLocalBounds(
919 integer_trail_, unique_lit, copy, &state.lbs);
920 };
921
922 // We always start with the first node in this case.
923 if (special_node >= 0) {
924 states.back().node_set = 1;
925 states.back().last_nodes_set = 1;
926
927 // If there is a unique way to enter (resp. leave) that subset, we can start
928 // with tighter bounds!
929 const UniqueArc unique_arc = use_forward_direction
930 ? incoming_arcs_from_outside_[subset[0]]
931 : outgoing_arcs_to_outside_[subset[0]];
932 if (!update_using_unique_arc(unique_arc, states.back())) {
933 return false; // We can't even start!
934 }
935 }
936
937 while (!states.empty()) {
938 if (allocated_memory_estimate > 1e7) {
939 ++num_full_dp_work_abort_;
940 return true; // Abort.
941 }
942 const State from_state = std::move(states.back());
943 states.pop_back();
944
945 // The number of routes is the hamming weight of from_state.last_nodes_set.
946 const int num_routes = absl::popcount(from_state.last_nodes_set);
947
948 // We start by choosing the first k starts (in increasing order).
949 // For that we only add after the maximum position already chosen.
950 if (num_routes < k) {
951 const int num_extra = k - num_routes - 1;
952 for (int i = 0; i + num_extra < size; ++i) {
953 if (from_state.node_set >> i) continue;
954 if (!can_be_first_set[subset[i]]) continue;
955
956 // All "initial-state" start with empty hash-map that correspond to
957 // the level zero bounds.
958 State to_state;
959 const uint32_t head_mask = (1 << i);
960 to_state.node_set = from_state.node_set | head_mask;
961 to_state.last_nodes_set = from_state.last_nodes_set | head_mask;
962 to_state.lbs = from_state.lbs;
963 to_state.sum_of_reduced_costs = from_state.sum_of_reduced_costs;
964
965 const UniqueArc unique_arc =
966 use_forward_direction ? incoming_arcs_from_outside_[subset[i]]
967 : outgoing_arcs_to_outside_[subset[i]];
968 if (!update_using_unique_arc(unique_arc, to_state)) {
969 continue;
970 }
971
972 if (to_state.node_set == final_mask) {
973 ++num_full_dp_early_abort_;
974 return true; // All served!
975 }
976 states.push_back(std::move(to_state));
977 }
978 continue;
979 }
980
981 // We have k routes, extend one of the last nodes.
982 for (int i = 0; i < size; ++i) {
983 const uint32_t tail_mask = 1 << i;
984 if ((from_state.last_nodes_set & tail_mask) == 0) continue;
985 const int tail = subset[i];
986 for (const auto [head, literal] : adjacency[i]) {
987 const uint32_t head_mask = (1 << index_in_subset_[head]);
988 if (from_state.node_set & head_mask) continue;
989
990 State to_state;
991
992 // Use reduced costs to exclude solutions that cannot improve our
993 // current best solution.
994 if (manager != nullptr) {
995 DCHECK_EQ(helper, nullptr);
996 to_state.sum_of_reduced_costs =
997 from_state.sum_of_reduced_costs +
998 manager->GetLiteralReducedCost(literal);
999 if (to_state.sum_of_reduced_costs > manager->ReducedCostsGap()) {
1000 ++num_full_dp_rc_skip_;
1001 continue;
1002 }
1003 }
1004
1005 // We start from the old bounds and update them.
1006 to_state.lbs = from_state.lbs;
1007 if (helper != nullptr) {
1009 integer_trail_, tail, head, from_state.lbs, &to_state.lbs)) {
1010 continue;
1011 }
1012 } else {
1013 if (!binary_relation_repository_.PropagateLocalBounds(
1014 integer_trail_, literal, from_state.lbs, &to_state.lbs)) {
1015 continue;
1016 }
1017 }
1018
1019 to_state.node_set = from_state.node_set | head_mask;
1020 to_state.last_nodes_set = from_state.last_nodes_set | head_mask;
1021 to_state.last_nodes_set ^= tail_mask;
1022 allocated_memory_estimate += to_state.lbs.size();
1023 if (to_state.node_set == final_mask) {
1024 // One of the last node has no arc to outside, this is not possible.
1025 if (to_state.last_nodes_set & ~can_be_last_mask) continue;
1026
1027 // Add the constraints implied by each last node that has an unique
1028 // way to leave the subset.
1029 bool infeasible = false;
1030 int last_mask = to_state.last_nodes_set;
1031 for (int j = 0; last_mask; j++, last_mask /= 2) {
1032 if ((last_mask & 1) == 0) continue;
1033
1034 const UniqueArc unique_arc =
1035 use_forward_direction ? outgoing_arcs_to_outside_[subset[j]]
1036 : incoming_arcs_from_outside_[subset[j]];
1037 if (!update_using_unique_arc(unique_arc, to_state)) {
1038 infeasible = true;
1039 break;
1040 }
1041 }
1042 if (infeasible) {
1043 continue;
1044 }
1045
1046 ++num_full_dp_early_abort_;
1047 return true;
1048 }
1049 states.push_back(std::move(to_state));
1050 }
1051 }
1052 }
1053
1054 // We explored everything, no way to serve this with only k routes!
1055 return false;
1056}
1057
1058namespace {
1059
1060// Represents the relation tail_coeff.X + head_coeff.Y \in [lhs, rhs] between
1061// the X and Y expressions associated with the tail and head of the given arc,
1062// respectively, and the given dimension (head_coeff is always positive).
1063//
1064// An "empty" struct with all fields set to 0 is used to represent no relation.
1065struct LocalRelation {
1066 IntegerValue tail_coeff = 0;
1067 IntegerValue head_coeff = 0;
1068 IntegerValue lhs;
1069 IntegerValue rhs;
1070
1071 bool empty() const { return tail_coeff == 0 && head_coeff == 0; }
1072 bool operator==(const LocalRelation& r) const {
1073 return tail_coeff == r.tail_coeff && head_coeff == r.head_coeff &&
1074 lhs == r.lhs && rhs == r.rhs;
1075 }
1076};
1077
1078IntegerVariable UniqueSharedVariable(const sat::Relation& r1,
1079 const sat::Relation& r2) {
1080 DCHECK_NE(r1.a.var, r1.b.var);
1081 DCHECK_NE(r2.a.var, r2.b.var);
1082 if (r1.a.var == r2.a.var && r1.b.var != r2.b.var) return r1.a.var;
1083 if (r1.a.var == r2.b.var && r1.b.var != r2.a.var) return r1.a.var;
1084 if (r1.b.var == r2.a.var && r1.a.var != r2.b.var) return r1.b.var;
1085 if (r1.b.var == r2.b.var && r1.a.var != r2.a.var) return r1.b.var;
1086 return kNoIntegerVariable;
1087}
1088
1089class RouteRelationsBuilder {
1090 public:
1091 using HeadMinusTailBounds = RouteRelationsHelper::HeadMinusTailBounds;
1092
1093 RouteRelationsBuilder(
1094 int num_nodes, absl::Span<const int> tails, absl::Span<const int> heads,
1095 absl::Span<const Literal> literals,
1096 absl::Span<const AffineExpression> flat_node_dim_expressions,
1097 const BinaryRelationRepository& binary_relation_repository)
1098 : num_nodes_(num_nodes),
1099 num_arcs_(tails.size()),
1100 tails_(tails),
1101 heads_(heads),
1102 literals_(literals),
1103 binary_relation_repository_(binary_relation_repository) {
1104 if (!flat_node_dim_expressions.empty()) {
1105 DCHECK_EQ(flat_node_dim_expressions.size() % num_nodes, 0);
1106 num_dimensions_ = flat_node_dim_expressions.size() / num_nodes;
1107 flat_node_dim_expressions_.assign(flat_node_dim_expressions.begin(),
1108 flat_node_dim_expressions.end());
1109 }
1110 }
1111
1112 int num_dimensions() const { return num_dimensions_; }
1113
1114 std::vector<AffineExpression>& flat_node_dim_expressions() {
1115 return flat_node_dim_expressions_;
1116 }
1117
1118 std::vector<HeadMinusTailBounds>& flat_arc_dim_relations() {
1119 return flat_arc_dim_relations_;
1120 }
1121
1122 std::vector<IntegerValue>& flat_shortest_path_lbs() {
1123 return flat_shortest_path_lbs_;
1124 }
1125
1126 bool BuildNodeExpressions() {
1127 if (flat_node_dim_expressions_.empty()) {
1128 // Step 1: find the number of dimensions (as the number of connected
1129 // components in the graph of binary relations), and find to which
1130 // dimension each variable belongs.
1131 ComputeDimensionOfEachVariable();
1132 if (num_dimensions_ == 0) return false;
1133
1134 // Step 2: find the variables which can be unambiguously associated with a
1135 // node and dimension.
1136 // - compute the indices of the binary relations which can be
1137 // unambiguously associated with the incoming and outgoing arcs of each
1138 // node, per dimension.
1139 ComputeAdjacentRelationsPerNodeAndDimension();
1140 // - find variable associations by using variables which are uniquely
1141 // shared by two adjacent relations of a node.
1142 std::queue<std::pair<int, int>> node_dim_pairs =
1143 ComputeVarAssociationsFromSharedVariableOfAdjacentRelations();
1144 if (node_dim_pairs.empty()) return false;
1145 // - find more variable associations by using arcs from nodes with
1146 // an associated variable, whose other end has no associated variable, and
1147 // where only one variable can be associated with it.
1148 ComputeVarAssociationsFromRelationsWithSingleFreeVar(node_dim_pairs);
1149 }
1150 return true;
1151 }
1152
1153 void BuildArcRelations(Model* model) {
1154 // Step 3: compute the relation for each arc and dimension, now that the
1155 // variables associated with each node and dimension are known.
1156 ComputeArcRelations(model);
1157 }
1158
1159 // We only consider lower bound for now.
1160 //
1161 // TODO(user): can we be smarter and compute a "scheduling" like lower bound
1162 // that exploit the [lb, ub] of the node expression even more?
1163 void BuildShortestPaths(Model* model) {
1164 if (num_nodes_ >= 100) return;
1165
1166 auto& integer_trail = *model->GetOrCreate<IntegerTrail>();
1167 std::vector<IntegerValue> flat_trivial_lbs(num_nodes_ * num_nodes_);
1168 const auto trivial = [this, &flat_trivial_lbs](int i,
1169 int j) -> IntegerValue& {
1170 CHECK_NE(i, j);
1171 CHECK_NE(i, 0);
1172 CHECK_NE(j, 0);
1173 return flat_trivial_lbs[i * num_nodes_ + j];
1174 };
1175 flat_shortest_path_lbs_.resize(num_nodes_ * num_nodes_ * num_dimensions_,
1176 std::numeric_limits<IntegerValue>::max());
1177
1178 for (int dim = 0; dim < num_dimensions_; ++dim) {
1179 // We can never beat trivial bound relations, starts by filling them.
1180 for (int i = 1; i < num_nodes_; ++i) {
1181 const AffineExpression& i_expr = node_expression(i, dim);
1182 const IntegerValue i_ub = integer_trail.LevelZeroUpperBound(i_expr);
1183 for (int j = 1; j < num_nodes_; ++j) {
1184 if (i == j) continue;
1185 const AffineExpression& j_expr = node_expression(j, dim);
1186 trivial(i, j) = integer_trail.LevelZeroLowerBound(j_expr) - i_ub;
1187 }
1188 }
1189
1190 // Starts by the known arc relations.
1191 const auto path = [this](int dim, int i, int j) -> IntegerValue& {
1192 return flat_shortest_path_lbs_[dim * num_nodes_ * num_nodes_ +
1193 i * num_nodes_ + j];
1194 };
1195 for (int arc = 0; arc < tails_.size(); ++arc) {
1196 const int tail = tails_[arc];
1197 const int head = heads_[arc];
1198 if (tail == 0 || head == 0 || tail == head) continue;
1199 path(dim, tail, head) =
1200 std::max(trivial(tail, head),
1201 flat_arc_dim_relations_[arc * num_dimensions_ + dim].lhs);
1202 }
1203
1204 // We do floyd-warshall to complete them into shortest path.
1205 for (int k = 1; k < num_nodes_; ++k) {
1206 for (int i = 1; i < num_nodes_; ++i) {
1207 if (i == k) continue;
1208 for (int j = 1; j < num_nodes_; ++j) {
1209 if (j == k || j == i) continue;
1210 path(dim, i, j) =
1211 std::min(path(dim, i, j),
1212 std::max(trivial(i, j),
1213 CapAddI(path(dim, i, k), path(dim, k, j))));
1214 }
1215 }
1216 }
1217 }
1218 }
1219
1220 private:
1221 AffineExpression& node_expression(int node, int dimension) {
1222 return flat_node_dim_expressions_[node * num_dimensions_ + dimension];
1223 };
1224
1225 void ProcessNewArcRelation(int arc_index, int dimension, LocalRelation r) {
1226 if (r.empty()) return;
1227 if (r.head_coeff < 0) {
1228 r.tail_coeff = -r.tail_coeff;
1229 r.head_coeff = -r.head_coeff;
1230 r.lhs = -r.lhs;
1231 r.rhs = -r.rhs;
1232 std::swap(r.lhs, r.rhs);
1233 }
1234 if (r.head_coeff != 1 || r.tail_coeff != -1) return;
1235
1236 // Combine the relation information.
1237 HeadMinusTailBounds& relation =
1238 flat_arc_dim_relations_[arc_index * num_dimensions_ + dimension];
1239 relation.lhs = std::max(relation.lhs, r.lhs);
1240 relation.rhs = std::min(relation.rhs, r.rhs);
1241 }
1242
1243 void ComputeDimensionOfEachVariable() {
1244 // Step 1: find the number of dimensions (as the number of connected
1245 // components in the graph of binary relations).
1246 // TODO(user): see if we can use a shared
1247 // DenseConnectedComponentsFinder with one node per variable of the whole
1248 // model instead.
1249 ConnectedComponentsFinder<IntegerVariable> cc_finder;
1250 for (int i = 0; i < num_arcs_; ++i) {
1251 if (tails_[i] == heads_[i]) continue;
1252 num_arcs_per_literal_[literals_[i]]++;
1253 for (const int relation_index :
1254 binary_relation_repository_.IndicesOfRelationsEnforcedBy(
1255 literals_[i])) {
1256 const auto& r = binary_relation_repository_.relation(relation_index);
1257 if (r.a.var == kNoIntegerVariable || r.b.var == kNoIntegerVariable) {
1258 continue;
1259 }
1260 cc_finder.AddEdge(r.a.var, r.b.var);
1261 }
1262 }
1263 const std::vector<std::vector<IntegerVariable>> connected_components =
1264 cc_finder.FindConnectedComponents();
1265 for (int i = 0; i < connected_components.size(); ++i) {
1266 for (const IntegerVariable var : connected_components[i]) {
1267 dimension_by_var_[var] = i;
1268 }
1269 }
1270 num_dimensions_ = connected_components.size();
1271 }
1272
1273 void ComputeAdjacentRelationsPerNodeAndDimension() {
1274 adjacent_relation_indices_ = std::vector<std::vector<std::vector<int>>>(
1275 num_dimensions_, std::vector<std::vector<int>>(num_nodes_));
1276 for (int i = 0; i < num_arcs_; ++i) {
1277 if (tails_[i] == heads_[i]) continue;
1278 // If a literal is associated with more than one arc, a relation
1279 // associated with this literal cannot be unambiguously associated with an
1280 // arc.
1281 if (num_arcs_per_literal_[literals_[i]] > 1) continue;
1282 for (const int relation_index :
1283 binary_relation_repository_.IndicesOfRelationsEnforcedBy(
1284 literals_[i])) {
1285 const auto& r = binary_relation_repository_.relation(relation_index);
1286 if (r.a.var == kNoIntegerVariable || r.b.var == kNoIntegerVariable) {
1287 continue;
1288 }
1289 const int dimension = dimension_by_var_[r.a.var];
1290 adjacent_relation_indices_[dimension][tails_[i]].push_back(
1291 relation_index);
1292 adjacent_relation_indices_[dimension][heads_[i]].push_back(
1293 relation_index);
1294 }
1295 }
1296 }
1297
1298 // Returns the (node, dimension) pairs for which a variable association has
1299 // been found.
1300 std::queue<std::pair<int, int>>
1301 ComputeVarAssociationsFromSharedVariableOfAdjacentRelations() {
1302 flat_node_dim_expressions_ = std::vector<AffineExpression>(
1303 num_nodes_ * num_dimensions_, AffineExpression());
1304 std::queue<std::pair<int, int>> result;
1305 for (int n = 0; n < num_nodes_; ++n) {
1306 for (int d = 0; d < num_dimensions_; ++d) {
1307 // If two relations on incoming or outgoing arcs of n have a unique
1308 // shared variable, such as in the case of l <-X,Y-> n <-Y,Z-> m (i.e. a
1309 // relation between X and Y on the (l,n) arc, and a relation between Y
1310 // and Z on the (n,m) arc), then this variable is necessarily associated
1311 // with n.
1312 for (const int r1_index : adjacent_relation_indices_[d][n]) {
1313 const auto& r1 = binary_relation_repository_.relation(r1_index);
1314 for (const int r2_index : adjacent_relation_indices_[d][n]) {
1315 if (r1_index == r2_index) continue;
1316 const auto& r2 = binary_relation_repository_.relation(r2_index);
1317 const IntegerVariable shared_var = UniqueSharedVariable(r1, r2);
1318 if (shared_var == kNoIntegerVariable) continue;
1319 DCHECK_EQ(dimension_by_var_[shared_var], d);
1320 AffineExpression& node_expr = node_expression(n, d);
1321 if (node_expr.IsConstant()) {
1322 result.push({n, d});
1323 } else if (node_expr.var != shared_var) {
1324 VLOG(2) << "Several vars per node and dimension in route with "
1325 << num_nodes_ << " nodes and " << num_arcs_ << " arcs";
1326 return {};
1327 }
1328 node_expr = shared_var;
1329 }
1330 }
1331 }
1332 }
1333 return result;
1334 }
1335
1336 void ComputeVarAssociationsFromRelationsWithSingleFreeVar(
1337 std::queue<std::pair<int, int>>& node_dim_pairs_to_process) {
1338 std::vector<std::vector<int>> adjacent_arcs_per_node(num_nodes_);
1339 for (int i = 0; i < num_arcs_; ++i) {
1340 if (tails_[i] == heads_[i]) continue;
1341 adjacent_arcs_per_node[tails_[i]].push_back(i);
1342 adjacent_arcs_per_node[heads_[i]].push_back(i);
1343 }
1344 while (!node_dim_pairs_to_process.empty()) {
1345 const auto [node, dimension] = node_dim_pairs_to_process.front();
1346 const AffineExpression node_expr = node_expression(node, dimension);
1347 DCHECK(!node_expr.IsConstant());
1348 node_dim_pairs_to_process.pop();
1349 for (const int arc_index : adjacent_arcs_per_node[node]) {
1350 const int tail = tails_[arc_index];
1351 const int head = heads_[arc_index];
1352 DCHECK(node == tail || node == head);
1353 int other_node = node == tail ? head : tail;
1354 if (!node_expression(other_node, dimension).IsConstant()) {
1355 continue;
1356 }
1357 IntegerVariable candidate_var = kNoIntegerVariable;
1358 bool candidate_var_is_unique = true;
1359 for (const int relation_index :
1360 binary_relation_repository_.IndicesOfRelationsEnforcedBy(
1361 literals_[arc_index])) {
1362 const auto& r = binary_relation_repository_.relation(relation_index);
1363 if (r.a.var == kNoIntegerVariable || r.b.var == kNoIntegerVariable) {
1364 continue;
1365 }
1366 if (r.a.var == node_expr.var) {
1367 if (candidate_var != kNoIntegerVariable &&
1368 candidate_var != r.b.var) {
1369 candidate_var_is_unique = false;
1370 break;
1371 }
1372 candidate_var = r.b.var;
1373 }
1374 if (r.b.var == node_expr.var) {
1375 if (candidate_var != kNoIntegerVariable &&
1376 candidate_var != r.a.var) {
1377 candidate_var_is_unique = false;
1378 break;
1379 }
1380 candidate_var = r.a.var;
1381 }
1382 }
1383 if (candidate_var != kNoIntegerVariable && candidate_var_is_unique) {
1384 node_expression(other_node, dimension) = candidate_var;
1385 node_dim_pairs_to_process.push({other_node, dimension});
1386 }
1387 }
1388 }
1389 }
1390
1391 void ComputeArcRelations(Model* model) {
1392 auto& binary_implication_graph =
1393 *model->GetOrCreate<BinaryImplicationGraph>();
1394 const auto& integer_encoder = *model->GetOrCreate<IntegerEncoder>();
1395 const auto& trail = *model->GetOrCreate<Trail>();
1396 const auto& integer_trail = *model->GetOrCreate<IntegerTrail>();
1397 DCHECK_EQ(trail.CurrentDecisionLevel(), 0);
1398
1399 flat_arc_dim_relations_ = std::vector<HeadMinusTailBounds>(
1400 num_arcs_ * num_dimensions_, HeadMinusTailBounds());
1401 binary_implication_graph.ResetWorkDone();
1402 for (int i = 0; i < num_arcs_; ++i) {
1403 const int tail = tails_[i];
1404 const int head = heads_[i];
1405 if (tail == head) continue;
1406 // The literals directly or indirectly implied by the arc literal,
1407 // including the arc literal itself.
1408 const Literal arc_literal_singleton[1] = {literals_[i]};
1409 absl::Span<const Literal> implied_literals =
1410 binary_implication_graph.WorkDone() < 1e8
1411 ? binary_implication_graph.GetAllImpliedLiterals(literals_[i])
1412 : absl::MakeSpan(arc_literal_singleton);
1413 // The integer view of the implied literals (resp. of their negation), and
1414 // the variable lower bounds implied directly or indirectly by the arc
1415 // literal.
1416 absl::flat_hash_set<IntegerVariable> implied_views;
1417 absl::flat_hash_set<IntegerVariable> negated_implied_views;
1418 absl::flat_hash_map<IntegerVariable, IntegerValue> implied_lower_bounds;
1419 for (const Literal implied : implied_literals) {
1420 implied_views.insert(integer_encoder.GetLiteralView(implied));
1421 negated_implied_views.insert(
1422 integer_encoder.GetLiteralView(implied.Negated()));
1423 for (const auto& [var, lb] :
1424 integer_encoder.GetIntegerLiterals(implied)) {
1425 auto it = implied_lower_bounds.find(var);
1426 if (it == implied_lower_bounds.end()) {
1427 implied_lower_bounds[var] = lb;
1428 } else {
1429 it->second = std::max(it->second, lb);
1430 }
1431 }
1432 }
1433 // Returns the bounds of the given expression, assuming that all the
1434 // literals implied by the arc literal are 1.
1435 auto get_bounds = [&](const NodeExpression& expr) {
1436 if (expr.var != kNoIntegerVariable) {
1437 if (implied_views.contains(expr.var)) {
1438 IntegerValue constant_value = expr.coeff + expr.offset;
1439 return std::make_pair(constant_value, constant_value);
1440 }
1441 if (implied_views.contains(NegationOf(expr.var))) {
1442 IntegerValue constant_value = -expr.coeff + expr.offset;
1443 return std::make_pair(constant_value, constant_value);
1444 }
1445 if (negated_implied_views.contains(expr.var) ||
1446 negated_implied_views.contains(NegationOf(expr.var))) {
1447 return std::make_pair(expr.offset, expr.offset);
1448 }
1449 }
1450 const AffineExpression e(expr.var, expr.coeff, expr.offset);
1451 IntegerValue lb = integer_trail.LevelZeroLowerBound(e);
1452 auto it = implied_lower_bounds.find(e.var);
1453 if (it != implied_lower_bounds.end()) {
1454 lb = std::max(lb, e.ValueAt(it->second));
1455 }
1456 IntegerValue ub = integer_trail.LevelZeroUpperBound(e);
1457 it = implied_lower_bounds.find(NegationOf(e.var));
1458 if (it != implied_lower_bounds.end()) {
1459 ub = std::min(ub, e.ValueAt(-it->second));
1460 }
1461 return std::make_pair(lb, ub);
1462 };
1463 // Changes `expr` to a constant expression if possible, and returns true.
1464 // Otherwise, returns false.
1465 auto to_constant = [&](NodeExpression& expr) {
1466 if (expr.var == kNoIntegerVariable) return true;
1467 if (integer_trail.IsFixed(expr.var)) {
1468 expr.offset += integer_trail.FixedValue(expr.var) * expr.coeff;
1469 expr.var = kNoIntegerVariable;
1470 expr.coeff = 0;
1471 return true;
1472 }
1473 const auto [min, max] = get_bounds(expr);
1474 if (min == max) {
1475 expr = NodeExpression(kNoIntegerVariable, 0, min);
1476 return true;
1477 }
1478 return false;
1479 };
1480 for (int dimension = 0; dimension < num_dimensions_; ++dimension) {
1481 NodeExpression tail_expr(node_expression(tail, dimension));
1482 NodeExpression head_expr(node_expression(head, dimension));
1483 for (const Literal implied_lit : implied_literals) {
1484 for (const int relation_index :
1485 binary_relation_repository_.IndicesOfRelationsEnforcedBy(
1486 implied_lit)) {
1487 auto r = binary_relation_repository_.relation(relation_index);
1488 // Try to match the relation variables with the node expression
1489 // variables. First swap the relation terms if needed (this does not
1490 // change the relation bounds).
1491 if ((r.a.var != kNoIntegerVariable && r.a.var == head_expr.var) ||
1492 (r.b.var != kNoIntegerVariable && r.b.var == tail_expr.var)) {
1493 std::swap(r.a, r.b);
1494 }
1495 // If the relation has only one term, try to remove the variable
1496 // in the node expression corresponding to the missing term.
1497 if (r.a.var == kNoIntegerVariable) {
1498 if (!to_constant(tail_expr)) continue;
1499 } else if (r.b.var == kNoIntegerVariable) {
1500 if (!to_constant(head_expr)) continue;
1501 }
1502 // If the relation and node expression variables do not match, we
1503 // cannot use this relation for this arc.
1504 if (!((tail_expr.var == r.a.var && head_expr.var == r.b.var) ||
1505 (tail_expr.var == r.b.var && head_expr.var == r.a.var))) {
1506 continue;
1507 }
1508 ComputeArcRelation(i, dimension, tail_expr, head_expr, r,
1509 integer_trail);
1510 }
1511 }
1512
1513 // Check if we can use non-enforced relations to improve the relations.
1514 if (!tail_expr.IsEmpty() && !head_expr.IsEmpty()) {
1515 for (const int relation_index :
1516 binary_relation_repository_.IndicesOfRelationsBetween(
1517 tail_expr.var, head_expr.var)) {
1518 ComputeArcRelation(
1519 i, dimension, tail_expr, head_expr,
1520 binary_relation_repository_.relation(relation_index),
1521 integer_trail);
1522 }
1523 }
1524
1525 // Check if we can use the default relation to improve the relations.
1526 // Compute the bounds of X_head - X_tail as [lb(X_head) - ub(X_tail),
1527 // ub(X_head) - lb(X_tail)], which is always true. Note that by our
1528 // overflow precondition, the difference of any two variable bounds
1529 // should fit on an int64_t.
1530 std::pair<IntegerValue, IntegerValue> bounds;
1531 if (head_expr.var == tail_expr.var) {
1532 const NodeExpression offset_expr(head_expr.var,
1533 head_expr.coeff - tail_expr.coeff,
1534 head_expr.offset - tail_expr.offset);
1535 bounds = get_bounds(offset_expr);
1536 } else {
1537 const auto [tail_min, tail_max] = get_bounds(tail_expr);
1538 const auto [head_min, head_max] = get_bounds(head_expr);
1539 bounds = {head_min - tail_max, head_max - tail_min};
1540 }
1541 ProcessNewArcRelation(i, dimension,
1542 {.tail_coeff = -1,
1543 .head_coeff = 1,
1544 .lhs = bounds.first,
1545 .rhs = bounds.second});
1546 }
1547 }
1548 }
1549
1550 // Infers a relation between two given expressions which is equivalent to a
1551 // given relation `r` between the variables used in these expressions.
1552 void ComputeArcRelation(int arc_index, int dimension,
1553 const NodeExpression& tail_expr,
1554 const NodeExpression& head_expr, sat::Relation r,
1555 const IntegerTrail& integer_trail) {
1556 DCHECK((r.a.var == tail_expr.var && r.b.var == head_expr.var) ||
1557 (r.a.var == head_expr.var && r.b.var == tail_expr.var));
1558 if (r.a.var != tail_expr.var) std::swap(r.a, r.b);
1559 if (r.a.coeff == 0 || tail_expr.coeff == 0) {
1560 LocalRelation result = ComputeArcUnaryRelation(head_expr, tail_expr,
1561 r.b.coeff, r.lhs, r.rhs);
1562 std::swap(result.tail_coeff, result.head_coeff);
1563 ProcessNewArcRelation(arc_index, dimension, result);
1564 return;
1565 }
1566 if (r.b.coeff == 0 || head_expr.coeff == 0) {
1567 ProcessNewArcRelation(arc_index, dimension,
1568 ComputeArcUnaryRelation(tail_expr, head_expr,
1569 r.a.coeff, r.lhs, r.rhs));
1570 return;
1571 }
1572 const auto [lhs, rhs] =
1573 GetDifferenceBounds(tail_expr, head_expr, r,
1574 {integer_trail.LowerBound(tail_expr.var),
1575 integer_trail.UpperBound(tail_expr.var)},
1576 {integer_trail.LowerBound(head_expr.var),
1577 integer_trail.UpperBound(head_expr.var)});
1578 ProcessNewArcRelation(
1579 arc_index, dimension,
1580 {.tail_coeff = -1, .head_coeff = 1, .lhs = lhs, .rhs = rhs});
1581 }
1582
1583 // Returns a relation between two given expressions which is equivalent to a
1584 // given relation "coeff * X \in [lhs, rhs]", where X is the variable used in
1585 // `tail_expression` (or an empty relation if this is not possible).
1586 LocalRelation ComputeArcUnaryRelation(const NodeExpression& tail_expr,
1587 const NodeExpression& head_expr,
1588 IntegerValue coeff, IntegerValue lhs,
1589 IntegerValue rhs) {
1590 DCHECK_NE(coeff, 0);
1591 // A relation a * X \in [lhs, rhs] is equivalent to a (k.a/A) * (A.X+α) + k'
1592 // * (B.Y+β) \in [k.lhs+ɣ, k.rhs+ɣ] relation between the A.X+α and B.Y+β
1593 // terms if A != 0, if the division is exact, and if k' = 0 or B = 0, where
1594 // ɣ=(k.a/A)*α+k'*β. The smallest k > 0 such that k.a/A is integer is
1595 // |A|/gcd(a, A).
1596 if (tail_expr.coeff == 0) return {};
1597 const int64_t k = std::abs(tail_expr.coeff.value()) /
1598 std::gcd(tail_expr.coeff.value(), coeff.value());
1599 // TODO(user): do not add the relation in case of overflow (this can
1600 // happen if the expressions are provided by the user in the model proto).
1601 // If B = 0 we can use any value for k' = head_coeff. We use the opposite of
1602 // tail_coeff to get a relation with +1/-1 coefficients if possible.
1603 const IntegerValue tail_coeff = (k * coeff) / tail_expr.coeff;
1604 const IntegerValue head_coeff = head_expr.coeff != 0 ? 0 : -tail_coeff;
1605 const IntegerValue domain_offset =
1606 tail_coeff * tail_expr.offset + head_coeff * head_expr.offset;
1607 return {
1608 tail_coeff,
1609 head_coeff,
1610 k * lhs + domain_offset,
1611 k * rhs + domain_offset,
1612 };
1613 }
1614
1615 const int num_nodes_;
1616 const int num_arcs_;
1617 absl::Span<const int> tails_;
1618 absl::Span<const int> heads_;
1619 absl::Span<const Literal> literals_;
1620 const BinaryRelationRepository& binary_relation_repository_;
1621
1622 int num_dimensions_;
1623 absl::flat_hash_map<IntegerVariable, int> dimension_by_var_;
1624 absl::flat_hash_map<Literal, int> num_arcs_per_literal_;
1625
1626 // The indices of the binary relations associated with the incoming and
1627 // outgoing arcs of each node, per dimension.
1628 std::vector<std::vector<std::vector<int>>> adjacent_relation_indices_;
1629
1630 // See comments for the same fields in RouteRelationsHelper() that this
1631 // builder creates.
1632 std::vector<AffineExpression> flat_node_dim_expressions_;
1633 std::vector<HeadMinusTailBounds> flat_arc_dim_relations_;
1634 std::vector<IntegerValue> flat_shortest_path_lbs_;
1635};
1636
1637} // namespace
1638
1640 int num_nodes, absl::Span<const int> tails, absl::Span<const int> heads,
1641 absl::Span<const Literal> literals,
1642 const BinaryRelationRepository& binary_relation_repository) {
1644 RouteRelationsBuilder builder(num_nodes, tails, heads, literals, {},
1645 binary_relation_repository);
1646 if (!builder.BuildNodeExpressions()) return result;
1647 result.num_dimensions = builder.num_dimensions();
1649 std::move(builder.flat_node_dim_expressions());
1650 return result;
1651}
1652
1653namespace {
1654// Returns a lower bound of y_expr - x_expr.
1655IntegerValue GetDifferenceLowerBound(
1656 const NodeExpression& x_expr, const NodeExpression& y_expr,
1657 const sat::Relation& r,
1658 const std::pair<IntegerValue, IntegerValue>& x_var_bounds,
1659 const std::pair<IntegerValue, IntegerValue>& y_var_bounds) {
1660 // Let's note x_expr = A.x+α, y_expr = B.x+β, and r = "b.y-a.x in [lhs, rhs]".
1661 // We have x in [lb(x), ub(x)] and y in [lb(y), ub(y)], and we want to find
1662 // the lower bound of y_expr - x_expr = lb(B.x - A.y) + β - α. We have:
1663 // B.y - A.x = [(B-k.b).y + k.b.y] - [(A-k.a).x + k.a.x]
1664 // = (B+k.b).y + (k.a-A).x + k.(b.y - a.x)
1665 // for any k. A lower bound on the first term is (B-k.b).lb(y) if B-k.b >= 0,
1666 // and (B-k.b).ub(y) otherwise. This can be written as:
1667 // (B-k.b).y >= max(0, B-k.b).lb(y) + min(0, B-k.b).ub(y)
1668 // Similarly:
1669 // (k.a-A).x >= max(0, k.a-A).lb(x) + min(0, k.a-A).ub(x)
1670 // And:
1671 // k.(b.y - a.x) >= max(0, k).lhs + min(0, k).rhs
1672 // Hence we get:
1673 // B.y - A.x >= max(0, B-k.b).lb(y) + min(0, B-k.b).ub(y) +
1674 // max(0, k.a-A).lb(x) + min(0, k.a-A).ub(x) +
1675 // max(0, k).lhs + min(0, k).rhs
1676 // The derivative of this expression with respect to k is piecewise constant
1677 // and can only change value around 0, A/a, and B/b. Hence we can compute an
1678 // "interesting" lower bound by taking the maximum of this expression around
1679 // these points, instead of over all k values.
1680 // TODO(user): overflows could happen if the node expressions are
1681 // provided by the user in the model proto.
1682 auto lower_bound = [&](IntegerValue k) {
1683 const IntegerValue y_coeff = y_expr.coeff - k * r.b.coeff;
1684 const IntegerValue x_coeff = k * (-r.a.coeff) - x_expr.coeff;
1685 return y_coeff * (y_coeff >= 0 ? y_var_bounds.first : y_var_bounds.second) +
1686 x_coeff * (x_coeff >= 0 ? x_var_bounds.first : x_var_bounds.second) +
1687 k * (k >= 0 ? r.lhs : r.rhs);
1688 };
1689 const IntegerValue k_x = MathUtil::FloorOfRatio(x_expr.coeff, -r.a.coeff);
1690 const IntegerValue k_y = MathUtil::FloorOfRatio(y_expr.coeff, r.b.coeff);
1691 IntegerValue result = lower_bound(0);
1692 result = std::max(result, lower_bound(k_x));
1693 result = std::max(result, lower_bound(k_x + 1));
1694 result = std::max(result, lower_bound(k_y));
1695 result = std::max(result, lower_bound(k_y + 1));
1696 return result + y_expr.offset - x_expr.offset;
1697}
1698} // namespace
1699
1700std::pair<IntegerValue, IntegerValue> GetDifferenceBounds(
1701 const NodeExpression& x_expr, const NodeExpression& y_expr,
1702 const sat::Relation& r,
1703 const std::pair<IntegerValue, IntegerValue>& x_var_bounds,
1704 const std::pair<IntegerValue, IntegerValue>& y_var_bounds) {
1705 DCHECK_EQ(x_expr.var, r.a.var);
1706 DCHECK_EQ(y_expr.var, r.b.var);
1707 DCHECK_NE(x_expr.var, kNoIntegerVariable);
1708 DCHECK_NE(y_expr.var, kNoIntegerVariable);
1709 DCHECK_NE(x_expr.coeff, 0);
1710 DCHECK_NE(y_expr.coeff, 0);
1711 DCHECK_NE(r.a.coeff, 0);
1712 DCHECK_NE(r.b.coeff, 0);
1713 const IntegerValue lb =
1714 GetDifferenceLowerBound(x_expr, y_expr, r, x_var_bounds, y_var_bounds);
1715 const IntegerValue ub = -GetDifferenceLowerBound(
1716 x_expr.Negated(), y_expr.Negated(), r, x_var_bounds, y_var_bounds);
1717 return {lb, ub};
1718}
1719
1720std::unique_ptr<RouteRelationsHelper> RouteRelationsHelper::Create(
1721 int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
1722 const std::vector<Literal>& literals,
1723 absl::Span<const AffineExpression> flat_node_dim_expressions,
1724 const BinaryRelationRepository& binary_relation_repository, Model* model) {
1725 CHECK(model != nullptr);
1726 if (flat_node_dim_expressions.empty()) return nullptr;
1727 RouteRelationsBuilder builder(num_nodes, tails, heads, literals,
1728 flat_node_dim_expressions,
1729 binary_relation_repository);
1730 builder.BuildArcRelations(model);
1731 builder.BuildShortestPaths(model);
1732 std::unique_ptr<RouteRelationsHelper> helper(new RouteRelationsHelper(
1733 builder.num_dimensions(), std::move(builder.flat_node_dim_expressions()),
1734 std::move(builder.flat_arc_dim_relations()),
1735 std::move(builder.flat_shortest_path_lbs())));
1736 if (VLOG_IS_ON(2)) helper->LogStats();
1737 return helper;
1738}
1739
1741 int num_dimensions, std::vector<AffineExpression> flat_node_dim_expressions,
1742 std::vector<HeadMinusTailBounds> flat_arc_dim_relations,
1743 std::vector<IntegerValue> flat_shortest_path_lbs)
1744 : num_nodes_(flat_node_dim_expressions.size() / num_dimensions),
1745 num_dimensions_(num_dimensions),
1746 flat_node_dim_expressions_(std::move(flat_node_dim_expressions)),
1747 flat_arc_dim_relations_(std::move(flat_arc_dim_relations)),
1748 flat_shortest_path_lbs_(std::move(flat_shortest_path_lbs)) {
1749 DCHECK_GE(num_dimensions_, 1);
1750 DCHECK_EQ(flat_node_dim_expressions_.size(), num_nodes_ * num_dimensions_);
1751}
1752
1754 int arc, int dimension, bool negate_expressions) const {
1755 // TODO(user): If level zero bounds got tighter since we precomputed the
1756 // relation, we might want to recompute it again.
1757 const auto& r = GetArcRelation(arc, dimension);
1758 return negate_expressions ? -r.rhs : r.lhs;
1759}
1760
1762 absl::Span<const int> sorted_arc_indices) {
1763 int new_size = 0;
1764 const int num_arcs = this->num_arcs();
1765 for (int i = 0; i < num_arcs; ++i) {
1766 if (!sorted_arc_indices.empty() && sorted_arc_indices.front() == i) {
1767 sorted_arc_indices.remove_prefix(1);
1768 continue;
1769 }
1770 for (int d = 0; d < num_dimensions_; ++d) {
1771 flat_arc_dim_relations_[new_size++] =
1772 flat_arc_dim_relations_[i * num_dimensions_ + d];
1773 }
1774 }
1775 flat_arc_dim_relations_.resize(new_size);
1776}
1777
1778void RouteRelationsHelper::LogStats() const {
1779 const int num_nodes = this->num_nodes();
1780 const int num_arcs = this->num_arcs();
1781 LOG(INFO) << "Route with " << num_nodes << " nodes and " << num_arcs
1782 << " arcs";
1783 for (int d = 0; d < num_dimensions_; ++d) {
1784 int num_vars = 0;
1785 int num_relations = 0;
1786 for (int i = 0; i < num_nodes; ++i) {
1787 if (!GetNodeExpression(i, d).IsConstant()) ++num_vars;
1788 }
1789 for (int i = 0; i < num_arcs; ++i) {
1790 if (GetArcRelation(i, d).lhs != kMinIntegerValue ||
1791 GetArcRelation(i, d).rhs != kMaxIntegerValue) {
1792 ++num_relations;
1793 }
1794 }
1795 LOG(INFO) << "dimension " << d << ": " << num_vars << " vars and "
1796 << num_relations << " relations";
1797 }
1798}
1799
1800namespace {
1801// Converts a literal index in a model proto to a Literal.
1802Literal ToLiteral(int lit) {
1803 return Literal(BooleanVariable(PositiveRef(lit)), RefIsPositive(lit));
1804}
1805
1806// Converts a variable index in a NodeVariables proto to an IntegerVariable.
1807IntegerVariable ToPositiveIntegerVariable(int i) {
1808 return IntegerVariable(PositiveRef(i) << 1);
1809}
1810
1811// Converts an IntegerVariable to variable indices in a NodeVariables proto.
1812int ToNodeVariableIndex(IntegerVariable var) {
1813 DCHECK(VariableIsPositive(var));
1814 return var.value() >> 1;
1815}
1816
1817// Returns a repository containing a partial view (i.e. without coefficients or
1818// domains) of the enforced linear constraints (of size 2 only) in `model`. This
1819// is the only information needed to infer the mapping from variables to nodes
1820// in routes constraints.
1821BinaryRelationRepository ComputePartialBinaryRelationRepository(
1822 const CpModelProto& model) {
1823 BinaryRelationRepository repository;
1824 for (const ConstraintProto& ct : model.constraints()) {
1825 if (ct.constraint_case() != ConstraintProto::kLinear) continue;
1826 const absl::Span<const int> vars = ct.linear().vars();
1827 if (ct.enforcement_literal().size() != 1 || vars.size() != 2) continue;
1828 if (vars[0] == vars[1]) continue;
1829 repository.AddPartialRelation(ToLiteral(ct.enforcement_literal(0)),
1830 ToPositiveIntegerVariable(vars[0]),
1831 ToPositiveIntegerVariable(vars[1]));
1832 }
1833 repository.Build();
1834 return repository;
1835}
1836
1837// Returns the number of dimensions added to the constraint.
1838int MaybeFillRoutesConstraintNodeExpressions(
1839 RoutesConstraintProto& routes, const BinaryRelationRepository& repository) {
1840 int max_node = 0;
1841 for (const int node : routes.tails()) {
1842 max_node = std::max(max_node, node);
1843 }
1844 for (const int node : routes.heads()) {
1845 max_node = std::max(max_node, node);
1846 }
1847 const int num_nodes = max_node + 1;
1848 std::vector<int> tails(routes.tails().begin(), routes.tails().end());
1849 std::vector<int> heads(routes.heads().begin(), routes.heads().end());
1850 std::vector<Literal> literals;
1851 literals.reserve(routes.literals_size());
1852 for (int lit : routes.literals()) {
1853 literals.push_back(ToLiteral(lit));
1854 }
1855
1857 num_nodes, tails, heads, literals, repository);
1858 if (cumuls.num_dimensions == 0) return 0;
1859
1860 for (int d = 0; d < cumuls.num_dimensions; ++d) {
1861 RoutesConstraintProto::NodeExpressions& dimension =
1862 *routes.add_dimensions();
1863 for (int n = 0; n < num_nodes; ++n) {
1864 const AffineExpression expr = cumuls.GetNodeExpression(n, d);
1865 LinearExpressionProto& node_expr = *dimension.add_exprs();
1866 if (expr.var != kNoIntegerVariable) {
1867 node_expr.add_vars(ToNodeVariableIndex(PositiveVariable(expr.var)));
1868 node_expr.add_coeffs(VariableIsPositive(expr.var)
1869 ? expr.coeff.value()
1870 : -expr.coeff.value());
1871 }
1872 node_expr.set_offset(expr.constant.value());
1873 }
1874 }
1875 return cumuls.num_dimensions;
1876}
1877
1878} // namespace
1879
1881 const CpModelProto& input_model, CpModelProto& output_model) {
1882 std::vector<RoutesConstraintProto*> routes_to_fill;
1883 for (ConstraintProto& ct : *output_model.mutable_constraints()) {
1884 if (ct.constraint_case() != ConstraintProto::kRoutes) continue;
1885 if (!ct.routes().dimensions().empty()) continue;
1886 routes_to_fill.push_back(ct.mutable_routes());
1887 }
1888 if (routes_to_fill.empty()) return {0, 0};
1889
1890 int total_num_dimensions = 0;
1891 const BinaryRelationRepository partial_repository =
1892 ComputePartialBinaryRelationRepository(input_model);
1893 for (RoutesConstraintProto* routes : routes_to_fill) {
1894 total_num_dimensions +=
1895 MaybeFillRoutesConstraintNodeExpressions(*routes, partial_repository);
1896 }
1897 return {static_cast<int>(routes_to_fill.size()), total_num_dimensions};
1898}
1899
1900namespace {
1901
1902class RoutingCutHelper {
1903 public:
1904 RoutingCutHelper(int num_nodes, bool is_route_constraint,
1905 absl::Span<const AffineExpression> flat_node_dim_expressions,
1906 absl::Span<const int> tails, absl::Span<const int> heads,
1907 absl::Span<const Literal> literals, Model* model)
1908 : num_nodes_(num_nodes),
1909 is_route_constraint_(is_route_constraint),
1910 tails_(tails.begin(), tails.end()),
1911 heads_(heads.begin(), heads.end()),
1912 literals_(literals.begin(), literals.end()),
1913 params_(*model->GetOrCreate<SatParameters>()),
1914 trail_(*model->GetOrCreate<Trail>()),
1915 integer_trail_(*model->GetOrCreate<IntegerTrail>()),
1916 binary_relation_repository_(
1917 *model->GetOrCreate<BinaryRelationRepository>()),
1918 random_(model->GetOrCreate<ModelRandomGenerator>()),
1919 encoder_(model->GetOrCreate<IntegerEncoder>()),
1920 in_subset_(num_nodes, false),
1921 self_arc_literal_(num_nodes_),
1922 self_arc_lp_value_(num_nodes_),
1923 nodes_incoming_weight_(num_nodes_),
1924 nodes_outgoing_weight_(num_nodes_),
1925 min_outgoing_flow_helper_(num_nodes, tails_, heads_, literals_, model),
1926 route_relations_helper_(RouteRelationsHelper::Create(
1927 num_nodes, tails_, heads_, literals_, flat_node_dim_expressions,
1928 *model->GetOrCreate<BinaryRelationRepository>(), model)) {}
1929
1930 int num_nodes() const { return num_nodes_; }
1931
1932 bool is_route_constraint() const { return is_route_constraint_; }
1933
1934 // Returns the arcs computed by InitializeForNewLpSolution().
1935 absl::Span<const ArcWithLpValue> relevant_arcs() const {
1936 return relevant_arcs_;
1937 }
1938
1939 // Returns the arcs computed in InitializeForNewLpSolution(), sorted by
1940 // decreasing lp_value.
1941 //
1942 // Note: If is_route_constraint() is true, we do not include arcs from/to
1943 // the depot here.
1944 absl::Span<const std::pair<int, int>> ordered_arcs() const {
1945 return ordered_arcs_;
1946 }
1947
1948 // This should be called each time a new LP solution is available, before
1949 // using the other methods.
1950 void InitializeForNewLpSolution(LinearConstraintManager* manager);
1951
1952 // Merge the relevant arcs (tail, head) and (head, tail) into a single (tail,
1953 // head) arc, with the sum of their lp_values.
1954 absl::Span<const ArcWithLpValue> SymmetrizedRelevantArcs() {
1955 symmetrized_relevant_arcs_ = relevant_arcs_;
1956 SymmetrizeArcs(&symmetrized_relevant_arcs_);
1957
1958 // We exclude arcs from/to the depot when we have a route constraint.
1959 if (is_route_constraint_) {
1960 int new_size = 0;
1961 for (const ArcWithLpValue& arc : symmetrized_relevant_arcs_) {
1962 if (arc.tail == 0 || arc.head == 0) continue;
1963 symmetrized_relevant_arcs_[new_size++] = arc;
1964 }
1965 symmetrized_relevant_arcs_.resize(new_size);
1966 }
1967
1968 return symmetrized_relevant_arcs_;
1969 }
1970
1971 // Try to add an outgoing cut from the given subset.
1972 bool TrySubsetCut(int known_bound, std::string name,
1973 absl::Span<const int> subset,
1974 LinearConstraintManager* manager);
1975
1976 // Returns a bound on the number of vehicle that is valid for this subset and
1977 // all superset. This takes a current valid bound. It might do nothing
1978 // depening on the parameters and just return the initial bound.
1979 int ShortestPathBound(int bound, absl::Span<const int> subset);
1980
1981 // If we look at the symmetrized version (tail <-> head = tail->head +
1982 // head->tail) and we split all the edges between a subset of nodes S and the
1983 // outside into a set A and the other d(S)\A, and |A| is odd, we have a
1984 // constraint of the form:
1985 // "all edge of A at 1" => sum other edges >= 1.
1986 // This is because a cycle or multiple-cycle must go in/out an even number
1987 // of time. This enforced constraint simply linearize to:
1988 // sum_d(S)\A x_e + sum_A (1 - x_e) >= 1.
1989 //
1990 // Given a subset of nodes, it is easy to identify the best subset A of edge
1991 // to consider.
1992 bool TryBlossomSubsetCut(std::string name,
1993 absl::Span<const ArcWithLpValue> symmetrized_edges,
1994 absl::Span<const int> subset,
1995 LinearConstraintManager* manager);
1996
1997 // Try adding cuts to exclude paths in the residual graph corresponding to the
1998 // LP solution which are actually infeasible.
1999 void TryInfeasiblePathCuts(LinearConstraintManager* manager);
2000
2001 private:
2002 // If bool is true it is a "incoming" cut, otherwise "outgoing" one.
2003 struct BestChoice {
2004 int node;
2005 bool incoming;
2006 double weight;
2007 };
2008 BestChoice PickBestNode(absl::Span<const int> cannot_be_first,
2009 absl::Span<const int> cannot_be_last);
2010
2011 // Removes the arcs with a literal fixed at false at level zero. This is
2012 // especially useful to remove fixed self loop.
2013 void FilterFalseArcsAtLevelZero();
2014
2015 // Add a cut of the form Sum_{outgoing arcs from S} lp >= rhs_lower_bound.
2016 //
2017 // Note that we used to also add the same cut for the incoming arcs, but
2018 // because of flow conservation on these problems, the outgoing flow is always
2019 // the same as the incoming flow, so adding this extra cut doesn't seem
2020 // relevant.
2021 bool AddOutgoingCut(LinearConstraintManager* manager, std::string name,
2022 int subset_size, const std::vector<bool>& in_subset,
2023 int64_t rhs_lower_bound);
2024
2025 bool MaybeAddStrongerFlowCut(LinearConstraintManager* manager,
2026 std::string name, int k,
2027 absl::Span<const int> cannot_be_first,
2028 absl::Span<const int> cannot_be_last,
2029 const std::vector<bool>& in_subset);
2030
2031 void GenerateCutsForInfeasiblePaths(
2032 int start_node, int max_path_length,
2033 const CompactVectorVector<int, std::pair<int, double>>&
2034 outgoing_arcs_with_lp_values,
2035 LinearConstraintManager* manager, TopNCuts& top_n_cuts);
2036
2037 const int num_nodes_;
2038 const bool is_route_constraint_;
2039 std::vector<int> tails_;
2040 std::vector<int> heads_;
2041 std::vector<Literal> literals_;
2042 std::vector<ArcWithLpValue> relevant_arcs_;
2043 std::vector<std::pair<int, double>> relevant_arc_indices_;
2044 std::vector<ArcWithLpValue> symmetrized_relevant_arcs_;
2045 std::vector<std::pair<int, int>> ordered_arcs_;
2046
2047 const SatParameters& params_;
2048 const Trail& trail_;
2049 const IntegerTrail& integer_trail_;
2050 const BinaryRelationRepository& binary_relation_repository_;
2051 ModelRandomGenerator* random_;
2052 IntegerEncoder* encoder_;
2053
2054 std::vector<bool> in_subset_;
2055
2056 // Self-arc information, indexed in [0, num_nodes_)
2057 std::vector<int> nodes_with_self_arc_;
2058 std::vector<Literal> self_arc_literal_;
2059 std::vector<double> self_arc_lp_value_;
2060
2061 // Initialized by TrySubsetCut() for each new subset.
2062 std::vector<double> nodes_incoming_weight_;
2063 std::vector<double> nodes_outgoing_weight_;
2064
2065 // Helper to compute bounds on the minimum number of vehicles to serve a
2066 // subset.
2067 MaxBoundedSubsetSum max_bounded_subset_sum_;
2068 MaxBoundedSubsetSumExact max_bounded_subset_sum_exact_;
2069 MinOutgoingFlowHelper min_outgoing_flow_helper_;
2070 std::unique_ptr<RouteRelationsHelper> route_relations_helper_;
2071};
2072
2073void RoutingCutHelper::FilterFalseArcsAtLevelZero() {
2074 if (trail_.CurrentDecisionLevel() != 0) return;
2075
2076 int new_size = 0;
2077 const int size = static_cast<int>(tails_.size());
2078 const VariablesAssignment& assignment = trail_.Assignment();
2079 std::vector<int> removed_arcs;
2080 for (int i = 0; i < size; ++i) {
2081 if (assignment.LiteralIsFalse(literals_[i])) {
2082 removed_arcs.push_back(i);
2083 continue;
2084 }
2085 tails_[new_size] = tails_[i];
2086 heads_[new_size] = heads_[i];
2087 literals_[new_size] = literals_[i];
2088 ++new_size;
2089 }
2090 if (new_size < size) {
2091 tails_.resize(new_size);
2092 heads_.resize(new_size);
2093 literals_.resize(new_size);
2094 if (route_relations_helper_ != nullptr) {
2095 route_relations_helper_->RemoveArcs(removed_arcs);
2096 }
2097 }
2098}
2099
2100void RoutingCutHelper::InitializeForNewLpSolution(
2101 LinearConstraintManager* manager) {
2102 FilterFalseArcsAtLevelZero();
2103
2104 // We will collect only the arcs with a positive lp_values to speed up some
2105 // computation below.
2106 relevant_arcs_.clear();
2107 relevant_arc_indices_.clear();
2108 nodes_with_self_arc_.clear();
2109
2110 // Sort the arcs by non-increasing lp_values.
2111 const auto& lp_values = manager->LpValues();
2112 std::vector<std::pair<double, int>> relevant_arc_by_decreasing_lp_values;
2113 for (int i = 0; i < literals_.size(); ++i) {
2114 const IntegerVariable direct_view = encoder_->GetLiteralView(literals_[i]);
2115 const double lp_value =
2116 direct_view != kNoIntegerVariable
2117 ? lp_values[direct_view]
2118 : 1.0 - lp_values[encoder_->GetLiteralView(literals_[i].Negated())];
2119
2120 // We treat self-edge separately.
2121 // Note also that we do not need to include them in relevant_arcs_.
2122 //
2123 // TODO(user): If there are multiple self-arc, the code should still
2124 // work, but is not ideal.
2125 if (tails_[i] == heads_[i]) {
2126 const int node = tails_[i];
2127 nodes_with_self_arc_.push_back(node);
2128 self_arc_lp_value_[node] = lp_value;
2129 self_arc_literal_[node] = literals_[i];
2130 continue;
2131 }
2132
2133 if (lp_value < 1e-6) continue;
2134 relevant_arcs_.push_back({tails_[i], heads_[i], lp_value});
2135 relevant_arc_indices_.push_back({i, lp_value});
2136 relevant_arc_by_decreasing_lp_values.push_back({lp_value, i});
2137 }
2138 std::sort(relevant_arc_by_decreasing_lp_values.begin(),
2139 relevant_arc_by_decreasing_lp_values.end(),
2140 std::greater<std::pair<double, int>>());
2141
2142 gtl::STLSortAndRemoveDuplicates(&nodes_with_self_arc_);
2143
2144 ordered_arcs_.clear();
2145 for (const auto& [score, arc] : relevant_arc_by_decreasing_lp_values) {
2146 if (is_route_constraint_) {
2147 if (tails_[arc] == 0 || heads_[arc] == 0) continue;
2148 }
2149 ordered_arcs_.push_back({tails_[arc], heads_[arc]});
2150 }
2151
2152 if (absl::GetFlag(FLAGS_cp_model_dump_routes_support_graphs) &&
2153 trail_.CurrentDecisionLevel() == 0) {
2154 static std::atomic<int> counter = 0;
2155 const std::string name =
2156 absl::StrCat(absl::GetFlag(FLAGS_cp_model_dump_prefix),
2157 "support_graph_", counter++, ".pb.txt");
2158 LOG(INFO) << "Dumping routes support graph to '" << name << "'.";
2159 RoutesSupportGraphProto support_graph_proto;
2160 for (auto& [tail, head, lp_value] : relevant_arcs_) {
2161 auto* arc = support_graph_proto.add_arc_lp_values();
2162 arc->set_tail(tail);
2163 arc->set_head(head);
2164 arc->set_lp_value(lp_value);
2165 }
2166 CHECK(WriteModelProtoToFile(support_graph_proto, name));
2167 }
2168}
2169
2170namespace {
2171
2172// Compute the current outgoing/incoming flow out of the subset.
2173// In many cases this will be the same, but not with outside_node_to_ignore
2174// or in case our LP does not contain all the constraints.
2175//
2176// Looping over all arcs can take a significant portion of the running time,
2177// it is why it is faster to do it only on arcs with non-zero lp values which
2178// should be in linear number rather than the total number of arc which can be
2179// quadratic.
2180//
2181// TODO(user): For the symmetric case there is an even faster algo. See if
2182// it can be generalized to the asymmetric one if become needed.
2183// Reference is algo 6.4 of the "The Traveling Salesman Problem" book
2184// mentioned above.
2185std::pair<double, double> GetIncomingAndOutgoingLpFlow(
2186 absl::Span<const ArcWithLpValue> relevant_arcs,
2187 const std::vector<bool>& in_subset, int outside_node_to_ignore = -1) {
2188 double outgoing_flow = 0.0;
2189 double incoming_flow = 0.0;
2190 for (const auto arc : relevant_arcs) {
2191 const bool tail_in = in_subset[arc.tail];
2192 const bool head_in = in_subset[arc.head];
2193 if (tail_in && !head_in) {
2194 if (arc.head == outside_node_to_ignore) continue;
2195 outgoing_flow += arc.lp_value;
2196 } else if (!tail_in && head_in) {
2197 if (arc.tail == outside_node_to_ignore) continue;
2198 incoming_flow += arc.lp_value;
2199 }
2200 }
2201 return {incoming_flow, outgoing_flow};
2202}
2203
2204} // namespace
2205
2206RoutingCutHelper::BestChoice RoutingCutHelper::PickBestNode(
2207 absl::Span<const int> cannot_be_first,
2208 absl::Span<const int> cannot_be_last) {
2209 // Pick the set of candidates with the highest lp weight.
2210 std::vector<BestChoice> bests;
2211 double best_weight = 0.0;
2212 for (const int n : cannot_be_first) {
2213 const double weight = nodes_incoming_weight_[n];
2214 if (bests.empty() || weight > best_weight) {
2215 bests.clear();
2216 bests.push_back({n, true, weight});
2217 best_weight = weight;
2218 } else if (weight == best_weight) {
2219 bests.push_back({n, true, weight});
2220 }
2221 }
2222 for (const int n : cannot_be_last) {
2223 const double weight = nodes_outgoing_weight_[n];
2224 if (bests.empty() || weight > best_weight) {
2225 bests.clear();
2226 bests.push_back({n, false, weight});
2227 best_weight = weight;
2228 } else if (weight == best_weight) {
2229 bests.push_back({n, false, weight});
2230 }
2231 }
2232 if (bests.empty()) return {-1, true};
2233 return bests.size() == 1
2234 ? bests[0]
2235 : bests[absl::Uniform<int>(*random_, 0, bests.size())];
2236}
2237
2238bool RoutingCutHelper::MaybeAddStrongerFlowCut(
2239 LinearConstraintManager* manager, std::string name, int k,
2240 absl::Span<const int> cannot_be_first, absl::Span<const int> cannot_be_last,
2241 const std::vector<bool>& in_subset) {
2242 if (cannot_be_first.empty() && cannot_be_last.empty()) return false;
2243
2244 // We pick the best node out of cannot_be_first/cannot_be_last to derive
2245 // a stronger cut.
2246 const BestChoice best_choice = PickBestNode(cannot_be_first, cannot_be_last);
2247 const int best_node = best_choice.node;
2248 const bool incoming = best_choice.incoming;
2249 CHECK_GE(best_node, 0);
2250
2251 const auto consider_arc = [incoming, &in_subset](int t, int h) {
2252 return incoming ? !in_subset[t] && in_subset[h]
2253 : in_subset[t] && !in_subset[h];
2254 };
2255
2256 // We consider the incoming (resp. outgoing_flow) not arriving at (resp. not
2257 // leaving from) best_node.
2258 //
2259 // This flow must be >= k.
2260 double flow = 0.0;
2261 for (const auto arc : relevant_arcs_) {
2262 if (!consider_arc(arc.tail, arc.head)) continue;
2263 if (arc.tail != best_node && arc.head != best_node) {
2264 flow += arc.lp_value;
2265 }
2266 }
2267
2268 // Add cut flow >= k.
2269 const double kTolerance = 1e-4;
2270 if (flow < static_cast<double>(k) - kTolerance) {
2271 LinearConstraintBuilder cut_builder(encoder_, IntegerValue(k),
2272 kMaxIntegerValue);
2273 for (int i = 0; i < tails_.size(); ++i) {
2274 if (!consider_arc(tails_[i], heads_[i])) continue;
2275 if (tails_[i] != best_node && heads_[i] != best_node) {
2276 CHECK(cut_builder.AddLiteralTerm(literals_[i], IntegerValue(1)));
2277 }
2278 }
2279 return manager->AddCut(cut_builder.Build(), name);
2280 }
2281 return false;
2282}
2283
2284bool RoutingCutHelper::AddOutgoingCut(LinearConstraintManager* manager,
2285 std::string name, int subset_size,
2286 const std::vector<bool>& in_subset,
2287 int64_t rhs_lower_bound) {
2288 // Skip cut if it is not violated.
2289 const auto [in_flow, out_flow] =
2290 GetIncomingAndOutgoingLpFlow(relevant_arcs_, in_subset);
2291 const double out_violation = static_cast<double>(rhs_lower_bound) - out_flow;
2292 const double in_violation = static_cast<double>(rhs_lower_bound) - in_flow;
2293 if (out_violation <= 1e-3 && in_violation <= 1e-3) return false;
2294
2295 // We create the cut and rely on AddCut() for computing its efficacy and
2296 // rejecting it if it is bad.
2297 LinearConstraintBuilder outgoing(encoder_, IntegerValue(rhs_lower_bound),
2298 kMaxIntegerValue);
2299 LinearConstraintBuilder incoming(encoder_, IntegerValue(rhs_lower_bound),
2300 kMaxIntegerValue);
2301
2302 // Rather than doing two loops, we initialize the cuts right away, even if
2303 // only one of them will be used.
2304 for (int i = 0; i < tails_.size(); ++i) {
2305 const bool tail_in = in_subset[tails_[i]];
2306 const bool head_in = in_subset[heads_[i]];
2307 if (tail_in && !head_in) {
2308 CHECK(outgoing.AddLiteralTerm(literals_[i], IntegerValue(1)));
2309 } else if (!tail_in && head_in) {
2310 CHECK(incoming.AddLiteralTerm(literals_[i], IntegerValue(1)));
2311 }
2312 }
2313
2314 // As arcs get fixed (this happens a lot in LNS subproblems), even if the
2315 // incoming flow is the same as the outgoing flow, the number of incoming arcs
2316 // might be widely different from the one of outgoing arcs. We prefer to pick
2317 // the sparser cut.
2318 const double out_efficacy = out_violation / std::sqrt(outgoing.NumTerms());
2319 const double in_efficacy = in_violation / std::sqrt(incoming.NumTerms());
2320
2321 // Select the best option between outgoing and incoming.
2322 LinearConstraintBuilder& cut_builder =
2323 (out_efficacy >= in_efficacy) ? outgoing : incoming;
2324
2325 // A node is said to be optional if it can be excluded from the subcircuit,
2326 // in which case there is a self-loop on that node.
2327 // If there are optional nodes, use extended formula:
2328 // sum(cut) >= 1 - optional_loop_in - optional_loop_out
2329 // where optional_loop_in's node is in subset, optional_loop_out's is out.
2330 // TODO(user): Favor optional loops fixed to zero at root.
2331 int num_optional_nodes_in = 0;
2332 int num_optional_nodes_out = 0;
2333 int optional_loop_in = -1;
2334 int optional_loop_out = -1;
2335 for (const int n : nodes_with_self_arc_) {
2336 if (in_subset[n]) {
2337 num_optional_nodes_in++;
2338 if (optional_loop_in == -1 ||
2339 self_arc_lp_value_[n] < self_arc_lp_value_[optional_loop_in]) {
2340 optional_loop_in = n;
2341 }
2342 } else {
2343 num_optional_nodes_out++;
2344 if (optional_loop_out == -1 ||
2345 self_arc_lp_value_[n] < self_arc_lp_value_[optional_loop_out]) {
2346 optional_loop_out = n;
2347 }
2348 }
2349 }
2350
2351 // This just makes sure we don't call this with a bound > 1 if there is
2352 // optional node inside the subset.
2353 CHECK(rhs_lower_bound == 1 || num_optional_nodes_in == 0);
2354
2355 // Support optional nodes if any.
2356 if (num_optional_nodes_in + num_optional_nodes_out > 0) {
2357 // When all optionals of one side are excluded in lp solution, no cut.
2358 if (num_optional_nodes_in == subset_size &&
2359 (optional_loop_in == -1 ||
2360 self_arc_lp_value_[optional_loop_in] > 1.0 - 1e-6)) {
2361 return false;
2362 }
2363 if (num_optional_nodes_out == num_nodes_ - subset_size &&
2364 (optional_loop_out == -1 ||
2365 self_arc_lp_value_[optional_loop_out] > 1.0 - 1e-6)) {
2366 return false;
2367 }
2368
2369 // There is no mandatory node in subset, add optional_loop_in.
2370 if (num_optional_nodes_in == subset_size) {
2371 CHECK_EQ(rhs_lower_bound, 1);
2372 CHECK(cut_builder.AddLiteralTerm(self_arc_literal_[optional_loop_in],
2373 IntegerValue(1)));
2374 }
2375
2376 // There is no mandatory node out of subset, add optional_loop_out.
2377 if (num_optional_nodes_out == num_nodes_ - subset_size) {
2378 CHECK_EQ(rhs_lower_bound, 1);
2379 CHECK(cut_builder.AddLiteralTerm(self_arc_literal_[optional_loop_out],
2380 IntegerValue(1)));
2381 }
2382 }
2383
2384 return manager->AddCut(cut_builder.Build(), name);
2385}
2386
2387int RoutingCutHelper::ShortestPathBound(int bound,
2388 absl::Span<const int> subset) {
2389 if (!is_route_constraint_) return bound;
2390 if (!nodes_with_self_arc_.empty()) return bound;
2391 if (route_relations_helper_ == nullptr) return bound;
2392 if (!route_relations_helper_->HasShortestPathsInformation()) return bound;
2393 if (subset.size() >
2394 params_.routing_cut_subset_size_for_shortest_paths_bound()) {
2395 return bound;
2396 }
2397
2398 while (bound < subset.size()) {
2399 if (min_outgoing_flow_helper_.SubsetMightBeServedWithKRoutes(
2400 bound, subset, route_relations_helper_.get())) {
2401 break;
2402 }
2403 bound += 1;
2404 }
2405
2406 return bound;
2407}
2408
2409bool RoutingCutHelper::TrySubsetCut(int known_bound, std::string name,
2410 absl::Span<const int> subset,
2411 LinearConstraintManager* manager) {
2412 DCHECK_GE(subset.size(), 1);
2413 DCHECK_LT(subset.size(), num_nodes_);
2414
2415 // Do some initialization.
2416 in_subset_.assign(num_nodes_, false);
2417 for (const int n : subset) {
2418 in_subset_[n] = true;
2419
2420 // We should aways generate subset without the depot for the route
2421 // constraint. This is because we remove arcs from/to the depot in the
2422 // tree based heuristics to generate all the subsets we try.
2423 if (is_route_constraint_) CHECK_NE(n, 0);
2424 }
2425
2426 // For now we can only apply fancy route cuts if all nodes in subset are
2427 // mandatory.
2428 bool all_subset_nodes_are_mandatory = true;
2429 if (is_route_constraint_) {
2430 for (const int n : nodes_with_self_arc_) {
2431 if (in_subset_[n]) {
2432 all_subset_nodes_are_mandatory = false;
2433 break;
2434 }
2435 }
2436 }
2437
2438 // The TSP case is "easy".
2439 //
2440 // TODO(user): Turn on some of the automatic detection for circuit constraint.
2441 // Even if we are looking for a full circuit of the mandatory nodes, some
2442 // side-constraint might require to go in and out of a subset more than once.
2443 //
2444 // TODO(user): deal with non-mandatory node in the route constraint?
2445 if (!is_route_constraint_ || !all_subset_nodes_are_mandatory) {
2446 return AddOutgoingCut(manager, name, subset.size(), in_subset_,
2447 /*rhs_lower_bound=*/1);
2448 }
2449
2450 // Compute the LP weight from/to the subset, and also the LP weight incoming
2451 // from outside to a given node in a subset or from the subset to a given node
2452 // outside (similarly for the outgoing case).
2453 double total_outgoing_weight = 0.0;
2454 double total_incoming_weight = 0.0;
2455 nodes_incoming_weight_.assign(num_nodes_, 0);
2456 nodes_outgoing_weight_.assign(num_nodes_, 0);
2457 for (const auto arc : relevant_arcs_) {
2458 if (in_subset_[arc.tail] != in_subset_[arc.head]) {
2459 if (in_subset_[arc.tail]) {
2460 total_outgoing_weight += arc.lp_value;
2461 } else {
2462 total_incoming_weight += arc.lp_value;
2463 }
2464 nodes_outgoing_weight_[arc.tail] += arc.lp_value;
2465 nodes_incoming_weight_[arc.head] += arc.lp_value;
2466 }
2467 }
2468
2469 BestBoundHelper best_bound(name);
2470 best_bound.Update(1, "");
2471 best_bound.Update(known_bound, "Hierarchical");
2472
2473 // Bounds inferred automatically from the enforced binary relation of the
2474 // model.
2475 //
2476 // TODO(user): This is still not as good as the "capacity" bounds below in
2477 // some cases. Fix! we should be able to use the same relation to infer the
2478 // capacity bounds somehow.
2479 if (route_relations_helper_ != nullptr) {
2480 min_outgoing_flow_helper_.ComputeDimensionBasedMinOutgoingFlow(
2481 subset, *route_relations_helper_, &best_bound);
2482 }
2483 if (subset.size() <
2484 params_.routing_cut_subset_size_for_tight_binary_relation_bound()) {
2485 const int bound =
2486 min_outgoing_flow_helper_.ComputeTightMinOutgoingFlow(subset);
2487 best_bound.Update(bound, "AutomaticTight");
2488 } else if (subset.size() <
2489 params_.routing_cut_subset_size_for_binary_relation_bound()) {
2490 const int bound = min_outgoing_flow_helper_.ComputeMinOutgoingFlow(subset);
2491 best_bound.Update(bound, "Automatic");
2492 }
2493
2494 if (subset.size() <=
2495 params_.routing_cut_subset_size_for_exact_binary_relation_bound()) {
2496 // Note that we only look for violated cuts, but also in order to not try
2497 // all nodes from the subset, we only look at large enough incoming/outgoing
2498 // weight.
2499 //
2500 // TODO(user): We could easily look for an arc that cannot be used to enter
2501 // or leave a subset rather than a node. This should lead to tighter bound,
2502 // and more cuts (that would just ignore that arc rather than all the arcs
2503 // entering/leaving from a node).
2504 const double threshold = static_cast<double>(best_bound.bound()) - 1e-4;
2505 for (const int n : subset) {
2506 if (nodes_incoming_weight_[n] > 0.1 &&
2507 total_incoming_weight - nodes_incoming_weight_[n] < threshold &&
2508 !min_outgoing_flow_helper_.SubsetMightBeServedWithKRoutes(
2509 best_bound.bound(), subset, nullptr, manager, /*special_node=*/n,
2510 /*use_forward_direction=*/true)) {
2511 best_bound.Update(best_bound.bound(), "", {n}, {});
2512 }
2513 if (nodes_outgoing_weight_[n] > 0.1 &&
2514 total_outgoing_weight - nodes_outgoing_weight_[n] < threshold &&
2515 !min_outgoing_flow_helper_.SubsetMightBeServedWithKRoutes(
2516 best_bound.bound(), subset, nullptr, manager, /*special_node=*/n,
2517 /*use_forward_direction=*/false)) {
2518 best_bound.Update(best_bound.bound(), "", {}, {n});
2519 }
2520 }
2521 }
2522
2523 if (MaybeAddStrongerFlowCut(manager,
2524 absl::StrCat(best_bound.name(), "Lifted"),
2525 best_bound.bound(), best_bound.CannotBeFirst(),
2526 best_bound.CannotBeLast(), in_subset_)) {
2527 return true;
2528 }
2529
2530 return AddOutgoingCut(manager, best_bound.name(), subset.size(), in_subset_,
2531 /*rhs_lower_bound=*/best_bound.bound());
2532}
2533
2534bool RoutingCutHelper::TryBlossomSubsetCut(
2535 std::string name, absl::Span<const ArcWithLpValue> symmetrized_edges,
2536 absl::Span<const int> subset, LinearConstraintManager* manager) {
2537 DCHECK_GE(subset.size(), 1);
2538 DCHECK_LT(subset.size(), num_nodes_);
2539
2540 // Initialize "in_subset".
2541 for (const int n : subset) in_subset_[n] = true;
2542 auto cleanup = ::absl::MakeCleanup([subset, this]() {
2543 for (const int n : subset) in_subset_[n] = false;
2544 });
2545
2546 // The heuristic assumes non-duplicates arcs, otherwise they are all bundled
2547 // together in the same symmetric edge, and the result is probably wrong.
2548 absl::flat_hash_set<std::pair<int, int>> special_edges;
2549 int num_inverted = 0;
2550 double sum_inverted = 0.0;
2551 double sum = 0.0;
2552 double best_change = 1.0;
2553 ArcWithLpValue const* best_swap = nullptr;
2554 for (const ArcWithLpValue& arc : symmetrized_edges) {
2555 if (in_subset_[arc.tail] == in_subset_[arc.head]) continue;
2556
2557 if (arc.lp_value > 0.5) {
2558 ++num_inverted;
2559 special_edges.insert({arc.tail, arc.head});
2560 sum_inverted += 1.0 - arc.lp_value;
2561 } else {
2562 sum += arc.lp_value;
2563 }
2564
2565 const double change = std::abs(2 * arc.lp_value - 1.0);
2566 if (change < best_change) {
2567 best_change = change;
2568 best_swap = &arc;
2569 }
2570 }
2571
2572 // If the we don't have an odd number, we move the best edge from one set to
2573 // the other.
2574 if (num_inverted % 2 == 0) {
2575 if (best_swap == nullptr) return false;
2576 if (special_edges.contains({best_swap->tail, best_swap->head})) {
2577 --num_inverted;
2578 special_edges.erase({best_swap->tail, best_swap->head});
2579 sum_inverted -= (1.0 - best_swap->lp_value);
2580 sum += best_swap->lp_value;
2581 } else {
2582 ++num_inverted;
2583 special_edges.insert({best_swap->tail, best_swap->head});
2584 sum_inverted += (1.0 - best_swap->lp_value);
2585 sum -= best_swap->lp_value;
2586 }
2587 }
2588 if (sum + sum_inverted > 0.99) return false;
2589
2590 // For the route constraint, it is actually allowed to have circuit of size
2591 // 2, so the reasoning is wrong if one of the edges touches the depot.
2592 if (!is_route_constraint_) {
2593 for (const auto [tail, head] : special_edges) {
2594 if (tail == 0) return false;
2595 }
2596 }
2597
2598 // If there is just one special edge, and all other node can be ignored, then
2599 // the reasonning is wrong too since we can have a 2-cycle. In that case
2600 // we enforce the constraint when an extra self-loop literal is at zero.
2601 int best_optional_index = -1;
2602 if (special_edges.size() == 1) {
2603 int num_other_optional = 0;
2604 const auto [special_tail, special_head] = *special_edges.begin();
2605 for (const int n : nodes_with_self_arc_) {
2606 if (n != special_head && n != special_tail) {
2607 ++num_other_optional;
2608 if (best_optional_index == -1 ||
2609 self_arc_lp_value_[n] < self_arc_lp_value_[best_optional_index]) {
2610 best_optional_index = n;
2611 }
2612 }
2613 }
2614 if (num_other_optional + 2 < num_nodes_) best_optional_index = -1;
2615 }
2616
2617 // Try to generate the cut.
2618 //
2619 // We deal with the corner case with duplicate arcs, or just one side of a
2620 // "symmetric" edge present.
2621 int num_actual_inverted = 0;
2622 absl::flat_hash_set<std::pair<int, int>> processed_arcs;
2623 LinearConstraintBuilder builder(encoder_, IntegerValue(1), kMaxIntegerValue);
2624
2625 // Add extra self-loop at zero enforcement if needed.
2626 // TODO(user): Can we deal with the 2-cycle case in a better way?
2627 if (best_optional_index != -1) {
2628 absl::StrAppend(&name, "_opt");
2629
2630 // This is tricky: The normal cut assume x_e <= 1, but in case of a single
2631 // 2 cycle, x_e can be equal to 2. So we need a coeff of 2 to disable that
2632 // cut.
2633 CHECK(builder.AddLiteralTerm(self_arc_literal_[best_optional_index],
2634 IntegerValue(2)));
2635 }
2636
2637 for (int i = 0; i < tails_.size(); ++i) {
2638 if (tails_[i] == heads_[i]) continue;
2639 if (in_subset_[tails_[i]] == in_subset_[heads_[i]]) continue;
2640
2641 const std::pair<int, int> key = {tails_[i], heads_[i]};
2642 const std::pair<int, int> r_key = {heads_[i], tails_[i]};
2643 const std::pair<int, int> s_key = std::min(key, r_key);
2644 if (special_edges.contains(s_key) && !processed_arcs.contains(key)) {
2645 processed_arcs.insert(key);
2646 CHECK(builder.AddLiteralTerm(literals_[i], IntegerValue(-1)));
2647 if (!processed_arcs.contains(r_key)) {
2648 ++num_actual_inverted;
2649 }
2650 continue;
2651 }
2652
2653 // Normal edge.
2654 CHECK(builder.AddLiteralTerm(literals_[i], IntegerValue(1)));
2655 }
2656 builder.AddConstant(IntegerValue(num_actual_inverted));
2657 if (num_actual_inverted % 2 == 0) return false;
2658
2659 return manager->AddCut(builder.Build(), name);
2660}
2661
2662void RoutingCutHelper::TryInfeasiblePathCuts(LinearConstraintManager* manager) {
2663 const int max_path_length = params_.routing_cut_max_infeasible_path_length();
2664 if (max_path_length <= 0) return;
2665
2666 // The outgoing arcs of the residual graph, per tail node, except self arcs
2667 // and arcs to or from the depot.
2668 std::vector<int> keys;
2669 std::vector<std::pair<int, double>> values;
2670 for (const auto [arc_index, lp_value] : relevant_arc_indices_) {
2671 const int tail = tails_[arc_index];
2672 const int head = heads_[arc_index];
2673
2674 // Self arcs are filtered out in InitializeForNewLpSolution().
2675 DCHECK_NE(tail, head);
2676 if (tail == 0 || head == 0) continue;
2677 keys.push_back(tail);
2678 values.push_back({arc_index, lp_value});
2679 }
2680 CompactVectorVector<int, std::pair<int, double>> outgoing_arcs_with_lp_values;
2681 outgoing_arcs_with_lp_values.ResetFromFlatMapping(keys, values, num_nodes_);
2682
2683 TopNCuts top_n_cuts(5);
2684 for (int i = 1; i < num_nodes_; ++i) {
2685 GenerateCutsForInfeasiblePaths(
2686 i, max_path_length, outgoing_arcs_with_lp_values, manager, top_n_cuts);
2687 }
2688 top_n_cuts.TransferToManager(manager);
2689}
2690
2691// Note that it makes little sense to use reduced cost here as the LP solution
2692// should likely satisfy the current LP optimality equation.
2693void RoutingCutHelper::GenerateCutsForInfeasiblePaths(
2694 int start_node, int max_path_length,
2695 const CompactVectorVector<int, std::pair<int, double>>&
2696 outgoing_arcs_with_lp_values,
2697 LinearConstraintManager* manager, TopNCuts& top_n_cuts) {
2698 // Performs a DFS on the residual graph, starting from the given node. To save
2699 // stack space, we use a manual stack on the heap instead of using recursion.
2700 // Each stack element corresponds to a feasible path from `start_node` to
2701 // `last_node`, included.
2702 struct State {
2703 // The last node of the path.
2704 int last_node;
2705
2706 // The sum of the lp values of the path's arcs.
2707 double sum_of_lp_values = 0.0;
2708
2709 // Variable lower bounds that can be inferred by assuming that the path's
2710 // arc literals are true (with the binary relation repository).
2711 absl::flat_hash_map<IntegerVariable, IntegerValue> bounds;
2712
2713 // The index of the next outgoing arc of `last_node` to explore (in
2714 // outgoing_arcs_with_lp_values[last_node]).
2715 int iteration = 0;
2716 };
2717 std::stack<State> states;
2718 // Whether each node is on the current path. The current path is the one
2719 // corresponding to the top stack element.
2720 std::vector<bool> path_nodes(num_nodes_, false);
2721 // The literals corresponding to the current path's arcs.
2722 std::vector<Literal> path_literals;
2723
2724 path_nodes[start_node] = true;
2725 states.push({start_node});
2726 while (!states.empty()) {
2727 State& state = states.top();
2728 DCHECK_EQ(states.size(), path_literals.size() + 1);
2729 // The number of arcs in the current path.
2730 const int path_length = static_cast<int>(path_literals.size());
2731 const auto& outgoing_arcs = outgoing_arcs_with_lp_values[state.last_node];
2732
2733 bool new_state_pushed = false;
2734 while (state.iteration < outgoing_arcs.size()) {
2735 const auto [arc_index, lp_value] = outgoing_arcs[state.iteration++];
2736 State next_state;
2737 next_state.last_node = heads_[arc_index];
2738 next_state.sum_of_lp_values = state.sum_of_lp_values + lp_value;
2739 next_state.iteration = 0;
2740 // We only consider simple paths.
2741 if (path_nodes[next_state.last_node]) continue;
2742 // If a path of length l is infeasible, we can exclude it by adding a cut
2743 // saying the sum of its arcs literals must be strictly less than l. But
2744 // this is only beneficial if the sum of the corresponding LP values is
2745 // currently greater than l. If it is not, and since it cannot become
2746 // greater for some extension of the path (the LP sum increases at most by
2747 // 1 for each arc, but l increases exactly by 1), we can stop the search
2748 // here.
2749 // Note: "<= path_length" is equivalent to "< l" (l = path_length + 1).
2750 if (next_state.sum_of_lp_values <=
2751 static_cast<double>(path_length) + 1e-4) {
2752 continue;
2753 }
2754
2755 const Literal next_literal = literals_[arc_index];
2756 next_state.bounds = state.bounds;
2757 if (binary_relation_repository_.PropagateLocalBounds(
2758 integer_trail_, next_literal, state.bounds, &next_state.bounds)) {
2759 // Do not explore "long" paths to keep the search time bounded.
2760 if (path_length < max_path_length) {
2761 path_nodes[next_state.last_node] = true;
2762 path_literals.push_back(literals_[arc_index]);
2763 states.push(std::move(next_state));
2764 new_state_pushed = true;
2765 break;
2766 }
2767 } else if (path_length > 0) {
2768 LinearConstraintBuilder cut_builder(encoder_, 0, path_length);
2769 for (const Literal literal : path_literals) {
2770 CHECK(cut_builder.AddLiteralTerm(literal, IntegerValue(1)));
2771 }
2772 CHECK(cut_builder.AddLiteralTerm(next_literal, IntegerValue(1)));
2773 top_n_cuts.AddCut(cut_builder.Build(), "CircuitInfeasiblePath",
2774 manager->LpValues());
2775 }
2776 }
2777 if (!new_state_pushed) {
2778 if (!path_literals.empty()) {
2779 path_literals.pop_back();
2780 }
2781 path_nodes[state.last_node] = false;
2782 states.pop();
2783 }
2784 }
2785}
2786
2787} // namespace
2788
2790 absl::Span<const std::pair<int, int>> arcs,
2791 int stop_at_num_components,
2792 std::vector<int>* subset_data,
2793 std::vector<absl::Span<const int>>* subsets) {
2794 // We will do a union-find by adding one by one the arc of the lp solution
2795 // in the order above. Every intermediate set during this construction will
2796 // be a candidate for a cut.
2797 //
2798 // In parallel to the union-find, to efficiently reconstruct these sets (at
2799 // most num_nodes), we construct a "decomposition forest" of the different
2800 // connected components. Note that we don't exploit any asymmetric nature of
2801 // the graph here. This is exactly the algo 6.3 in the book above.
2802 int num_components = num_nodes;
2803 std::vector<int> parent(num_nodes);
2804 std::vector<int> root(num_nodes);
2805 for (int i = 0; i < num_nodes; ++i) {
2806 parent[i] = i;
2807 root[i] = i;
2808 }
2809 auto get_root_and_compress_path = [&root](int node) {
2810 int r = node;
2811 while (root[r] != r) r = root[r];
2812 while (root[node] != r) {
2813 const int next = root[node];
2814 root[node] = r;
2815 node = next;
2816 }
2817 return r;
2818 };
2819 for (const auto& [initial_tail, initial_head] : arcs) {
2820 if (num_components <= stop_at_num_components) break;
2821 const int tail = get_root_and_compress_path(initial_tail);
2822 const int head = get_root_and_compress_path(initial_head);
2823 if (tail != head) {
2824 // Update the decomposition forest, note that the number of nodes is
2825 // growing.
2826 const int new_node = parent.size();
2827 parent.push_back(new_node);
2828 parent[head] = new_node;
2829 parent[tail] = new_node;
2830 --num_components;
2831
2832 // It is important that the union-find representative is the same node.
2833 root.push_back(new_node);
2834 root[head] = new_node;
2835 root[tail] = new_node;
2836 }
2837 }
2838
2839 // For each node in the decomposition forest, try to add a cut for the set
2840 // formed by the nodes and its children. To do that efficiently, we first
2841 // order the nodes so that for each node in a tree, the set of children forms
2842 // a consecutive span in the subset_data vector. This vector just lists the
2843 // nodes in the "pre-order" graph traversal order. The Spans will point inside
2844 // the subset_data vector, it is why we initialize it once and for all.
2845 ExtractAllSubsetsFromForest(parent, subset_data, subsets,
2846 /*node_limit=*/num_nodes);
2847}
2848
2849void ExtractAllSubsetsFromForest(absl::Span<const int> parent,
2850 std::vector<int>* subset_data,
2851 std::vector<absl::Span<const int>>* subsets,
2852 int node_limit) {
2853 // To not reallocate memory since we need the span to point inside this
2854 // vector, we resize subset_data right away.
2855 int out_index = 0;
2856 const int num_nodes = parent.size();
2857 subset_data->resize(std::min(num_nodes, node_limit));
2858 subsets->clear();
2859
2860 // Starts by creating the corresponding graph and find the root.
2861 util::StaticGraph<int> graph(num_nodes, num_nodes - 1);
2862 for (int i = 0; i < num_nodes; ++i) {
2863 if (parent[i] != i) {
2864 graph.AddArc(parent[i], i);
2865 }
2866 }
2867 graph.Build();
2868
2869 // Perform a dfs on the rooted tree.
2870 // The subset_data will just be the node in post-order.
2871 std::vector<int> subtree_starts(num_nodes, -1);
2872 std::vector<int> stack;
2873 stack.reserve(num_nodes);
2874 for (int i = 0; i < parent.size(); ++i) {
2875 if (parent[i] != i) continue;
2876
2877 stack.push_back(i); // root.
2878 while (!stack.empty()) {
2879 const int node = stack.back();
2880
2881 // The node was already explored, output its subtree and pop it.
2882 if (subtree_starts[node] >= 0) {
2883 stack.pop_back();
2884 if (node < node_limit) {
2885 (*subset_data)[out_index++] = node;
2886 }
2887 const int start = subtree_starts[node];
2888 const int size = out_index - start;
2889 subsets->push_back(absl::MakeSpan(&(*subset_data)[start], size));
2890 continue;
2891 }
2892
2893 // Explore.
2894 subtree_starts[node] = out_index;
2895 for (const int child : graph[node]) {
2896 stack.push_back(child);
2897 }
2898 }
2899 }
2900}
2901
2902std::vector<int> ComputeGomoryHuTree(
2903 int num_nodes, absl::Span<const ArcWithLpValue> relevant_arcs) {
2904 // Initialize the graph. Note that we use only arcs with a relevant lp
2905 // value, so this should be small in practice.
2906 SimpleMaxFlow max_flow;
2907 for (const auto& [tail, head, lp_value] : relevant_arcs) {
2908 max_flow.AddArcWithCapacity(tail, head, std::round(1.0e6 * lp_value));
2909 max_flow.AddArcWithCapacity(head, tail, std::round(1.0e6 * lp_value));
2910 }
2911
2912 // Compute an equivalent max-flow tree, according to the paper.
2913 // This version should actually produce a Gomory-Hu cut tree.
2914 std::vector<int> min_cut_subset;
2915 std::vector<int> parent(num_nodes, 0);
2916 for (int s = 1; s < num_nodes; ++s) {
2917 const int t = parent[s];
2918 if (max_flow.Solve(s, t) != SimpleMaxFlow::OPTIMAL) break;
2919 max_flow.GetSourceSideMinCut(&min_cut_subset);
2920 bool parent_of_t_in_subset = false;
2921 for (const int i : min_cut_subset) {
2922 if (i == parent[t]) parent_of_t_in_subset = true;
2923 if (i != s && parent[i] == t) parent[i] = s;
2924 }
2925 if (parent_of_t_in_subset) {
2926 parent[s] = parent[t];
2927 parent[t] = s;
2928 }
2929 }
2930
2931 return parent;
2932}
2933
2934void SymmetrizeArcs(std::vector<ArcWithLpValue>* arcs) {
2935 for (ArcWithLpValue& arc : *arcs) {
2936 if (arc.tail <= arc.head) continue;
2937 std::swap(arc.tail, arc.head);
2938 }
2939 std::sort(arcs->begin(), arcs->end(),
2940 [](const ArcWithLpValue& a, const ArcWithLpValue& b) {
2941 return std::tie(a.tail, a.head) < std::tie(b.tail, b.head);
2942 });
2943
2944 int new_size = 0;
2945 int last_tail = -1;
2946 int last_head = -1;
2947 for (const ArcWithLpValue& arc : *arcs) {
2948 if (arc.tail == last_tail && arc.head == last_head) {
2949 (*arcs)[new_size - 1].lp_value += arc.lp_value;
2950 continue;
2951 }
2952 (*arcs)[new_size++] = arc;
2953 last_tail = arc.tail;
2954 last_head = arc.head;
2955 }
2956 arcs->resize(new_size);
2957}
2958
2959// Processes each subsets and add any violated cut.
2960// Returns the number of added cuts.
2961int TryAllSubsets(std::string cut_name, absl::Span<const int> subset_data,
2962 std::vector<absl::Span<const int>> subsets,
2963 RoutingCutHelper& helper, LinearConstraintManager* manager) {
2964 const int num_nodes = subset_data.size();
2965
2966 // This exploit the fact that all subsets point into subset_data of size
2967 // num_nodes. Moreover, if S1 is included in S2, we will always process S1
2968 // before S2.
2969 //
2970 // TODO(user): There is probably a better way than doing a linear scan per
2971 // subset.
2972 std::vector<bool> cut_was_added(num_nodes, false);
2973 std::vector<int> shortest_path_lb(num_nodes, 0);
2974
2975 int num_added = 0;
2976 for (const absl::Span<const int> subset : subsets) {
2977 if (subset.size() <= 1) continue;
2978 if (subset.size() == num_nodes) continue;
2979
2980 // we exploit the tree structure of the subsets to not add a cut
2981 // for a larger subset if we added a cut from one included in it. Also,
2982 // since the "shortest path" lower bound is valid for any superset, we use
2983 // it to have a starting lb for the number of vehicles.
2984 //
2985 // TODO(user): Currently if we add too many not so relevant cuts, our
2986 // generic MIP cut heuritic are way too slow on TSP/VRP problems.
2987 int lb_for_that_subset = 1;
2988 bool included_cut_was_added = false;
2989 const int start = static_cast<int>(subset.data() - subset_data.data());
2990 for (int i = start; i < subset.size(); ++i) {
2991 lb_for_that_subset = std::max(lb_for_that_subset, shortest_path_lb[i]);
2992 if (cut_was_added[i]) {
2993 included_cut_was_added = true;
2994 }
2995 }
2996
2997 // If the subset is small enough and the parameters ask for it, compute
2998 // a lower bound on the number of vehicle for that subset. This uses
2999 // "shortest path bounds" and thus that bounds will also be valid for
3000 // any superset !
3001 lb_for_that_subset = helper.ShortestPathBound(lb_for_that_subset, subset);
3002 shortest_path_lb[start] = lb_for_that_subset;
3003
3004 // Note that we still try the num_nodes - 1 subset cut as that gives
3005 // directly a lower bound on the number of vehicles.
3006 if (included_cut_was_added && subset.size() + 1 < num_nodes) continue;
3007 if (helper.TrySubsetCut(lb_for_that_subset, cut_name, subset, manager)) {
3008 ++num_added;
3009 cut_was_added[start] = true;
3010 }
3011 }
3012
3013 return num_added;
3014}
3015
3016// We roughly follow the algorithm described in section 6 of "The Traveling
3017// Salesman Problem, A computational Study", David L. Applegate, Robert E.
3018// Bixby, Vasek Chvatal, William J. Cook.
3019//
3020// Note that this is mainly a "symmetric" case algo, but it does still work for
3021// the asymmetric case.
3022void SeparateSubtourInequalities(RoutingCutHelper& helper,
3023 LinearConstraintManager* manager) {
3024 const int num_nodes = helper.num_nodes();
3025 if (num_nodes <= 2) return;
3026
3027 std::vector<int> subset_data;
3028 std::vector<absl::Span<const int>> subsets;
3029 GenerateInterestingSubsets(num_nodes, helper.ordered_arcs(),
3030 /*stop_at_num_components=*/2, &subset_data,
3031 &subsets);
3032
3033 // For a routing problem, we always try with all nodes but the root as this
3034 // gives a global lower bound on the number of vehicles. Note that usually
3035 // the arcs with non-zero lp values should connect everything, but that only
3036 // happen after many cuts on large problems.
3037 if (helper.is_route_constraint()) {
3038 subsets.push_back(absl::MakeSpan(&subset_data[1], num_nodes - 1));
3039 }
3040
3041 int num_added =
3042 TryAllSubsets("Circuit", subset_data, subsets, helper, manager);
3043
3044 // If there were no cut added by the heuristic above, we try exact separation.
3045 //
3046 // With n-1 max_flow from a source to all destination, we can get the global
3047 // min-cut. Here, we use a slightly more advanced algorithm that will find a
3048 // min-cut for all possible pair of nodes. This is achieved by computing a
3049 // Gomory-Hu tree, still with n-1 max flow call.
3050 //
3051 // Note(user): Compared to any min-cut, these cut have some nice properties
3052 // since they are "included" in each other. This might help with combining
3053 // them within our generic IP cuts framework.
3054 //
3055 // TODO(user): I had an older version that tried the n-cuts generated during
3056 // the course of the algorithm. This could also be interesting. But it is
3057 // hard to tell with our current benchmark setup.
3058 if (num_added != 0) return;
3059 absl::Span<const ArcWithLpValue> symmetrized_relevant_arcs =
3060 helper.SymmetrizedRelevantArcs();
3061 std::vector<int> parent =
3062 ComputeGomoryHuTree(num_nodes, symmetrized_relevant_arcs);
3063
3064 // Try all interesting subset from the Gomory-Hu tree.
3065 ExtractAllSubsetsFromForest(parent, &subset_data, &subsets);
3066 num_added +=
3067 TryAllSubsets("CircuitExact", subset_data, subsets, helper, manager);
3068
3069 // Exact separation of symmetric Blossom cut. We use the algorithm in the
3070 // paper: "A Faster Exact Separation Algorithm for Blossom Inequalities", Adam
3071 // N. Letchford, Gerhard Reinelt, Dirk Oliver Theis, 2004.
3072 if (num_added != 0) return;
3073 if (num_nodes <= 2) return;
3074 std::vector<ArcWithLpValue> for_blossom;
3075 for_blossom.reserve(symmetrized_relevant_arcs.size());
3076 for (ArcWithLpValue arc : symmetrized_relevant_arcs) {
3077 if (arc.lp_value > 0.5) arc.lp_value = 1.0 - arc.lp_value;
3078 if (arc.lp_value < 1e-6) continue;
3079 for_blossom.push_back(arc);
3080 }
3081 parent = ComputeGomoryHuTree(num_nodes, for_blossom);
3082 ExtractAllSubsetsFromForest(parent, &subset_data, &subsets);
3083 int last_added_start = -1;
3084 for (const absl::Span<const int> subset : subsets) {
3085 if (subset.size() <= 1) continue;
3086 if (subset.size() == num_nodes) continue;
3087 const int start = static_cast<int>(subset.data() - subset_data.data());
3088 if (start <= last_added_start) continue;
3089 if (helper.TryBlossomSubsetCut("CircuitBlossom", symmetrized_relevant_arcs,
3090 subset, manager)) {
3091 ++num_added;
3092 last_added_start = start;
3093 }
3094 }
3095}
3096
3097namespace {
3098
3099// Returns for each literal its integer view, or the view of its negation.
3100std::vector<IntegerVariable> GetAssociatedVariables(
3101 absl::Span<const Literal> literals, Model* model) {
3102 auto* encoder = model->GetOrCreate<IntegerEncoder>();
3103 std::vector<IntegerVariable> result;
3104 for (const Literal l : literals) {
3105 const IntegerVariable direct_view = encoder->GetLiteralView(l);
3106 if (direct_view != kNoIntegerVariable) {
3107 result.push_back(direct_view);
3108 } else {
3109 result.push_back(encoder->GetLiteralView(l.Negated()));
3110 DCHECK_NE(result.back(), kNoIntegerVariable);
3111 }
3112 }
3113 return result;
3114}
3115
3116} // namespace
3117
3118// We use a basic algorithm to detect components that are not connected to the
3119// rest of the graph in the LP solution, and add cuts to force some arcs to
3120// enter and leave this component from outside.
3122 int num_nodes, absl::Span<const int> tails, absl::Span<const int> heads,
3123 absl::Span<const Literal> literals, Model* model) {
3124 auto helper = std::make_unique<RoutingCutHelper>(
3125 num_nodes, /*is_route_constraint=*/false,
3126 /*flat_node_dim_expressions=*/
3127 absl::Span<const AffineExpression>{}, tails, heads, literals, model);
3128 CutGenerator result;
3129 result.vars = GetAssociatedVariables(literals, model);
3130 result.generate_cuts =
3131 [helper = std::move(helper)](LinearConstraintManager* manager) {
3132 helper->InitializeForNewLpSolution(manager);
3133 SeparateSubtourInequalities(*helper, manager);
3134 return true;
3135 };
3136 return result;
3137}
3138
3140 int num_nodes, absl::Span<const int> tails, absl::Span<const int> heads,
3141 absl::Span<const Literal> literals,
3142 absl::Span<const AffineExpression> flat_node_dim_expressions,
3143 Model* model) {
3144 auto helper = std::make_unique<RoutingCutHelper>(
3145 num_nodes, /*is_route_constraint=*/true, flat_node_dim_expressions, tails,
3146 heads, literals, model);
3147 CutGenerator result;
3148 result.vars = GetAssociatedVariables(literals, model);
3149 result.generate_cuts =
3150 [helper = std::move(helper)](LinearConstraintManager* manager) {
3151 helper->InitializeForNewLpSolution(manager);
3152 SeparateSubtourInequalities(*helper, manager);
3153 helper->TryInfeasiblePathCuts(manager);
3154 return true;
3155 };
3156 return result;
3157}
3158
3159// This is really similar to SeparateSubtourInequalities, see the reference
3160// there.
3162 int num_nodes, absl::Span<const int> tails, absl::Span<const int> heads,
3163 absl::Span<const AffineExpression> arc_capacities,
3164 std::function<void(const std::vector<bool>& in_subset,
3165 IntegerValue* min_incoming_flow,
3166 IntegerValue* min_outgoing_flow)>
3167 get_flows,
3169 LinearConstraintManager* manager, Model* model) {
3170 // We will collect only the arcs with a positive lp capacity value to speed up
3171 // some computation below.
3172 struct Arc {
3173 int tail;
3174 int head;
3175 double lp_value;
3176 IntegerValue offset;
3177 };
3178 std::vector<Arc> relevant_arcs;
3179
3180 // Often capacities have a coeff > 1.
3181 // We currently exploit this if all coeff have a gcd > 1.
3182 int64_t gcd = 0;
3183
3184 // Sort the arcs by non-increasing lp_values.
3185 std::vector<std::pair<double, int>> arc_by_decreasing_lp_values;
3186 for (int i = 0; i < arc_capacities.size(); ++i) {
3187 const double lp_value = arc_capacities[i].LpValue(lp_values);
3188 if (!arc_capacities[i].IsConstant()) {
3189 gcd = std::gcd(gcd, std::abs(arc_capacities[i].coeff.value()));
3190 }
3191 if (lp_value < 1e-6 && arc_capacities[i].constant == 0) continue;
3192 relevant_arcs.push_back(
3193 {tails[i], heads[i], lp_value, arc_capacities[i].constant});
3194 arc_by_decreasing_lp_values.push_back({lp_value, i});
3195 }
3196 if (gcd == 0) return;
3197 std::sort(arc_by_decreasing_lp_values.begin(),
3198 arc_by_decreasing_lp_values.end(),
3199 std::greater<std::pair<double, int>>());
3200
3201 std::vector<std::pair<int, int>> ordered_arcs;
3202 for (const auto& [score, arc] : arc_by_decreasing_lp_values) {
3203 if (tails[arc] == -1) continue;
3204 if (heads[arc] == -1) continue;
3205 ordered_arcs.push_back({tails[arc], heads[arc]});
3206 }
3207 std::vector<int> subset_data;
3208 std::vector<absl::Span<const int>> subsets;
3209 GenerateInterestingSubsets(num_nodes, ordered_arcs,
3210 /*stop_at_num_components=*/1, &subset_data,
3211 &subsets);
3212
3213 // Process each subsets and add any violated cut.
3214 std::vector<bool> in_subset(num_nodes, false);
3215 for (const absl::Span<const int> subset : subsets) {
3216 DCHECK(!subset.empty());
3217 DCHECK_LE(subset.size(), num_nodes);
3218
3219 // Initialize "in_subset" and the subset demands.
3220 for (const int n : subset) in_subset[n] = true;
3221
3222 IntegerValue min_incoming_flow;
3223 IntegerValue min_outgoing_flow;
3224 get_flows(in_subset, &min_incoming_flow, &min_outgoing_flow);
3225
3226 // We will sum the offset of all incoming/outgoing arc capacities.
3227 // Note that all arcs with a non-zero offset are part of relevant_arcs.
3228 IntegerValue incoming_offset(0);
3229 IntegerValue outgoing_offset(0);
3230
3231 // Compute the current flow in and out of the subset.
3232 //
3233 // This can take a significant portion of the running time, it is why it is
3234 // faster to do it only on arcs with non-zero lp values which should be in
3235 // linear number rather than the total number of arc which can be quadratic.
3236 double lp_outgoing_flow = 0.0;
3237 double lp_incoming_flow = 0.0;
3238 for (const auto arc : relevant_arcs) {
3239 if (arc.tail != -1 && in_subset[arc.tail]) {
3240 if (arc.head == -1 || !in_subset[arc.head]) {
3241 incoming_offset += arc.offset;
3242 lp_outgoing_flow += arc.lp_value;
3243 }
3244 } else {
3245 if (arc.head != -1 && in_subset[arc.head]) {
3246 outgoing_offset += arc.offset;
3247 lp_incoming_flow += arc.lp_value;
3248 }
3249 }
3250 }
3251
3252 // If the gcd is greater than one, because all variables are integer we
3253 // can round the flow lower bound to the next multiple of the gcd.
3254 //
3255 // TODO(user): Alternatively, try MIR heuristics if the coefficients in
3256 // the capacities are not all the same.
3257 if (gcd > 1) {
3258 const IntegerValue test_incoming = min_incoming_flow - incoming_offset;
3259 const IntegerValue new_incoming =
3260 CeilRatio(test_incoming, IntegerValue(gcd)) * IntegerValue(gcd);
3261 const IntegerValue incoming_delta = new_incoming - test_incoming;
3262 if (incoming_delta > 0) min_incoming_flow += incoming_delta;
3263 }
3264 if (gcd > 1) {
3265 const IntegerValue test_outgoing = min_outgoing_flow - outgoing_offset;
3266 const IntegerValue new_outgoing =
3267 CeilRatio(test_outgoing, IntegerValue(gcd)) * IntegerValue(gcd);
3268 const IntegerValue outgoing_delta = new_outgoing - test_outgoing;
3269 if (outgoing_delta > 0) min_outgoing_flow += outgoing_delta;
3270 }
3271
3272 if (lp_incoming_flow < ToDouble(min_incoming_flow) - 1e-6) {
3273 VLOG(2) << "INCOMING CUT " << lp_incoming_flow
3274 << " >= " << min_incoming_flow << " size " << subset.size()
3275 << " offset " << incoming_offset << " gcd " << gcd;
3276 LinearConstraintBuilder cut(model, min_incoming_flow, kMaxIntegerValue);
3277 for (int i = 0; i < tails.size(); ++i) {
3278 if ((tails[i] == -1 || !in_subset[tails[i]]) &&
3279 (heads[i] != -1 && in_subset[heads[i]])) {
3280 cut.AddTerm(arc_capacities[i], 1.0);
3281 }
3282 }
3283 manager->AddCut(cut.Build(), "IncomingFlow");
3284 }
3285
3286 if (lp_outgoing_flow < ToDouble(min_outgoing_flow) - 1e-6) {
3287 VLOG(2) << "OUGOING CUT " << lp_outgoing_flow
3288 << " >= " << min_outgoing_flow << " size " << subset.size()
3289 << " offset " << outgoing_offset << " gcd " << gcd;
3290 LinearConstraintBuilder cut(model, min_outgoing_flow, kMaxIntegerValue);
3291 for (int i = 0; i < tails.size(); ++i) {
3292 if ((tails[i] != -1 && in_subset[tails[i]]) &&
3293 (heads[i] == -1 || !in_subset[heads[i]])) {
3294 cut.AddTerm(arc_capacities[i], 1.0);
3295 }
3296 }
3297 manager->AddCut(cut.Build(), "OutgoingFlow");
3298 }
3299
3300 // Sparse clean up.
3301 for (const int n : subset) in_subset[n] = false;
3302 }
3303}
3304
3306 int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
3307 const std::vector<AffineExpression>& arc_capacities,
3308 std::function<void(const std::vector<bool>& in_subset,
3309 IntegerValue* min_incoming_flow,
3310 IntegerValue* min_outgoing_flow)>
3311 get_flows,
3312 Model* model) {
3313 CutGenerator result;
3314 for (const AffineExpression expr : arc_capacities) {
3315 if (!expr.IsConstant()) result.vars.push_back(expr.var);
3316 }
3317 result.generate_cuts = [=](LinearConstraintManager* manager) {
3318 SeparateFlowInequalities(num_nodes, tails, heads, arc_capacities, get_flows,
3319 manager->LpValues(), manager, model);
3320 return true;
3321 };
3322 return result;
3323}
3324
3325} // namespace sat
3326} // namespace operations_research
std::vector< std::vector< T > > FindConnectedComponents()
bool AddEdge(T node1, T node2)
Definition model.h:341
static IntegralType CeilOfRatio(IntegralType numerator, IntegralType denominator)
Definition mathutil.h:39
static IntegralType FloorOfRatio(IntegralType numerator, IntegralType denominator)
Definition mathutil.h:53
ArcIndex AddArcWithCapacity(NodeIndex tail, NodeIndex head, FlowQuantity capacity)
Definition max_flow.cc:27
Status Solve(NodeIndex source, NodeIndex sink)
Definition max_flow.cc:60
void GetSourceSideMinCut(std::vector< NodeIndex > *result)
Definition max_flow.cc:119
Keep the best min outgoing/incoming flow out of a subset.
void Update(int new_bound, std::string name, absl::Span< const int > cannot_be_first={}, absl::Span< const int > cannot_be_last={})
bool PropagateLocalBounds(const IntegerTrail &integer_trail, Literal lit, const absl::flat_hash_map< IntegerVariable, IntegerValue > &input, absl::flat_hash_map< IntegerVariable, IntegerValue > *output) const
void reserve(int size)
Reserve memory if it is already known or tightly estimated.
Definition util.h:89
void AppendToLastVector(const V &value)
Definition util.h:764
int Add(absl::Span< const V > values)
Definition util.h:755
IntegerVariable GetLiteralView(Literal lit) const
Definition integer.h:249
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition integer.h:1419
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Returns globally valid lower/upper bound on the given integer variable.
Definition integer.h:1412
void AddTerm(IntegerVariable var, IntegerValue coeff)
bool AddCut(LinearConstraint ct, std::string type_name, std::string extra_info="")
int ComputeTightMinOutgoingFlow(absl::Span< const int > subset)
int ComputeMinOutgoingFlow(absl::Span< const int > subset)
bool SubsetMightBeServedWithKRoutes(int k, absl::Span< const int > subset, RouteRelationsHelper *helper=nullptr, LinearConstraintManager *manager=nullptr, int special_node=-1, bool use_forward_direction=true)
int ComputeDimensionBasedMinOutgoingFlow(absl::Span< const int > subset, const RouteRelationsHelper &helper, BestBoundHelper *best_bound)
MinOutgoingFlowHelper(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, Model *model)
bool PathExists(int tail, int head) const
const HeadMinusTailBounds & GetArcRelation(int arc, int dimension) const
IntegerValue GetArcOffsetLowerBound(int arc, int dimension, bool negate_expressions) const
const AffineExpression & GetNodeExpression(int node, int dimension) const
Returns the expression associated with the given node and dimension.
static std::unique_ptr< RouteRelationsHelper > Create(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, absl::Span< const AffineExpression > flat_node_dim_expressions, const BinaryRelationRepository &binary_relation_repository, Model *model)
bool PropagateLocalBoundsUsingShortestPaths(const IntegerTrail &integer_trail, int tail, int head, const absl::flat_hash_map< IntegerVariable, IntegerValue > &input, absl::flat_hash_map< IntegerVariable, IntegerValue > *output) const
int num_dimensions() const
Returns the number of "dimensions", such as time or vehicle load.
void RemoveArcs(absl::Span< const int > sorted_arc_indices)
RouteRelationsHelper()=default
Visible for testing.
Simple class to add statistics by name and print them at the end.
int ComputeMinNumberOfBins(absl::Span< ItemOrBin > objects, std::vector< int > &objects_that_cannot_be_bin_and_reach_minimum, std::string &info)
bool UseDpToTightenCapacities(absl::Span< ItemOrBin > objects)
bool GreedyPackingWorks(int num_bins, absl::Span< const ItemOrBin > objects)
const VariablesAssignment & Assignment() const
Definition sat_base.h:462
ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)
Definition graph.h:1419
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition stl_util.h:55
H AbslHashValue(H h, std::shared_ptr< Variable > i)
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
CutGenerator CreateStronglyConnectedGraphCutGenerator(int num_nodes, absl::Span< const int > tails, absl::Span< const int > heads, absl::Span< const Literal > literals, Model *model)
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
std::vector< int > ComputeGomoryHuTree(int num_nodes, absl::Span< const ArcWithLpValue > relevant_arcs)
const LiteralIndex kNoLiteralIndex(-1)
std::pair< int, int > MaybeFillMissingRoutesConstraintNodeExpressions(const CpModelProto &input_model, CpModelProto &output_model)
int TryAllSubsets(std::string cut_name, absl::Span< const int > subset_data, std::vector< absl::Span< const int > > subsets, RoutingCutHelper &helper, LinearConstraintManager *manager)
std::vector< IntegerVariable > NegationOf(absl::Span< const IntegerVariable > vars)
Returns the vector of the negated variables.
Definition integer.cc:52
bool WriteModelProtoToFile(const M &proto, absl::string_view filename)
RoutingCumulExpressions DetectDimensionsAndCumulExpressions(int num_nodes, absl::Span< const int > tails, absl::Span< const int > heads, absl::Span< const Literal > literals, const BinaryRelationRepository &binary_relation_repository)
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue.value())
std::pair< IntegerValue, IntegerValue > GetDifferenceBounds(const NodeExpression &x_expr, const NodeExpression &y_expr, const sat::Relation &r, const std::pair< IntegerValue, IntegerValue > &x_var_bounds, const std::pair< IntegerValue, IntegerValue > &y_var_bounds)
const IntegerVariable kNoIntegerVariable(-1)
IntegerVariable PositiveVariable(IntegerVariable i)
CutGenerator CreateFlowCutGenerator(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< AffineExpression > &arc_capacities, std::function< void(const std::vector< bool > &in_subset, IntegerValue *min_incoming_flow, IntegerValue *min_outgoing_flow)> get_flows, Model *model)
IntegerValue CapAddI(IntegerValue a, IntegerValue b)
void SeparateSubtourInequalities(RoutingCutHelper &helper, LinearConstraintManager *manager)
CutGenerator CreateCVRPCutGenerator(int num_nodes, absl::Span< const int > tails, absl::Span< const int > heads, absl::Span< const Literal > literals, absl::Span< const AffineExpression > flat_node_dim_expressions, Model *model)
void SeparateFlowInequalities(int num_nodes, absl::Span< const int > tails, absl::Span< const int > heads, absl::Span< const AffineExpression > arc_capacities, std::function< void(const std::vector< bool > &in_subset, IntegerValue *min_incoming_flow, IntegerValue *min_outgoing_flow)> get_flows, const util_intops::StrongVector< IntegerVariable, double > &lp_values, LinearConstraintManager *manager, Model *model)
void SymmetrizeArcs(std::vector< ArcWithLpValue > *arcs)
void ExtractAllSubsetsFromForest(absl::Span< const int > parent, std::vector< int > *subset_data, std::vector< absl::Span< const int > > *subsets, int node_limit)
bool VariableIsPositive(IntegerVariable i)
bool AtMinOrMaxInt64I(IntegerValue t)
void GenerateInterestingSubsets(int num_nodes, absl::Span< const std::pair< int, int > > arcs, int stop_at_num_components, std::vector< int > *subset_data, std::vector< absl::Span< const int > > *subsets)
double ToDouble(IntegerValue value)
In SWIG mode, we don't want anything besides these top-level includes.
LinearRange operator==(const LinearExpr &lhs, const LinearExpr &rhs)
ClosedInterval::Iterator end(ClosedInterval interval)
ClosedInterval::Iterator begin(ClosedInterval interval)
STL namespace.
ABSL_FLAG(bool, cp_model_dump_routes_support_graphs, false, "DEBUG ONLY. When set to true, SolveCpModel() dumps the arcs with " "non-zero LP values of the routes constraints, at decision level 0, " "which are used to subsequently generate cuts. The values are " "written as a SupportGraphProto in text format to " "'FLAGS_cp_model_dump_prefix'support_graph_{counter}.pb.txt.")
absl::AnyInvocable< bool(LinearConstraintManager *manager)> generate_cuts
Definition cuts.h:59
std::vector< IntegerVariable > vars
Definition cuts.h:58
std::vector< AffineExpression > flat_node_dim_expressions