1 // Copyright 2010-2021 Google LLC
2 // Licensed under the Apache License, Version 2.0 (the "License");
3 // you may not use this file except in compliance with the License.
4 // You may obtain a copy of the License at
5 //
6 //     http://www.apache.org/licenses/LICENSE-2.0
7 //
8 // Unless required by applicable law or agreed to in writing, software
9 // distributed under the License is distributed on an "AS IS" BASIS,
10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 // See the License for the specific language governing permissions and
12 // limitations under the License.
13 
14 // Implementation of local search filters for routing models.
15 
16 #include "ortools/constraint_solver/routing_filters.h"
17 
18 #include <stddef.h>
19 
20 #include <algorithm>
21 #include <cstdint>
22 #include <deque>
23 #include <functional>
24 #include <iterator>
25 #include <limits>
26 #include <map>
27 #include <memory>
28 #include <numeric>
29 #include <set>
30 #include <string>
31 #include <utility>
32 #include <vector>
33 
34 #include "absl/algorithm/container.h"
35 #include "absl/container/flat_hash_map.h"
36 #include "absl/container/flat_hash_set.h"
37 #include "absl/flags/flag.h"
38 #include "absl/memory/memory.h"
39 #include "absl/strings/string_view.h"
40 #include "ortools/base/int_type.h"
41 #include "ortools/base/integral_types.h"
42 #include "ortools/base/logging.h"
43 #include "ortools/base/map_util.h"
44 #include "ortools/base/small_map.h"
45 #include "ortools/base/small_ordered_set.h"
46 #include "ortools/base/strong_vector.h"
47 #include "ortools/constraint_solver/constraint_solver.h"
48 #include "ortools/constraint_solver/constraint_solveri.h"
49 #include "ortools/constraint_solver/routing.h"
50 #include "ortools/constraint_solver/routing_lp_scheduling.h"
51 #include "ortools/constraint_solver/routing_parameters.pb.h"
52 #include "ortools/graph/min_cost_flow.h"
53 #include "ortools/util/bitset.h"
54 #include "ortools/util/piecewise_linear_function.h"
55 #include "ortools/util/saturated_arithmetic.h"
56 #include "ortools/util/sorted_interval_list.h"
57 
58 ABSL_FLAG(bool, routing_strong_debug_checks, false,
59           "Run stronger checks in debug; these stronger tests might change "
60           "the complexity of the code in particular.");
61 
62 namespace operations_research {
63 
64 namespace {
65 
66 // Max active vehicles filter.
67 
68 class MaxActiveVehiclesFilter : public IntVarLocalSearchFilter {
69  public:
MaxActiveVehiclesFilter(const RoutingModel & routing_model)70   explicit MaxActiveVehiclesFilter(const RoutingModel& routing_model)
71       : IntVarLocalSearchFilter(routing_model.Nexts()),
72         routing_model_(routing_model),
73         is_active_(routing_model.vehicles(), false),
74         active_vehicles_(0) {}
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)75   bool Accept(const Assignment* delta, const Assignment* deltadelta,
76               int64_t objective_min, int64_t objective_max) override {
77     const int64_t kUnassigned = -1;
78     const Assignment::IntContainer& container = delta->IntVarContainer();
79     const int delta_size = container.Size();
80     int current_active_vehicles = active_vehicles_;
81     for (int i = 0; i < delta_size; ++i) {
82       const IntVarElement& new_element = container.Element(i);
83       IntVar* const var = new_element.Var();
84       int64_t index = kUnassigned;
85       if (FindIndex(var, &index) && routing_model_.IsStart(index)) {
86         if (new_element.Min() != new_element.Max()) {
87           // LNS detected.
88           return true;
89         }
90         const int vehicle = routing_model_.VehicleIndex(index);
91         const bool is_active =
92             (new_element.Min() != routing_model_.End(vehicle));
93         if (is_active && !is_active_[vehicle]) {
94           ++current_active_vehicles;
95         } else if (!is_active && is_active_[vehicle]) {
96           --current_active_vehicles;
97         }
98       }
99     }
100     return current_active_vehicles <=
101            routing_model_.GetMaximumNumberOfActiveVehicles();
102   }
103 
104  private:
OnSynchronize(const Assignment * delta)105   void OnSynchronize(const Assignment* delta) override {
106     active_vehicles_ = 0;
107     for (int i = 0; i < routing_model_.vehicles(); ++i) {
108       const int index = routing_model_.Start(i);
109       if (IsVarSynced(index) && Value(index) != routing_model_.End(i)) {
110         is_active_[i] = true;
111         ++active_vehicles_;
112       } else {
113         is_active_[i] = false;
114       }
115     }
116   }
117 
118   const RoutingModel& routing_model_;
119   std::vector<bool> is_active_;
120   int active_vehicles_;
121 };
122 }  // namespace
123 
MakeMaxActiveVehiclesFilter(const RoutingModel & routing_model)124 IntVarLocalSearchFilter* MakeMaxActiveVehiclesFilter(
125     const RoutingModel& routing_model) {
126   return routing_model.solver()->RevAlloc(
127       new MaxActiveVehiclesFilter(routing_model));
128 }
129 
130 namespace {
131 
132 // Node disjunction filter class.
133 class NodeDisjunctionFilter : public IntVarLocalSearchFilter {
134  public:
NodeDisjunctionFilter(const RoutingModel & routing_model,bool filter_cost)135   explicit NodeDisjunctionFilter(const RoutingModel& routing_model,
136                                  bool filter_cost)
137       : IntVarLocalSearchFilter(routing_model.Nexts()),
138         routing_model_(routing_model),
139         active_per_disjunction_(routing_model.GetNumberOfDisjunctions(), 0),
140         inactive_per_disjunction_(routing_model.GetNumberOfDisjunctions(), 0),
141         synchronized_objective_value_(std::numeric_limits<int64_t>::min()),
142         accepted_objective_value_(std::numeric_limits<int64_t>::min()),
143         filter_cost_(filter_cost),
144         has_mandatory_disjunctions_(routing_model.HasMandatoryDisjunctions()) {}
145 
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)146   bool Accept(const Assignment* delta, const Assignment* deltadelta,
147               int64_t objective_min, int64_t objective_max) override {
148     const int64_t kUnassigned = -1;
149     const Assignment::IntContainer& container = delta->IntVarContainer();
150     const int delta_size = container.Size();
151     gtl::small_map<absl::flat_hash_map<RoutingModel::DisjunctionIndex, int>>
152         disjunction_active_deltas;
153     gtl::small_map<absl::flat_hash_map<RoutingModel::DisjunctionIndex, int>>
154         disjunction_inactive_deltas;
155     bool lns_detected = false;
156     // Update active/inactive count per disjunction for each element of delta.
157     for (int i = 0; i < delta_size; ++i) {
158       const IntVarElement& new_element = container.Element(i);
159       IntVar* const var = new_element.Var();
160       int64_t index = kUnassigned;
161       if (FindIndex(var, &index)) {
162         const bool is_inactive =
163             (new_element.Min() <= index && new_element.Max() >= index);
164         if (new_element.Min() != new_element.Max()) {
165           lns_detected = true;
166         }
167         for (const RoutingModel::DisjunctionIndex disjunction_index :
168              routing_model_.GetDisjunctionIndices(index)) {
169           const bool is_var_synced = IsVarSynced(index);
170           if (!is_var_synced || (Value(index) == index) != is_inactive) {
171             ++gtl::LookupOrInsert(is_inactive ? &disjunction_inactive_deltas
172                                               : &disjunction_active_deltas,
173                                   disjunction_index, 0);
174             if (is_var_synced) {
175               --gtl::LookupOrInsert(is_inactive ? &disjunction_active_deltas
176                                                 : &disjunction_inactive_deltas,
177                                     disjunction_index, 0);
178             }
179           }
180         }
181       }
182     }
183     // Check if any disjunction has too many active nodes.
184     for (const auto [disjunction_index, active_nodes] :
185          disjunction_active_deltas) {
186       // Too many active nodes.
187       if (active_per_disjunction_[disjunction_index] + active_nodes >
188           routing_model_.GetDisjunctionMaxCardinality(disjunction_index)) {
189         return false;
190       }
191     }
192     if (lns_detected || (!filter_cost_ && !has_mandatory_disjunctions_)) {
193       accepted_objective_value_ = 0;
194       return true;
195     }
196     // Update penalty costs for disjunctions.
197     accepted_objective_value_ = synchronized_objective_value_;
198     for (const auto [disjunction_index, inactive_nodes] :
199          disjunction_inactive_deltas) {
200       const int64_t penalty =
201           routing_model_.GetDisjunctionPenalty(disjunction_index);
202       if (penalty == 0) continue;
203       const int current_inactive_nodes =
204           inactive_per_disjunction_[disjunction_index];
205       const int max_inactive_cardinality =
206           routing_model_.GetDisjunctionNodeIndices(disjunction_index).size() -
207           routing_model_.GetDisjunctionMaxCardinality(disjunction_index);
208       // Too many inactive nodes.
209       if (current_inactive_nodes + inactive_nodes > max_inactive_cardinality) {
210         if (penalty < 0) {
211           // Nodes are mandatory, i.e. exactly max_cardinality nodes must be
212           // performed, so the move is not acceptable.
213           return false;
214         } else if (current_inactive_nodes <= max_inactive_cardinality) {
215           // Add penalty if there were not too many inactive nodes before the
216           // move.
217           accepted_objective_value_ =
218               CapAdd(accepted_objective_value_, penalty);
219         }
220       } else if (current_inactive_nodes > max_inactive_cardinality) {
221         // Remove penalty if there were too many inactive nodes before the
222         // move and there are not too many after the move.
223         accepted_objective_value_ = CapSub(accepted_objective_value_, penalty);
224       }
225     }
226     // Only compare to max as a cost lower bound is computed.
227     return accepted_objective_value_ <= objective_max;
228   }
DebugString() const229   std::string DebugString() const override { return "NodeDisjunctionFilter"; }
GetSynchronizedObjectiveValue() const230   int64_t GetSynchronizedObjectiveValue() const override {
231     return synchronized_objective_value_;
232   }
GetAcceptedObjectiveValue() const233   int64_t GetAcceptedObjectiveValue() const override {
234     return accepted_objective_value_;
235   }
236 
237  private:
OnSynchronize(const Assignment * delta)238   void OnSynchronize(const Assignment* delta) override {
239     synchronized_objective_value_ = 0;
240     for (RoutingModel::DisjunctionIndex i(0);
241          i < active_per_disjunction_.size(); ++i) {
242       active_per_disjunction_[i] = 0;
243       inactive_per_disjunction_[i] = 0;
244       const std::vector<int64_t>& disjunction_indices =
245           routing_model_.GetDisjunctionNodeIndices(i);
246       for (const int64_t index : disjunction_indices) {
247         if (IsVarSynced(index)) {
248           if (Value(index) != index) {
249             ++active_per_disjunction_[i];
250           } else {
251             ++inactive_per_disjunction_[i];
252           }
253         }
254       }
255       if (!filter_cost_) continue;
256       const int64_t penalty = routing_model_.GetDisjunctionPenalty(i);
257       const int max_cardinality =
258           routing_model_.GetDisjunctionMaxCardinality(i);
259       if (inactive_per_disjunction_[i] >
260               disjunction_indices.size() - max_cardinality &&
261           penalty > 0) {
262         synchronized_objective_value_ =
263             CapAdd(synchronized_objective_value_, penalty);
264       }
265     }
266   }
267 
268   const RoutingModel& routing_model_;
269 
270   absl::StrongVector<RoutingModel::DisjunctionIndex, int>
271       active_per_disjunction_;
272   absl::StrongVector<RoutingModel::DisjunctionIndex, int>
273       inactive_per_disjunction_;
274   int64_t synchronized_objective_value_;
275   int64_t accepted_objective_value_;
276   const bool filter_cost_;
277   const bool has_mandatory_disjunctions_;
278 };
279 }  // namespace
280 
MakeNodeDisjunctionFilter(const RoutingModel & routing_model,bool filter_cost)281 IntVarLocalSearchFilter* MakeNodeDisjunctionFilter(
282     const RoutingModel& routing_model, bool filter_cost) {
283   return routing_model.solver()->RevAlloc(
284       new NodeDisjunctionFilter(routing_model, filter_cost));
285 }
286 
287 const int64_t BasePathFilter::kUnassigned = -1;
288 
BasePathFilter(const std::vector<IntVar * > & nexts,int next_domain_size)289 BasePathFilter::BasePathFilter(const std::vector<IntVar*>& nexts,
290                                int next_domain_size)
291     : IntVarLocalSearchFilter(nexts),
292       node_path_starts_(next_domain_size, kUnassigned),
293       paths_(nexts.size(), -1),
294       new_synchronized_unperformed_nodes_(nexts.size()),
295       new_nexts_(nexts.size(), kUnassigned),
296       touched_paths_(nexts.size()),
297       touched_path_chain_start_ends_(nexts.size(), {kUnassigned, kUnassigned}),
298       ranks_(next_domain_size, -1),
299       status_(BasePathFilter::UNKNOWN) {}
300 
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)301 bool BasePathFilter::Accept(const Assignment* delta,
302                             const Assignment* deltadelta, int64_t objective_min,
303                             int64_t objective_max) {
304   if (IsDisabled()) return true;
305   for (const int touched : delta_touched_) {
306     new_nexts_[touched] = kUnassigned;
307   }
308   delta_touched_.clear();
309   const Assignment::IntContainer& container = delta->IntVarContainer();
310   const int delta_size = container.Size();
311   delta_touched_.reserve(delta_size);
312   // Determining touched paths and their touched chain start and ends (a node is
313   // touched if it corresponds to an element of delta or that an element of
314   // delta points to it).
315   // The start and end of a touched path subchain will have remained on the same
316   // path and will correspond to the min and max ranks of touched nodes in the
317   // current assignment.
318   for (int64_t touched_path : touched_paths_.PositionsSetAtLeastOnce()) {
319     touched_path_chain_start_ends_[touched_path] = {kUnassigned, kUnassigned};
320   }
321   touched_paths_.SparseClearAll();
322 
323   const auto update_touched_path_chain_start_end = [this](int64_t index) {
324     const int64_t start = node_path_starts_[index];
325     if (start == kUnassigned) return;
326     touched_paths_.Set(start);
327 
328     int64_t& chain_start = touched_path_chain_start_ends_[start].first;
329     if (chain_start == kUnassigned || ranks_[index] < ranks_[chain_start]) {
330       chain_start = index;
331     }
332 
333     int64_t& chain_end = touched_path_chain_start_ends_[start].second;
334     if (chain_end == kUnassigned || ranks_[index] > ranks_[chain_end]) {
335       chain_end = index;
336     }
337   };
338 
339   for (int i = 0; i < delta_size; ++i) {
340     const IntVarElement& new_element = container.Element(i);
341     IntVar* const var = new_element.Var();
342     int64_t index = kUnassigned;
343     if (FindIndex(var, &index)) {
344       if (!new_element.Bound()) {
345         // LNS detected
346         return true;
347       }
348       new_nexts_[index] = new_element.Value();
349       delta_touched_.push_back(index);
350       update_touched_path_chain_start_end(index);
351       update_touched_path_chain_start_end(new_nexts_[index]);
352     }
353   }
354   // Checking feasibility of touched paths.
355   if (!InitializeAcceptPath()) return false;
356   for (const int64_t touched_start : touched_paths_.PositionsSetAtLeastOnce()) {
357     const std::pair<int64_t, int64_t> start_end =
358         touched_path_chain_start_ends_[touched_start];
359     if (!AcceptPath(touched_start, start_end.first, start_end.second)) {
360       return false;
361     }
362   }
363   // NOTE: FinalizeAcceptPath() is only called if InitializeAcceptPath() is true
364   // and all paths are accepted.
365   return FinalizeAcceptPath(objective_min, objective_max);
366 }
367 
ComputePathStarts(std::vector<int64_t> * path_starts,std::vector<int> * index_to_path)368 void BasePathFilter::ComputePathStarts(std::vector<int64_t>* path_starts,
369                                        std::vector<int>* index_to_path) {
370   path_starts->clear();
371   const int nexts_size = Size();
372   index_to_path->assign(nexts_size, kUnassigned);
373   Bitset64<> has_prevs(nexts_size);
374   for (int i = 0; i < nexts_size; ++i) {
375     if (!IsVarSynced(i)) {
376       has_prevs.Set(i);
377     } else {
378       const int next = Value(i);
379       if (next < nexts_size) {
380         has_prevs.Set(next);
381       }
382     }
383   }
384   for (int i = 0; i < nexts_size; ++i) {
385     if (!has_prevs[i]) {
386       (*index_to_path)[i] = path_starts->size();
387       path_starts->push_back(i);
388     }
389   }
390 }
391 
HavePathsChanged()392 bool BasePathFilter::HavePathsChanged() {
393   std::vector<int64_t> path_starts;
394   std::vector<int> index_to_path(Size(), kUnassigned);
395   ComputePathStarts(&path_starts, &index_to_path);
396   if (path_starts.size() != starts_.size()) {
397     return true;
398   }
399   for (int i = 0; i < path_starts.size(); ++i) {
400     if (path_starts[i] != starts_[i]) {
401       return true;
402     }
403   }
404   for (int i = 0; i < Size(); ++i) {
405     if (index_to_path[i] != paths_[i]) {
406       return true;
407     }
408   }
409   return false;
410 }
411 
SynchronizeFullAssignment()412 void BasePathFilter::SynchronizeFullAssignment() {
413   // Subclasses of BasePathFilter might not propagate injected objective values
414   // so making sure it is done here (can be done again by the subclass if
415   // needed).
416   ComputePathStarts(&starts_, &paths_);
417   for (int64_t index = 0; index < Size(); index++) {
418     if (IsVarSynced(index) && Value(index) == index &&
419         node_path_starts_[index] != kUnassigned) {
420       // index was performed before and is now unperformed.
421       new_synchronized_unperformed_nodes_.Set(index);
422     }
423   }
424   // Marking unactive nodes (which are not on a path).
425   node_path_starts_.assign(node_path_starts_.size(), kUnassigned);
426   // Marking nodes on a path and storing next values.
427   const int nexts_size = Size();
428   for (const int64_t start : starts_) {
429     int node = start;
430     node_path_starts_[node] = start;
431     DCHECK(IsVarSynced(node));
432     int next = Value(node);
433     while (next < nexts_size) {
434       node = next;
435       node_path_starts_[node] = start;
436       DCHECK(IsVarSynced(node));
437       next = Value(node);
438     }
439     node_path_starts_[next] = start;
440   }
441   OnBeforeSynchronizePaths();
442   UpdateAllRanks();
443   OnAfterSynchronizePaths();
444 }
445 
OnSynchronize(const Assignment * delta)446 void BasePathFilter::OnSynchronize(const Assignment* delta) {
447   if (status_ == BasePathFilter::UNKNOWN) {
448     status_ =
449         DisableFiltering() ? BasePathFilter::DISABLED : BasePathFilter::ENABLED;
450   }
451   if (IsDisabled()) return;
452   new_synchronized_unperformed_nodes_.ClearAll();
453   if (delta == nullptr || delta->Empty() || starts_.empty()) {
454     SynchronizeFullAssignment();
455     return;
456   }
457   // Subclasses of BasePathFilter might not propagate injected objective values
458   // so making sure it is done here (can be done again by the subclass if
459   // needed).
460   // This code supposes that path starts didn't change.
461   DCHECK(!absl::GetFlag(FLAGS_routing_strong_debug_checks) ||
462          !HavePathsChanged());
463   const Assignment::IntContainer& container = delta->IntVarContainer();
464   touched_paths_.SparseClearAll();
465   for (int i = 0; i < container.Size(); ++i) {
466     const IntVarElement& new_element = container.Element(i);
467     int64_t index = kUnassigned;
468     if (FindIndex(new_element.Var(), &index)) {
469       const int64_t start = node_path_starts_[index];
470       if (start != kUnassigned) {
471         touched_paths_.Set(start);
472         if (Value(index) == index) {
473           // New unperformed node (its previous start isn't unassigned).
474           DCHECK_LT(index, new_nexts_.size());
475           new_synchronized_unperformed_nodes_.Set(index);
476           node_path_starts_[index] = kUnassigned;
477         }
478       }
479     }
480   }
481   OnBeforeSynchronizePaths();
482   for (const int64_t touched_start : touched_paths_.PositionsSetAtLeastOnce()) {
483     int64_t node = touched_start;
484     while (node < Size()) {
485       node_path_starts_[node] = touched_start;
486       node = Value(node);
487     }
488     node_path_starts_[node] = touched_start;
489     UpdatePathRanksFromStart(touched_start);
490     OnSynchronizePathFromStart(touched_start);
491   }
492   OnAfterSynchronizePaths();
493 }
494 
UpdateAllRanks()495 void BasePathFilter::UpdateAllRanks() {
496   for (int i = 0; i < ranks_.size(); ++i) {
497     ranks_[i] = kUnassigned;
498   }
499   for (int r = 0; r < NumPaths(); ++r) {
500     UpdatePathRanksFromStart(Start(r));
501     OnSynchronizePathFromStart(Start(r));
502   }
503 }
504 
UpdatePathRanksFromStart(int start)505 void BasePathFilter::UpdatePathRanksFromStart(int start) {
506   int rank = 0;
507   int64_t node = start;
508   while (node < Size()) {
509     ranks_[node] = rank;
510     rank++;
511     node = Value(node);
512   }
513   ranks_[node] = rank;
514 }
515 
516 namespace {
517 
518 class VehicleAmortizedCostFilter : public BasePathFilter {
519  public:
520   explicit VehicleAmortizedCostFilter(const RoutingModel& routing_model);
~VehicleAmortizedCostFilter()521   ~VehicleAmortizedCostFilter() override {}
DebugString() const522   std::string DebugString() const override {
523     return "VehicleAmortizedCostFilter";
524   }
GetSynchronizedObjectiveValue() const525   int64_t GetSynchronizedObjectiveValue() const override {
526     return current_vehicle_cost_;
527   }
GetAcceptedObjectiveValue() const528   int64_t GetAcceptedObjectiveValue() const override {
529     return delta_vehicle_cost_;
530   }
531 
532  private:
533   void OnSynchronizePathFromStart(int64_t start) override;
534   void OnAfterSynchronizePaths() override;
535   bool InitializeAcceptPath() override;
536   bool AcceptPath(int64_t path_start, int64_t chain_start,
537                   int64_t chain_end) override;
538   bool FinalizeAcceptPath(int64_t objective_min,
539                           int64_t objective_max) override;
540 
541   int64_t current_vehicle_cost_;
542   int64_t delta_vehicle_cost_;
543   std::vector<int> current_route_lengths_;
544   std::vector<int64_t> start_to_end_;
545   std::vector<int> start_to_vehicle_;
546   std::vector<int64_t> vehicle_to_start_;
547   const std::vector<int64_t>& linear_cost_factor_of_vehicle_;
548   const std::vector<int64_t>& quadratic_cost_factor_of_vehicle_;
549 };
550 
VehicleAmortizedCostFilter(const RoutingModel & routing_model)551 VehicleAmortizedCostFilter::VehicleAmortizedCostFilter(
552     const RoutingModel& routing_model)
553     : BasePathFilter(routing_model.Nexts(),
554                      routing_model.Size() + routing_model.vehicles()),
555       current_vehicle_cost_(0),
556       delta_vehicle_cost_(0),
557       current_route_lengths_(Size(), -1),
558       linear_cost_factor_of_vehicle_(
559           routing_model.GetAmortizedLinearCostFactorOfVehicles()),
560       quadratic_cost_factor_of_vehicle_(
561           routing_model.GetAmortizedQuadraticCostFactorOfVehicles()) {
562   start_to_end_.resize(Size(), -1);
563   start_to_vehicle_.resize(Size(), -1);
564   vehicle_to_start_.resize(routing_model.vehicles());
565   for (int v = 0; v < routing_model.vehicles(); v++) {
566     const int64_t start = routing_model.Start(v);
567     start_to_vehicle_[start] = v;
568     start_to_end_[start] = routing_model.End(v);
569     vehicle_to_start_[v] = start;
570   }
571 }
572 
OnSynchronizePathFromStart(int64_t start)573 void VehicleAmortizedCostFilter::OnSynchronizePathFromStart(int64_t start) {
574   const int64_t end = start_to_end_[start];
575   CHECK_GE(end, 0);
576   const int route_length = Rank(end) - 1;
577   CHECK_GE(route_length, 0);
578   current_route_lengths_[start] = route_length;
579 }
580 
OnAfterSynchronizePaths()581 void VehicleAmortizedCostFilter::OnAfterSynchronizePaths() {
582   current_vehicle_cost_ = 0;
583   for (int vehicle = 0; vehicle < vehicle_to_start_.size(); vehicle++) {
584     const int64_t start = vehicle_to_start_[vehicle];
585     DCHECK_EQ(vehicle, start_to_vehicle_[start]);
586 
587     const int route_length = current_route_lengths_[start];
588     DCHECK_GE(route_length, 0);
589 
590     if (route_length == 0) {
591       // The path is empty.
592       continue;
593     }
594 
595     const int64_t linear_cost_factor = linear_cost_factor_of_vehicle_[vehicle];
596     const int64_t route_length_cost =
597         CapProd(quadratic_cost_factor_of_vehicle_[vehicle],
598                 route_length * route_length);
599 
600     current_vehicle_cost_ = CapAdd(
601         current_vehicle_cost_, CapSub(linear_cost_factor, route_length_cost));
602   }
603 }
604 
InitializeAcceptPath()605 bool VehicleAmortizedCostFilter::InitializeAcceptPath() {
606   delta_vehicle_cost_ = current_vehicle_cost_;
607   return true;
608 }
609 
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)610 bool VehicleAmortizedCostFilter::AcceptPath(int64_t path_start,
611                                             int64_t chain_start,
612                                             int64_t chain_end) {
613   // Number of nodes previously between chain_start and chain_end
614   const int previous_chain_nodes = Rank(chain_end) - 1 - Rank(chain_start);
615   CHECK_GE(previous_chain_nodes, 0);
616   int new_chain_nodes = 0;
617   int64_t node = GetNext(chain_start);
618   while (node != chain_end) {
619     new_chain_nodes++;
620     node = GetNext(node);
621   }
622 
623   const int previous_route_length = current_route_lengths_[path_start];
624   CHECK_GE(previous_route_length, 0);
625   const int new_route_length =
626       previous_route_length - previous_chain_nodes + new_chain_nodes;
627 
628   const int vehicle = start_to_vehicle_[path_start];
629   CHECK_GE(vehicle, 0);
630   DCHECK_EQ(path_start, vehicle_to_start_[vehicle]);
631 
632   // Update the cost related to used vehicles.
633   // TODO(user): Handle possible overflows.
634   if (previous_route_length == 0) {
635     // The route was empty before, it is no longer the case (changed path).
636     CHECK_GT(new_route_length, 0);
637     delta_vehicle_cost_ =
638         CapAdd(delta_vehicle_cost_, linear_cost_factor_of_vehicle_[vehicle]);
639   } else if (new_route_length == 0) {
640     // The route is now empty.
641     delta_vehicle_cost_ =
642         CapSub(delta_vehicle_cost_, linear_cost_factor_of_vehicle_[vehicle]);
643   }
644 
645   // Update the cost related to the sum of the squares of the route lengths.
646   const int64_t quadratic_cost_factor =
647       quadratic_cost_factor_of_vehicle_[vehicle];
648   delta_vehicle_cost_ =
649       CapAdd(delta_vehicle_cost_,
650              CapProd(quadratic_cost_factor,
651                      previous_route_length * previous_route_length));
652   delta_vehicle_cost_ = CapSub(
653       delta_vehicle_cost_,
654       CapProd(quadratic_cost_factor, new_route_length * new_route_length));
655 
656   return true;
657 }
658 
FinalizeAcceptPath(int64_t objective_min,int64_t objective_max)659 bool VehicleAmortizedCostFilter::FinalizeAcceptPath(int64_t objective_min,
660                                                     int64_t objective_max) {
661   return delta_vehicle_cost_ <= objective_max;
662 }
663 
664 }  // namespace
665 
MakeVehicleAmortizedCostFilter(const RoutingModel & routing_model)666 IntVarLocalSearchFilter* MakeVehicleAmortizedCostFilter(
667     const RoutingModel& routing_model) {
668   return routing_model.solver()->RevAlloc(
669       new VehicleAmortizedCostFilter(routing_model));
670 }
671 
672 namespace {
673 
674 class TypeRegulationsFilter : public BasePathFilter {
675  public:
676   explicit TypeRegulationsFilter(const RoutingModel& model);
~TypeRegulationsFilter()677   ~TypeRegulationsFilter() override {}
DebugString() const678   std::string DebugString() const override { return "TypeRegulationsFilter"; }
679 
680  private:
681   void OnSynchronizePathFromStart(int64_t start) override;
682   bool AcceptPath(int64_t path_start, int64_t chain_start,
683                   int64_t chain_end) override;
684 
685   bool HardIncompatibilitiesRespected(int vehicle, int64_t chain_start,
686                                       int64_t chain_end);
687 
688   const RoutingModel& routing_model_;
689   std::vector<int> start_to_vehicle_;
690   // The following vector is used to keep track of the type counts for hard
691   // incompatibilities.
692   std::vector<std::vector<int>> hard_incompatibility_type_counts_per_vehicle_;
693   // Used to verify the temporal incompatibilities and requirements.
694   TypeIncompatibilityChecker temporal_incompatibility_checker_;
695   TypeRequirementChecker requirement_checker_;
696 };
697 
TypeRegulationsFilter(const RoutingModel & model)698 TypeRegulationsFilter::TypeRegulationsFilter(const RoutingModel& model)
699     : BasePathFilter(model.Nexts(), model.Size() + model.vehicles()),
700       routing_model_(model),
701       start_to_vehicle_(model.Size(), -1),
702       temporal_incompatibility_checker_(model,
703                                         /*check_hard_incompatibilities*/ false),
704       requirement_checker_(model) {
705   const int num_vehicles = model.vehicles();
706   const bool has_hard_type_incompatibilities =
707       model.HasHardTypeIncompatibilities();
708   if (has_hard_type_incompatibilities) {
709     hard_incompatibility_type_counts_per_vehicle_.resize(num_vehicles);
710   }
711   const int num_visit_types = model.GetNumberOfVisitTypes();
712   for (int vehicle = 0; vehicle < num_vehicles; vehicle++) {
713     const int64_t start = model.Start(vehicle);
714     start_to_vehicle_[start] = vehicle;
715     if (has_hard_type_incompatibilities) {
716       hard_incompatibility_type_counts_per_vehicle_[vehicle].resize(
717           num_visit_types, 0);
718     }
719   }
720 }
721 
OnSynchronizePathFromStart(int64_t start)722 void TypeRegulationsFilter::OnSynchronizePathFromStart(int64_t start) {
723   if (!routing_model_.HasHardTypeIncompatibilities()) return;
724 
725   const int vehicle = start_to_vehicle_[start];
726   CHECK_GE(vehicle, 0);
727   std::vector<int>& type_counts =
728       hard_incompatibility_type_counts_per_vehicle_[vehicle];
729   std::fill(type_counts.begin(), type_counts.end(), 0);
730   const int num_types = type_counts.size();
731 
732   int64_t node = start;
733   while (node < Size()) {
734     DCHECK(IsVarSynced(node));
735     const int type = routing_model_.GetVisitType(node);
736     if (type >= 0 && routing_model_.GetVisitTypePolicy(node) !=
737                          RoutingModel::ADDED_TYPE_REMOVED_FROM_VEHICLE) {
738       CHECK_LT(type, num_types);
739       type_counts[type]++;
740     }
741     node = Value(node);
742   }
743 }
744 
HardIncompatibilitiesRespected(int vehicle,int64_t chain_start,int64_t chain_end)745 bool TypeRegulationsFilter::HardIncompatibilitiesRespected(int vehicle,
746                                                            int64_t chain_start,
747                                                            int64_t chain_end) {
748   if (!routing_model_.HasHardTypeIncompatibilities()) return true;
749 
750   const std::vector<int>& previous_type_counts =
751       hard_incompatibility_type_counts_per_vehicle_[vehicle];
752 
753   absl::flat_hash_map</*type*/ int, /*new_count*/ int> new_type_counts;
754   absl::flat_hash_set<int> types_to_check;
755 
756   // Go through the new nodes on the path and increment their type counts.
757   int64_t node = GetNext(chain_start);
758   while (node != chain_end) {
759     const int type = routing_model_.GetVisitType(node);
760     if (type >= 0 && routing_model_.GetVisitTypePolicy(node) !=
761                          RoutingModel::ADDED_TYPE_REMOVED_FROM_VEHICLE) {
762       DCHECK_LT(type, previous_type_counts.size());
763       int& type_count = gtl::LookupOrInsert(&new_type_counts, type,
764                                             previous_type_counts[type]);
765       if (type_count++ == 0) {
766         // New type on the route, mark to check its incompatibilities.
767         types_to_check.insert(type);
768       }
769     }
770     node = GetNext(node);
771   }
772 
773   // Update new_type_counts by decrementing the occurrence of the types of the
774   // nodes no longer on the route.
775   node = Value(chain_start);
776   while (node != chain_end) {
777     const int type = routing_model_.GetVisitType(node);
778     if (type >= 0 && routing_model_.GetVisitTypePolicy(node) !=
779                          RoutingModel::ADDED_TYPE_REMOVED_FROM_VEHICLE) {
780       DCHECK_LT(type, previous_type_counts.size());
781       int& type_count = gtl::LookupOrInsert(&new_type_counts, type,
782                                             previous_type_counts[type]);
783       CHECK_GE(type_count, 1);
784       type_count--;
785     }
786     node = Value(node);
787   }
788 
789   // Check the incompatibilities for types in types_to_check.
790   for (int type : types_to_check) {
791     for (int incompatible_type :
792          routing_model_.GetHardTypeIncompatibilitiesOfType(type)) {
793       if (gtl::FindWithDefault(new_type_counts, incompatible_type,
794                                previous_type_counts[incompatible_type]) > 0) {
795         return false;
796       }
797     }
798   }
799   return true;
800 }
801 
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)802 bool TypeRegulationsFilter::AcceptPath(int64_t path_start, int64_t chain_start,
803                                        int64_t chain_end) {
804   const int vehicle = start_to_vehicle_[path_start];
805   CHECK_GE(vehicle, 0);
806   const auto next_accessor = [this](int64_t node) { return GetNext(node); };
807   return HardIncompatibilitiesRespected(vehicle, chain_start, chain_end) &&
808          temporal_incompatibility_checker_.CheckVehicle(vehicle,
809                                                         next_accessor) &&
810          requirement_checker_.CheckVehicle(vehicle, next_accessor);
811 }
812 
813 }  // namespace
814 
MakeTypeRegulationsFilter(const RoutingModel & routing_model)815 IntVarLocalSearchFilter* MakeTypeRegulationsFilter(
816     const RoutingModel& routing_model) {
817   return routing_model.solver()->RevAlloc(
818       new TypeRegulationsFilter(routing_model));
819 }
820 
821 namespace {
822 
823 // ChainCumul filter. Version of dimension path filter which is O(delta) rather
824 // than O(length of touched paths). Currently only supports dimensions without
825 // costs (global and local span cost, soft bounds) and with unconstrained
826 // cumul variables except overall capacity and cumul variables of path ends.
827 
828 class ChainCumulFilter : public BasePathFilter {
829  public:
830   ChainCumulFilter(const RoutingModel& routing_model,
831                    const RoutingDimension& dimension);
~ChainCumulFilter()832   ~ChainCumulFilter() override {}
DebugString() const833   std::string DebugString() const override {
834     return "ChainCumulFilter(" + name_ + ")";
835   }
836 
837  private:
838   void OnSynchronizePathFromStart(int64_t start) override;
839   bool AcceptPath(int64_t path_start, int64_t chain_start,
840                   int64_t chain_end) override;
841 
842   const std::vector<IntVar*> cumuls_;
843   std::vector<int64_t> start_to_vehicle_;
844   std::vector<int64_t> start_to_end_;
845   std::vector<const RoutingModel::TransitCallback2*> evaluators_;
846   const std::vector<int64_t> vehicle_capacities_;
847   std::vector<int64_t> current_path_cumul_mins_;
848   std::vector<int64_t> current_max_of_path_end_cumul_mins_;
849   std::vector<int64_t> old_nexts_;
850   std::vector<int> old_vehicles_;
851   std::vector<int64_t> current_transits_;
852   const std::string name_;
853 };
854 
ChainCumulFilter(const RoutingModel & routing_model,const RoutingDimension & dimension)855 ChainCumulFilter::ChainCumulFilter(const RoutingModel& routing_model,
856                                    const RoutingDimension& dimension)
857     : BasePathFilter(routing_model.Nexts(), dimension.cumuls().size()),
858       cumuls_(dimension.cumuls()),
859       evaluators_(routing_model.vehicles(), nullptr),
860       vehicle_capacities_(dimension.vehicle_capacities()),
861       current_path_cumul_mins_(dimension.cumuls().size(), 0),
862       current_max_of_path_end_cumul_mins_(dimension.cumuls().size(), 0),
863       old_nexts_(routing_model.Size(), kUnassigned),
864       old_vehicles_(routing_model.Size(), kUnassigned),
865       current_transits_(routing_model.Size(), 0),
866       name_(dimension.name()) {
867   start_to_vehicle_.resize(Size(), -1);
868   start_to_end_.resize(Size(), -1);
869   for (int i = 0; i < routing_model.vehicles(); ++i) {
870     start_to_vehicle_[routing_model.Start(i)] = i;
871     start_to_end_[routing_model.Start(i)] = routing_model.End(i);
872     evaluators_[i] = &dimension.transit_evaluator(i);
873   }
874 }
875 
876 // On synchronization, maintain "propagated" cumul mins and max level of cumul
877 // from each node to the end of the path; to be used by AcceptPath to
878 // incrementally check feasibility.
OnSynchronizePathFromStart(int64_t start)879 void ChainCumulFilter::OnSynchronizePathFromStart(int64_t start) {
880   const int vehicle = start_to_vehicle_[start];
881   std::vector<int64_t> path_nodes;
882   int64_t node = start;
883   int64_t cumul = cumuls_[node]->Min();
884   while (node < Size()) {
885     path_nodes.push_back(node);
886     current_path_cumul_mins_[node] = cumul;
887     const int64_t next = Value(node);
888     if (next != old_nexts_[node] || vehicle != old_vehicles_[node]) {
889       old_nexts_[node] = next;
890       old_vehicles_[node] = vehicle;
891       current_transits_[node] = (*evaluators_[vehicle])(node, next);
892     }
893     cumul = CapAdd(cumul, current_transits_[node]);
894     cumul = std::max(cumuls_[next]->Min(), cumul);
895     node = next;
896   }
897   path_nodes.push_back(node);
898   current_path_cumul_mins_[node] = cumul;
899   int64_t max_cumuls = cumul;
900   for (int i = path_nodes.size() - 1; i >= 0; --i) {
901     const int64_t node = path_nodes[i];
902     max_cumuls = std::max(max_cumuls, current_path_cumul_mins_[node]);
903     current_max_of_path_end_cumul_mins_[node] = max_cumuls;
904   }
905 }
906 
907 // The complexity of the method is O(size of chain (chain_start...chain_end).
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)908 bool ChainCumulFilter::AcceptPath(int64_t path_start, int64_t chain_start,
909                                   int64_t chain_end) {
910   const int vehicle = start_to_vehicle_[path_start];
911   const int64_t capacity = vehicle_capacities_[vehicle];
912   int64_t node = chain_start;
913   int64_t cumul = current_path_cumul_mins_[node];
914   while (node != chain_end) {
915     const int64_t next = GetNext(node);
916     if (IsVarSynced(node) && next == Value(node) &&
917         vehicle == old_vehicles_[node]) {
918       cumul = CapAdd(cumul, current_transits_[node]);
919     } else {
920       cumul = CapAdd(cumul, (*evaluators_[vehicle])(node, next));
921     }
922     cumul = std::max(cumuls_[next]->Min(), cumul);
923     if (cumul > capacity) return false;
924     node = next;
925   }
926   const int64_t end = start_to_end_[path_start];
927   const int64_t end_cumul_delta =
928       CapSub(current_path_cumul_mins_[end], current_path_cumul_mins_[node]);
929   const int64_t after_chain_cumul_delta =
930       CapSub(current_max_of_path_end_cumul_mins_[node],
931              current_path_cumul_mins_[node]);
932   return CapAdd(cumul, after_chain_cumul_delta) <= capacity &&
933          CapAdd(cumul, end_cumul_delta) <= cumuls_[end]->Max();
934 }
935 
936 // PathCumul filter.
937 
938 class PathCumulFilter : public BasePathFilter {
939  public:
940   PathCumulFilter(const RoutingModel& routing_model,
941                   const RoutingDimension& dimension,
942                   const RoutingSearchParameters& parameters,
943                   bool propagate_own_objective_value,
944                   bool filter_objective_cost, bool can_use_lp);
~PathCumulFilter()945   ~PathCumulFilter() override {}
DebugString() const946   std::string DebugString() const override {
947     return "PathCumulFilter(" + name_ + ")";
948   }
GetSynchronizedObjectiveValue() const949   int64_t GetSynchronizedObjectiveValue() const override {
950     return propagate_own_objective_value_ ? synchronized_objective_value_ : 0;
951   }
GetAcceptedObjectiveValue() const952   int64_t GetAcceptedObjectiveValue() const override {
953     return propagate_own_objective_value_ ? accepted_objective_value_ : 0;
954   }
955 
956  private:
957   // This structure stores the "best" path cumul value for a solution, the path
958   // supporting this value, and the corresponding path cumul values for all
959   // paths.
960   struct SupportedPathCumul {
SupportedPathCumuloperations_research::__anonbeef3bbe0711::PathCumulFilter::SupportedPathCumul961     SupportedPathCumul() : cumul_value(0), cumul_value_support(0) {}
962     int64_t cumul_value;
963     int cumul_value_support;
964     std::vector<int64_t> path_values;
965   };
966 
967   struct SoftBound {
SoftBoundoperations_research::__anonbeef3bbe0711::PathCumulFilter::SoftBound968     SoftBound() : bound(-1), coefficient(0) {}
969     int64_t bound;
970     int64_t coefficient;
971   };
972 
973   // This class caches transit values between nodes of paths. Transit and path
974   // nodes are to be added in the order in which they appear on a path.
975   class PathTransits {
976    public:
Clear()977     void Clear() {
978       paths_.clear();
979       transits_.clear();
980     }
ClearPath(int path)981     void ClearPath(int path) {
982       paths_[path].clear();
983       transits_[path].clear();
984     }
AddPaths(int num_paths)985     int AddPaths(int num_paths) {
986       const int first_path = paths_.size();
987       paths_.resize(first_path + num_paths);
988       transits_.resize(first_path + num_paths);
989       return first_path;
990     }
ReserveTransits(int path,int number_of_route_arcs)991     void ReserveTransits(int path, int number_of_route_arcs) {
992       transits_[path].reserve(number_of_route_arcs);
993       paths_[path].reserve(number_of_route_arcs + 1);
994     }
995     // Stores the transit between node and next on path. For a given non-empty
996     // path, node must correspond to next in the previous call to PushTransit.
PushTransit(int path,int node,int next,int64_t transit)997     void PushTransit(int path, int node, int next, int64_t transit) {
998       transits_[path].push_back(transit);
999       if (paths_[path].empty()) {
1000         paths_[path].push_back(node);
1001       }
1002       DCHECK_EQ(paths_[path].back(), node);
1003       paths_[path].push_back(next);
1004     }
NumPaths() const1005     int NumPaths() const { return paths_.size(); }
PathSize(int path) const1006     int PathSize(int path) const { return paths_[path].size(); }
Node(int path,int position) const1007     int Node(int path, int position) const { return paths_[path][position]; }
Transit(int path,int position) const1008     int64_t Transit(int path, int position) const {
1009       return transits_[path][position];
1010     }
1011 
1012    private:
1013     // paths_[r][i] is the ith node on path r.
1014     std::vector<std::vector<int64_t>> paths_;
1015     // transits_[r][i] is the transit value between nodes path_[i] and
1016     // path_[i+1] on path r.
1017     std::vector<std::vector<int64_t>> transits_;
1018   };
1019 
InitializeAcceptPath()1020   bool InitializeAcceptPath() override {
1021     cumul_cost_delta_ = total_current_cumul_cost_value_;
1022     node_with_precedence_to_delta_min_max_cumuls_.clear();
1023     // Cleaning up for the new delta.
1024     delta_max_end_cumul_ = std::numeric_limits<int64_t>::min();
1025     delta_paths_.clear();
1026     delta_path_transits_.Clear();
1027     lns_detected_ = false;
1028     delta_nodes_with_precedences_and_changed_cumul_.ClearAll();
1029     return true;
1030   }
1031   bool AcceptPath(int64_t path_start, int64_t chain_start,
1032                   int64_t chain_end) override;
1033   bool FinalizeAcceptPath(int64_t objective_min,
1034                           int64_t objective_max) override;
1035   void OnBeforeSynchronizePaths() override;
1036 
FilterSpanCost() const1037   bool FilterSpanCost() const { return global_span_cost_coefficient_ != 0; }
1038 
FilterSlackCost() const1039   bool FilterSlackCost() const {
1040     return has_nonzero_vehicle_span_cost_coefficients_ ||
1041            has_vehicle_span_upper_bounds_;
1042   }
1043 
FilterBreakCost(int vehicle) const1044   bool FilterBreakCost(int vehicle) const {
1045     return dimension_.HasBreakConstraints() &&
1046            !dimension_.GetBreakIntervalsOfVehicle(vehicle).empty();
1047   }
1048 
FilterCumulSoftBounds() const1049   bool FilterCumulSoftBounds() const { return !cumul_soft_bounds_.empty(); }
1050 
1051   int64_t GetCumulSoftCost(int64_t node, int64_t cumul_value) const;
1052 
FilterCumulPiecewiseLinearCosts() const1053   bool FilterCumulPiecewiseLinearCosts() const {
1054     return !cumul_piecewise_linear_costs_.empty();
1055   }
1056 
FilterWithDimensionCumulOptimizerForVehicle(int vehicle) const1057   bool FilterWithDimensionCumulOptimizerForVehicle(int vehicle) const {
1058     if (!can_use_lp_ || FilterCumulPiecewiseLinearCosts()) {
1059       return false;
1060     }
1061 
1062     int num_linear_constraints = 0;
1063     if (dimension_.GetSpanCostCoefficientForVehicle(vehicle) > 0)
1064       ++num_linear_constraints;
1065     if (FilterSoftSpanCost(vehicle)) ++num_linear_constraints;
1066     if (FilterCumulSoftLowerBounds()) ++num_linear_constraints;
1067     if (FilterCumulSoftBounds()) ++num_linear_constraints;
1068     if (vehicle_span_upper_bounds_[vehicle] <
1069         std::numeric_limits<int64_t>::max()) {
1070       ++num_linear_constraints;
1071     }
1072     const bool has_breaks = FilterBreakCost(vehicle);
1073     if (has_breaks) ++num_linear_constraints;
1074 
1075     // The DimensionCumulOptimizer is used to compute a more precise value of
1076     // the cost related to the cumul values (soft bounds and span costs).
1077     // It is also used to garantee feasibility with complex mixes of constraints
1078     // and in particular in the presence of break requests along other
1079     // constraints.
1080     // Therefore, without breaks, we only use the optimizer when the costs are
1081     // actually used to filter the solutions, i.e. when filter_objective_cost_
1082     // is true.
1083     return num_linear_constraints >= 2 &&
1084            (has_breaks || filter_objective_cost_);
1085   }
1086 
FilterDimensionForbiddenIntervals() const1087   bool FilterDimensionForbiddenIntervals() const {
1088     for (const SortedDisjointIntervalList& intervals :
1089          dimension_.forbidden_intervals()) {
1090       // TODO(user): Change the following test to check intervals within
1091       // the domain of the corresponding variables.
1092       if (intervals.NumIntervals() > 0) {
1093         return true;
1094       }
1095     }
1096     return false;
1097   }
1098 
1099   int64_t GetCumulPiecewiseLinearCost(int64_t node, int64_t cumul_value) const;
1100 
FilterCumulSoftLowerBounds() const1101   bool FilterCumulSoftLowerBounds() const {
1102     return !cumul_soft_lower_bounds_.empty();
1103   }
1104 
FilterPrecedences() const1105   bool FilterPrecedences() const { return !node_index_to_precedences_.empty(); }
1106 
FilterSoftSpanCost() const1107   bool FilterSoftSpanCost() const {
1108     return dimension_.HasSoftSpanUpperBounds();
1109   }
FilterSoftSpanCost(int vehicle) const1110   bool FilterSoftSpanCost(int vehicle) const {
1111     return dimension_.HasSoftSpanUpperBounds() &&
1112            dimension_.GetSoftSpanUpperBoundForVehicle(vehicle).cost > 0;
1113   }
FilterSoftSpanQuadraticCost() const1114   bool FilterSoftSpanQuadraticCost() const {
1115     return dimension_.HasQuadraticCostSoftSpanUpperBounds();
1116   }
FilterSoftSpanQuadraticCost(int vehicle) const1117   bool FilterSoftSpanQuadraticCost(int vehicle) const {
1118     return dimension_.HasQuadraticCostSoftSpanUpperBounds() &&
1119            dimension_.GetQuadraticCostSoftSpanUpperBoundForVehicle(vehicle)
1120                    .cost > 0;
1121   }
1122 
1123   int64_t GetCumulSoftLowerBoundCost(int64_t node, int64_t cumul_value) const;
1124 
1125   int64_t GetPathCumulSoftLowerBoundCost(const PathTransits& path_transits,
1126                                          int path) const;
1127 
1128   void InitializeSupportedPathCumul(SupportedPathCumul* supported_cumul,
1129                                     int64_t default_value);
1130 
1131   // Given the vector of minimum cumuls on the path, determines if the pickup to
1132   // delivery limits for this dimension (if there are any) can be respected by
1133   // this path.
1134   // Returns true if for every pickup/delivery nodes visited on this path,
1135   // min_cumul_value(delivery) - max_cumul_value(pickup) is less than the limit
1136   // set for this pickup to delivery.
1137   // TODO(user): Verify if we should filter the pickup/delivery limits using
1138   // the LP, for a perfect filtering.
1139   bool PickupToDeliveryLimitsRespected(
1140       const PathTransits& path_transits, int path,
1141       const std::vector<int64_t>& min_path_cumuls) const;
1142 
1143   // Computes the maximum cumul value of nodes along the path using
1144   // [current|delta]_path_transits_, and stores the min/max cumul
1145   // related to each node in the corresponding vector
1146   // [current|delta]_[min|max]_node_cumuls_.
1147   // The boolean is_delta indicates if the computations should take place on the
1148   // "delta" or "current" members. When true, the nodes for which the min/max
1149   // cumul has changed from the current value are marked in
1150   // delta_nodes_with_precedences_and_changed_cumul_.
1151   void StoreMinMaxCumulOfNodesOnPath(
1152       int path, const std::vector<int64_t>& min_path_cumuls, bool is_delta);
1153 
1154   // Compute the max start cumul value for a given path and a given minimal end
1155   // cumul value.
1156   // NOTE: Since this function is used to compute a lower bound on the span of
1157   // the routes, we don't "jump" over the forbidden intervals with this min end
1158   // cumul value. We do however concurrently compute the max possible start
1159   // given the max end cumul, for which we can "jump" over forbidden intervals,
1160   // and return the minimum of the two.
1161   int64_t ComputePathMaxStartFromEndCumul(const PathTransits& path_transits,
1162                                           int path, int64_t path_start,
1163                                           int64_t min_end_cumul) const;
1164 
1165   const RoutingModel& routing_model_;
1166   const RoutingDimension& dimension_;
1167   const std::vector<IntVar*> cumuls_;
1168   const std::vector<IntVar*> slacks_;
1169   std::vector<int64_t> start_to_vehicle_;
1170   std::vector<const RoutingModel::TransitCallback2*> evaluators_;
1171   std::vector<int64_t> vehicle_span_upper_bounds_;
1172   bool has_vehicle_span_upper_bounds_;
1173   int64_t total_current_cumul_cost_value_;
1174   int64_t synchronized_objective_value_;
1175   int64_t accepted_objective_value_;
1176   // Map between paths and path soft cumul bound costs. The paths are indexed
1177   // by the index of the start node of the path.
1178   absl::flat_hash_map<int64_t, int64_t> current_cumul_cost_values_;
1179   int64_t cumul_cost_delta_;
1180   // Cumul cost values for paths in delta, indexed by vehicle.
1181   std::vector<int64_t> delta_path_cumul_cost_values_;
1182   const int64_t global_span_cost_coefficient_;
1183   std::vector<SoftBound> cumul_soft_bounds_;
1184   std::vector<SoftBound> cumul_soft_lower_bounds_;
1185   std::vector<const PiecewiseLinearFunction*> cumul_piecewise_linear_costs_;
1186   std::vector<int64_t> vehicle_span_cost_coefficients_;
1187   bool has_nonzero_vehicle_span_cost_coefficients_;
1188   const std::vector<int64_t> vehicle_capacities_;
1189   // node_index_to_precedences_[node_index] contains all NodePrecedence elements
1190   // with node_index as either "first_node" or "second_node".
1191   // This vector is empty if there are no precedences on the dimension_.
1192   std::vector<std::vector<RoutingDimension::NodePrecedence>>
1193       node_index_to_precedences_;
1194   // Data reflecting information on paths and cumul variables for the solution
1195   // to which the filter was synchronized.
1196   SupportedPathCumul current_min_start_;
1197   SupportedPathCumul current_max_end_;
1198   PathTransits current_path_transits_;
1199   // Current min/max cumul values, indexed by node.
1200   std::vector<std::pair<int64_t, int64_t>> current_min_max_node_cumuls_;
1201   // Data reflecting information on paths and cumul variables for the "delta"
1202   // solution (aka neighbor solution) being examined.
1203   PathTransits delta_path_transits_;
1204   int64_t delta_max_end_cumul_;
1205   SparseBitset<int64_t> delta_nodes_with_precedences_and_changed_cumul_;
1206   absl::flat_hash_map<int64_t, std::pair<int64_t, int64_t>>
1207       node_with_precedence_to_delta_min_max_cumuls_;
1208   // Note: small_ordered_set only support non-hash sets.
1209   gtl::small_ordered_set<std::set<int>> delta_paths_;
1210   const std::string name_;
1211 
1212   LocalDimensionCumulOptimizer* optimizer_;
1213   LocalDimensionCumulOptimizer* mp_optimizer_;
1214   const bool filter_objective_cost_;
1215   // This boolean indicates if the LP optimizer can be used if necessary to
1216   // optimize the dimension cumuls, and is only used for testing purposes.
1217   const bool can_use_lp_;
1218   const bool propagate_own_objective_value_;
1219 
1220   // Used to do span lower bounding in presence of vehicle breaks.
1221   DisjunctivePropagator disjunctive_propagator_;
1222   DisjunctivePropagator::Tasks tasks_;
1223   TravelBounds travel_bounds_;
1224   std::vector<int64_t> current_path_;
1225 
1226   bool lns_detected_;
1227 };
1228 
PathCumulFilter(const RoutingModel & routing_model,const RoutingDimension & dimension,const RoutingSearchParameters & parameters,bool propagate_own_objective_value,bool filter_objective_cost,bool can_use_lp)1229 PathCumulFilter::PathCumulFilter(const RoutingModel& routing_model,
1230                                  const RoutingDimension& dimension,
1231                                  const RoutingSearchParameters& parameters,
1232                                  bool propagate_own_objective_value,
1233                                  bool filter_objective_cost, bool can_use_lp)
1234     : BasePathFilter(routing_model.Nexts(), dimension.cumuls().size()),
1235       routing_model_(routing_model),
1236       dimension_(dimension),
1237       cumuls_(dimension.cumuls()),
1238       slacks_(dimension.slacks()),
1239       evaluators_(routing_model.vehicles(), nullptr),
1240       vehicle_span_upper_bounds_(dimension.vehicle_span_upper_bounds()),
1241       has_vehicle_span_upper_bounds_(false),
1242       total_current_cumul_cost_value_(0),
1243       synchronized_objective_value_(0),
1244       accepted_objective_value_(0),
1245       current_cumul_cost_values_(),
1246       cumul_cost_delta_(0),
1247       delta_path_cumul_cost_values_(routing_model.vehicles(),
1248                                     std::numeric_limits<int64_t>::min()),
1249       global_span_cost_coefficient_(dimension.global_span_cost_coefficient()),
1250       vehicle_span_cost_coefficients_(
1251           dimension.vehicle_span_cost_coefficients()),
1252       has_nonzero_vehicle_span_cost_coefficients_(false),
1253       vehicle_capacities_(dimension.vehicle_capacities()),
1254       delta_max_end_cumul_(0),
1255       delta_nodes_with_precedences_and_changed_cumul_(routing_model.Size()),
1256       name_(dimension.name()),
1257       optimizer_(routing_model.GetMutableLocalCumulOptimizer(dimension)),
1258       mp_optimizer_(routing_model.GetMutableLocalCumulMPOptimizer(dimension)),
1259       filter_objective_cost_(filter_objective_cost),
1260       can_use_lp_(can_use_lp),
1261       propagate_own_objective_value_(propagate_own_objective_value),
1262       lns_detected_(false) {
1263   for (const int64_t upper_bound : vehicle_span_upper_bounds_) {
1264     if (upper_bound != std::numeric_limits<int64_t>::max()) {
1265       has_vehicle_span_upper_bounds_ = true;
1266       break;
1267     }
1268   }
1269   for (const int64_t coefficient : vehicle_span_cost_coefficients_) {
1270     if (coefficient != 0) {
1271       has_nonzero_vehicle_span_cost_coefficients_ = true;
1272       break;
1273     }
1274   }
1275   cumul_soft_bounds_.resize(cumuls_.size());
1276   cumul_soft_lower_bounds_.resize(cumuls_.size());
1277   cumul_piecewise_linear_costs_.resize(cumuls_.size());
1278   bool has_cumul_soft_bounds = false;
1279   bool has_cumul_soft_lower_bounds = false;
1280   bool has_cumul_piecewise_linear_costs = false;
1281   bool has_cumul_hard_bounds = false;
1282   for (const IntVar* const slack : slacks_) {
1283     if (slack->Min() > 0) {
1284       has_cumul_hard_bounds = true;
1285       break;
1286     }
1287   }
1288   for (int i = 0; i < cumuls_.size(); ++i) {
1289     if (dimension.HasCumulVarSoftUpperBound(i)) {
1290       has_cumul_soft_bounds = true;
1291       cumul_soft_bounds_[i].bound = dimension.GetCumulVarSoftUpperBound(i);
1292       cumul_soft_bounds_[i].coefficient =
1293           dimension.GetCumulVarSoftUpperBoundCoefficient(i);
1294     }
1295     if (dimension.HasCumulVarSoftLowerBound(i)) {
1296       has_cumul_soft_lower_bounds = true;
1297       cumul_soft_lower_bounds_[i].bound =
1298           dimension.GetCumulVarSoftLowerBound(i);
1299       cumul_soft_lower_bounds_[i].coefficient =
1300           dimension.GetCumulVarSoftLowerBoundCoefficient(i);
1301     }
1302     if (dimension.HasCumulVarPiecewiseLinearCost(i)) {
1303       has_cumul_piecewise_linear_costs = true;
1304       cumul_piecewise_linear_costs_[i] =
1305           dimension.GetCumulVarPiecewiseLinearCost(i);
1306     }
1307     IntVar* const cumul_var = cumuls_[i];
1308     if (cumul_var->Min() > 0 ||
1309         cumul_var->Max() < std::numeric_limits<int64_t>::max()) {
1310       has_cumul_hard_bounds = true;
1311     }
1312   }
1313   if (!has_cumul_soft_bounds) {
1314     cumul_soft_bounds_.clear();
1315   }
1316   if (!has_cumul_soft_lower_bounds) {
1317     cumul_soft_lower_bounds_.clear();
1318   }
1319   if (!has_cumul_piecewise_linear_costs) {
1320     cumul_piecewise_linear_costs_.clear();
1321   }
1322   if (!has_cumul_hard_bounds) {
1323     // Slacks don't need to be constrained if the cumuls don't have hard bounds;
1324     // therefore we can ignore the vehicle span cost coefficient (note that the
1325     // transit part is already handled by the arc cost filters).
1326     // This doesn't concern the global span filter though.
1327     vehicle_span_cost_coefficients_.assign(routing_model.vehicles(), 0);
1328     has_nonzero_vehicle_span_cost_coefficients_ = false;
1329   }
1330   start_to_vehicle_.resize(Size(), -1);
1331   for (int i = 0; i < routing_model.vehicles(); ++i) {
1332     start_to_vehicle_[routing_model.Start(i)] = i;
1333     evaluators_[i] = &dimension.transit_evaluator(i);
1334   }
1335 
1336   const std::vector<RoutingDimension::NodePrecedence>& node_precedences =
1337       dimension.GetNodePrecedences();
1338   if (!node_precedences.empty()) {
1339     current_min_max_node_cumuls_.resize(cumuls_.size(), {-1, -1});
1340     node_index_to_precedences_.resize(cumuls_.size());
1341     for (const auto& node_precedence : node_precedences) {
1342       node_index_to_precedences_[node_precedence.first_node].push_back(
1343           node_precedence);
1344       node_index_to_precedences_[node_precedence.second_node].push_back(
1345           node_precedence);
1346     }
1347   }
1348 
1349 #ifndef NDEBUG
1350   for (int vehicle = 0; vehicle < routing_model.vehicles(); vehicle++) {
1351     if (FilterWithDimensionCumulOptimizerForVehicle(vehicle)) {
1352       DCHECK_NE(optimizer_, nullptr);
1353       DCHECK_NE(mp_optimizer_, nullptr);
1354     }
1355   }
1356 #endif  // NDEBUG
1357 }
1358 
GetCumulSoftCost(int64_t node,int64_t cumul_value) const1359 int64_t PathCumulFilter::GetCumulSoftCost(int64_t node,
1360                                           int64_t cumul_value) const {
1361   if (node < cumul_soft_bounds_.size()) {
1362     const int64_t bound = cumul_soft_bounds_[node].bound;
1363     const int64_t coefficient = cumul_soft_bounds_[node].coefficient;
1364     if (coefficient > 0 && bound < cumul_value) {
1365       return CapProd(CapSub(cumul_value, bound), coefficient);
1366     }
1367   }
1368   return 0;
1369 }
1370 
GetCumulPiecewiseLinearCost(int64_t node,int64_t cumul_value) const1371 int64_t PathCumulFilter::GetCumulPiecewiseLinearCost(
1372     int64_t node, int64_t cumul_value) const {
1373   if (node < cumul_piecewise_linear_costs_.size()) {
1374     const PiecewiseLinearFunction* cost = cumul_piecewise_linear_costs_[node];
1375     if (cost != nullptr) {
1376       return cost->Value(cumul_value);
1377     }
1378   }
1379   return 0;
1380 }
1381 
GetCumulSoftLowerBoundCost(int64_t node,int64_t cumul_value) const1382 int64_t PathCumulFilter::GetCumulSoftLowerBoundCost(int64_t node,
1383                                                     int64_t cumul_value) const {
1384   if (node < cumul_soft_lower_bounds_.size()) {
1385     const int64_t bound = cumul_soft_lower_bounds_[node].bound;
1386     const int64_t coefficient = cumul_soft_lower_bounds_[node].coefficient;
1387     if (coefficient > 0 && bound > cumul_value) {
1388       return CapProd(CapSub(bound, cumul_value), coefficient);
1389     }
1390   }
1391   return 0;
1392 }
1393 
GetPathCumulSoftLowerBoundCost(const PathTransits & path_transits,int path) const1394 int64_t PathCumulFilter::GetPathCumulSoftLowerBoundCost(
1395     const PathTransits& path_transits, int path) const {
1396   int64_t node = path_transits.Node(path, path_transits.PathSize(path) - 1);
1397   int64_t cumul = cumuls_[node]->Max();
1398   int64_t current_cumul_cost_value = GetCumulSoftLowerBoundCost(node, cumul);
1399   for (int i = path_transits.PathSize(path) - 2; i >= 0; --i) {
1400     node = path_transits.Node(path, i);
1401     cumul = CapSub(cumul, path_transits.Transit(path, i));
1402     cumul = std::min(cumuls_[node]->Max(), cumul);
1403     current_cumul_cost_value = CapAdd(current_cumul_cost_value,
1404                                       GetCumulSoftLowerBoundCost(node, cumul));
1405   }
1406   return current_cumul_cost_value;
1407 }
1408 
OnBeforeSynchronizePaths()1409 void PathCumulFilter::OnBeforeSynchronizePaths() {
1410   total_current_cumul_cost_value_ = 0;
1411   cumul_cost_delta_ = 0;
1412   current_cumul_cost_values_.clear();
1413   if (NumPaths() > 0 &&
1414       (FilterSpanCost() || FilterCumulSoftBounds() || FilterSlackCost() ||
1415        FilterCumulSoftLowerBounds() || FilterCumulPiecewiseLinearCosts() ||
1416        FilterPrecedences() || FilterSoftSpanCost() ||
1417        FilterSoftSpanQuadraticCost())) {
1418     InitializeSupportedPathCumul(&current_min_start_,
1419                                  std::numeric_limits<int64_t>::max());
1420     InitializeSupportedPathCumul(&current_max_end_,
1421                                  std::numeric_limits<int64_t>::min());
1422     current_path_transits_.Clear();
1423     current_path_transits_.AddPaths(NumPaths());
1424     // For each path, compute the minimum end cumul and store the max of these.
1425     for (int r = 0; r < NumPaths(); ++r) {
1426       int64_t node = Start(r);
1427       const int vehicle = start_to_vehicle_[Start(r)];
1428       // First pass: evaluating route length to reserve memory to store route
1429       // information.
1430       int number_of_route_arcs = 0;
1431       while (node < Size()) {
1432         ++number_of_route_arcs;
1433         node = Value(node);
1434       }
1435       current_path_transits_.ReserveTransits(r, number_of_route_arcs);
1436       // Second pass: update cumul, transit and cost values.
1437       node = Start(r);
1438       int64_t cumul = cumuls_[node]->Min();
1439       std::vector<int64_t> min_path_cumuls;
1440       min_path_cumuls.reserve(number_of_route_arcs + 1);
1441       min_path_cumuls.push_back(cumul);
1442 
1443       int64_t current_cumul_cost_value = GetCumulSoftCost(node, cumul);
1444       current_cumul_cost_value = CapAdd(
1445           current_cumul_cost_value, GetCumulPiecewiseLinearCost(node, cumul));
1446 
1447       int64_t total_transit = 0;
1448       while (node < Size()) {
1449         const int64_t next = Value(node);
1450         const int64_t transit = (*evaluators_[vehicle])(node, next);
1451         total_transit = CapAdd(total_transit, transit);
1452         const int64_t transit_slack = CapAdd(transit, slacks_[node]->Min());
1453         current_path_transits_.PushTransit(r, node, next, transit_slack);
1454         cumul = CapAdd(cumul, transit_slack);
1455         cumul =
1456             dimension_.GetFirstPossibleGreaterOrEqualValueForNode(next, cumul);
1457         cumul = std::max(cumuls_[next]->Min(), cumul);
1458         min_path_cumuls.push_back(cumul);
1459         node = next;
1460         current_cumul_cost_value =
1461             CapAdd(current_cumul_cost_value, GetCumulSoftCost(node, cumul));
1462         current_cumul_cost_value = CapAdd(
1463             current_cumul_cost_value, GetCumulPiecewiseLinearCost(node, cumul));
1464       }
1465       if (FilterPrecedences()) {
1466         StoreMinMaxCumulOfNodesOnPath(/*path=*/r, min_path_cumuls,
1467                                       /*is_delta=*/false);
1468       }
1469       if (number_of_route_arcs == 1 &&
1470           !routing_model_.IsVehicleUsedWhenEmpty(vehicle)) {
1471         // This is an empty route (single start->end arc) which we don't take
1472         // into account for costs.
1473         current_cumul_cost_values_[Start(r)] = 0;
1474         current_path_transits_.ClearPath(r);
1475         continue;
1476       }
1477       if (FilterSlackCost() || FilterSoftSpanCost() ||
1478           FilterSoftSpanQuadraticCost()) {
1479         const int64_t start = ComputePathMaxStartFromEndCumul(
1480             current_path_transits_, r, Start(r), cumul);
1481         const int64_t span_lower_bound = CapSub(cumul, start);
1482         if (FilterSlackCost()) {
1483           current_cumul_cost_value =
1484               CapAdd(current_cumul_cost_value,
1485                      CapProd(vehicle_span_cost_coefficients_[vehicle],
1486                              CapSub(span_lower_bound, total_transit)));
1487         }
1488         if (FilterSoftSpanCost()) {
1489           const SimpleBoundCosts::BoundCost bound_cost =
1490               dimension_.GetSoftSpanUpperBoundForVehicle(vehicle);
1491           if (bound_cost.bound < span_lower_bound) {
1492             const int64_t violation =
1493                 CapSub(span_lower_bound, bound_cost.bound);
1494             current_cumul_cost_value = CapAdd(
1495                 current_cumul_cost_value, CapProd(bound_cost.cost, violation));
1496           }
1497         }
1498         if (FilterSoftSpanQuadraticCost()) {
1499           const SimpleBoundCosts::BoundCost bound_cost =
1500               dimension_.GetQuadraticCostSoftSpanUpperBoundForVehicle(vehicle);
1501           if (bound_cost.bound < span_lower_bound) {
1502             const int64_t violation =
1503                 CapSub(span_lower_bound, bound_cost.bound);
1504             current_cumul_cost_value =
1505                 CapAdd(current_cumul_cost_value,
1506                        CapProd(bound_cost.cost, CapProd(violation, violation)));
1507           }
1508         }
1509       }
1510       if (FilterCumulSoftLowerBounds()) {
1511         current_cumul_cost_value =
1512             CapAdd(current_cumul_cost_value,
1513                    GetPathCumulSoftLowerBoundCost(current_path_transits_, r));
1514       }
1515       if (FilterWithDimensionCumulOptimizerForVehicle(vehicle)) {
1516         // TODO(user): Return a status from the optimizer to detect failures
1517         // The only admissible failures here are because of LP timeout.
1518         int64_t lp_cumul_cost_value = 0;
1519         LocalDimensionCumulOptimizer* const optimizer =
1520             FilterBreakCost(vehicle) ? mp_optimizer_ : optimizer_;
1521         DCHECK(optimizer != nullptr);
1522         const DimensionSchedulingStatus status =
1523             optimizer->ComputeRouteCumulCostWithoutFixedTransits(
1524                 vehicle, [this](int64_t node) { return Value(node); },
1525                 &lp_cumul_cost_value);
1526         switch (status) {
1527           case DimensionSchedulingStatus::INFEASIBLE:
1528             lp_cumul_cost_value = 0;
1529             break;
1530           case DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY:
1531             DCHECK(mp_optimizer_ != nullptr);
1532             if (mp_optimizer_->ComputeRouteCumulCostWithoutFixedTransits(
1533                     vehicle, [this](int64_t node) { return Value(node); },
1534                     &lp_cumul_cost_value) ==
1535                 DimensionSchedulingStatus::INFEASIBLE) {
1536               lp_cumul_cost_value = 0;
1537             }
1538             break;
1539           default:
1540             DCHECK(status == DimensionSchedulingStatus::OPTIMAL);
1541         }
1542         current_cumul_cost_value =
1543             std::max(current_cumul_cost_value, lp_cumul_cost_value);
1544       }
1545       current_cumul_cost_values_[Start(r)] = current_cumul_cost_value;
1546       current_max_end_.path_values[r] = cumul;
1547       if (current_max_end_.cumul_value < cumul) {
1548         current_max_end_.cumul_value = cumul;
1549         current_max_end_.cumul_value_support = r;
1550       }
1551       total_current_cumul_cost_value_ =
1552           CapAdd(total_current_cumul_cost_value_, current_cumul_cost_value);
1553     }
1554     if (FilterPrecedences()) {
1555       // Update the min/max node cumuls of new unperformed nodes.
1556       for (int64_t node : GetNewSynchronizedUnperformedNodes()) {
1557         current_min_max_node_cumuls_[node] = {-1, -1};
1558       }
1559     }
1560     // Use the max of the path end cumul mins to compute the corresponding
1561     // maximum start cumul of each path; store the minimum of these.
1562     for (int r = 0; r < NumPaths(); ++r) {
1563       const int64_t start = ComputePathMaxStartFromEndCumul(
1564           current_path_transits_, r, Start(r), current_max_end_.cumul_value);
1565       current_min_start_.path_values[r] = start;
1566       if (current_min_start_.cumul_value > start) {
1567         current_min_start_.cumul_value = start;
1568         current_min_start_.cumul_value_support = r;
1569       }
1570     }
1571   }
1572   // Initialize this before considering any deltas (neighbor).
1573   delta_max_end_cumul_ = std::numeric_limits<int64_t>::min();
1574   lns_detected_ = false;
1575 
1576   DCHECK(global_span_cost_coefficient_ == 0 ||
1577          current_min_start_.cumul_value <= current_max_end_.cumul_value);
1578   synchronized_objective_value_ =
1579       CapAdd(total_current_cumul_cost_value_,
1580              CapProd(global_span_cost_coefficient_,
1581                      CapSub(current_max_end_.cumul_value,
1582                             current_min_start_.cumul_value)));
1583 }
1584 
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)1585 bool PathCumulFilter::AcceptPath(int64_t path_start, int64_t chain_start,
1586                                  int64_t chain_end) {
1587   int64_t node = path_start;
1588   int64_t cumul = cumuls_[node]->Min();
1589   int64_t cumul_cost_delta = 0;
1590   int64_t total_transit = 0;
1591   const int path = delta_path_transits_.AddPaths(1);
1592   const int vehicle = start_to_vehicle_[path_start];
1593   const int64_t capacity = vehicle_capacities_[vehicle];
1594   const bool filter_vehicle_costs =
1595       !routing_model_.IsEnd(GetNext(node)) ||
1596       routing_model_.IsVehicleUsedWhenEmpty(vehicle);
1597   if (filter_vehicle_costs) {
1598     cumul_cost_delta = CapAdd(GetCumulSoftCost(node, cumul),
1599                               GetCumulPiecewiseLinearCost(node, cumul));
1600   }
1601   // Evaluating route length to reserve memory to store transit information.
1602   int number_of_route_arcs = 0;
1603   while (node < Size()) {
1604     const int64_t next = GetNext(node);
1605     // TODO(user): This shouldn't be needed anymore as such deltas should
1606     // have been filtered already.
1607     if (next == kUnassigned) {
1608       // LNS detected, return true since other paths were ok up to now.
1609       lns_detected_ = true;
1610       return true;
1611     }
1612     ++number_of_route_arcs;
1613     node = next;
1614   }
1615   delta_path_transits_.ReserveTransits(path, number_of_route_arcs);
1616   std::vector<int64_t> min_path_cumuls;
1617   min_path_cumuls.reserve(number_of_route_arcs + 1);
1618   min_path_cumuls.push_back(cumul);
1619   // Check that the path is feasible with regards to cumul bounds, scanning
1620   // the paths from start to end (caching path node sequences and transits
1621   // for further span cost filtering).
1622   node = path_start;
1623   while (node < Size()) {
1624     const int64_t next = GetNext(node);
1625     const int64_t transit = (*evaluators_[vehicle])(node, next);
1626     total_transit = CapAdd(total_transit, transit);
1627     const int64_t transit_slack = CapAdd(transit, slacks_[node]->Min());
1628     delta_path_transits_.PushTransit(path, node, next, transit_slack);
1629     cumul = CapAdd(cumul, transit_slack);
1630     cumul = dimension_.GetFirstPossibleGreaterOrEqualValueForNode(next, cumul);
1631     if (cumul > std::min(capacity, cumuls_[next]->Max())) {
1632       return false;
1633     }
1634     cumul = std::max(cumuls_[next]->Min(), cumul);
1635     min_path_cumuls.push_back(cumul);
1636     node = next;
1637     if (filter_vehicle_costs) {
1638       cumul_cost_delta =
1639           CapAdd(cumul_cost_delta, GetCumulSoftCost(node, cumul));
1640       cumul_cost_delta =
1641           CapAdd(cumul_cost_delta, GetCumulPiecewiseLinearCost(node, cumul));
1642     }
1643   }
1644   const int64_t min_end = cumul;
1645 
1646   if (!PickupToDeliveryLimitsRespected(delta_path_transits_, path,
1647                                        min_path_cumuls)) {
1648     return false;
1649   }
1650   if (FilterSlackCost() || FilterBreakCost(vehicle) ||
1651       FilterSoftSpanCost(vehicle) || FilterSoftSpanQuadraticCost(vehicle)) {
1652     int64_t slack_max = std::numeric_limits<int64_t>::max();
1653     if (vehicle_span_upper_bounds_[vehicle] <
1654         std::numeric_limits<int64_t>::max()) {
1655       const int64_t span_max = vehicle_span_upper_bounds_[vehicle];
1656       slack_max = std::min(slack_max, CapSub(span_max, total_transit));
1657     }
1658     const int64_t max_start_from_min_end = ComputePathMaxStartFromEndCumul(
1659         delta_path_transits_, path, path_start, min_end);
1660     const int64_t span_lb = CapSub(min_end, max_start_from_min_end);
1661     int64_t min_total_slack = CapSub(span_lb, total_transit);
1662     if (min_total_slack > slack_max) return false;
1663 
1664     if (dimension_.HasBreakConstraints()) {
1665       for (const auto [limit, min_break_duration] :
1666            dimension_.GetBreakDistanceDurationOfVehicle(vehicle)) {
1667         // Minimal number of breaks depends on total transit:
1668         // 0 breaks for 0 <= total transit <= limit,
1669         // 1 break for limit + 1 <= total transit <= 2 * limit,
1670         // i breaks for i * limit + 1 <= total transit <= (i+1) * limit, ...
1671         if (limit == 0 || total_transit == 0) continue;
1672         const int num_breaks_lb = (total_transit - 1) / limit;
1673         const int64_t slack_lb = CapProd(num_breaks_lb, min_break_duration);
1674         if (slack_lb > slack_max) return false;
1675         min_total_slack = std::max(min_total_slack, slack_lb);
1676       }
1677       // Compute a lower bound of the amount of break that must be made inside
1678       // the route. We compute a mandatory interval (might be empty)
1679       // [max_start, min_end[ during which the route will have to happen,
1680       // then the duration of break that must happen during this interval.
1681       int64_t min_total_break = 0;
1682       int64_t max_path_end = cumuls_[routing_model_.End(vehicle)]->Max();
1683       const int64_t max_start = ComputePathMaxStartFromEndCumul(
1684           delta_path_transits_, path, path_start, max_path_end);
1685       for (const IntervalVar* br :
1686            dimension_.GetBreakIntervalsOfVehicle(vehicle)) {
1687         if (!br->MustBePerformed()) continue;
1688         if (max_start < br->EndMin() && br->StartMax() < min_end) {
1689           min_total_break = CapAdd(min_total_break, br->DurationMin());
1690         }
1691       }
1692       if (min_total_break > slack_max) return false;
1693       min_total_slack = std::max(min_total_slack, min_total_break);
1694     }
1695     if (filter_vehicle_costs) {
1696       cumul_cost_delta = CapAdd(
1697           cumul_cost_delta,
1698           CapProd(vehicle_span_cost_coefficients_[vehicle], min_total_slack));
1699       const int64_t span_lower_bound = CapAdd(total_transit, min_total_slack);
1700       if (FilterSoftSpanCost()) {
1701         const SimpleBoundCosts::BoundCost bound_cost =
1702             dimension_.GetSoftSpanUpperBoundForVehicle(vehicle);
1703         if (bound_cost.bound < span_lower_bound) {
1704           const int64_t violation = CapSub(span_lower_bound, bound_cost.bound);
1705           cumul_cost_delta =
1706               CapAdd(cumul_cost_delta, CapProd(bound_cost.cost, violation));
1707         }
1708       }
1709       if (FilterSoftSpanQuadraticCost()) {
1710         const SimpleBoundCosts::BoundCost bound_cost =
1711             dimension_.GetQuadraticCostSoftSpanUpperBoundForVehicle(vehicle);
1712         if (bound_cost.bound < span_lower_bound) {
1713           const int64_t violation = CapSub(span_lower_bound, bound_cost.bound);
1714           cumul_cost_delta =
1715               CapAdd(cumul_cost_delta,
1716                      CapProd(bound_cost.cost, CapProd(violation, violation)));
1717         }
1718       }
1719     }
1720     if (CapAdd(total_transit, min_total_slack) >
1721         vehicle_span_upper_bounds_[vehicle]) {
1722       return false;
1723     }
1724   }
1725   if (FilterCumulSoftLowerBounds() && filter_vehicle_costs) {
1726     cumul_cost_delta =
1727         CapAdd(cumul_cost_delta,
1728                GetPathCumulSoftLowerBoundCost(delta_path_transits_, path));
1729   }
1730   if (FilterPrecedences()) {
1731     StoreMinMaxCumulOfNodesOnPath(path, min_path_cumuls, /*is_delta=*/true);
1732   }
1733   if (!filter_vehicle_costs) {
1734     // If this route's costs should't be taken into account, reset the
1735     // cumul_cost_delta and delta_path_transits_ for this path.
1736     cumul_cost_delta = 0;
1737     delta_path_transits_.ClearPath(path);
1738   }
1739   if (FilterSpanCost() || FilterCumulSoftBounds() || FilterSlackCost() ||
1740       FilterCumulSoftLowerBounds() || FilterCumulPiecewiseLinearCosts() ||
1741       FilterSoftSpanCost(vehicle) || FilterSoftSpanQuadraticCost(vehicle)) {
1742     delta_paths_.insert(GetPath(path_start));
1743     delta_path_cumul_cost_values_[vehicle] = cumul_cost_delta;
1744     cumul_cost_delta =
1745         CapSub(cumul_cost_delta, current_cumul_cost_values_[path_start]);
1746     if (filter_vehicle_costs) {
1747       delta_max_end_cumul_ = std::max(delta_max_end_cumul_, min_end);
1748     }
1749   }
1750   cumul_cost_delta_ = CapAdd(cumul_cost_delta_, cumul_cost_delta);
1751   return true;
1752 }
1753 
FinalizeAcceptPath(int64_t objective_min,int64_t objective_max)1754 bool PathCumulFilter::FinalizeAcceptPath(int64_t objective_min,
1755                                          int64_t objective_max) {
1756   if ((!FilterSpanCost() && !FilterCumulSoftBounds() && !FilterSlackCost() &&
1757        !FilterCumulSoftLowerBounds() && !FilterCumulPiecewiseLinearCosts() &&
1758        !FilterPrecedences() && !FilterSoftSpanCost() &&
1759        !FilterSoftSpanQuadraticCost()) ||
1760       lns_detected_) {
1761     return true;
1762   }
1763   if (FilterPrecedences()) {
1764     for (int64_t node : delta_nodes_with_precedences_and_changed_cumul_
1765                             .PositionsSetAtLeastOnce()) {
1766       const std::pair<int64_t, int64_t> node_min_max_cumul_in_delta =
1767           gtl::FindWithDefault(node_with_precedence_to_delta_min_max_cumuls_,
1768                                node, {-1, -1});
1769       // NOTE: This node was seen in delta, so its delta min/max cumul should be
1770       // stored in the map.
1771       DCHECK(node_min_max_cumul_in_delta.first >= 0 &&
1772              node_min_max_cumul_in_delta.second >= 0);
1773       for (const RoutingDimension::NodePrecedence& precedence :
1774            node_index_to_precedences_[node]) {
1775         const bool node_is_first = (precedence.first_node == node);
1776         const int64_t other_node =
1777             node_is_first ? precedence.second_node : precedence.first_node;
1778         if (GetNext(other_node) == kUnassigned ||
1779             GetNext(other_node) == other_node) {
1780           // The other node is unperformed, so the precedence constraint is
1781           // inactive.
1782           continue;
1783         }
1784         // max_cumul[second_node] should be greater or equal than
1785         // min_cumul[first_node] + offset.
1786         const std::pair<int64_t, int64_t>& other_min_max_cumul_in_delta =
1787             gtl::FindWithDefault(node_with_precedence_to_delta_min_max_cumuls_,
1788                                  other_node,
1789                                  current_min_max_node_cumuls_[other_node]);
1790 
1791         const int64_t first_min_cumul =
1792             node_is_first ? node_min_max_cumul_in_delta.first
1793                           : other_min_max_cumul_in_delta.first;
1794         const int64_t second_max_cumul =
1795             node_is_first ? other_min_max_cumul_in_delta.second
1796                           : node_min_max_cumul_in_delta.second;
1797 
1798         if (second_max_cumul < first_min_cumul + precedence.offset) {
1799           return false;
1800         }
1801       }
1802     }
1803   }
1804   int64_t new_max_end = delta_max_end_cumul_;
1805   int64_t new_min_start = std::numeric_limits<int64_t>::max();
1806   if (FilterSpanCost()) {
1807     if (new_max_end < current_max_end_.cumul_value) {
1808       // Delta max end is lower than the current solution one.
1809       // If the path supporting the current max end has been modified, we need
1810       // to check all paths to find the largest max end.
1811       if (!gtl::ContainsKey(delta_paths_,
1812                             current_max_end_.cumul_value_support)) {
1813         new_max_end = current_max_end_.cumul_value;
1814       } else {
1815         for (int i = 0; i < current_max_end_.path_values.size(); ++i) {
1816           if (current_max_end_.path_values[i] > new_max_end &&
1817               !gtl::ContainsKey(delta_paths_, i)) {
1818             new_max_end = current_max_end_.path_values[i];
1819           }
1820         }
1821       }
1822     }
1823     // Now that the max end cumul has been found, compute the corresponding
1824     // min start cumul, first from the delta, then if the max end cumul has
1825     // changed, from the unchanged paths as well.
1826     for (int r = 0; r < delta_path_transits_.NumPaths(); ++r) {
1827       new_min_start =
1828           std::min(ComputePathMaxStartFromEndCumul(delta_path_transits_, r,
1829                                                    Start(r), new_max_end),
1830                    new_min_start);
1831     }
1832     if (new_max_end != current_max_end_.cumul_value) {
1833       for (int r = 0; r < NumPaths(); ++r) {
1834         if (gtl::ContainsKey(delta_paths_, r)) {
1835           continue;
1836         }
1837         new_min_start = std::min(new_min_start, ComputePathMaxStartFromEndCumul(
1838                                                     current_path_transits_, r,
1839                                                     Start(r), new_max_end));
1840       }
1841     } else if (new_min_start > current_min_start_.cumul_value) {
1842       // Delta min start is greater than the current solution one.
1843       // If the path supporting the current min start has been modified, we need
1844       // to check all paths to find the smallest min start.
1845       if (!gtl::ContainsKey(delta_paths_,
1846                             current_min_start_.cumul_value_support)) {
1847         new_min_start = current_min_start_.cumul_value;
1848       } else {
1849         for (int i = 0; i < current_min_start_.path_values.size(); ++i) {
1850           if (current_min_start_.path_values[i] < new_min_start &&
1851               !gtl::ContainsKey(delta_paths_, i)) {
1852             new_min_start = current_min_start_.path_values[i];
1853           }
1854         }
1855       }
1856     }
1857   }
1858 
1859   // Filtering on objective value, calling LPs and MIPs if needed..
1860   accepted_objective_value_ =
1861       CapAdd(cumul_cost_delta_, CapProd(global_span_cost_coefficient_,
1862                                         CapSub(new_max_end, new_min_start)));
1863 
1864   if (can_use_lp_ && optimizer_ != nullptr &&
1865       accepted_objective_value_ <= objective_max) {
1866     const size_t num_touched_paths = GetTouchedPathStarts().size();
1867     std::vector<int64_t> path_delta_cost_values(num_touched_paths, 0);
1868     std::vector<bool> requires_mp(num_touched_paths, false);
1869     for (int i = 0; i < num_touched_paths; ++i) {
1870       const int64_t start = GetTouchedPathStarts()[i];
1871       const int vehicle = start_to_vehicle_[start];
1872       if (!FilterWithDimensionCumulOptimizerForVehicle(vehicle)) {
1873         continue;
1874       }
1875       int64_t path_delta_cost_with_lp = 0;
1876       const DimensionSchedulingStatus status =
1877           optimizer_->ComputeRouteCumulCostWithoutFixedTransits(
1878               vehicle, [this](int64_t node) { return GetNext(node); },
1879               &path_delta_cost_with_lp);
1880       if (status == DimensionSchedulingStatus::INFEASIBLE) {
1881         return false;
1882       }
1883       DCHECK(gtl::ContainsKey(delta_paths_, GetPath(start)));
1884       const int64_t path_cost_diff_with_lp = CapSub(
1885           path_delta_cost_with_lp, delta_path_cumul_cost_values_[vehicle]);
1886       if (path_cost_diff_with_lp > 0) {
1887         path_delta_cost_values[i] = path_delta_cost_with_lp;
1888         accepted_objective_value_ =
1889             CapAdd(accepted_objective_value_, path_cost_diff_with_lp);
1890         if (accepted_objective_value_ > objective_max) {
1891           return false;
1892         }
1893       } else {
1894         path_delta_cost_values[i] = delta_path_cumul_cost_values_[vehicle];
1895       }
1896       DCHECK_NE(mp_optimizer_, nullptr);
1897       requires_mp[i] =
1898           FilterBreakCost(vehicle) ||
1899           (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY);
1900     }
1901 
1902     DCHECK_LE(accepted_objective_value_, objective_max);
1903 
1904     for (int i = 0; i < num_touched_paths; ++i) {
1905       if (!requires_mp[i]) {
1906         continue;
1907       }
1908       const int64_t start = GetTouchedPathStarts()[i];
1909       const int vehicle = start_to_vehicle_[start];
1910       int64_t path_delta_cost_with_mp = 0;
1911       if (mp_optimizer_->ComputeRouteCumulCostWithoutFixedTransits(
1912               vehicle, [this](int64_t node) { return GetNext(node); },
1913               &path_delta_cost_with_mp) ==
1914           DimensionSchedulingStatus::INFEASIBLE) {
1915         return false;
1916       }
1917       DCHECK(gtl::ContainsKey(delta_paths_, GetPath(start)));
1918       const int64_t path_cost_diff_with_mp =
1919           CapSub(path_delta_cost_with_mp, path_delta_cost_values[i]);
1920       if (path_cost_diff_with_mp > 0) {
1921         accepted_objective_value_ =
1922             CapAdd(accepted_objective_value_, path_cost_diff_with_mp);
1923         if (accepted_objective_value_ > objective_max) {
1924           return false;
1925         }
1926       }
1927     }
1928   }
1929 
1930   return accepted_objective_value_ <= objective_max;
1931 }
1932 
InitializeSupportedPathCumul(SupportedPathCumul * supported_cumul,int64_t default_value)1933 void PathCumulFilter::InitializeSupportedPathCumul(
1934     SupportedPathCumul* supported_cumul, int64_t default_value) {
1935   supported_cumul->cumul_value = default_value;
1936   supported_cumul->cumul_value_support = -1;
1937   supported_cumul->path_values.resize(NumPaths(), default_value);
1938 }
1939 
PickupToDeliveryLimitsRespected(const PathTransits & path_transits,int path,const std::vector<int64_t> & min_path_cumuls) const1940 bool PathCumulFilter::PickupToDeliveryLimitsRespected(
1941     const PathTransits& path_transits, int path,
1942     const std::vector<int64_t>& min_path_cumuls) const {
1943   if (!dimension_.HasPickupToDeliveryLimits()) {
1944     return true;
1945   }
1946   const int num_pairs = routing_model_.GetPickupAndDeliveryPairs().size();
1947   DCHECK_GT(num_pairs, 0);
1948   std::vector<std::pair<int, int64_t>> visited_delivery_and_min_cumul_per_pair(
1949       num_pairs, {-1, -1});
1950 
1951   const int path_size = path_transits.PathSize(path);
1952   CHECK_EQ(min_path_cumuls.size(), path_size);
1953 
1954   int64_t max_cumul = min_path_cumuls.back();
1955   for (int i = path_transits.PathSize(path) - 2; i >= 0; i--) {
1956     const int node_index = path_transits.Node(path, i);
1957     max_cumul = CapSub(max_cumul, path_transits.Transit(path, i));
1958     max_cumul = std::min(cumuls_[node_index]->Max(), max_cumul);
1959 
1960     const std::vector<std::pair<int, int>>& pickup_index_pairs =
1961         routing_model_.GetPickupIndexPairs(node_index);
1962     const std::vector<std::pair<int, int>>& delivery_index_pairs =
1963         routing_model_.GetDeliveryIndexPairs(node_index);
1964     if (!pickup_index_pairs.empty()) {
1965       // The node is a pickup. Check that it is not a delivery and that it
1966       // appears in a single pickup/delivery pair (as required when limits are
1967       // set on dimension cumuls for pickup and deliveries).
1968       DCHECK(delivery_index_pairs.empty());
1969       DCHECK_EQ(pickup_index_pairs.size(), 1);
1970       const int pair_index = pickup_index_pairs[0].first;
1971       // Get the delivery visited for this pair.
1972       const int delivery_index =
1973           visited_delivery_and_min_cumul_per_pair[pair_index].first;
1974       if (delivery_index < 0) {
1975         // No delivery visited after this pickup for this pickup/delivery pair.
1976         continue;
1977       }
1978       const int64_t cumul_diff_limit =
1979           dimension_.GetPickupToDeliveryLimitForPair(
1980               pair_index, pickup_index_pairs[0].second, delivery_index);
1981       if (CapSub(visited_delivery_and_min_cumul_per_pair[pair_index].second,
1982                  max_cumul) > cumul_diff_limit) {
1983         return false;
1984       }
1985     }
1986     if (!delivery_index_pairs.empty()) {
1987       // The node is a delivery. Check that it's not a pickup and it belongs to
1988       // a single pair.
1989       DCHECK(pickup_index_pairs.empty());
1990       DCHECK_EQ(delivery_index_pairs.size(), 1);
1991       const int pair_index = delivery_index_pairs[0].first;
1992       std::pair<int, int64_t>& delivery_index_and_cumul =
1993           visited_delivery_and_min_cumul_per_pair[pair_index];
1994       int& delivery_index = delivery_index_and_cumul.first;
1995       DCHECK_EQ(delivery_index, -1);
1996       delivery_index = delivery_index_pairs[0].second;
1997       delivery_index_and_cumul.second = min_path_cumuls[i];
1998     }
1999   }
2000   return true;
2001 }
2002 
StoreMinMaxCumulOfNodesOnPath(int path,const std::vector<int64_t> & min_path_cumuls,bool is_delta)2003 void PathCumulFilter::StoreMinMaxCumulOfNodesOnPath(
2004     int path, const std::vector<int64_t>& min_path_cumuls, bool is_delta) {
2005   const PathTransits& path_transits =
2006       is_delta ? delta_path_transits_ : current_path_transits_;
2007 
2008   const int path_size = path_transits.PathSize(path);
2009   DCHECK_EQ(min_path_cumuls.size(), path_size);
2010 
2011   int64_t max_cumul = cumuls_[path_transits.Node(path, path_size - 1)]->Max();
2012   for (int i = path_size - 1; i >= 0; i--) {
2013     const int node_index = path_transits.Node(path, i);
2014 
2015     if (i < path_size - 1) {
2016       max_cumul = CapSub(max_cumul, path_transits.Transit(path, i));
2017       max_cumul = std::min(cumuls_[node_index]->Max(), max_cumul);
2018     }
2019 
2020     if (is_delta && node_index_to_precedences_[node_index].empty()) {
2021       // No need to update the delta cumul map for nodes without precedences.
2022       continue;
2023     }
2024 
2025     std::pair<int64_t, int64_t>& min_max_cumuls =
2026         is_delta ? node_with_precedence_to_delta_min_max_cumuls_[node_index]
2027                  : current_min_max_node_cumuls_[node_index];
2028     min_max_cumuls.first = min_path_cumuls[i];
2029     min_max_cumuls.second = max_cumul;
2030 
2031     if (is_delta && !routing_model_.IsEnd(node_index) &&
2032         (min_max_cumuls.first !=
2033              current_min_max_node_cumuls_[node_index].first ||
2034          max_cumul != current_min_max_node_cumuls_[node_index].second)) {
2035       delta_nodes_with_precedences_and_changed_cumul_.Set(node_index);
2036     }
2037   }
2038 }
2039 
ComputePathMaxStartFromEndCumul(const PathTransits & path_transits,int path,int64_t path_start,int64_t min_end_cumul) const2040 int64_t PathCumulFilter::ComputePathMaxStartFromEndCumul(
2041     const PathTransits& path_transits, int path, int64_t path_start,
2042     int64_t min_end_cumul) const {
2043   int64_t cumul_from_min_end = min_end_cumul;
2044   int64_t cumul_from_max_end =
2045       cumuls_[routing_model_.End(start_to_vehicle_[path_start])]->Max();
2046   for (int i = path_transits.PathSize(path) - 2; i >= 0; --i) {
2047     const int64_t transit = path_transits.Transit(path, i);
2048     const int64_t node = path_transits.Node(path, i);
2049     cumul_from_min_end =
2050         std::min(cumuls_[node]->Max(), CapSub(cumul_from_min_end, transit));
2051     cumul_from_max_end = dimension_.GetLastPossibleLessOrEqualValueForNode(
2052         node, CapSub(cumul_from_max_end, transit));
2053   }
2054   return std::min(cumul_from_min_end, cumul_from_max_end);
2055 }
2056 
2057 }  // namespace
2058 
MakePathCumulFilter(const RoutingDimension & dimension,const RoutingSearchParameters & parameters,bool propagate_own_objective_value,bool filter_objective_cost,bool can_use_lp)2059 IntVarLocalSearchFilter* MakePathCumulFilter(
2060     const RoutingDimension& dimension,
2061     const RoutingSearchParameters& parameters,
2062     bool propagate_own_objective_value, bool filter_objective_cost,
2063     bool can_use_lp) {
2064   RoutingModel& model = *dimension.model();
2065   return model.solver()->RevAlloc(new PathCumulFilter(
2066       model, dimension, parameters, propagate_own_objective_value,
2067       filter_objective_cost, can_use_lp));
2068 }
2069 
2070 namespace {
2071 
DimensionHasCumulCost(const RoutingDimension & dimension)2072 bool DimensionHasCumulCost(const RoutingDimension& dimension) {
2073   if (dimension.global_span_cost_coefficient() != 0) return true;
2074   if (dimension.HasSoftSpanUpperBounds()) return true;
2075   if (dimension.HasQuadraticCostSoftSpanUpperBounds()) return true;
2076   for (const int64_t coefficient : dimension.vehicle_span_cost_coefficients()) {
2077     if (coefficient != 0) return true;
2078   }
2079   for (int i = 0; i < dimension.cumuls().size(); ++i) {
2080     if (dimension.HasCumulVarSoftUpperBound(i)) return true;
2081     if (dimension.HasCumulVarSoftLowerBound(i)) return true;
2082     if (dimension.HasCumulVarPiecewiseLinearCost(i)) return true;
2083   }
2084   return false;
2085 }
2086 
DimensionHasPathCumulConstraint(const RoutingDimension & dimension)2087 bool DimensionHasPathCumulConstraint(const RoutingDimension& dimension) {
2088   if (dimension.HasBreakConstraints()) return true;
2089   if (dimension.HasPickupToDeliveryLimits()) return true;
2090   for (const int64_t upper_bound : dimension.vehicle_span_upper_bounds()) {
2091     if (upper_bound != std::numeric_limits<int64_t>::max()) return true;
2092   }
2093   for (const IntVar* const slack : dimension.slacks()) {
2094     if (slack->Min() > 0) return true;
2095   }
2096   const std::vector<IntVar*>& cumuls = dimension.cumuls();
2097   for (int i = 0; i < cumuls.size(); ++i) {
2098     IntVar* const cumul_var = cumuls[i];
2099     if (cumul_var->Min() > 0 &&
2100         cumul_var->Max() < std::numeric_limits<int64_t>::max() &&
2101         !dimension.model()->IsEnd(i)) {
2102       return true;
2103     }
2104     if (dimension.forbidden_intervals()[i].NumIntervals() > 0) return true;
2105   }
2106   return false;
2107 }
2108 
2109 }  // namespace
2110 
AppendLightWeightDimensionFilters(const PathState * path_state,const std::vector<RoutingDimension * > & dimensions,std::vector<LocalSearchFilterManager::FilterEvent> * filters)2111 void AppendLightWeightDimensionFilters(
2112     const PathState* path_state,
2113     const std::vector<RoutingDimension*>& dimensions,
2114     std::vector<LocalSearchFilterManager::FilterEvent>* filters) {
2115   // For every dimension that fits, add a UnaryDimensionChecker.
2116   for (const RoutingDimension* dimension : dimensions) {
2117     // Skip dimension if not unary.
2118     if (dimension->GetUnaryTransitEvaluator(0) == nullptr) continue;
2119 
2120     using Intervals = std::vector<UnaryDimensionChecker::Interval>;
2121     // Fill path capacities and classes.
2122     const int num_vehicles = dimension->model()->vehicles();
2123     Intervals path_capacity(num_vehicles);
2124     std::vector<int> path_class(num_vehicles);
2125     for (int v = 0; v < num_vehicles; ++v) {
2126       const auto& vehicle_capacities = dimension->vehicle_capacities();
2127       path_capacity[v] = {0, vehicle_capacities[v]};
2128       path_class[v] = dimension->vehicle_to_class(v);
2129     }
2130     // For each class, retrieve the demands of each node.
2131     // Dimension store evaluators with a double indirection for compacity:
2132     // vehicle -> vehicle_class -> evaluator_index.
2133     // We replicate this in UnaryDimensionChecker,
2134     // except we expand evaluator_index to an array of values for all nodes.
2135     const int num_vehicle_classes =
2136         1 + *std::max_element(path_class.begin(), path_class.end());
2137     std::vector<Intervals> demands(num_vehicle_classes);
2138     const int num_cumuls = dimension->cumuls().size();
2139     const int num_slacks = dimension->slacks().size();
2140     for (int vehicle = 0; vehicle < num_vehicles; ++vehicle) {
2141       const int vehicle_class = path_class[vehicle];
2142       if (!demands[vehicle_class].empty()) continue;
2143       const auto& evaluator = dimension->GetUnaryTransitEvaluator(vehicle);
2144       Intervals class_demands(num_cumuls);
2145       for (int node = 0; node < num_cumuls; ++node) {
2146         if (node < num_slacks) {
2147           const int64_t demand_min = evaluator(node);
2148           const int64_t slack_max = dimension->SlackVar(node)->Max();
2149           class_demands[node] = {demand_min, CapAdd(demand_min, slack_max)};
2150         } else {
2151           class_demands[node] = {0, 0};
2152         }
2153       }
2154       demands[vehicle_class] = std::move(class_demands);
2155     }
2156     // Fill node capacities.
2157     Intervals node_capacity(num_cumuls);
2158     for (int node = 0; node < num_cumuls; ++node) {
2159       const IntVar* cumul = dimension->CumulVar(node);
2160       node_capacity[node] = {cumul->Min(), cumul->Max()};
2161     }
2162     // Make the dimension checker and pass ownership to the filter.
2163     auto checker = absl::make_unique<UnaryDimensionChecker>(
2164         path_state, std::move(path_capacity), std::move(path_class),
2165         std::move(demands), std::move(node_capacity));
2166     const auto kAccept = LocalSearchFilterManager::FilterEventType::kAccept;
2167     LocalSearchFilter* filter = MakeUnaryDimensionFilter(
2168         dimension->model()->solver(), std::move(checker), dimension->name());
2169     filters->push_back({filter, kAccept});
2170   }
2171 }
2172 
AppendDimensionCumulFilters(const std::vector<RoutingDimension * > & dimensions,const RoutingSearchParameters & parameters,bool filter_objective_cost,bool filter_light_weight_unary_dimensions,std::vector<LocalSearchFilterManager::FilterEvent> * filters)2173 void AppendDimensionCumulFilters(
2174     const std::vector<RoutingDimension*>& dimensions,
2175     const RoutingSearchParameters& parameters, bool filter_objective_cost,
2176     bool filter_light_weight_unary_dimensions,
2177     std::vector<LocalSearchFilterManager::FilterEvent>* filters) {
2178   const auto kAccept = LocalSearchFilterManager::FilterEventType::kAccept;
2179   // NOTE: We first sort the dimensions by increasing complexity of filtering:
2180   // - Dimensions without any cumul-related costs or constraints will have a
2181   //   ChainCumulFilter.
2182   // - Dimensions with cumul costs or constraints, but no global span cost
2183   //   and/or precedences will have a PathCumulFilter.
2184   // - Dimensions with a global span cost coefficient and/or precedences will
2185   //   have a global LP filter.
2186   const int num_dimensions = dimensions.size();
2187 
2188   std::vector<bool> use_path_cumul_filter(num_dimensions);
2189   std::vector<bool> use_cumul_bounds_propagator_filter(num_dimensions);
2190   std::vector<bool> use_global_lp_filter(num_dimensions);
2191   std::vector<bool> use_resource_assignment_filter(num_dimensions);
2192   std::vector<int> filtering_difficulty(num_dimensions);
2193   for (int d = 0; d < num_dimensions; d++) {
2194     const RoutingDimension& dimension = *dimensions[d];
2195     const bool has_cumul_cost = DimensionHasCumulCost(dimension);
2196     use_path_cumul_filter[d] =
2197         has_cumul_cost || DimensionHasPathCumulConstraint(dimension);
2198 
2199     const int num_dimension_resource_groups =
2200         dimension.model()->GetDimensionResourceGroupIndices(&dimension).size();
2201     const bool can_use_cumul_bounds_propagator_filter =
2202         !dimension.HasBreakConstraints() &&
2203         num_dimension_resource_groups == 0 &&
2204         (!filter_objective_cost || !has_cumul_cost);
2205     const bool has_precedences = !dimension.GetNodePrecedences().empty();
2206     use_global_lp_filter[d] =
2207         (has_precedences && !can_use_cumul_bounds_propagator_filter) ||
2208         (filter_objective_cost &&
2209          dimension.global_span_cost_coefficient() > 0) ||
2210         num_dimension_resource_groups > 1;
2211 
2212     use_cumul_bounds_propagator_filter[d] =
2213         has_precedences && !use_global_lp_filter[d];
2214 
2215     use_resource_assignment_filter[d] = num_dimension_resource_groups > 0;
2216 
2217     filtering_difficulty[d] =
2218         8 * use_global_lp_filter[d] + 4 * use_resource_assignment_filter[d] +
2219         2 * use_cumul_bounds_propagator_filter[d] + use_path_cumul_filter[d];
2220   }
2221 
2222   std::vector<int> sorted_dimension_indices(num_dimensions);
2223   std::iota(sorted_dimension_indices.begin(), sorted_dimension_indices.end(),
2224             0);
2225   std::sort(sorted_dimension_indices.begin(), sorted_dimension_indices.end(),
2226             [&filtering_difficulty](int d1, int d2) {
2227               return filtering_difficulty[d1] < filtering_difficulty[d2];
2228             });
2229 
2230   for (const int d : sorted_dimension_indices) {
2231     const RoutingDimension& dimension = *dimensions[d];
2232     const RoutingModel& model = *dimension.model();
2233     // NOTE: We always add the [Chain|Path]CumulFilter to filter each route's
2234     // feasibility separately to try and cut bad decisions earlier in the
2235     // search, but we don't propagate the computed cost if the LPCumulFilter is
2236     // already doing it.
2237     const bool use_global_lp = use_global_lp_filter[d];
2238     const bool filter_resource_assignment = use_resource_assignment_filter[d];
2239     if (use_path_cumul_filter[d]) {
2240       filters->push_back(
2241           {MakePathCumulFilter(dimension, parameters,
2242                                /*propagate_own_objective_value*/
2243                                !use_global_lp && !filter_resource_assignment,
2244                                filter_objective_cost),
2245            kAccept});
2246     } else if (filter_light_weight_unary_dimensions ||
2247                dimension.GetUnaryTransitEvaluator(0) == nullptr) {
2248       filters->push_back(
2249           {model.solver()->RevAlloc(new ChainCumulFilter(model, dimension)),
2250            kAccept});
2251     }
2252 
2253     if (use_cumul_bounds_propagator_filter[d]) {
2254       DCHECK(!use_global_lp);
2255       DCHECK(!filter_resource_assignment);
2256       filters->push_back({MakeCumulBoundsPropagatorFilter(dimension), kAccept});
2257     }
2258 
2259     if (filter_resource_assignment) {
2260       filters->push_back({MakeResourceAssignmentFilter(
2261           model.GetMutableLocalCumulOptimizer(dimension),
2262           model.GetMutableLocalCumulMPOptimizer(dimension),
2263           /*propagate_own_objective_value*/ !use_global_lp,
2264           filter_objective_cost)});
2265     }
2266 
2267     if (use_global_lp) {
2268       DCHECK(model.GetMutableGlobalCumulOptimizer(dimension) != nullptr);
2269       filters->push_back({MakeGlobalLPCumulFilter(
2270                               model.GetMutableGlobalCumulOptimizer(dimension),
2271                               model.GetMutableGlobalCumulMPOptimizer(dimension),
2272                               filter_objective_cost),
2273                           kAccept});
2274     }
2275   }
2276 }
2277 
2278 namespace {
2279 
2280 // Filter for pickup/delivery precedences.
2281 class PickupDeliveryFilter : public BasePathFilter {
2282  public:
2283   PickupDeliveryFilter(const std::vector<IntVar*>& nexts, int next_domain_size,
2284                        const RoutingModel::IndexPairs& pairs,
2285                        const std::vector<RoutingModel::PickupAndDeliveryPolicy>&
2286                            vehicle_policies);
~PickupDeliveryFilter()2287   ~PickupDeliveryFilter() override {}
2288   bool AcceptPath(int64_t path_start, int64_t chain_start,
2289                   int64_t chain_end) override;
DebugString() const2290   std::string DebugString() const override { return "PickupDeliveryFilter"; }
2291 
2292  private:
2293   bool AcceptPathDefault(int64_t path_start);
2294   template <bool lifo>
2295   bool AcceptPathOrdered(int64_t path_start);
2296 
2297   std::vector<int> pair_firsts_;
2298   std::vector<int> pair_seconds_;
2299   const RoutingModel::IndexPairs pairs_;
2300   SparseBitset<> visited_;
2301   std::deque<int> visited_deque_;
2302   const std::vector<RoutingModel::PickupAndDeliveryPolicy> vehicle_policies_;
2303 };
2304 
PickupDeliveryFilter(const std::vector<IntVar * > & nexts,int next_domain_size,const RoutingModel::IndexPairs & pairs,const std::vector<RoutingModel::PickupAndDeliveryPolicy> & vehicle_policies)2305 PickupDeliveryFilter::PickupDeliveryFilter(
2306     const std::vector<IntVar*>& nexts, int next_domain_size,
2307     const RoutingModel::IndexPairs& pairs,
2308     const std::vector<RoutingModel::PickupAndDeliveryPolicy>& vehicle_policies)
2309     : BasePathFilter(nexts, next_domain_size),
2310       pair_firsts_(next_domain_size, kUnassigned),
2311       pair_seconds_(next_domain_size, kUnassigned),
2312       pairs_(pairs),
2313       visited_(Size()),
2314       vehicle_policies_(vehicle_policies) {
2315   for (int i = 0; i < pairs.size(); ++i) {
2316     const auto& index_pair = pairs[i];
2317     for (int first : index_pair.first) {
2318       pair_firsts_[first] = i;
2319     }
2320     for (int second : index_pair.second) {
2321       pair_seconds_[second] = i;
2322     }
2323   }
2324 }
2325 
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)2326 bool PickupDeliveryFilter::AcceptPath(int64_t path_start, int64_t chain_start,
2327                                       int64_t chain_end) {
2328   switch (vehicle_policies_[GetPath(path_start)]) {
2329     case RoutingModel::PICKUP_AND_DELIVERY_NO_ORDER:
2330       return AcceptPathDefault(path_start);
2331     case RoutingModel::PICKUP_AND_DELIVERY_LIFO:
2332       return AcceptPathOrdered<true>(path_start);
2333     case RoutingModel::PICKUP_AND_DELIVERY_FIFO:
2334       return AcceptPathOrdered<false>(path_start);
2335     default:
2336       return true;
2337   }
2338 }
2339 
AcceptPathDefault(int64_t path_start)2340 bool PickupDeliveryFilter::AcceptPathDefault(int64_t path_start) {
2341   visited_.ClearAll();
2342   int64_t node = path_start;
2343   int64_t path_length = 1;
2344   while (node < Size()) {
2345     // Detect sub-cycles (path is longer than longest possible path).
2346     if (path_length > Size()) {
2347       return false;
2348     }
2349     if (pair_firsts_[node] != kUnassigned) {
2350       // Checking on pair firsts is not actually necessary (inconsistencies
2351       // will get caught when checking pair seconds); doing it anyway to
2352       // cut checks early.
2353       for (int second : pairs_[pair_firsts_[node]].second) {
2354         if (visited_[second]) {
2355           return false;
2356         }
2357       }
2358     }
2359     if (pair_seconds_[node] != kUnassigned) {
2360       bool found_first = false;
2361       bool some_synced = false;
2362       for (int first : pairs_[pair_seconds_[node]].first) {
2363         if (visited_[first]) {
2364           found_first = true;
2365           break;
2366         }
2367         if (IsVarSynced(first)) {
2368           some_synced = true;
2369         }
2370       }
2371       if (!found_first && some_synced) {
2372         return false;
2373       }
2374     }
2375     visited_.Set(node);
2376     const int64_t next = GetNext(node);
2377     if (next == kUnassigned) {
2378       // LNS detected, return true since path was ok up to now.
2379       return true;
2380     }
2381     node = next;
2382     ++path_length;
2383   }
2384   for (const int64_t node : visited_.PositionsSetAtLeastOnce()) {
2385     if (pair_firsts_[node] != kUnassigned) {
2386       bool found_second = false;
2387       bool some_synced = false;
2388       for (int second : pairs_[pair_firsts_[node]].second) {
2389         if (visited_[second]) {
2390           found_second = true;
2391           break;
2392         }
2393         if (IsVarSynced(second)) {
2394           some_synced = true;
2395         }
2396       }
2397       if (!found_second && some_synced) {
2398         return false;
2399       }
2400     }
2401   }
2402   return true;
2403 }
2404 
2405 template <bool lifo>
AcceptPathOrdered(int64_t path_start)2406 bool PickupDeliveryFilter::AcceptPathOrdered(int64_t path_start) {
2407   visited_deque_.clear();
2408   int64_t node = path_start;
2409   int64_t path_length = 1;
2410   while (node < Size()) {
2411     // Detect sub-cycles (path is longer than longest possible path).
2412     if (path_length > Size()) {
2413       return false;
2414     }
2415     if (pair_firsts_[node] != kUnassigned) {
2416       if (lifo) {
2417         visited_deque_.push_back(node);
2418       } else {
2419         visited_deque_.push_front(node);
2420       }
2421     }
2422     if (pair_seconds_[node] != kUnassigned) {
2423       bool found_first = false;
2424       bool some_synced = false;
2425       for (int first : pairs_[pair_seconds_[node]].first) {
2426         if (!visited_deque_.empty() && visited_deque_.back() == first) {
2427           found_first = true;
2428           break;
2429         }
2430         if (IsVarSynced(first)) {
2431           some_synced = true;
2432         }
2433       }
2434       if (!found_first && some_synced) {
2435         return false;
2436       } else if (!visited_deque_.empty()) {
2437         visited_deque_.pop_back();
2438       }
2439     }
2440     const int64_t next = GetNext(node);
2441     if (next == kUnassigned) {
2442       // LNS detected, return true since path was ok up to now.
2443       return true;
2444     }
2445     node = next;
2446     ++path_length;
2447   }
2448   while (!visited_deque_.empty()) {
2449     for (int second : pairs_[pair_firsts_[visited_deque_.back()]].second) {
2450       if (IsVarSynced(second)) {
2451         return false;
2452       }
2453     }
2454     visited_deque_.pop_back();
2455   }
2456   return true;
2457 }
2458 
2459 }  // namespace
2460 
MakePickupDeliveryFilter(const RoutingModel & routing_model,const RoutingModel::IndexPairs & pairs,const std::vector<RoutingModel::PickupAndDeliveryPolicy> & vehicle_policies)2461 IntVarLocalSearchFilter* MakePickupDeliveryFilter(
2462     const RoutingModel& routing_model, const RoutingModel::IndexPairs& pairs,
2463     const std::vector<RoutingModel::PickupAndDeliveryPolicy>&
2464         vehicle_policies) {
2465   return routing_model.solver()->RevAlloc(new PickupDeliveryFilter(
2466       routing_model.Nexts(), routing_model.Size() + routing_model.vehicles(),
2467       pairs, vehicle_policies));
2468 }
2469 
2470 namespace {
2471 
2472 // Vehicle variable filter
2473 class VehicleVarFilter : public BasePathFilter {
2474  public:
2475   explicit VehicleVarFilter(const RoutingModel& routing_model);
~VehicleVarFilter()2476   ~VehicleVarFilter() override {}
2477   bool AcceptPath(int64_t path_start, int64_t chain_start,
2478                   int64_t chain_end) override;
DebugString() const2479   std::string DebugString() const override { return "VehicleVariableFilter"; }
2480 
2481  private:
2482   bool DisableFiltering() const override;
2483   bool IsVehicleVariableConstrained(int index) const;
2484 
2485   std::vector<int64_t> start_to_vehicle_;
2486   std::vector<IntVar*> vehicle_vars_;
2487   const int64_t unconstrained_vehicle_var_domain_size_;
2488 };
2489 
VehicleVarFilter(const RoutingModel & routing_model)2490 VehicleVarFilter::VehicleVarFilter(const RoutingModel& routing_model)
2491     : BasePathFilter(routing_model.Nexts(),
2492                      routing_model.Size() + routing_model.vehicles()),
2493       vehicle_vars_(routing_model.VehicleVars()),
2494       unconstrained_vehicle_var_domain_size_(routing_model.vehicles()) {
2495   start_to_vehicle_.resize(Size(), -1);
2496   for (int i = 0; i < routing_model.vehicles(); ++i) {
2497     start_to_vehicle_[routing_model.Start(i)] = i;
2498   }
2499 }
2500 
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)2501 bool VehicleVarFilter::AcceptPath(int64_t path_start, int64_t chain_start,
2502                                   int64_t chain_end) {
2503   const int64_t vehicle = start_to_vehicle_[path_start];
2504   int64_t node = chain_start;
2505   while (node != chain_end) {
2506     if (!vehicle_vars_[node]->Contains(vehicle)) {
2507       return false;
2508     }
2509     node = GetNext(node);
2510   }
2511   return vehicle_vars_[node]->Contains(vehicle);
2512 }
2513 
DisableFiltering() const2514 bool VehicleVarFilter::DisableFiltering() const {
2515   for (int i = 0; i < vehicle_vars_.size(); ++i) {
2516     if (IsVehicleVariableConstrained(i)) return false;
2517   }
2518   return true;
2519 }
2520 
IsVehicleVariableConstrained(int index) const2521 bool VehicleVarFilter::IsVehicleVariableConstrained(int index) const {
2522   const IntVar* const vehicle_var = vehicle_vars_[index];
2523   // If vehicle variable contains -1 (optional node), then we need to
2524   // add it to the "unconstrained" domain. Impact we don't filter mandatory
2525   // nodes made inactive here, but it is covered by other filters.
2526   const int adjusted_unconstrained_vehicle_var_domain_size =
2527       vehicle_var->Min() >= 0 ? unconstrained_vehicle_var_domain_size_
2528                               : unconstrained_vehicle_var_domain_size_ + 1;
2529   return vehicle_var->Size() != adjusted_unconstrained_vehicle_var_domain_size;
2530 }
2531 
2532 }  // namespace
2533 
MakeVehicleVarFilter(const RoutingModel & routing_model)2534 IntVarLocalSearchFilter* MakeVehicleVarFilter(
2535     const RoutingModel& routing_model) {
2536   return routing_model.solver()->RevAlloc(new VehicleVarFilter(routing_model));
2537 }
2538 
2539 namespace {
2540 
2541 class CumulBoundsPropagatorFilter : public IntVarLocalSearchFilter {
2542  public:
2543   explicit CumulBoundsPropagatorFilter(const RoutingDimension& dimension);
2544   bool Accept(const Assignment* delta, const Assignment* deltadelta,
2545               int64_t objective_min, int64_t objective_max) override;
DebugString() const2546   std::string DebugString() const override {
2547     return "CumulBoundsPropagatorFilter(" + propagator_.dimension().name() +
2548            ")";
2549   }
2550 
2551  private:
2552   CumulBoundsPropagator propagator_;
2553   const int64_t cumul_offset_;
2554   SparseBitset<int64_t> delta_touched_;
2555   std::vector<int64_t> delta_nexts_;
2556 };
2557 
CumulBoundsPropagatorFilter(const RoutingDimension & dimension)2558 CumulBoundsPropagatorFilter::CumulBoundsPropagatorFilter(
2559     const RoutingDimension& dimension)
2560     : IntVarLocalSearchFilter(dimension.model()->Nexts()),
2561       propagator_(&dimension),
2562       cumul_offset_(dimension.GetGlobalOptimizerOffset()),
2563       delta_touched_(Size()),
2564       delta_nexts_(Size()) {}
2565 
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)2566 bool CumulBoundsPropagatorFilter::Accept(const Assignment* delta,
2567                                          const Assignment* deltadelta,
2568                                          int64_t objective_min,
2569                                          int64_t objective_max) {
2570   delta_touched_.ClearAll();
2571   for (const IntVarElement& delta_element :
2572        delta->IntVarContainer().elements()) {
2573     int64_t index = -1;
2574     if (FindIndex(delta_element.Var(), &index)) {
2575       if (!delta_element.Bound()) {
2576         // LNS detected
2577         return true;
2578       }
2579       delta_touched_.Set(index);
2580       delta_nexts_[index] = delta_element.Value();
2581     }
2582   }
2583   const auto& next_accessor = [this](int64_t index) {
2584     return delta_touched_[index] ? delta_nexts_[index] : Value(index);
2585   };
2586 
2587   return propagator_.PropagateCumulBounds(next_accessor, cumul_offset_);
2588 }
2589 
2590 }  // namespace
2591 
MakeCumulBoundsPropagatorFilter(const RoutingDimension & dimension)2592 IntVarLocalSearchFilter* MakeCumulBoundsPropagatorFilter(
2593     const RoutingDimension& dimension) {
2594   return dimension.model()->solver()->RevAlloc(
2595       new CumulBoundsPropagatorFilter(dimension));
2596 }
2597 
2598 namespace {
2599 
2600 class LPCumulFilter : public IntVarLocalSearchFilter {
2601  public:
2602   LPCumulFilter(const std::vector<IntVar*>& nexts,
2603                 GlobalDimensionCumulOptimizer* optimizer,
2604                 GlobalDimensionCumulOptimizer* mp_optimizer,
2605                 bool filter_objective_cost);
2606   bool Accept(const Assignment* delta, const Assignment* deltadelta,
2607               int64_t objective_min, int64_t objective_max) override;
2608   int64_t GetAcceptedObjectiveValue() const override;
2609   void OnSynchronize(const Assignment* delta) override;
2610   int64_t GetSynchronizedObjectiveValue() const override;
DebugString() const2611   std::string DebugString() const override {
2612     return "LPCumulFilter(" + optimizer_.dimension()->name() + ")";
2613   }
2614 
2615  private:
2616   GlobalDimensionCumulOptimizer& optimizer_;
2617   GlobalDimensionCumulOptimizer& mp_optimizer_;
2618   const bool filter_objective_cost_;
2619   int64_t synchronized_cost_without_transit_;
2620   int64_t delta_cost_without_transit_;
2621   SparseBitset<int64_t> delta_touched_;
2622   std::vector<int64_t> delta_nexts_;
2623 };
2624 
LPCumulFilter(const std::vector<IntVar * > & nexts,GlobalDimensionCumulOptimizer * optimizer,GlobalDimensionCumulOptimizer * mp_optimizer,bool filter_objective_cost)2625 LPCumulFilter::LPCumulFilter(const std::vector<IntVar*>& nexts,
2626                              GlobalDimensionCumulOptimizer* optimizer,
2627                              GlobalDimensionCumulOptimizer* mp_optimizer,
2628                              bool filter_objective_cost)
2629     : IntVarLocalSearchFilter(nexts),
2630       optimizer_(*optimizer),
2631       mp_optimizer_(*mp_optimizer),
2632       filter_objective_cost_(filter_objective_cost),
2633       synchronized_cost_without_transit_(-1),
2634       delta_cost_without_transit_(-1),
2635       delta_touched_(Size()),
2636       delta_nexts_(Size()) {}
2637 
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)2638 bool LPCumulFilter::Accept(const Assignment* delta,
2639                            const Assignment* deltadelta, int64_t objective_min,
2640                            int64_t objective_max) {
2641   delta_touched_.ClearAll();
2642   for (const IntVarElement& delta_element :
2643        delta->IntVarContainer().elements()) {
2644     int64_t index = -1;
2645     if (FindIndex(delta_element.Var(), &index)) {
2646       if (!delta_element.Bound()) {
2647         // LNS detected
2648         return true;
2649       }
2650       delta_touched_.Set(index);
2651       delta_nexts_[index] = delta_element.Value();
2652     }
2653   }
2654   const auto& next_accessor = [this](int64_t index) {
2655     return delta_touched_[index] ? delta_nexts_[index] : Value(index);
2656   };
2657 
2658   if (!filter_objective_cost_) {
2659     // No need to compute the cost of the LP, only verify its feasibility.
2660     delta_cost_without_transit_ = 0;
2661     const DimensionSchedulingStatus status =
2662         optimizer_.ComputeCumuls(next_accessor, nullptr, nullptr, nullptr);
2663     if (status == DimensionSchedulingStatus::OPTIMAL) return true;
2664     if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY &&
2665         mp_optimizer_.ComputeCumuls(next_accessor, nullptr, nullptr, nullptr) ==
2666             DimensionSchedulingStatus::OPTIMAL) {
2667       return true;
2668     }
2669     return false;
2670   }
2671 
2672   const DimensionSchedulingStatus status =
2673       optimizer_.ComputeCumulCostWithoutFixedTransits(
2674           next_accessor, &delta_cost_without_transit_);
2675   if (status == DimensionSchedulingStatus::INFEASIBLE) {
2676     delta_cost_without_transit_ = std::numeric_limits<int64_t>::max();
2677     return false;
2678   }
2679   if (delta_cost_without_transit_ > objective_max) return false;
2680 
2681   if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY &&
2682       mp_optimizer_.ComputeCumulCostWithoutFixedTransits(
2683           next_accessor, &delta_cost_without_transit_) !=
2684           DimensionSchedulingStatus::OPTIMAL) {
2685     delta_cost_without_transit_ = std::numeric_limits<int64_t>::max();
2686     return false;
2687   }
2688   return delta_cost_without_transit_ <= objective_max;
2689 }
2690 
GetAcceptedObjectiveValue() const2691 int64_t LPCumulFilter::GetAcceptedObjectiveValue() const {
2692   return delta_cost_without_transit_;
2693 }
2694 
OnSynchronize(const Assignment * delta)2695 void LPCumulFilter::OnSynchronize(const Assignment* delta) {
2696   // TODO(user): Try to optimize this so the LP is not called when the last
2697   // computed delta cost corresponds to the solution being synchronized.
2698   const RoutingModel& model = *optimizer_.dimension()->model();
2699   const auto& next_accessor = [this, &model](int64_t index) {
2700     return IsVarSynced(index)     ? Value(index)
2701            : model.IsStart(index) ? model.End(model.VehicleIndex(index))
2702                                   : index;
2703   };
2704 
2705   const DimensionSchedulingStatus status =
2706       optimizer_.ComputeCumulCostWithoutFixedTransits(
2707           next_accessor, &synchronized_cost_without_transit_);
2708   if (status == DimensionSchedulingStatus::INFEASIBLE) {
2709     // TODO(user): This should only happen if the LP solver times out.
2710     // DCHECK the fail wasn't due to an infeasible model.
2711     synchronized_cost_without_transit_ = 0;
2712   }
2713   if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY &&
2714       mp_optimizer_.ComputeCumulCostWithoutFixedTransits(
2715           next_accessor, &synchronized_cost_without_transit_) !=
2716           DimensionSchedulingStatus::OPTIMAL) {
2717     // TODO(user): This should only happen if the MP solver times out.
2718     // DCHECK the fail wasn't due to an infeasible model.
2719     synchronized_cost_without_transit_ = 0;
2720   }
2721 }
2722 
GetSynchronizedObjectiveValue() const2723 int64_t LPCumulFilter::GetSynchronizedObjectiveValue() const {
2724   return synchronized_cost_without_transit_;
2725 }
2726 
2727 }  // namespace
2728 
MakeGlobalLPCumulFilter(GlobalDimensionCumulOptimizer * optimizer,GlobalDimensionCumulOptimizer * mp_optimizer,bool filter_objective_cost)2729 IntVarLocalSearchFilter* MakeGlobalLPCumulFilter(
2730     GlobalDimensionCumulOptimizer* optimizer,
2731     GlobalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost) {
2732   const RoutingModel& model = *optimizer->dimension()->model();
2733   return model.solver()->RevAlloc(new LPCumulFilter(
2734       model.Nexts(), optimizer, mp_optimizer, filter_objective_cost));
2735 }
2736 
2737 namespace {
2738 
2739 using ResourceGroup = RoutingModel::ResourceGroup;
2740 
2741 class ResourceGroupAssignmentFilter : public BasePathFilter {
2742  public:
2743   ResourceGroupAssignmentFilter(const std::vector<IntVar*>& nexts,
2744                                 const ResourceGroup* resource_group,
2745                                 LocalDimensionCumulOptimizer* optimizer,
2746                                 LocalDimensionCumulOptimizer* mp_optimizer,
2747                                 bool filter_objective_cost);
2748   bool InitializeAcceptPath() override;
2749   bool AcceptPath(int64_t path_start, int64_t chain_start,
2750                   int64_t chain_end) override;
2751   bool FinalizeAcceptPath(int64_t objective_min,
2752                           int64_t objective_max) override;
2753   void OnSynchronizePathFromStart(int64_t start) override;
2754   void OnAfterSynchronizePaths() override;
2755 
GetAcceptedObjectiveValue() const2756   int64_t GetAcceptedObjectiveValue() const override {
2757     return delta_cost_without_transit_;
2758   }
GetSynchronizedObjectiveValue() const2759   int64_t GetSynchronizedObjectiveValue() const override {
2760     return synchronized_cost_without_transit_;
2761   }
DebugString() const2762   std::string DebugString() const override {
2763     return "ResourceGroupAssignmentFilter(" +
2764            resource_assignment_optimizer_.dimension()->name() + ")";
2765   }
2766 
2767  private:
2768   ResourceAssignmentOptimizer resource_assignment_optimizer_;
2769   const RoutingModel& model_;
2770   const ResourceGroup& resource_group_;
2771   const bool filter_objective_cost_;
2772   bool synch_timed_out_;
2773   int64_t synchronized_cost_without_transit_;
2774   int64_t delta_cost_without_transit_;
2775   std::vector<std::vector<int64_t>> vehicle_to_resource_assignment_costs_;
2776   std::vector<std::vector<int64_t>> delta_vehicle_to_resource_assignment_costs_;
2777 };
2778 
ResourceGroupAssignmentFilter(const std::vector<IntVar * > & nexts,const ResourceGroup * resource_group,LocalDimensionCumulOptimizer * optimizer,LocalDimensionCumulOptimizer * mp_optimizer,bool filter_objective_cost)2779 ResourceGroupAssignmentFilter::ResourceGroupAssignmentFilter(
2780     const std::vector<IntVar*>& nexts, const ResourceGroup* resource_group,
2781     LocalDimensionCumulOptimizer* optimizer,
2782     LocalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost)
2783     : BasePathFilter(nexts, optimizer->dimension()->cumuls().size()),
2784       resource_assignment_optimizer_(resource_group, optimizer, mp_optimizer),
2785       model_(*optimizer->dimension()->model()),
2786       resource_group_(*resource_group),
2787       filter_objective_cost_(filter_objective_cost),
2788       synch_timed_out_(false),
2789       synchronized_cost_without_transit_(-1),
2790       delta_cost_without_transit_(-1) {
2791   vehicle_to_resource_assignment_costs_.resize(model_.vehicles());
2792   delta_vehicle_to_resource_assignment_costs_.resize(model_.vehicles());
2793 }
2794 
InitializeAcceptPath()2795 bool ResourceGroupAssignmentFilter::InitializeAcceptPath() {
2796   delta_vehicle_to_resource_assignment_costs_.assign(model_.vehicles(), {});
2797   // TODO(user): Keep track of num_used_vehicles internally and compute its
2798   // new value here by only going through the touched_paths_.
2799   int num_used_vehicles = 0;
2800   const int num_resources = resource_group_.Size();
2801   for (int v = 0; v < model_.vehicles(); v++) {
2802     if (GetNext(model_.Start(v)) != model_.End(v) ||
2803         model_.IsVehicleUsedWhenEmpty(v)) {
2804       if (++num_used_vehicles > num_resources) {
2805         return false;
2806       }
2807     }
2808   }
2809   return true;
2810 }
2811 
AcceptPath(int64_t path_start,int64_t chain_start,int64_t chain_end)2812 bool ResourceGroupAssignmentFilter::AcceptPath(int64_t path_start,
2813                                                int64_t chain_start,
2814                                                int64_t chain_end) {
2815   const int vehicle = model_.VehicleIndex(path_start);
2816   return resource_assignment_optimizer_.ComputeAssignmentCostsForVehicle(
2817       vehicle, [this](int64_t index) { return GetNext(index); },
2818       filter_objective_cost_,
2819       &delta_vehicle_to_resource_assignment_costs_[vehicle], nullptr, nullptr);
2820 }
2821 
FinalizeAcceptPath(int64_t objective_min,int64_t objective_max)2822 bool ResourceGroupAssignmentFilter::FinalizeAcceptPath(int64_t objective_min,
2823                                                        int64_t objective_max) {
2824   delta_cost_without_transit_ =
2825       resource_assignment_optimizer_.ComputeBestAssignmentCost(
2826           delta_vehicle_to_resource_assignment_costs_,
2827           vehicle_to_resource_assignment_costs_,
2828           [this](int v) { return PathStartTouched(model_.Start(v)); }, nullptr);
2829   return delta_cost_without_transit_ >= 0 &&
2830          delta_cost_without_transit_ <= objective_max;
2831 }
2832 
OnSynchronizePathFromStart(int64_t start)2833 void ResourceGroupAssignmentFilter::OnSynchronizePathFromStart(int64_t start) {
2834   if (synch_timed_out_) return;
2835   // NOTE(user): Even if filter_objective_cost_ is false, we still need to
2836   // call ComputeAssignmentCostsForVehicle() for every vehicle to keep track
2837   // of whether or not a given vehicle-to-resource assignment is possible by
2838   // storing 0 or -1 in vehicle_to_resource_assignment_costs_.
2839   const auto& next_accessor = [this](int64_t index) {
2840     return IsVarSynced(index)      ? Value(index)
2841            : model_.IsStart(index) ? model_.End(model_.VehicleIndex(index))
2842                                    : index;
2843   };
2844   const int v = model_.VehicleIndex(start);
2845   if (!resource_assignment_optimizer_.ComputeAssignmentCostsForVehicle(
2846           v, next_accessor, filter_objective_cost_,
2847           &vehicle_to_resource_assignment_costs_[v], nullptr, nullptr)) {
2848     // A timeout was reached.
2849     synch_timed_out_ = true;
2850   }
2851 }
2852 
OnAfterSynchronizePaths()2853 void ResourceGroupAssignmentFilter::OnAfterSynchronizePaths() {
2854   synchronized_cost_without_transit_ =
2855       (synch_timed_out_ || !filter_objective_cost_)
2856           ? 0
2857           : resource_assignment_optimizer_.ComputeBestAssignmentCost(
2858                 vehicle_to_resource_assignment_costs_,
2859                 vehicle_to_resource_assignment_costs_, [](int) { return true; },
2860                 nullptr);
2861   DCHECK_GE(synchronized_cost_without_transit_, 0);
2862 }
2863 
2864 // ResourceAssignmentFilter
2865 class ResourceAssignmentFilter : public LocalSearchFilter {
2866  public:
2867   ResourceAssignmentFilter(const std::vector<IntVar*>& nexts,
2868                            LocalDimensionCumulOptimizer* optimizer,
2869                            LocalDimensionCumulOptimizer* mp_optimizer,
2870                            bool propagate_own_objective_value,
2871                            bool filter_objective_cost);
2872   bool Accept(const Assignment* delta, const Assignment* deltadelta,
2873               int64_t objective_min, int64_t objective_max) override;
2874   void Synchronize(const Assignment* assignment,
2875                    const Assignment* delta) override;
2876 
GetAcceptedObjectiveValue() const2877   int64_t GetAcceptedObjectiveValue() const override {
2878     return propagate_own_objective_value_ ? delta_cost_ : 0;
2879   }
GetSynchronizedObjectiveValue() const2880   int64_t GetSynchronizedObjectiveValue() const override {
2881     return propagate_own_objective_value_ ? synchronized_cost_ : 0;
2882   }
DebugString() const2883   std::string DebugString() const override {
2884     return "ResourceAssignmentFilter(" + dimension_name_ + ")";
2885   }
2886 
2887  private:
2888   std::vector<IntVarLocalSearchFilter*> resource_group_assignment_filters_;
2889   int64_t synchronized_cost_;
2890   int64_t delta_cost_;
2891   const bool propagate_own_objective_value_;
2892   const std::string dimension_name_;
2893 };
2894 
ResourceAssignmentFilter(const std::vector<IntVar * > & nexts,LocalDimensionCumulOptimizer * optimizer,LocalDimensionCumulOptimizer * mp_optimizer,bool propagate_own_objective_value,bool filter_objective_cost)2895 ResourceAssignmentFilter::ResourceAssignmentFilter(
2896     const std::vector<IntVar*>& nexts, LocalDimensionCumulOptimizer* optimizer,
2897     LocalDimensionCumulOptimizer* mp_optimizer,
2898     bool propagate_own_objective_value, bool filter_objective_cost)
2899     : propagate_own_objective_value_(propagate_own_objective_value),
2900       dimension_name_(optimizer->dimension()->name()) {
2901   const RoutingModel& model = *optimizer->dimension()->model();
2902   for (const auto& resource_group : model.GetResourceGroups()) {
2903     resource_group_assignment_filters_.push_back(
2904         model.solver()->RevAlloc(new ResourceGroupAssignmentFilter(
2905             nexts, resource_group.get(), optimizer, mp_optimizer,
2906             filter_objective_cost)));
2907   }
2908 }
2909 
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)2910 bool ResourceAssignmentFilter::Accept(const Assignment* delta,
2911                                       const Assignment* deltadelta,
2912                                       int64_t objective_min,
2913                                       int64_t objective_max) {
2914   delta_cost_ = 0;
2915   for (LocalSearchFilter* group_filter : resource_group_assignment_filters_) {
2916     if (!group_filter->Accept(delta, deltadelta, objective_min,
2917                               objective_max)) {
2918       return false;
2919     }
2920     delta_cost_ =
2921         std::max(delta_cost_, group_filter->GetAcceptedObjectiveValue());
2922     DCHECK_LE(delta_cost_, objective_max)
2923         << "ResourceGroupAssignmentFilter should return false when the "
2924            "objective_max is exceeded.";
2925   }
2926   return true;
2927 }
2928 
Synchronize(const Assignment * assignment,const Assignment * delta)2929 void ResourceAssignmentFilter::Synchronize(const Assignment* assignment,
2930                                            const Assignment* delta) {
2931   synchronized_cost_ = 0;
2932   for (LocalSearchFilter* group_filter : resource_group_assignment_filters_) {
2933     group_filter->Synchronize(assignment, delta);
2934     synchronized_cost_ = std::max(
2935         synchronized_cost_, group_filter->GetSynchronizedObjectiveValue());
2936   }
2937 }
2938 
2939 }  // namespace
2940 
MakeResourceAssignmentFilter(LocalDimensionCumulOptimizer * optimizer,LocalDimensionCumulOptimizer * mp_optimizer,bool propagate_own_objective_value,bool filter_objective_cost)2941 LocalSearchFilter* MakeResourceAssignmentFilter(
2942     LocalDimensionCumulOptimizer* optimizer,
2943     LocalDimensionCumulOptimizer* mp_optimizer,
2944     bool propagate_own_objective_value, bool filter_objective_cost) {
2945   const RoutingModel& model = *optimizer->dimension()->model();
2946   return model.solver()->RevAlloc(new ResourceAssignmentFilter(
2947       model.Nexts(), optimizer, mp_optimizer, propagate_own_objective_value,
2948       filter_objective_cost));
2949 }
2950 
2951 namespace {
2952 
2953 // This filter accepts deltas for which the assignment satisfies the
2954 // constraints of the Solver. This is verified by keeping an internal copy of
2955 // the assignment with all Next vars and their updated values, and calling
2956 // RestoreAssignment() on the assignment+delta.
2957 // TODO(user): Also call the solution finalizer on variables, with the
2958 // exception of Next Vars (woud fail on large instances).
2959 // WARNING: In the case of mandatory nodes, when all vehicles are currently
2960 // being used in the solution but uninserted nodes still remain, this filter
2961 // will reject the solution, even if the node could be inserted on one of these
2962 // routes, because all Next vars of vehicle starts are already instantiated.
2963 // TODO(user): Avoid such false negatives.
2964 
2965 class CPFeasibilityFilter : public IntVarLocalSearchFilter {
2966  public:
2967   explicit CPFeasibilityFilter(RoutingModel* routing_model);
~CPFeasibilityFilter()2968   ~CPFeasibilityFilter() override {}
DebugString() const2969   std::string DebugString() const override { return "CPFeasibilityFilter"; }
2970   bool Accept(const Assignment* delta, const Assignment* deltadelta,
2971               int64_t objective_min, int64_t objective_max) override;
2972   void OnSynchronize(const Assignment* delta) override;
2973 
2974  private:
2975   void AddDeltaToAssignment(const Assignment* delta, Assignment* assignment);
2976 
2977   static const int64_t kUnassigned;
2978   const RoutingModel* const model_;
2979   Solver* const solver_;
2980   Assignment* const assignment_;
2981   Assignment* const temp_assignment_;
2982   DecisionBuilder* const restore_;
2983   SearchLimit* const limit_;
2984 };
2985 
2986 const int64_t CPFeasibilityFilter::kUnassigned = -1;
2987 
CPFeasibilityFilter(RoutingModel * routing_model)2988 CPFeasibilityFilter::CPFeasibilityFilter(RoutingModel* routing_model)
2989     : IntVarLocalSearchFilter(routing_model->Nexts()),
2990       model_(routing_model),
2991       solver_(routing_model->solver()),
2992       assignment_(solver_->MakeAssignment()),
2993       temp_assignment_(solver_->MakeAssignment()),
2994       restore_(solver_->MakeRestoreAssignment(temp_assignment_)),
2995       limit_(solver_->MakeCustomLimit(
2996           [routing_model]() { return routing_model->CheckLimit(); })) {
2997   assignment_->Add(routing_model->Nexts());
2998 }
2999 
Accept(const Assignment * delta,const Assignment * deltadelta,int64_t objective_min,int64_t objective_max)3000 bool CPFeasibilityFilter::Accept(const Assignment* delta,
3001                                  const Assignment* deltadelta,
3002                                  int64_t objective_min, int64_t objective_max) {
3003   temp_assignment_->Copy(assignment_);
3004   AddDeltaToAssignment(delta, temp_assignment_);
3005 
3006   return solver_->Solve(restore_, limit_);
3007 }
3008 
OnSynchronize(const Assignment * delta)3009 void CPFeasibilityFilter::OnSynchronize(const Assignment* delta) {
3010   AddDeltaToAssignment(delta, assignment_);
3011 }
3012 
AddDeltaToAssignment(const Assignment * delta,Assignment * assignment)3013 void CPFeasibilityFilter::AddDeltaToAssignment(const Assignment* delta,
3014                                                Assignment* assignment) {
3015   if (delta == nullptr) {
3016     return;
3017   }
3018   Assignment::IntContainer* const container =
3019       assignment->MutableIntVarContainer();
3020   const Assignment::IntContainer& delta_container = delta->IntVarContainer();
3021   const int delta_size = delta_container.Size();
3022 
3023   for (int i = 0; i < delta_size; i++) {
3024     const IntVarElement& delta_element = delta_container.Element(i);
3025     IntVar* const var = delta_element.Var();
3026     int64_t index = kUnassigned;
3027     CHECK(FindIndex(var, &index));
3028     DCHECK_EQ(var, Var(index));
3029     const int64_t value = delta_element.Value();
3030 
3031     container->AddAtPosition(var, index)->SetValue(value);
3032     if (model_->IsStart(index)) {
3033       if (model_->IsEnd(value)) {
3034         // Do not restore unused routes.
3035         container->MutableElement(index)->Deactivate();
3036       } else {
3037         // Re-activate the route's start in case it was deactivated before.
3038         container->MutableElement(index)->Activate();
3039       }
3040     }
3041   }
3042 }
3043 
3044 }  // namespace
3045 
MakeCPFeasibilityFilter(RoutingModel * routing_model)3046 IntVarLocalSearchFilter* MakeCPFeasibilityFilter(RoutingModel* routing_model) {
3047   return routing_model->solver()->RevAlloc(
3048       new CPFeasibilityFilter(routing_model));
3049 }
3050 
3051 // TODO(user): Implement same-vehicle filter. Could be merged with node
3052 // precedence filter.
3053 
3054 }  // namespace operations_research
3055