diff --git a/ortools/constraint_solver/BUILD.bazel b/ortools/constraint_solver/BUILD.bazel index dc2e10904f..0b76a700c0 100644 --- a/ortools/constraint_solver/BUILD.bazel +++ b/ortools/constraint_solver/BUILD.bazel @@ -333,25 +333,17 @@ cc_library( ":routing_parameters_cc_proto", ":routing_types", "//ortools/base", - "@com_google_protobuf//:protobuf", - "@com_google_absl//absl/container:btree", - "@com_google_absl//absl/container:flat_hash_map", - "@com_google_absl//absl/container:flat_hash_set", - "@com_google_absl//absl/functional:bind_front", - "@com_google_absl//absl/hash", - "@com_google_absl//absl/memory", - "@com_google_absl//absl/strings", - "@com_google_absl//absl/strings:str_format", - "@com_google_absl//absl/time", - # "//kdtree", - "//ortools/graph", - "//ortools/base:strong_vector", + "//ortools/base:adjustable_priority_queue", + "//ortools/base:dump_vars", + "//ortools/base:hash", "//ortools/base:map_util", + "//ortools/base:murmur", + "//ortools/base:protoutil", "//ortools/base:small_map", "//ortools/base:stl_util", - "//ortools/base:hash", - "//ortools/base:murmur", + "//ortools/base:strong_vector", "//ortools/glop:lp_solver", + "//ortools/graph", "//ortools/graph:christofides", "//ortools/graph:connected_components", "//ortools/graph:linear_assignment", @@ -372,11 +364,21 @@ cc_library( "//ortools/sat:optimization", "//ortools/sat:theta_tree", "//ortools/util:bitset", + "//ortools/util:flat_matrix", "//ortools/util:optional_boolean_cc_proto", "//ortools/util:range_query_function", "//ortools/util:saturated_arithmetic", "//ortools/util:sorted_interval_list", - "//ortools/base:adjustable_priority_queue", - "//ortools/base:protoutil", + "@com_google_absl//absl/container:btree", + "@com_google_absl//absl/container:flat_hash_map", + "@com_google_absl//absl/container:flat_hash_set", + "@com_google_absl//absl/container:inlined_vector", + "@com_google_absl//absl/functional:bind_front", + "@com_google_absl//absl/hash", + "@com_google_absl//absl/memory", + "@com_google_absl//absl/strings", + "@com_google_absl//absl/strings:str_format", + "@com_google_absl//absl/time", + "@com_google_protobuf//:protobuf", ], ) diff --git a/ortools/constraint_solver/assignment.cc b/ortools/constraint_solver/assignment.cc index 6808a70a71..df77d2eb59 100644 --- a/ortools/constraint_solver/assignment.cc +++ b/ortools/constraint_solver/assignment.cc @@ -875,14 +875,6 @@ void Assignment::SetUnperformed(const SequenceVar* const var, // ----- Objective ----- -void Assignment::AddObjective(IntVar* const v) { - // Check if adding twice an objective to the solution. - CHECK(!HasObjective()); - objective_element_.Reset(v); -} - -IntVar* Assignment::Objective() const { return objective_element_.Var(); } - int64_t Assignment::ObjectiveMin() const { if (HasObjective()) { return objective_element_.Min(); diff --git a/ortools/constraint_solver/constraint_solver.cc b/ortools/constraint_solver/constraint_solver.cc index f2a9ce757b..4aa8eff994 100644 --- a/ortools/constraint_solver/constraint_solver.cc +++ b/ortools/constraint_solver/constraint_solver.cc @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -118,6 +119,13 @@ void ForAll(const std::vector& objects, MethodPointer method, (object->*method)(args...); } } + +// Converts a scoped enum to its underlying type. +template +constexpr typename std::underlying_type::type to_underlying(E e) { + return static_cast::type>(e); +} + } // namespace // ----- ConstraintSolverParameters ----- @@ -707,7 +715,7 @@ class CompressedTrail { } } - std::unique_ptr > packer_; + std::unique_ptr> packer_; const int block_size_; Block* blocks_; Block* free_blocks_; @@ -962,6 +970,7 @@ class Search { explicit Search(Solver* const s) : solver_(s), marker_stack_(), + monitor_event_listeners_(to_underlying(Solver::MonitorEvent::kLast)), fail_buffer_(), solution_counter_(0), unchecked_solution_counter_(0), @@ -981,6 +990,7 @@ class Search { Search(Solver* const s, int /* dummy_argument */) : solver_(s), marker_stack_(), + monitor_event_listeners_(to_underlying(Solver::MonitorEvent::kLast)), fail_buffer_(), solution_counter_(0), unchecked_solution_counter_(0), @@ -1019,7 +1029,15 @@ class Search { void PeriodicCheck(); int ProgressPercent(); void Accept(ModelVisitor* const visitor) const; - void push_monitor(SearchMonitor* const m); + void AddEventListener(Solver::MonitorEvent event, SearchMonitor* monitor) { + if (monitor != nullptr) { + monitor_event_listeners_[to_underlying(event)].push_back(monitor); + } + } + const std::vector& GetEventListeners( + Solver::MonitorEvent event) const { + return monitor_event_listeners_[to_underlying(event)]; + } void Clear(); void IncrementSolutionCounter() { ++solution_counter_; } int64_t solution_counter() const { return solution_counter_; } @@ -1075,7 +1093,7 @@ class Search { Solver* const solver_; std::vector marker_stack_; - std::vector monitors_; + std::vector> monitor_event_listeners_; jmp_buf fail_buffer_; int64_t solution_counter_; int64_t unchecked_solution_counter_; @@ -1192,20 +1210,20 @@ Solver::DecisionModification Search::ModifyDecision() { return Solver::NO_CHANGE; } -void Search::push_monitor(SearchMonitor* const m) { - if (m) { - monitors_.push_back(m); - } -} - void Search::Clear() { - monitors_.clear(); + for (auto& listeners : monitor_event_listeners_) listeners.clear(); search_depth_ = 0; left_search_depth_ = 0; selector_ = nullptr; backtrack_at_the_end_of_the_search_ = true; } +#define CALL_EVENT_LISTENERS(Event) \ + do { \ + ForAll(GetEventListeners(Solver::MonitorEvent::k##Event), \ + &SearchMonitor::Event); \ + } while (false) + void Search::EnterSearch() { // The solution counter is reset when entering search and not when // leaving search. This enables the information to persist outside of @@ -1213,58 +1231,62 @@ void Search::EnterSearch() { solution_counter_ = 0; unchecked_solution_counter_ = 0; - ForAll(monitors_, &SearchMonitor::EnterSearch); + CALL_EVENT_LISTENERS(EnterSearch); } void Search::ExitSearch() { // Backtrack to the correct state. - ForAll(monitors_, &SearchMonitor::ExitSearch); + CALL_EVENT_LISTENERS(ExitSearch); } -void Search::RestartSearch() { - ForAll(monitors_, &SearchMonitor::RestartSearch); -} +void Search::RestartSearch() { CALL_EVENT_LISTENERS(RestartSearch); } void Search::BeginNextDecision(DecisionBuilder* const db) { - ForAll(monitors_, &SearchMonitor::BeginNextDecision, db); + ForAll(GetEventListeners(Solver::MonitorEvent::kBeginNextDecision), + &SearchMonitor::BeginNextDecision, db); CheckFail(); } void Search::EndNextDecision(DecisionBuilder* const db, Decision* const d) { - ForAll(monitors_, &SearchMonitor::EndNextDecision, db, d); + ForAll(GetEventListeners(Solver::MonitorEvent::kEndNextDecision), + &SearchMonitor::EndNextDecision, db, d); CheckFail(); } void Search::ApplyDecision(Decision* const d) { - ForAll(monitors_, &SearchMonitor::ApplyDecision, d); + ForAll(GetEventListeners(Solver::MonitorEvent::kApplyDecision), + &SearchMonitor::ApplyDecision, d); CheckFail(); } void Search::AfterDecision(Decision* const d, bool apply) { - ForAll(monitors_, &SearchMonitor::AfterDecision, d, apply); + ForAll(GetEventListeners(Solver::MonitorEvent::kAfterDecision), + &SearchMonitor::AfterDecision, d, apply); CheckFail(); } void Search::RefuteDecision(Decision* const d) { - ForAll(monitors_, &SearchMonitor::RefuteDecision, d); + ForAll(GetEventListeners(Solver::MonitorEvent::kRefuteDecision), + &SearchMonitor::RefuteDecision, d); CheckFail(); } -void Search::BeginFail() { ForAll(monitors_, &SearchMonitor::BeginFail); } +void Search::BeginFail() { CALL_EVENT_LISTENERS(BeginFail); } -void Search::EndFail() { ForAll(monitors_, &SearchMonitor::EndFail); } +void Search::EndFail() { CALL_EVENT_LISTENERS(EndFail); } void Search::BeginInitialPropagation() { - ForAll(monitors_, &SearchMonitor::BeginInitialPropagation); + CALL_EVENT_LISTENERS(BeginInitialPropagation); } void Search::EndInitialPropagation() { - ForAll(monitors_, &SearchMonitor::EndInitialPropagation); + CALL_EVENT_LISTENERS(EndInitialPropagation); } bool Search::AcceptSolution() { bool valid = true; - for (SearchMonitor* const monitor : monitors_) { + for (SearchMonitor* const monitor : + GetEventListeners(Solver::MonitorEvent::kAcceptSolution)) { if (!monitor->AcceptSolution()) { // Even though we know the return value, we cannot return yet: this would // break the contract we have with solution monitors. They all deserve @@ -1277,7 +1299,8 @@ bool Search::AcceptSolution() { bool Search::AtSolution() { bool should_continue = false; - for (SearchMonitor* const monitor : monitors_) { + for (SearchMonitor* const monitor : + GetEventListeners(Solver::MonitorEvent::kAtSolution)) { if (monitor->AtSolution()) { // Even though we know the return value, we cannot return yet: this would // break the contract we have with solution monitors. They all deserve @@ -1288,23 +1311,23 @@ bool Search::AtSolution() { return should_continue; } -void Search::NoMoreSolutions() { - ForAll(monitors_, &SearchMonitor::NoMoreSolutions); -} +void Search::NoMoreSolutions() { CALL_EVENT_LISTENERS(NoMoreSolutions); } bool Search::LocalOptimum() { - bool res = false; - for (SearchMonitor* const monitor : monitors_) { + bool at_local_optimum = false; + for (SearchMonitor* const monitor : + GetEventListeners(Solver::MonitorEvent::kLocalOptimum)) { if (monitor->LocalOptimum()) { - res = true; + at_local_optimum = true; } } - return res; + return at_local_optimum; } bool Search::AcceptDelta(Assignment* delta, Assignment* deltadelta) { bool accept = true; - for (SearchMonitor* const monitor : monitors_) { + for (SearchMonitor* const monitor : + GetEventListeners(Solver::MonitorEvent::kAcceptDelta)) { if (!monitor->AcceptDelta(delta, deltadelta)) { accept = false; } @@ -1312,16 +1335,15 @@ bool Search::AcceptDelta(Assignment* delta, Assignment* deltadelta) { return accept; } -void Search::AcceptNeighbor() { - ForAll(monitors_, &SearchMonitor::AcceptNeighbor); -} +void Search::AcceptNeighbor() { CALL_EVENT_LISTENERS(AcceptNeighbor); } void Search::AcceptUncheckedNeighbor() { - ForAll(monitors_, &SearchMonitor::AcceptUncheckedNeighbor); + CALL_EVENT_LISTENERS(AcceptUncheckedNeighbor); } bool Search::IsUncheckedSolutionLimitReached() { - for (SearchMonitor* const monitor : monitors_) { + for (SearchMonitor* const monitor : GetEventListeners( + Solver::MonitorEvent::kIsUncheckedSolutionLimitReached)) { if (monitor->IsUncheckedSolutionLimitReached()) { return true; } @@ -1329,25 +1351,27 @@ bool Search::IsUncheckedSolutionLimitReached() { return false; } -void Search::PeriodicCheck() { - ForAll(monitors_, &SearchMonitor::PeriodicCheck); -} +void Search::PeriodicCheck() { CALL_EVENT_LISTENERS(PeriodicCheck); } int Search::ProgressPercent() { int progress = SearchMonitor::kNoProgress; - for (SearchMonitor* const monitor : monitors_) { + for (SearchMonitor* const monitor : + GetEventListeners(Solver::MonitorEvent::kProgressPercent)) { progress = std::max(progress, monitor->ProgressPercent()); } return progress; } void Search::Accept(ModelVisitor* const visitor) const { - ForAll(monitors_, &SearchMonitor::Accept, visitor); + ForAll(GetEventListeners(Solver::MonitorEvent::kAccept), + &SearchMonitor::Accept, visitor); if (decision_builder_ != nullptr) { decision_builder_->Accept(visitor); } } +#undef CALL_EVENT_LISTENERS + bool LocalOptimumReached(Search* const search) { return search->LocalOptimum(); } @@ -2587,6 +2611,7 @@ const char ModelVisitor::kDistribute[] = "Distribute"; const char ModelVisitor::kDivide[] = "Divide"; const char ModelVisitor::kDurationExpr[] = "DurationExpression"; const char ModelVisitor::kElement[] = "Element"; +const char ModelVisitor::kLightElementEqual[] = "LightElementEqual"; const char ModelVisitor::kElementEqual[] = "ElementEqual"; const char ModelVisitor::kEndExpr[] = "EndExpression"; const char ModelVisitor::kEquality[] = "Equal"; @@ -2903,7 +2928,14 @@ void SearchMonitor::PeriodicCheck() {} void SearchMonitor::Accept(ModelVisitor* const visitor) const {} // A search monitors adds itself on the active search. void SearchMonitor::Install() { - solver()->searches_.back()->push_monitor(this); + for (std::underlying_type::type event = 0; + event != to_underlying(Solver::MonitorEvent::kLast); ++event) { + ListenToEvent(static_cast(event)); + } +} + +void SearchMonitor::ListenToEvent(Solver::MonitorEvent event) { + solver()->searches_.back()->AddEventListener(event, this); } // ---------- Propagation Monitor ----------- diff --git a/ortools/constraint_solver/constraint_solver.h b/ortools/constraint_solver/constraint_solver.h index 2398e2387b..9ed0695176 100644 --- a/ortools/constraint_solver/constraint_solver.h +++ b/ortools/constraint_solver/constraint_solver.h @@ -156,6 +156,10 @@ struct StateInfo; struct Trail; template class SimpleRevFIFO; +template +class LightIntFunctionElementCt; +template +class LightIntIntFunctionElementCt; inline int64_t CpRandomSeed() { return absl::GetFlag(FLAGS_cp_random_seed) == -1 @@ -731,6 +735,37 @@ class Solver { /// Optimization directions. enum OptimizationDirection { NOT_SET, MAXIMIZATION, MINIMIZATION }; +#ifndef SWIG + /// Search monitor events. + enum class MonitorEvent : int { + kEnterSearch = 0, + kRestartSearch, + kExitSearch, + kBeginNextDecision, + kEndNextDecision, + kApplyDecision, + kRefuteDecision, + kAfterDecision, + kBeginFail, + kEndFail, + kBeginInitialPropagation, + kEndInitialPropagation, + kAcceptSolution, + kAtSolution, + kNoMoreSolutions, + kLocalOptimum, + kAcceptDelta, + kAcceptNeighbor, + kAcceptUncheckedNeighbor, + kIsUncheckedSolutionLimitReached, + kPeriodicCheck, + kProgressPercent, + kAccept, + // Dummy event whose underlying int is the number of MonitorEvent enums. + kLast, + }; +#endif // SWIG + /// Callback typedefs typedef std::function IndexEvaluator1; typedef std::function IndexEvaluator2; @@ -1154,6 +1189,39 @@ class Solver { int64_t range_end, IntVar* argument); #endif // SWIG + /// Light versions of function-based elements, in constraint version only, + /// well-suited for use within Local Search. + /// These constraints are "checking" constraints, only triggered on WhenBound + /// events. They provide very little (or no) domain filtering. + + /// Returns a light one-dimension function-based element constraint ensuring + /// var == values(index). + /// The constraint does not perform bound reduction of the resulting variable + /// until the index variable is bound. + /// If deep_serialize returns false, the model visitor will not extract all + /// possible values from the values function. + template + Constraint* MakeLightElement(F values, IntVar* const var, IntVar* const index, + std::function deep_serialize = nullptr) { + return RevAlloc(new LightIntFunctionElementCt( + this, var, index, std::move(values), std::move(deep_serialize))); + } + + /// Light two-dimension function-based element constraint ensuring + /// var == values(index1, index2). + /// The constraint does not perform bound reduction of the resulting variable + /// until the index variables are bound. + /// If deep_serialize returns false, the model visitor will not extract all + /// possible values from the values function. + template + Constraint* MakeLightElement(F values, IntVar* const var, + IntVar* const index1, IntVar* const index2, + std::function deep_serialize = nullptr) { + return RevAlloc(new LightIntIntFunctionElementCt( + this, var, index1, index2, std::move(values), + std::move(deep_serialize))); + } + /// Returns the expression expr such that vars[expr] == value. /// It assumes that vars are all different. IntExpr* MakeIndexExpression(const std::vector& vars, int64_t value); @@ -3399,6 +3467,7 @@ class ModelVisitor : public BaseObject { static const char kDivide[]; static const char kDurationExpr[]; static const char kElement[]; + static const char kLightElementEqual[]; static const char kElementEqual[]; static const char kEndExpr[]; static const char kEquality[]; @@ -3751,8 +3820,6 @@ class SearchMonitor : public BaseObject { /// unchecked solutions. virtual bool IsUncheckedSolutionLimitReached() { return false; } - Solver* solver() const { return solver_; } - /// Periodic call to check limits in long running methods. virtual void PeriodicCheck(); @@ -3764,9 +3831,15 @@ class SearchMonitor : public BaseObject { virtual void Accept(ModelVisitor* const visitor) const; /// Registers itself on the solver such that it gets notified of the search - /// and propagation events. + /// and propagation events. Override to incrementally install listeners for + /// specific events. virtual void Install(); + Solver* solver() const { return solver_; } + + protected: + void ListenToEvent(Solver::MonitorEvent event); + private: Solver* const solver_; DISALLOW_COPY_AND_ASSIGN(SearchMonitor); @@ -4150,6 +4223,7 @@ class SolutionCollector : public SearchMonitor { SolutionCollector(Solver* const solver, const Assignment* assignment); explicit SolutionCollector(Solver* const solver); ~SolutionCollector() override; + void Install() override; std::string DebugString() const override { return "SolutionCollector"; } /// Add API. @@ -4315,6 +4389,7 @@ class SearchLimit : public SearchMonitor { std::string DebugString() const override { return absl::StrFormat("SearchLimit(crossed = %i)", crossed_); } + void Install() override; private: void TopPeriodicCheck(); @@ -4351,6 +4426,7 @@ class RegularLimit : public SearchLimit { bool IsUncheckedSolutionLimitReached() override; int ProgressPercent() override; std::string DebugString() const override; + void Install() override; absl::Time AbsoluteSolverDeadline() const { return solver_time_at_limit_start_ + duration_limit_; @@ -4410,6 +4486,7 @@ class ImprovementSearchLimit : public SearchLimit { bool Check() override; bool AtSolution() override; void Init() override; + void Install() override; private: IntVar* objective_var_; @@ -5122,10 +5199,14 @@ class Assignment : public PropagationBaseObject { #endif // #if !defined(SWIG) void Save(AssignmentProto* const assignment_proto) const; - void AddObjective(IntVar* const v); + void AddObjective(IntVar* const v) { + // Objective can only set once. + DCHECK(!HasObjective()); + objective_element_.Reset(v); + } void ClearObjective() { objective_element_.Reset(nullptr); } - IntVar* Objective() const; - bool HasObjective() const { return (objective_element_.Var() != nullptr); } + IntVar* Objective() const { return objective_element_.Var(); } + bool HasObjective() const { return (Objective() != nullptr); } int64_t ObjectiveMin() const; int64_t ObjectiveMax() const; int64_t ObjectiveValue() const; diff --git a/ortools/constraint_solver/constraint_solveri.h b/ortools/constraint_solver/constraint_solveri.h index 3de82c1f15..97fd420fd1 100644 --- a/ortools/constraint_solver/constraint_solveri.h +++ b/ortools/constraint_solver/constraint_solveri.h @@ -770,6 +770,126 @@ Demon* MakeDelayedConstraintDemon2(Solver* const s, T* const ct, #endif // !defined(SWIG) +// ----- LightIntFunctionElementCt ----- + +template +class LightIntFunctionElementCt : public Constraint { + public: + LightIntFunctionElementCt(Solver* const solver, IntVar* const var, + IntVar* const index, F values, + std::function deep_serialize) + : Constraint(solver), + var_(var), + index_(index), + values_(std::move(values)), + deep_serialize_(std::move(deep_serialize)) {} + ~LightIntFunctionElementCt() override {} + + void Post() override { + Demon* demon = MakeConstraintDemon0( + solver(), this, &LightIntFunctionElementCt::IndexBound, "IndexBound"); + index_->WhenBound(demon); + } + + void InitialPropagate() override { + if (index_->Bound()) { + IndexBound(); + } + } + + std::string DebugString() const override { + return absl::StrFormat("LightIntFunctionElementCt(%s, %s)", + var_->DebugString(), index_->DebugString()); + } + + void Accept(ModelVisitor* const visitor) const override { + visitor->BeginVisitConstraint(ModelVisitor::kLightElementEqual, this); + visitor->VisitIntegerExpressionArgument(ModelVisitor::kTargetArgument, + var_); + visitor->VisitIntegerExpressionArgument(ModelVisitor::kIndexArgument, + index_); + // Warning: This will expand all values into a vector. + if (deep_serialize_ == nullptr || deep_serialize_()) { + visitor->VisitInt64ToInt64Extension(values_, index_->Min(), + index_->Max()); + } + visitor->EndVisitConstraint(ModelVisitor::kLightElementEqual, this); + } + + private: + void IndexBound() { var_->SetValue(values_(index_->Min())); } + + IntVar* const var_; + IntVar* const index_; + F values_; + std::function deep_serialize_; +}; + +// ----- LightIntIntFunctionElementCt ----- + +template +class LightIntIntFunctionElementCt : public Constraint { + public: + LightIntIntFunctionElementCt(Solver* const solver, IntVar* const var, + IntVar* const index1, IntVar* const index2, + F values, std::function deep_serialize) + : Constraint(solver), + var_(var), + index1_(index1), + index2_(index2), + values_(std::move(values)), + deep_serialize_(std::move(deep_serialize)) {} + ~LightIntIntFunctionElementCt() override {} + void Post() override { + Demon* demon = MakeConstraintDemon0( + solver(), this, &LightIntIntFunctionElementCt::IndexBound, + "IndexBound"); + index1_->WhenBound(demon); + index2_->WhenBound(demon); + } + void InitialPropagate() override { IndexBound(); } + + std::string DebugString() const override { + return "LightIntIntFunctionElementCt"; + } + + void Accept(ModelVisitor* const visitor) const override { + visitor->BeginVisitConstraint(ModelVisitor::kLightElementEqual, this); + visitor->VisitIntegerExpressionArgument(ModelVisitor::kTargetArgument, + var_); + visitor->VisitIntegerExpressionArgument(ModelVisitor::kIndexArgument, + index1_); + visitor->VisitIntegerExpressionArgument(ModelVisitor::kIndex2Argument, + index2_); + // Warning: This will expand all values into a vector. + const int64_t index1_min = index1_->Min(); + const int64_t index1_max = index1_->Max(); + visitor->VisitIntegerArgument(ModelVisitor::kMinArgument, index1_min); + visitor->VisitIntegerArgument(ModelVisitor::kMaxArgument, index1_max); + if (deep_serialize_ == nullptr || deep_serialize_()) { + for (int i = index1_min; i <= index1_max; ++i) { + visitor->VisitInt64ToInt64Extension( + [this, i](int64_t j) { return values_(i, j); }, index2_->Min(), + index2_->Max()); + } + } + visitor->EndVisitConstraint(ModelVisitor::kLightElementEqual, this); + } + + private: + void IndexBound() { + if (index1_->Bound() && index2_->Bound()) { + var_->SetValue(values_(index1_->Min(), index2_->Min())); + } + } + + IntVar* const var_; + IntVar* const index1_; + IntVar* const index2_; + Solver::IndexEvaluator2 values_; + std::function deep_serialize_; +}; + /// The base class for all local search operators. /// /// A local search operator is an object that defines the neighborhood of a @@ -803,146 +923,268 @@ class LocalSearchOperator : public BaseObject { virtual bool HoldsDelta() const { return false; } }; -/// Base operator class for operators manipulating variables. -template -class VarLocalSearchOperator : public LocalSearchOperator { +class LocalSearchOperatorState { public: - VarLocalSearchOperator() : activated_(), was_activated_(), cleared_(true) {} - explicit VarLocalSearchOperator(Handler var_handler) - : activated_(), - was_activated_(), - cleared_(true), - var_handler_(var_handler) {} - ~VarLocalSearchOperator() override {} + LocalSearchOperatorState() {} + + void SetCurrentDomainInjectiveAndKeepInverseValues(int max_value) { + max_inversible_index_ = candidate_values_.size(); + candidate_value_to_index_.resize(max_value + 1, -1); + committed_value_to_index_.resize(max_value + 1, -1); + } + + /// Returns the value in the current assignment of the variable of given + /// index. + int64_t CandidateValue(int64_t index) const { + DCHECK_LT(index, candidate_values_.size()); + return candidate_values_[index]; + } + int64_t CommittedValue(int64_t index) const { + return committed_values_[index]; + } + int64_t CheckPointValue(int64_t index) const { + return checkpoint_values_[index]; + } + void SetCandidateValue(int64_t index, int64_t value) { + candidate_values_[index] = value; + if (index < max_inversible_index_) { + candidate_value_to_index_[value] = index; + } + MarkChange(index); + } + + bool CandidateIsActive(int64_t index) const { + return candidate_is_active_[index]; + } + void SetCandidateActive(int64_t index, bool active) { + if (active) { + candidate_is_active_.Set(index); + } else { + candidate_is_active_.Clear(index); + } + MarkChange(index); + } + + void Commit() { + for (const int64_t index : changes_.PositionsSetAtLeastOnce()) { + const int64_t value = candidate_values_[index]; + committed_values_[index] = value; + if (index < max_inversible_index_) { + committed_value_to_index_[value] = index; + } + committed_is_active_.CopyBucket(candidate_is_active_, index); + } + changes_.SparseClearAll(); + incremental_changes_.SparseClearAll(); + } + + void CheckPoint() { checkpoint_values_ = committed_values_; } + + void Revert(bool only_incremental) { + incremental_changes_.SparseClearAll(); + if (only_incremental) return; + + for (const int64_t index : changes_.PositionsSetAtLeastOnce()) { + const int64_t committed_value = committed_values_[index]; + candidate_values_[index] = committed_value; + if (index < max_inversible_index_) { + candidate_value_to_index_[committed_value] = index; + } + candidate_is_active_.CopyBucket(committed_is_active_, index); + } + changes_.SparseClearAll(); + } + + const std::vector& CandidateIndicesChanged() const { + return changes_.PositionsSetAtLeastOnce(); + } + const std::vector& IncrementalIndicesChanged() const { + return incremental_changes_.PositionsSetAtLeastOnce(); + } + + void Resize(int size) { + candidate_values_.resize(size); + committed_values_.resize(size); + checkpoint_values_.resize(size); + candidate_is_active_.Resize(size); + committed_is_active_.Resize(size); + changes_.ClearAndResize(size); + incremental_changes_.ClearAndResize(size); + } + + int64_t CandidateInverseValue(int64_t value) const { + return candidate_value_to_index_[value]; + } + int64_t CommittedInverseValue(int64_t value) const { + return committed_value_to_index_[value]; + } + + private: + void MarkChange(int64_t index) { + incremental_changes_.Set(index); + changes_.Set(index); + } + + std::vector candidate_values_; + std::vector committed_values_; + std::vector checkpoint_values_; + + Bitset64<> candidate_is_active_; + Bitset64<> committed_is_active_; + + SparseBitset<> changes_; + SparseBitset<> incremental_changes_; + + int64_t max_inversible_index_ = -1; + std::vector candidate_value_to_index_; + std::vector committed_value_to_index_; +}; + +/// Specialization of LocalSearchOperator built from an array of IntVars +/// which specifies the scope of the operator. +/// This class also takes care of storing current variable values in Start(), +/// keeps track of changes done by the operator and builds the delta. +/// The Deactivate() method can be used to perform Large Neighborhood Search. +class IntVarLocalSearchOperator : public LocalSearchOperator { + public: + // If keep_inverse_values is true, assumes that vars models an injective + // function f with domain [0, vars.size()) in which case the operator will + // maintain the inverse function. + explicit IntVarLocalSearchOperator(const std::vector& vars, + bool keep_inverse_values = false) { + AddVars(vars); + if (keep_inverse_values) { + int64_t max_value = -1; + for (const IntVar* const var : vars) { + max_value = std::max(max_value, var->Max()); + } + state_.SetCurrentDomainInjectiveAndKeepInverseValues(max_value); + } + } + ~IntVarLocalSearchOperator() override {} + bool HoldsDelta() const override { return true; } /// This method should not be overridden. Override OnStart() instead which is /// called before exiting this method. void Start(const Assignment* assignment) override { + state_.CheckPoint(); + RevertChanges(false); const int size = Size(); CHECK_LE(size, assignment->Size()) << "Assignment contains fewer variables than operator"; + const Assignment::IntContainer& container = assignment->IntVarContainer(); for (int i = 0; i < size; ++i) { - activated_.Set(i, var_handler_.ValueFromAssignment(*assignment, vars_[i], - i, &values_[i])); + const IntVarElement* element = &(container.Element(i)); + if (element->Var() != vars_[i]) { + CHECK(container.Contains(vars_[i])) + << "Assignment does not contain operator variable " << vars_[i]; + element = &(container.Element(vars_[i])); + } + state_.SetCandidateValue(i, element->Value()); + state_.SetCandidateActive(i, element->Activated()); } - prev_values_ = old_values_; - old_values_ = values_; - was_activated_.SetContentFromBitsetOfSameSize(activated_); + state_.Commit(); OnStart(); } virtual bool IsIncremental() const { return false; } + int Size() const { return vars_.size(); } /// Returns the value in the current assignment of the variable of given /// index. - const Val& Value(int64_t index) const { + int64_t Value(int64_t index) const { DCHECK_LT(index, vars_.size()); - return values_[index]; + return state_.CandidateValue(index); } /// Returns the variable of given index. - V* Var(int64_t index) const { return vars_[index]; } + IntVar* Var(int64_t index) const { return vars_[index]; } virtual bool SkipUnchanged(int index) const { return false; } - const Val& OldValue(int64_t index) const { return old_values_[index]; } - void SetValue(int64_t index, const Val& value) { - values_[index] = value; - MarkChange(index); + int64_t OldValue(int64_t index) const { return state_.CommittedValue(index); } + int64_t PrevValue(int64_t index) const { + return state_.CheckPointValue(index); } - bool Activated(int64_t index) const { return activated_[index]; } - void Activate(int64_t index) { - activated_.Set(index); - MarkChange(index); + void SetValue(int64_t index, int64_t value) { + state_.SetCandidateValue(index, value); } - void Deactivate(int64_t index) { - activated_.Clear(index); - MarkChange(index); + bool Activated(int64_t index) const { + return state_.CandidateIsActive(index); } + void Activate(int64_t index) { state_.SetCandidateActive(index, true); } + void Deactivate(int64_t index) { state_.SetCandidateActive(index, false); } + bool ApplyChanges(Assignment* delta, Assignment* deltadelta) const { - if (IsIncremental() && !cleared_) { - for (const int64_t index : delta_changes_.PositionsSetAtLeastOnce()) { - V* var = Var(index); - const Val& value = Value(index); - const bool activated = activated_[index]; - var_handler_.AddToAssignment(var, value, activated, nullptr, index, - deltadelta); - var_handler_.AddToAssignment(var, value, activated, - &assignment_indices_, index, delta); + if (IsIncremental() && candidate_has_changes_) { + for (const int64_t index : state_.IncrementalIndicesChanged()) { + IntVar* var = Var(index); + const int64_t value = Value(index); + const bool activated = Activated(index); + AddToAssignment(var, value, activated, nullptr, index, deltadelta); + AddToAssignment(var, value, activated, &assignment_indices_, index, + delta); } } else { delta->Clear(); - for (const int64_t index : changes_.PositionsSetAtLeastOnce()) { - const Val& value = Value(index); - const bool activated = activated_[index]; + for (const int64_t index : state_.CandidateIndicesChanged()) { + const int64_t value = Value(index); + const bool activated = Activated(index); if (!activated || value != OldValue(index) || !SkipUnchanged(index)) { - var_handler_.AddToAssignment(Var(index), value, activated_[index], - &assignment_indices_, index, delta); + AddToAssignment(Var(index), value, activated, &assignment_indices_, + index, delta); } } } return true; } - void RevertChanges(bool incremental) { - cleared_ = false; - delta_changes_.SparseClearAll(); - if (incremental && IsIncremental()) return; - cleared_ = true; - for (const int64_t index : changes_.PositionsSetAtLeastOnce()) { - values_[index] = old_values_[index]; - var_handler_.OnRevertChanges(index, values_[index]); - activated_.CopyBucket(was_activated_, index); - assignment_indices_[index] = -1; + + void RevertChanges(bool change_was_incremental) { + candidate_has_changes_ = change_was_incremental && IsIncremental(); + + if (!candidate_has_changes_) { + for (const int64_t index : state_.CandidateIndicesChanged()) { + assignment_indices_[index] = -1; + } } - changes_.SparseClearAll(); + state_.Revert(candidate_has_changes_); } - void AddVars(const std::vector& vars) { + + void AddVars(const std::vector& vars) { if (!vars.empty()) { vars_.insert(vars_.end(), vars.begin(), vars.end()); const int64_t size = Size(); - values_.resize(size); - old_values_.resize(size); - prev_values_.resize(size); assignment_indices_.resize(size, -1); - activated_.Resize(size); - was_activated_.Resize(size); - changes_.ClearAndResize(size); - delta_changes_.ClearAndResize(size); - var_handler_.OnAddVars(); + state_.Resize(size); } } /// Called by Start() after synchronizing the operator with the current /// assignment. Should be overridden instead of Start() to avoid calling - /// VarLocalSearchOperator::Start explicitly. + /// IntVarLocalSearchOperator::Start explicitly. virtual void OnStart() {} /// OnStart() should really be protected, but then SWIG doesn't see it. So we /// make it public, but only subclasses should access to it (to override it). + + /// Redefines MakeNextNeighbor to export a simpler interface. The calls to + /// ApplyChanges() and RevertChanges() are factored in this method, hiding + /// both delta and deltadelta from subclasses which only need to override + /// MakeOneNeighbor(). + /// Therefore this method should not be overridden. Override MakeOneNeighbor() + /// instead. + bool MakeNextNeighbor(Assignment* delta, Assignment* deltadelta) override; + protected: - void MarkChange(int64_t index) { - delta_changes_.Set(index); - changes_.Set(index); + /// Creates a new neighbor. It returns false when the neighborhood is + /// completely explored. + // TODO(user): make it pure virtual, implies porting all apps overriding + /// MakeNextNeighbor() in a subclass of IntVarLocalSearchOperator. + virtual bool MakeOneNeighbor(); + + int64_t InverseValue(int64_t index) const { + return state_.CandidateInverseValue(index); + } + int64_t OldInverseValue(int64_t index) const { + return state_.CommittedInverseValue(index); } - std::vector vars_; - std::vector values_; - std::vector old_values_; - std::vector prev_values_; - mutable std::vector assignment_indices_; - Bitset64<> activated_; - Bitset64<> was_activated_; - SparseBitset<> changes_; - SparseBitset<> delta_changes_; - bool cleared_; - Handler var_handler_; -}; - -/// Base operator class for operators manipulating IntVars. -class IntVarLocalSearchOperator; - -class IntVarLocalSearchHandler { - public: - IntVarLocalSearchHandler() : op_(nullptr) {} - IntVarLocalSearchHandler(const IntVarLocalSearchHandler& other) - : op_(other.op_) {} - explicit IntVarLocalSearchHandler(IntVarLocalSearchOperator* op) : op_(op) {} void AddToAssignment(IntVar* var, int64_t value, bool active, std::vector* assignment_indices, int64_t index, Assignment* assignment) const { @@ -966,270 +1208,15 @@ class IntVarLocalSearchHandler { element->Deactivate(); } } - bool ValueFromAssignment(const Assignment& assignment, IntVar* var, - int64_t index, int64_t* value); - void OnRevertChanges(int64_t index, int64_t value); - void OnAddVars() {} private: - IntVarLocalSearchOperator* const op_; + std::vector vars_; + mutable std::vector assignment_indices_; + bool candidate_has_changes_ = false; + + LocalSearchOperatorState state_; }; -/// Specialization of LocalSearchOperator built from an array of IntVars -/// which specifies the scope of the operator. -/// This class also takes care of storing current variable values in Start(), -/// keeps track of changes done by the operator and builds the delta. -/// The Deactivate() method can be used to perform Large Neighborhood Search. - -#ifdef SWIG -/// Unfortunately, we must put this code here and not in -/// */constraint_solver.i, because it must be parsed by SWIG before the -/// derived C++ class. -// TODO(user): find a way to move this code back to the .i file, where it -/// belongs. -/// In python, we use an allow-list to expose the API. This list must also -/// be extended here. -#if defined(SWIGPYTHON) -// clang-format off -%unignore VarLocalSearchOperator::Size; -%unignore VarLocalSearchOperator::Value; -%unignore VarLocalSearchOperator::OldValue; -%unignore VarLocalSearchOperator::SetValue; -%feature("director") VarLocalSearchOperator::IsIncremental; -%feature("director") VarLocalSearchOperator::OnStart; -%unignore VarLocalSearchOperator::IsIncremental; -%unignore VarLocalSearchOperator::OnStart; -// clang-format on -#endif // SWIGPYTHON - -// clang-format off -%rename(IntVarLocalSearchOperatorTemplate) - VarLocalSearchOperator; -%template(IntVarLocalSearchOperatorTemplate) - VarLocalSearchOperator; -// clang-format on -#endif // SWIG - -class IntVarLocalSearchOperator - : public VarLocalSearchOperator { - public: - IntVarLocalSearchOperator() : max_inverse_value_(-1) {} - // If keep_inverse_values is true, assumes that vars models an injective - // function f with domain [0, vars.size()) in which case the operator will - // maintain the inverse function. - explicit IntVarLocalSearchOperator(const std::vector& vars, - bool keep_inverse_values = false) - : VarLocalSearchOperator( - IntVarLocalSearchHandler(this)), - max_inverse_value_(keep_inverse_values ? vars.size() - 1 : -1) { - AddVars(vars); - if (keep_inverse_values) { - int64_t max_value = -1; - for (const IntVar* const var : vars) { - max_value = std::max(max_value, var->Max()); - } - inverse_values_.resize(max_value + 1, -1); - old_inverse_values_.resize(max_value + 1, -1); - } - } - ~IntVarLocalSearchOperator() override {} - /// Redefines MakeNextNeighbor to export a simpler interface. The calls to - /// ApplyChanges() and RevertChanges() are factored in this method, hiding - /// both delta and deltadelta from subclasses which only need to override - /// MakeOneNeighbor(). - /// Therefore this method should not be overridden. Override MakeOneNeighbor() - /// instead. - bool MakeNextNeighbor(Assignment* delta, Assignment* deltadelta) override; - - protected: - friend class IntVarLocalSearchHandler; - - /// Creates a new neighbor. It returns false when the neighborhood is - /// completely explored. - // TODO(user): make it pure virtual, implies porting all apps overriding - /// MakeNextNeighbor() in a subclass of IntVarLocalSearchOperator. - virtual bool MakeOneNeighbor(); - - bool IsInverseValue(int64_t index) const { - DCHECK_GE(index, 0); - return index <= max_inverse_value_; - } - - int64_t InverseValue(int64_t index) const { return inverse_values_[index]; } - - int64_t OldInverseValue(int64_t index) const { - return old_inverse_values_[index]; - } - - void SetInverseValue(int64_t index, int64_t value) { - inverse_values_[index] = value; - } - - void SetOldInverseValue(int64_t index, int64_t value) { - old_inverse_values_[index] = value; - } - - private: - const int64_t max_inverse_value_; - std::vector old_inverse_values_; - std::vector inverse_values_; -}; - -inline bool IntVarLocalSearchHandler::ValueFromAssignment( - const Assignment& assignment, IntVar* var, int64_t index, int64_t* value) { - const Assignment::IntContainer& container = assignment.IntVarContainer(); - const IntVarElement* element = &(container.Element(index)); - if (element->Var() != var) { - CHECK(container.Contains(var)) - << "Assignment does not contain operator variable " << var; - element = &(container.Element(var)); - } - *value = element->Value(); - if (op_->IsInverseValue(index)) { - op_->SetInverseValue(*value, index); - op_->SetOldInverseValue(*value, index); - } - return element->Activated(); -} - -inline void IntVarLocalSearchHandler::OnRevertChanges(int64_t index, - int64_t value) { - if (op_->IsInverseValue(index)) { - op_->SetInverseValue(value, index); - } -} - -/// SequenceVarLocalSearchOperator -class SequenceVarLocalSearchOperator; - -class SequenceVarLocalSearchHandler { - public: - SequenceVarLocalSearchHandler() : op_(nullptr) {} - SequenceVarLocalSearchHandler(const SequenceVarLocalSearchHandler& other) - : op_(other.op_) {} - explicit SequenceVarLocalSearchHandler(SequenceVarLocalSearchOperator* op) - : op_(op) {} - void AddToAssignment(SequenceVar* var, const std::vector& value, - bool active, std::vector* assignment_indices, - int64_t index, Assignment* assignment) const; - bool ValueFromAssignment(const Assignment& assignment, SequenceVar* var, - int64_t index, std::vector* value); - void OnRevertChanges(int64_t index, const std::vector& value); - void OnAddVars(); - - private: - SequenceVarLocalSearchOperator* const op_; -}; - -#ifdef SWIG -/// Unfortunately, we must put this code here and not in -/// */constraint_solver.i, because it must be parsed by SWIG before the -/// derived C++ class. -// TODO(user): find a way to move this code back to the .i file, where it -/// belongs. -// clang-format off -%rename(SequenceVarLocalSearchOperatorTemplate) VarLocalSearchOperator< - SequenceVar, std::vector, SequenceVarLocalSearchHandler>; -%template(SequenceVarLocalSearchOperatorTemplate) VarLocalSearchOperator< - SequenceVar, std::vector, SequenceVarLocalSearchHandler>; -// clang-format on -#endif - -typedef VarLocalSearchOperator, - SequenceVarLocalSearchHandler> - SequenceVarLocalSearchOperatorTemplate; - -class SequenceVarLocalSearchOperator - : public SequenceVarLocalSearchOperatorTemplate { - public: - SequenceVarLocalSearchOperator() {} - explicit SequenceVarLocalSearchOperator(const std::vector& vars) - : SequenceVarLocalSearchOperatorTemplate( - SequenceVarLocalSearchHandler(this)) { - AddVars(vars); - } - ~SequenceVarLocalSearchOperator() override {} - /// Returns the value in the current assignment of the variable of given - /// index. - const std::vector& Sequence(int64_t index) const { return Value(index); } - const std::vector& OldSequence(int64_t index) const { - return OldValue(index); - } - void SetForwardSequence(int64_t index, const std::vector& value) { - SetValue(index, value); - } - void SetBackwardSequence(int64_t index, const std::vector& value) { - backward_values_[index] = value; - MarkChange(index); - } - - protected: - friend class SequenceVarLocalSearchHandler; - - std::vector> backward_values_; -}; - -inline void SequenceVarLocalSearchHandler::AddToAssignment( - SequenceVar* var, const std::vector& value, bool active, - std::vector* assignment_indices, int64_t index, - Assignment* assignment) const { - Assignment::SequenceContainer* const container = - assignment->MutableSequenceVarContainer(); - SequenceVarElement* element = nullptr; - if (assignment_indices != nullptr) { - if ((*assignment_indices)[index] == -1) { - (*assignment_indices)[index] = container->Size(); - element = assignment->FastAdd(var); - } else { - element = container->MutableElement((*assignment_indices)[index]); - } - } else { - element = assignment->FastAdd(var); - } - if (active) { - element->SetForwardSequence(value); - element->SetBackwardSequence(op_->backward_values_[index]); - element->Activate(); - } else { - element->Deactivate(); - } -} - -inline bool SequenceVarLocalSearchHandler::ValueFromAssignment( - const Assignment& assignment, SequenceVar* var, int64_t index, - std::vector* value) { - const Assignment::SequenceContainer& container = - assignment.SequenceVarContainer(); - const SequenceVarElement* element = &(container.Element(index)); - if (element->Var() != var) { - CHECK(container.Contains(var)) - << "Assignment does not contain operator variable " << var; - element = &(container.Element(var)); - } - const std::vector& element_value = element->ForwardSequence(); - CHECK_GE(var->size(), element_value.size()); - op_->backward_values_[index].clear(); - *value = element_value; - return element->Activated(); -} - -inline void SequenceVarLocalSearchHandler::OnRevertChanges( - int64_t index, const std::vector& value) { - op_->backward_values_[index].clear(); -} - -inline void SequenceVarLocalSearchHandler::OnAddVars() { - op_->backward_values_.resize(op_->Size()); -} - /// This is the base class for building an Lns operator. An Lns fragment is a /// collection of variables which will be relaxed. Fragments are built with /// NextFragment(), which returns false if there are no more fragments to build. @@ -1410,6 +1397,8 @@ class PathOperator : public IntVarLocalSearchOperator { } /// Returns the start node of the ith base node. int64_t StartNode(int i) const { return path_starts_[base_paths_[i]]; } + /// Returns the end node of the ith base node. + int64_t EndNode(int i) const { return path_ends_[base_paths_[i]]; } /// Returns the vector of path start nodes. const std::vector& path_starts() const { return path_starts_; } /// Returns the class of the path of the ith base node. @@ -1456,6 +1445,11 @@ class PathOperator : public IntVarLocalSearchOperator { return OldValue(node); } + int64_t PrevNext(int64_t node) const { + DCHECK(!IsPathEnd(node)); + return PrevValue(node); + } + int64_t OldPrev(int64_t node) const { DCHECK(!IsPathStart(node)); return OldInverseValue(node); @@ -1486,7 +1480,6 @@ class PathOperator : public IntVarLocalSearchOperator { void SetNext(int64_t from, int64_t to, int64_t path) { DCHECK_LT(from, number_of_nexts_); SetValue(from, to); - SetInverseValue(to, from); if (!ignore_path_vars_) { DCHECK_LT(from + number_of_nexts_, Size()); SetValue(from + number_of_nexts_, path); @@ -1602,6 +1595,7 @@ class PathOperator : public IntVarLocalSearchOperator { std::vector end_nodes_; std::vector base_paths_; std::vector path_starts_; + std::vector path_ends_; std::vector inactives_; bool just_started_; bool first_start_; @@ -3376,13 +3370,15 @@ class DimensionChecker { int64_t min; int64_t max; }; - + // TODO(user): the addition of kMinRangeSizeForRIQ slowed down Check(). + // See if using a template parameter makes it faster. DimensionChecker(const PathState* path_state, std::vector path_capacity, std::vector path_class, - std::vector> - min_max_demand_per_path_class, - std::vector node_capacity); + std::vector> + demand_per_path_class, + std::vector node_capacity, + int min_range_size_for_riq = kOptimalMinRangeSizeForRIQ); // Given the change made in PathState, checks that the dimension // constraint is still feasible. @@ -3392,6 +3388,8 @@ class DimensionChecker { // must be called before PathState::Commit(). void Commit(); + static constexpr int kOptimalMinRangeSizeForRIQ = 4; + private: // Range intersection query on backwards_demand_sums_. // The first_node and last_node MUST form a subpath in the committed state. @@ -3426,8 +3424,8 @@ class DimensionChecker { const PathState* const path_state_; const std::vector path_capacity_; const std::vector path_class_; - const std::vector> - min_max_demand_per_path_class_; + const std::vector> + demand_per_path_class_; std::vector cached_demand_; const std::vector node_capacity_; @@ -3449,6 +3447,8 @@ class DimensionChecker { // The incremental branch of Commit() may waste space in the layers of the // RIQ structure. This is the upper limit of a layer's size. const int maximum_riq_layer_size_; + // Range queries are used on a chain only if the range is larger than this. + const int min_range_size_for_riq_; // previous_nontrivial_index_[index_[node]] has the index of the previous // node on its committed path that has nonfixed demand or nontrivial node // capacity. This allows for O(1) queries that all nodes on a subpath diff --git a/ortools/constraint_solver/csharp/constraint_solver.i b/ortools/constraint_solver/csharp/constraint_solver.i index 83a66ccaeb..9c2103b51a 100644 --- a/ortools/constraint_solver/csharp/constraint_solver.i +++ b/ortools/constraint_solver/csharp/constraint_solver.i @@ -710,12 +710,6 @@ namespace operations_research { %unignore SearchLog::Maintain; %unignore SearchLog::OutputDecision; -// IntVarLocalSearchHandler -%ignore IntVarLocalSearchHandler; - -// SequenceVarLocalSearchHandler -%ignore SequenceVarLocalSearchHandler; - // LocalSearchOperator %feature("director") LocalSearchOperator; %unignore LocalSearchOperator; @@ -724,26 +718,6 @@ namespace operations_research { %unignore LocalSearchOperator::Reset; %unignore LocalSearchOperator::Start; -// VarLocalSearchOperator<> -%unignore VarLocalSearchOperator; -// Ignored: -%ignore VarLocalSearchOperator::Start; -%ignore VarLocalSearchOperator::ApplyChanges; -%ignore VarLocalSearchOperator::RevertChanges; -%ignore VarLocalSearchOperator::SkipUnchanged; -// Methods: -%unignore VarLocalSearchOperator::Size; -%unignore VarLocalSearchOperator::Value; -%unignore VarLocalSearchOperator::IsIncremental; -%unignore VarLocalSearchOperator::OnStart; -%unignore VarLocalSearchOperator::OldValue; -%unignore VarLocalSearchOperator::SetValue; -%unignore VarLocalSearchOperator::Var; -%unignore VarLocalSearchOperator::Activated; -%unignore VarLocalSearchOperator::Activate; -%unignore VarLocalSearchOperator::Deactivate; -%unignore VarLocalSearchOperator::AddVars; - // IntVarLocalSearchOperator %feature("director") IntVarLocalSearchOperator; %unignore IntVarLocalSearchOperator; @@ -782,17 +756,6 @@ namespace operations_research { // Methods: %unignore ChangeValue::ModifyValue; -// SequenceVarLocalSearchOperator -%feature("director") SequenceVarLocalSearchOperator; -%unignore SequenceVarLocalSearchOperator; -// Ignored: -%ignore SequenceVarLocalSearchOperator::SetBackwardSequence; -%ignore SequenceVarLocalSearchOperator::SetForwardSequence; -// Methods: -%unignore SequenceVarLocalSearchOperator::OldSequence; -%unignore SequenceVarLocalSearchOperator::Sequence; -%unignore SequenceVarLocalSearchOperator::Start; - // PathOperator %feature("director") PathOperator; %unignore PathOperator; diff --git a/ortools/constraint_solver/element.cc b/ortools/constraint_solver/element.cc index e5afa249d1..9275e13aa5 100644 --- a/ortools/constraint_solver/element.cc +++ b/ortools/constraint_solver/element.cc @@ -89,6 +89,27 @@ class BaseIntExprElement : public BaseIntExpr { private: void UpdateSupports() const; + template + void UpdateElementIndexBounds(T check_value) { + const int64_t emin = ExprMin(); + const int64_t emax = ExprMax(); + int64_t nmin = emin; + int64_t value = ElementValue(nmin); + while (nmin < emax && check_value(value)) { + nmin++; + value = ElementValue(nmin); + } + if (nmin == emax && check_value(value)) { + solver()->Fail(); + } + int64_t nmax = emax; + value = ElementValue(nmax); + while (nmax >= nmin && check_value(value)) { + nmax--; + value = ElementValue(nmax); + } + expr_->SetRange(nmin, nmax); + } mutable int64_t min_; mutable int min_support_; @@ -127,43 +148,22 @@ void BaseIntExprElement::Range(int64_t* mi, int64_t* ma) { *ma = max_; } -#define UPDATE_BASE_ELEMENT_INDEX_BOUNDS(test) \ - const int64_t emin = ExprMin(); \ - const int64_t emax = ExprMax(); \ - int64_t nmin = emin; \ - int64_t value = ElementValue(nmin); \ - while (nmin < emax && test) { \ - nmin++; \ - value = ElementValue(nmin); \ - } \ - if (nmin == emax && test) { \ - solver()->Fail(); \ - } \ - int64_t nmax = emax; \ - value = ElementValue(nmax); \ - while (nmax >= nmin && test) { \ - nmax--; \ - value = ElementValue(nmax); \ - } \ - expr_->SetRange(nmin, nmax); - void BaseIntExprElement::SetMin(int64_t m) { - UPDATE_BASE_ELEMENT_INDEX_BOUNDS(value < m); + UpdateElementIndexBounds([m](int64_t value) { return value < m; }); } void BaseIntExprElement::SetMax(int64_t m) { - UPDATE_BASE_ELEMENT_INDEX_BOUNDS(value > m); + UpdateElementIndexBounds([m](int64_t value) { return value > m; }); } void BaseIntExprElement::SetRange(int64_t mi, int64_t ma) { if (mi > ma) { solver()->Fail(); } - UPDATE_BASE_ELEMENT_INDEX_BOUNDS((value < mi || value > ma)); + UpdateElementIndexBounds( + [mi, ma](int64_t value) { return value < mi || value > ma; }); } -#undef UPDATE_BASE_ELEMENT_INDEX_BOUNDS - void BaseIntExprElement::UpdateSupports() const { if (initial_update_ || !expr_->Contains(min_support_) || !expr_->Contains(max_support_)) { diff --git a/ortools/constraint_solver/expr_cst.cc b/ortools/constraint_solver/expr_cst.cc index 81e63301ff..f01c2b3e73 100644 --- a/ortools/constraint_solver/expr_cst.cc +++ b/ortools/constraint_solver/expr_cst.cc @@ -360,9 +360,9 @@ void DiffCst::BoundPropagate() { if (var_min > value_ || var_max < value_) { demon_->inhibit(solver()); } else if (var_min == value_) { - var_->SetMin(value_ + 1); + var_->SetMin(CapAdd(value_, 1)); } else if (var_max == value_) { - var_->SetMax(value_ - 1); + var_->SetMax(CapSub(value_, 1)); } else if (!HasLargeDomain(var_)) { demon_->inhibit(solver()); var_->RemoveValue(value_); diff --git a/ortools/constraint_solver/expressions.cc b/ortools/constraint_solver/expressions.cc index 9d475ac21b..234d76850b 100644 --- a/ortools/constraint_solver/expressions.cc +++ b/ortools/constraint_solver/expressions.cc @@ -2466,15 +2466,18 @@ void DomainIntVar::Process() { } } -#define COND_REV_ALLOC(rev, alloc) rev ? solver()->RevAlloc(alloc) : alloc; +template +T* CondRevAlloc(Solver* solver, bool reversible, T* object) { + return reversible ? solver->RevAlloc(object) : object; +} IntVarIterator* DomainIntVar::MakeHoleIterator(bool reversible) const { - return COND_REV_ALLOC(reversible, new DomainIntVarHoleIterator(this)); + return CondRevAlloc(solver(), reversible, new DomainIntVarHoleIterator(this)); } IntVarIterator* DomainIntVar::MakeDomainIterator(bool reversible) const { - return COND_REV_ALLOC(reversible, - new DomainIntVarDomainIterator(this, reversible)); + return CondRevAlloc(solver(), reversible, + new DomainIntVarDomainIterator(this, reversible)); } std::string DomainIntVar::DebugString() const { @@ -2605,10 +2608,10 @@ class IntConst : public IntVar { uint64_t Size() const override { return 1; } bool Contains(int64_t v) const override { return (v == value_); } IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new EmptyIterator()); + return CondRevAlloc(solver(), reversible, new EmptyIterator()); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new RangeIterator(this)); + return CondRevAlloc(solver(), reversible, new RangeIterator(this)); } int64_t OldMin() const override { return value_; } int64_t OldMax() const override { return value_; } @@ -2768,12 +2771,14 @@ class PlusCstIntVar : public PlusCstVar { bool Contains(int64_t v) const override { return var_->Contains(v - cst_); } IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC( - reversible, new PlusCstIntVarIterator(var_, cst_, true, reversible)); + return CondRevAlloc( + solver(), reversible, + new PlusCstIntVarIterator(var_, cst_, true, reversible)); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC( - reversible, new PlusCstIntVarIterator(var_, cst_, false, reversible)); + return CondRevAlloc( + solver(), reversible, + new PlusCstIntVarIterator(var_, cst_, false, reversible)); } }; @@ -2816,12 +2821,14 @@ class PlusCstDomainIntVar : public PlusCstVar { } IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new PlusCstDomainIntVarIterator( - var_, cst_, true, reversible)); + return CondRevAlloc( + solver(), reversible, + new PlusCstDomainIntVarIterator(var_, cst_, true, reversible)); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new PlusCstDomainIntVarIterator( - var_, cst_, false, reversible)); + return CondRevAlloc( + solver(), reversible, + new PlusCstDomainIntVarIterator(var_, cst_, false, reversible)); } }; @@ -2910,12 +2917,13 @@ class SubCstIntVar : public IntVar { void WhenBound(Demon* d) override; void WhenDomain(Demon* d) override; IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC( - reversible, new SubCstIntVarIterator(var_, cst_, true, reversible)); + return CondRevAlloc(solver(), reversible, + new SubCstIntVarIterator(var_, cst_, true, reversible)); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC( - reversible, new SubCstIntVarIterator(var_, cst_, false, reversible)); + return CondRevAlloc( + solver(), reversible, + new SubCstIntVarIterator(var_, cst_, false, reversible)); } int64_t OldMin() const override { return CapSub(cst_, var_->OldMax()); } int64_t OldMax() const override { return CapSub(cst_, var_->OldMin()); } @@ -3043,12 +3051,12 @@ class OppIntVar : public IntVar { void WhenBound(Demon* d) override; void WhenDomain(Demon* d) override; IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, - new OppIntVarIterator(var_, true, reversible)); + return CondRevAlloc(solver(), reversible, + new OppIntVarIterator(var_, true, reversible)); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, - new OppIntVarIterator(var_, false, reversible)); + return CondRevAlloc(solver(), reversible, + new OppIntVarIterator(var_, false, reversible)); } int64_t OldMin() const override { return CapOpp(var_->OldMax()); } int64_t OldMax() const override { return CapOpp(var_->OldMin()); } @@ -3219,12 +3227,14 @@ class TimesPosCstIntVar : public TimesCstIntVar { void WhenBound(Demon* d) override; void WhenDomain(Demon* d) override; IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new TimesPosCstIntVarIterator( - var_, cst_, true, reversible)); + return CondRevAlloc( + solver(), reversible, + new TimesPosCstIntVarIterator(var_, cst_, true, reversible)); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new TimesPosCstIntVarIterator( - var_, cst_, false, reversible)); + return CondRevAlloc( + solver(), reversible, + new TimesPosCstIntVarIterator(var_, cst_, false, reversible)); } int64_t OldMin() const override { return CapProd(var_->OldMin(), cst_); } int64_t OldMax() const override { return CapProd(var_->OldMax(), cst_); } @@ -3332,11 +3342,11 @@ class TimesPosCstBoolVar : public TimesCstIntVar { void WhenBound(Demon* d) override; void WhenDomain(Demon* d) override; IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new EmptyIterator()); + return CondRevAlloc(solver(), reversible, new EmptyIterator()); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC( - reversible, + return CondRevAlloc( + solver(), reversible, new TimesPosCstBoolVarIterator(boolean_var(), cst_, false, reversible)); } int64_t OldMin() const override { return 0; } @@ -3483,12 +3493,14 @@ class TimesNegCstIntVar : public TimesCstIntVar { void WhenBound(Demon* d) override; void WhenDomain(Demon* d) override; IntVarIterator* MakeHoleIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new TimesNegCstIntVarIterator( - var_, cst_, true, reversible)); + return CondRevAlloc( + solver(), reversible, + new TimesNegCstIntVarIterator(var_, cst_, true, reversible)); } IntVarIterator* MakeDomainIterator(bool reversible) const override { - return COND_REV_ALLOC(reversible, new TimesNegCstIntVarIterator( - var_, cst_, false, reversible)); + return CondRevAlloc( + solver(), reversible, + new TimesNegCstIntVarIterator(var_, cst_, false, reversible)); } int64_t OldMin() const override { return CapProd(var_->OldMax(), cst_); } int64_t OldMax() const override { return CapProd(var_->OldMin(), cst_); } @@ -3518,7 +3530,8 @@ void TimesNegCstIntVar::SetMax(int64_t m) { } void TimesNegCstIntVar::SetRange(int64_t l, int64_t u) { - var_->SetRange(PosIntDivUp(-u, -cst_), PosIntDivDown(-l, -cst_)); + var_->SetRange(PosIntDivUp(CapOpp(u), CapOpp(cst_)), + PosIntDivDown(CapOpp(l), CapOpp(cst_))); } void TimesNegCstIntVar::SetValue(int64_t v) { @@ -6348,10 +6361,10 @@ class LinkExprAndDomainIntVar : public CastConstraint { // ----- Misc ----- IntVarIterator* BooleanVar::MakeHoleIterator(bool reversible) const { - return COND_REV_ALLOC(reversible, new EmptyIterator()); + return CondRevAlloc(solver(), reversible, new EmptyIterator()); } IntVarIterator* BooleanVar::MakeDomainIterator(bool reversible) const { - return COND_REV_ALLOC(reversible, new RangeIterator(this)); + return CondRevAlloc(solver(), reversible, new RangeIterator(this)); } // ----- API ----- @@ -7385,7 +7398,7 @@ void IntVar::SetValues(const std::vector& values) { // that uses a global, static shared (and locked) storage. // TODO(user): [optional] consider porting // STLSortAndRemoveDuplicates from ortools/base/stl_util.h to the - // existing base/stl_util.h and using it here. + // existing open_source/base/stl_util.h and using it here. // TODO(user): We could filter out values not in the var. std::vector& tmp = solver()->tmp_vector_; tmp.clear(); @@ -7509,6 +7522,4 @@ bool Solver::IsProduct(IntExpr* const expr, IntExpr** inner_expr, return false; } -#undef COND_REV_ALLOC - } // namespace operations_research diff --git a/ortools/constraint_solver/java/constraint_solver.i b/ortools/constraint_solver/java/constraint_solver.i index 29e90f2745..241ec13c89 100644 --- a/ortools/constraint_solver/java/constraint_solver.i +++ b/ortools/constraint_solver/java/constraint_solver.i @@ -1400,20 +1400,6 @@ import java.util.function.Supplier; %rename (setStartRange) PropagationMonitor::SetStartRange; %rename (startProcessingIntegerVariable) PropagationMonitor::StartProcessingIntegerVariable; -// IntVarLocalSearchHandler -%unignore IntVarLocalSearchHandler; -%rename (addToAssignment) IntVarLocalSearchHandler::AddToAssignment; -%rename (onAddVars) IntVarLocalSearchHandler::OnAddVars; -%rename (onRevertChanges) IntVarLocalSearchHandler::OnRevertChanges; -%rename (valueFromAssignent) IntVarLocalSearchHandler::ValueFromAssignent; - -// SequenceVarLocalSearchHandler -%unignore SequenceVarLocalSearchHandler; -%rename (addToAssignment) SequenceVarLocalSearchHandler::AddToAssignment; -%rename (onAddVars) SequenceVarLocalSearchHandler::OnAddVars; -%rename (onRevertChanges) SequenceVarLocalSearchHandler::OnRevertChanges; -%rename (valueFromAssignent) SequenceVarLocalSearchHandler::ValueFromAssignent; - // LocalSearchOperator %feature("director") LocalSearchOperator; %unignore LocalSearchOperator; @@ -1421,24 +1407,6 @@ import java.util.function.Supplier; %rename (reset) LocalSearchOperator::Reset; %rename (start) LocalSearchOperator::Start; -// VarLocalSearchOperator<> -%unignore VarLocalSearchOperator; -%ignore VarLocalSearchOperator::Start; -%ignore VarLocalSearchOperator::ApplyChanges; -%ignore VarLocalSearchOperator::RevertChanges; -%ignore VarLocalSearchOperator::SkipUnchanged; -%rename (size) VarLocalSearchOperator::Size; -%rename (value) VarLocalSearchOperator::Value; -%rename (isIncremental) VarLocalSearchOperator::IsIncremental; -%rename (onStart) VarLocalSearchOperator::OnStart; -%rename (oldValue) VarLocalSearchOperator::OldValue; -%rename (setValue) VarLocalSearchOperator::SetValue; -%rename (var) VarLocalSearchOperator::Var; -%rename (activated) VarLocalSearchOperator::Activated; -%rename (activate) VarLocalSearchOperator::Activate; -%rename (deactivate) VarLocalSearchOperator::Deactivate; -%rename (addVars) VarLocalSearchOperator::AddVars; - // IntVarLocalSearchOperator %feature("director") IntVarLocalSearchOperator; %unignore IntVarLocalSearchOperator; @@ -1473,15 +1441,6 @@ import java.util.function.Supplier; %unignore ChangeValue; %rename (modifyValue) ChangeValue::ModifyValue; -// SequenceVarLocalSearchOperator -%feature("director") SequenceVarLocalSearchOperator; -%unignore SequenceVarLocalSearchOperator; -%ignore SequenceVarLocalSearchOperator::OldSequence; -%ignore SequenceVarLocalSearchOperator::Sequence; -%ignore SequenceVarLocalSearchOperator::SetBackwardSequence; -%ignore SequenceVarLocalSearchOperator::SetForwardSequence; -%rename (start) SequenceVarLocalSearchOperator::Start; - // PathOperator %feature("director") PathOperator; %unignore PathOperator; diff --git a/ortools/constraint_solver/local_search.cc b/ortools/constraint_solver/local_search.cc index af9c639a3a..6c8ad971d4 100644 --- a/ortools/constraint_solver/local_search.cc +++ b/ortools/constraint_solver/local_search.cc @@ -666,7 +666,7 @@ void PathOperator::InitializePathStarts() { if (!has_prevs[i]) { int current = i; while (!IsPathEnd(current)) { - if ((OldNext(current) != prev_values_[current])) { + if ((OldNext(current) != PrevNext(current))) { for (int j = 0; j < num_paths_; ++j) { optimal_paths_[num_paths_ * start_to_path_[i] + j] = false; optimal_paths_[num_paths_ * j + start_to_path_[i]] = false; @@ -753,6 +753,15 @@ void PathOperator::InitializePathStarts() { } } path_starts_.swap(new_path_starts); + // For every base path, store the end corresponding to the path start. + // TODO(user): make this faster, maybe by pairing starts with ends. + path_ends_.clear(); + path_ends_.reserve(path_starts_.size()); + for (const int start_node : path_starts_) { + int64_t node = start_node; + while (!IsPathEnd(node)) node = OldNext(node); + path_ends_.push_back(node); + } } void PathOperator::InitializeInactives() { @@ -1062,17 +1071,39 @@ bool Cross::MakeNeighbor() { if (node0 == start0) return false; const int64_t node1 = BaseNode(1); if (node1 == start1) return false; - if (!IsPathEnd(node0) && !IsPathEnd(node1)) { - // If two paths are equivalent don't exchange them. - if (PathClass(0) == PathClass(1) && IsPathEnd(Next(node0)) && - IsPathEnd(Next(node1))) { + + bool moved = false; + if (start0 < start1) { + // Cross path starts. + // If two paths are equivalent don't exchange the full paths. + if (PathClass(0) == PathClass(1) && !IsPathEnd(node0) && + IsPathEnd(Next(node0)) && !IsPathEnd(node1) && IsPathEnd(Next(node1))) { return false; } - return MoveChain(start0, node0, start1) && MoveChain(node0, node1, start0); + + const int first1 = Next(start1); + if (!IsPathEnd(node0)) moved |= MoveChain(start0, node0, start1); + if (!IsPathEnd(node1)) moved |= MoveChain(Prev(first1), node1, start0); + } else { // start1 > start0. + // Cross path ends. + // If paths are equivalent, every end crossing has a corresponding start + // crossing, we don't generate those symmetric neighbors. + if (PathClass(0) == PathClass(1)) return false; + // Never exchange full paths, equivalent or not. + // Full path exchange is only performed when start0 < start1. + if (IsPathStart(Prev(node0)) && IsPathStart(Prev(node1))) { + return false; + } + + const int prev_end_node1 = Prev(EndNode(1)); + if (!IsPathEnd(node0)) { + moved |= MoveChain(Prev(node0), Prev(EndNode(0)), prev_end_node1); + } + if (!IsPathEnd(node1)) { + moved |= MoveChain(Prev(node1), prev_end_node1, Prev(EndNode(0))); + } } - if (!IsPathEnd(node0)) return MoveChain(start0, node0, start1); - if (!IsPathEnd(node1)) return MoveChain(start1, node1, start0); - return false; + return moved; } // ----- BaseInactiveNodeToPathOperator ----- @@ -2472,66 +2503,6 @@ LocalSearchOperator* Solver::MakeOperator( } namespace { -// Classes for Local Search Operation used in Local Search filters. - -class SumOperation { - public: - SumOperation() : value_(0) {} - void Init() { value_ = 0; } - void Update(int64_t update) { value_ = CapAdd(value_, update); } - void Remove(int64_t remove) { value_ = CapSub(value_, remove); } - int64_t value() const { return value_; } - void set_value(int64_t new_value) { value_ = new_value; } - - private: - int64_t value_; -}; - -class ProductOperation { - public: - ProductOperation() : value_(1) {} - void Init() { value_ = 1; } - void Update(int64_t update) { value_ *= update; } - void Remove(int64_t remove) { - if (remove != 0) { - value_ /= remove; - } - } - int64_t value() const { return value_; } - void set_value(int64_t new_value) { value_ = new_value; } - - private: - int64_t value_; -}; - -class MinOperation { - public: - void Init() { values_set_.clear(); } - void Update(int64_t update) { values_set_.insert(update); } - void Remove(int64_t remove) { values_set_.erase(remove); } - int64_t value() const { - return (!values_set_.empty()) ? *values_set_.begin() : 0; - } - void set_value(int64_t new_value) {} - - private: - std::set values_set_; -}; - -class MaxOperation { - public: - void Init() { values_set_.clear(); } - void Update(int64_t update) { values_set_.insert(update); } - void Remove(int64_t remove) { values_set_.erase(remove); } - int64_t value() const { - return (!values_set_.empty()) ? *values_set_.rbegin() : 0; - } - void set_value(int64_t new_value) {} - - private: - std::set values_set_; -}; - // Always accepts deltas, cost 0. class AcceptFilter : public LocalSearchFilter { public: @@ -3035,17 +3006,18 @@ Interval Delta(const Interval& from, const Interval& to) { DimensionChecker::DimensionChecker( const PathState* path_state, std::vector path_capacity, std::vector path_class, - std::vector> min_max_demand_per_path_class, - std::vector node_capacity) + std::vector> + demand_per_path_class, + std::vector node_capacity, int min_range_size_for_riq) : path_state_(path_state), path_capacity_(std::move(path_capacity)), path_class_(std::move(path_class)), - min_max_demand_per_path_class_(std::move(min_max_demand_per_path_class)), + demand_per_path_class_(std::move(demand_per_path_class)), node_capacity_(std::move(node_capacity)), index_(path_state_->NumNodes(), 0), - maximum_riq_layer_size_( - std::max(16, 4 * path_state_->NumNodes())) // 16 and 4 are arbitrary. -{ + maximum_riq_layer_size_(std::max( + 16, 4 * path_state_->NumNodes())), // 16 and 4 are arbitrary. + min_range_size_for_riq_(min_range_size_for_riq) { const int num_nodes = path_state_->NumNodes(); cached_demand_.resize(num_nodes); const int num_paths = path_state_->NumPaths(); @@ -3062,8 +3034,10 @@ bool DimensionChecker::Check() const { for (const int path : path_state_->ChangedPaths()) { const Interval path_capacity = path_capacity_[path]; const int path_class = path_class_[path]; - // Loop invariant: cumul is nonempty and within path_capacity. - Interval cumul = node_capacity_[path_state_->Start(path)]; + // Loop invariant: except for the first chain, cumul represents the cumul + // state of the last node of the previous chain, and it is nonempty. + int prev_node = path_state_->Start(path); + Interval cumul = node_capacity_[prev_node]; cumul = Intersect(cumul, path_capacity); if (IsEmpty(cumul)) return false; @@ -3071,21 +3045,33 @@ bool DimensionChecker::Check() const { const int first_node = chain.First(); const int last_node = chain.Last(); + if (prev_node != first_node) { + // Bring cumul state from last node of previous chain to first node of + // current chain. + const Interval demand = + demand_per_path_class_[path_class](prev_node, first_node); + cumul += demand; + cumul = Intersect(cumul, path_capacity); + cumul = Intersect(cumul, node_capacity_[first_node]); + if (IsEmpty(cumul)) return false; + prev_node = first_node; + } + + // Bring cumul state from first node to last node of the current chain. const int first_index = index_[first_node]; const int last_index = index_[last_node]; - const int chain_path = path_state_->Path(first_node); const int chain_path_class = chain_path == -1 ? -1 : path_class_[chain_path]; - // Call the RIQ if the chain size is large enough; - // the optimal value was found with the associated benchmark in tests, + // Use a RIQ if the chain size is large enough; + // the optimal size was found with the associated benchmark in tests, // in particular BM_DimensionChecker. - constexpr int kMinRangeSizeForRIQ = 4; const bool chain_is_cached = chain_path_class == path_class; - if (last_index - first_index > kMinRangeSizeForRIQ && chain_is_cached && + if (last_index - first_index > min_range_size_for_riq_ && + chain_is_cached && SubpathOnlyHasTrivialNodes(first_index, last_index)) { - // Feasible cumuls at chain end are constrained by: - // - feasible cumuls at chain start plus demands of the chain. + // Feasible cumuls at chain's last node are constrained by: + // - feasible cumuls at chain's first node plus demands of the chain. // - tightest demand backwards sum from the path capacity interval. const Interval tightest_sum = GetTightestBackwardsDemandSum(first_index, last_index); @@ -3094,24 +3080,22 @@ bool DimensionChecker::Check() const { cumul += Delta(last_sum, prev_sum); cumul = Intersect(cumul, path_capacity + Delta(last_sum, tightest_sum)); if (IsEmpty(cumul)) return false; - cumul += min_max_demand_per_path_class_[path_class](last_node); cumul = Intersect(cumul, path_capacity); + prev_node = chain.Last(); } else { - for (const int node : chain) { - cumul = Intersect(cumul, node_capacity_[node]); - if (IsEmpty(cumul)) return false; + for (const int node : chain.WithoutFirstNode()) { const Interval demand = chain_is_cached - ? cached_demand_[node] - : min_max_demand_per_path_class_[path_class](node); + ? cached_demand_[prev_node] + : demand_per_path_class_[path_class](prev_node, node); cumul += demand; + cumul = Intersect(cumul, node_capacity_[node]); cumul = Intersect(cumul, path_capacity); + if (IsEmpty(cumul)) return false; + prev_node = node; } - if (IsEmpty(cumul)) return false; } } - cumul = Intersect(cumul, path_capacity); - if (IsEmpty(cumul)) return false; } return true; } @@ -3157,11 +3141,16 @@ void DimensionChecker::AppendPathDemandsToSums(int path) { // Compute sum of all demands for this path. // TODO(user): backwards Nodes() iterator, to compute sum directly. Interval demand_sum = {0, 0}; - for (const int node : path_state_->Nodes(path)) { - const Interval demand = min_max_demand_per_path_class_[path_class](node); + int prev = path_state_->Start(path); + const auto node_range = path_state_->Nodes(path); + for (auto it = ++node_range.begin(); it != node_range.end(); ++it) { + const int node = *it; + const Interval demand = demand_per_path_class_[path_class](prev, node); demand_sum += demand; - cached_demand_[node] = demand; + cached_demand_[prev] = demand; + prev = node; } + cached_demand_[path_state_->End(path)] = {0, 0}; int previous_nontrivial_index = -1; int index = backwards_demand_sums_riq_[0].size(); // Value of backwards_demand_sums_riq_ at node_index must be the sum @@ -3186,21 +3175,19 @@ void DimensionChecker::AppendPathDemandsToSums(int path) { void DimensionChecker::UpdateRIQStructure(int begin_index, int end_index) { // The max layer is the one used by // GetMinMaxPartialDemandSum(begin_index, end_index - 1). - const int maximum_riq_exponent = - MostSignificantBitPosition32(end_index - 1 - begin_index); - for (int layer = 1, window_size = 1; layer <= maximum_riq_exponent; - ++layer, window_size *= 2) { + const int max_layer = + MostSignificantBitPosition32(end_index - begin_index - 1); + for (int layer = 1, window = 1; layer <= max_layer; ++layer, window *= 2) { backwards_demand_sums_riq_[layer].resize(end_index); - for (int i = begin_index; i < end_index - window_size; ++i) { - const Interval i1 = backwards_demand_sums_riq_[layer - 1][i]; - const Interval i2 = - backwards_demand_sums_riq_[layer - 1][i + window_size]; - backwards_demand_sums_riq_[layer][i] = Intersect(i1, i2); + for (int i = begin_index; i < end_index - window; ++i) { + backwards_demand_sums_riq_[layer][i] = + Intersect(backwards_demand_sums_riq_[layer - 1][i], + backwards_demand_sums_riq_[layer - 1][i + window]); } std::copy( - backwards_demand_sums_riq_[layer - 1].begin() + end_index - window_size, + backwards_demand_sums_riq_[layer - 1].begin() + end_index - window, backwards_demand_sums_riq_[layer - 1].begin() + end_index, - backwards_demand_sums_riq_[layer].begin() + end_index - window_size); + backwards_demand_sums_riq_[layer].begin() + end_index - window); } } @@ -3215,16 +3202,14 @@ DimensionChecker::Interval DimensionChecker::GetTightestBackwardsDemandSum( DCHECK_LE(0, first_node_index); DCHECK_LT(first_node_index, last_node_index); DCHECK_LT(last_node_index, backwards_demand_sums_riq_[0].size()); - // Find largest window_size = 2^layer such that - // first_node_index < last_node_index - window_size + 1. + // Find largest window = 2^layer such that + // first_node_index < last_node_index - window + 1. const int layer = MostSignificantBitPosition32(last_node_index - first_node_index); - const int window_size = 1 << layer; - // Adaptation of the conventional O(1) Range Max Query scheme. - Interval i1 = backwards_demand_sums_riq_[layer][first_node_index]; - const Interval i2 = - backwards_demand_sums_riq_[layer][last_node_index - window_size + 1]; - return Intersect(i1, i2); + const int window = 1 << layer; + return Intersect( + backwards_demand_sums_riq_[layer][first_node_index], + backwards_demand_sums_riq_[layer][last_node_index - window + 1]); } bool DimensionChecker::SubpathOnlyHasTrivialNodes(int first_node_index, diff --git a/ortools/constraint_solver/python/constraint_solver.i b/ortools/constraint_solver/python/constraint_solver.i index 8b416f6cd0..608faf1768 100644 --- a/ortools/constraint_solver/python/constraint_solver.i +++ b/ortools/constraint_solver/python/constraint_solver.i @@ -1918,8 +1918,6 @@ namespace operations_research { // Ignored top-level classes, enums and methods: // - BaseIntExpr // - VarTypes (enum) -// - IntVarLocalSearchHandler -// - SequenceVarLocalSearchHandler // - ChangeValue // - PathOperator // - MakeLocalSearchOperator() @@ -1967,29 +1965,6 @@ namespace operations_research { %rename (NextNeighbor) LocalSearchOperator::MakeNextNeighbor; %unignore LocalSearchOperator::Start; -// VarLocalSearchOperator<> -// Ignored: -// - VarLocalSearchOperator() -// - ~VarLocalSearchOperator() -// - Start() -// - Var() -// - SkipUnchanged() -// - Activated() -// - Activate() -// - Deactivate() -// - ApplyChanges() -// - RevertChanges() -// - AddVars() -%unignore VarLocalSearchOperator; -%unignore VarLocalSearchOperator::Size; -%unignore VarLocalSearchOperator::Value; -%unignore VarLocalSearchOperator::IsIncremental; -%unignore VarLocalSearchOperator::IsIncremental; -%unignore VarLocalSearchOperator::OnStart; -%unignore VarLocalSearchOperator::OnStart; -%unignore VarLocalSearchOperator::OldValue; -%unignore VarLocalSearchOperator::SetValue; - // IntVarLocalSearchOperator // Ignored: @@ -2024,17 +1999,6 @@ namespace operations_research { %unignore ChangeValue::~ChangeValue; %unignore ChangeValue::ModifyValue; -// SequenceVarLocalSearchOperator -// Ignored: -// - SequenceVarLocalSearchOperator() -// - ~SequenceVarLocalSearchOperator() -// - Sequence() -// - OldSequence() -// - SetForwardSequence() -// - SetBackwardSequence() -%unignore SequenceVarLocalSearchOperator; -%unignore SequenceVarLocalSearchOperator::Start; - // PathOperator // Ignored: // - PathOperator() diff --git a/ortools/constraint_solver/python/routing.i b/ortools/constraint_solver/python/routing.i index 14f32906c1..cbd9d84aea 100644 --- a/ortools/constraint_solver/python/routing.i +++ b/ortools/constraint_solver/python/routing.i @@ -52,6 +52,7 @@ DEFINE_INDEX_TYPE_TYPEDEF( operations_research::RoutingVehicleClassIndex, operations_research::RoutingModel::VehicleClassIndex); + %ignore operations_research::RoutingModel::RegisterStateDependentTransitCallback; %ignore operations_research::RoutingModel::StateDependentTransitCallback; %ignore operations_research::RoutingModel::MakeStateDependentTransit; diff --git a/ortools/constraint_solver/resource.cc b/ortools/constraint_solver/resource.cc index 2ba5b22381..02d0e87487 100644 --- a/ortools/constraint_solver/resource.cc +++ b/ortools/constraint_solver/resource.cc @@ -2194,7 +2194,7 @@ class CumulativeConstraint : public Constraint { // For the cumulative constraint, there are many propagators, and they // don't dominate each other. So the strongest propagation is obtained // by posting a bunch of different propagators. - const ConstraintSolverParameters& params = solver()->parameters(); + const ConstraintSolverParameters& params = solver()->const_parameters(); if (params.use_cumulative_time_table()) { if (params.use_cumulative_time_table_sync()) { PostOneSidedConstraint(false, false, true); @@ -2329,7 +2329,7 @@ class CumulativeConstraint : public Constraint { } else { Solver* const s = solver(); if (edge_finder) { - const ConstraintSolverParameters& params = solver()->parameters(); + const ConstraintSolverParameters& params = solver()->const_parameters(); return useful_tasks.size() < params.max_edge_finder_size() ? s->RevAlloc(new EdgeFinder(s, useful_tasks, capacity_)) @@ -2388,7 +2388,7 @@ class VariableDemandCumulativeConstraint : public Constraint { // For the cumulative constraint, there are many propagators, and they // don't dominate each other. So the strongest propagation is obtained // by posting a bunch of different propagators. - const ConstraintSolverParameters& params = solver()->parameters(); + const ConstraintSolverParameters& params = solver()->const_parameters(); if (params.use_cumulative_time_table()) { PostOneSidedConstraint(false, false, false); PostOneSidedConstraint(true, false, false); diff --git a/ortools/constraint_solver/routing.cc b/ortools/constraint_solver/routing.cc index f3dbf24ba7..68a525bd40 100644 --- a/ortools/constraint_solver/routing.cc +++ b/ortools/constraint_solver/routing.cc @@ -42,6 +42,7 @@ #include "absl/strings/str_format.h" #include "absl/strings/string_view.h" #include "absl/time/time.h" +#include "ortools/base/dump_vars.h" #include "ortools/base/int_type.h" #include "ortools/base/integral_types.h" #include "ortools/base/logging.h" @@ -97,6 +98,46 @@ class TwoOpt; namespace operations_research { +std::string RoutingModel::RouteDimensionTravelInfo::DebugString( + std::string line_prefix) const { + std::string s = absl::StrFormat("%stravel_cost_coefficient: %ld", line_prefix, + travel_cost_coefficient); + for (int i = 0; i < transition_info.size(); ++i) { + absl::StrAppendFormat(&s, "\ntransition[%d] {\n%s\n}\n", i, + transition_info[i].DebugString(line_prefix + "\t")); + } + return s; +} + +std::string RoutingModel::RouteDimensionTravelInfo::TransitionInfo::DebugString( + std::string line_prefix) const { + return absl::StrFormat( + R"({ +%spre: %ld +%spost: %ld +%slower_bound: %ld +%supper_bound: %ld +%stravel_value: %s +%scost: %s +})", + line_prefix, pre_travel_transit_value, line_prefix, + post_travel_transit_value, line_prefix, + compressed_travel_value_lower_bound, line_prefix, + travel_value_upper_bound, line_prefix, + travel_start_dependent_travel.DebugString(line_prefix + "\t"), + line_prefix, travel_compression_cost.DebugString(line_prefix + "\t")); +} + +std::string RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation::DebugString(std::string line_prefix) const { + if (x_anchors.size() <= 10) { + return "{ " + DUMP_VARS(x_anchors, y_anchors).str() + "}"; + } + return absl::StrFormat("{\n%s%s\n%s%s\n}", line_prefix, + DUMP_VARS(x_anchors).str(), line_prefix, + DUMP_VARS(y_anchors).str()); +} + namespace { using ResourceGroup = RoutingModel::ResourceGroup; @@ -252,11 +293,18 @@ class SetCumulsFromLocalDimensionCosts : public DecisionBuilder { SetCumulsFromLocalDimensionCosts( LocalDimensionCumulOptimizer* local_optimizer, LocalDimensionCumulOptimizer* local_mp_optimizer, SearchMonitor* monitor, - bool optimize_and_pack = false) + bool optimize_and_pack = false, + std::vector + dimension_travel_info_per_route = {}) : local_optimizer_(local_optimizer), local_mp_optimizer_(local_mp_optimizer), monitor_(monitor), - optimize_and_pack_(optimize_and_pack) { + optimize_and_pack_(optimize_and_pack), + dimension_travel_info_per_route_( + std::move(dimension_travel_info_per_route)) { + DCHECK(dimension_travel_info_per_route_.empty() || + dimension_travel_info_per_route_.size() == + local_optimizer_->dimension()->model()->vehicles()); const RoutingDimension* const dimension = local_optimizer->dimension(); const std::vector& resource_groups = dimension->model()->GetDimensionResourceGroupIndices(dimension); @@ -327,6 +375,7 @@ class SetCumulsFromLocalDimensionCosts : public DecisionBuilder { private: using Resource = RoutingModel::ResourceGroup::Resource; + using RouteDimensionTravelInfo = RoutingModel::RouteDimensionTravelInfo; DimensionSchedulingStatus ComputeCumulAndBreakValuesForVehicle( LocalDimensionCumulOptimizer* optimizer, int vehicle, @@ -336,6 +385,10 @@ class SetCumulsFromLocalDimensionCosts : public DecisionBuilder { break_start_end_values->clear(); RoutingModel* const model = optimizer->dimension()->model(); const auto next = [model](int64_t n) { return model->NextVar(n)->Value(); }; + const RouteDimensionTravelInfo& dimension_travel_info = + dimension_travel_info_per_route_.empty() + ? RouteDimensionTravelInfo() + : dimension_travel_info_per_route_[vehicle]; if (optimize_and_pack_) { const int resource_index = resource_group_index_ < 0 @@ -346,10 +399,12 @@ class SetCumulsFromLocalDimensionCosts : public DecisionBuilder { : &model->GetResourceGroup(resource_group_index_) ->GetResource(resource_index); return optimizer->ComputePackedRouteCumuls( - vehicle, next, resource, cumul_values, break_start_end_values); + vehicle, next, dimension_travel_info, resource, cumul_values, + break_start_end_values); } else { // TODO(user): Add the resource to the call in this case too! - return optimizer->ComputeRouteCumuls(vehicle, next, cumul_values, + return optimizer->ComputeRouteCumuls(vehicle, next, dimension_travel_info, + cumul_values, break_start_end_values); } } @@ -360,6 +415,7 @@ class SetCumulsFromLocalDimensionCosts : public DecisionBuilder { int resource_group_index_; SearchMonitor* const monitor_; const bool optimize_and_pack_; + const std::vector dimension_travel_info_per_route_; }; class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { @@ -367,11 +423,19 @@ class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { SetCumulsFromGlobalDimensionCosts( GlobalDimensionCumulOptimizer* global_optimizer, GlobalDimensionCumulOptimizer* global_mp_optimizer, - SearchMonitor* monitor, bool optimize_and_pack = false) + SearchMonitor* monitor, bool optimize_and_pack = false, + std::vector + dimension_travel_info_per_route = {}) : global_optimizer_(global_optimizer), global_mp_optimizer_(global_mp_optimizer), monitor_(monitor), - optimize_and_pack_(optimize_and_pack) {} + optimize_and_pack_(optimize_and_pack), + dimension_travel_info_per_route_( + std::move(dimension_travel_info_per_route)) { + DCHECK(dimension_travel_info_per_route_.empty() || + dimension_travel_info_per_route_.size() == + global_optimizer_->dimension()->model()->vehicles()); + } Decision* Next(Solver* const solver) override { // The following boolean variable indicates if the solver should fail, in @@ -473,31 +537,33 @@ class SetCumulsFromGlobalDimensionCosts : public DecisionBuilder { RoutingModel* const model = optimizer->dimension()->model(); const auto next = [model](int64_t n) { return model->NextVar(n)->Value(); }; return optimize_and_pack_ - ? optimizer->ComputePackedCumuls(next, cumul_values, - break_start_end_values, - resource_indices_per_group) - : optimizer->ComputeCumuls(next, cumul_values, - break_start_end_values, - resource_indices_per_group); + ? optimizer->ComputePackedCumuls( + next, dimension_travel_info_per_route_, cumul_values, + break_start_end_values, resource_indices_per_group) + : optimizer->ComputeCumuls( + next, dimension_travel_info_per_route_, cumul_values, + break_start_end_values, resource_indices_per_group); } GlobalDimensionCumulOptimizer* const global_optimizer_; GlobalDimensionCumulOptimizer* const global_mp_optimizer_; SearchMonitor* const monitor_; const bool optimize_and_pack_; + const std::vector + dimension_travel_info_per_route_; }; class SetCumulsFromResourceAssignmentCosts : public DecisionBuilder { public: SetCumulsFromResourceAssignmentCosts( - LocalDimensionCumulOptimizer* optimizer, + LocalDimensionCumulOptimizer* lp_optimizer, LocalDimensionCumulOptimizer* mp_optimizer, SearchMonitor* monitor) - : model_(*optimizer->dimension()->model()), - dimension_(*optimizer->dimension()), + : model_(*lp_optimizer->dimension()->model()), + dimension_(*lp_optimizer->dimension()), + lp_optimizer_(lp_optimizer), + mp_optimizer_(mp_optimizer), rg_index_(model_.GetDimensionResourceGroupIndex(&dimension_)), resource_group_(*model_.GetResourceGroup(rg_index_)), - resource_assignment_optimizer_(&resource_group_, optimizer, - mp_optimizer), monitor_(monitor) {} Decision* Next(Solver* const solver) override { @@ -514,20 +580,23 @@ class SetCumulsFromResourceAssignmentCosts : public DecisionBuilder { DCHECK(DimensionFixedTransitsEqualTransitEvaluators(dimension_)); for (int v : resource_group_.GetVehiclesRequiringAResource()) { - if (!resource_assignment_optimizer_.ComputeAssignmentCostsForVehicle( - v, next, dimension_.transit_evaluator(v), - /*optimize_vehicle_costs*/ true, &assignment_costs[v], - &cumul_values[v], &break_values[v])) { + if (!ComputeVehicleToResourcesAssignmentCosts( + v, resource_group_, next, dimension_.transit_evaluator(v), + /*optimize_vehicle_costs*/ true, lp_optimizer_, mp_optimizer_, + &assignment_costs[v], &cumul_values[v], &break_values[v])) { should_fail = true; break; } } - std::vector resource_indices; - should_fail = should_fail || - resource_assignment_optimizer_.ComputeBestAssignmentCost( - assignment_costs, assignment_costs, - [](int) { return true; }, &resource_indices) < 0; + std::vector resource_indices(num_vehicles); + should_fail = + should_fail || + ComputeBestVehicleToResourceAssignment( + resource_group_.GetVehiclesRequiringAResource(), + resource_group_.Size(), + [&assignment_costs](int v) { return &assignment_costs[v]; }, + &resource_indices) < 0; if (!should_fail) { DCHECK_EQ(resource_indices.size(), num_vehicles); @@ -577,9 +646,10 @@ class SetCumulsFromResourceAssignmentCosts : public DecisionBuilder { private: const RoutingModel& model_; const RoutingDimension& dimension_; + LocalDimensionCumulOptimizer* lp_optimizer_; + LocalDimensionCumulOptimizer* mp_optimizer_; const int rg_index_; const ResourceGroup& resource_group_; - ResourceAssignmentOptimizer resource_assignment_optimizer_; SearchMonitor* const monitor_; }; @@ -698,156 +768,6 @@ class DifferentFromValues : public Constraint { const std::vector values_; }; -// Set of "light" constraints, well-suited for use within Local Search. -// These constraints are "checking" constraints, only triggered on WhenBound -// events. The provide very little (or no) domain filtering. -// TODO(user): Move to core constraintsolver library. - -// Light one-dimension function-based element constraint ensuring: -// var == values(index). -// Doesn't perform bound reduction of the resulting variable until the index -// variable is bound. -// If deep_serialize returns false, the model visitor will not extract all -// possible values from the values function. -template -class LightFunctionElementConstraint : public Constraint { - public: - LightFunctionElementConstraint(Solver* const solver, IntVar* const var, - IntVar* const index, F values, - std::function deep_serialize) - : Constraint(solver), - var_(var), - index_(index), - values_(std::move(values)), - deep_serialize_(std::move(deep_serialize)) {} - ~LightFunctionElementConstraint() override {} - - void Post() override { - Demon* demon = MakeConstraintDemon0( - solver(), this, &LightFunctionElementConstraint::IndexBound, - "IndexBound"); - index_->WhenBound(demon); - } - - void InitialPropagate() override { - if (index_->Bound()) { - IndexBound(); - } - } - - std::string DebugString() const override { - return "LightFunctionElementConstraint"; - } - - void Accept(ModelVisitor* const visitor) const override { - visitor->BeginVisitConstraint(RoutingModelVisitor::kLightElement, this); - visitor->VisitIntegerExpressionArgument(ModelVisitor::kTargetArgument, - var_); - visitor->VisitIntegerExpressionArgument(ModelVisitor::kIndexArgument, - index_); - // Warning: This will expand all values into a vector. - if (deep_serialize_()) { - visitor->VisitInt64ToInt64Extension(values_, index_->Min(), - index_->Max()); - } - visitor->EndVisitConstraint(RoutingModelVisitor::kLightElement, this); - } - - private: - void IndexBound() { var_->SetValue(values_(index_->Min())); } - - IntVar* const var_; - IntVar* const index_; - F values_; - std::function deep_serialize_; -}; - -template -Constraint* MakeLightElement(Solver* const solver, IntVar* const var, - IntVar* const index, F values, - std::function deep_serialize) { - return solver->RevAlloc(new LightFunctionElementConstraint( - solver, var, index, std::move(values), std::move(deep_serialize))); -} - -// Light two-dimension function-based element constraint ensuring: -// var == values(index1, index2). -// Doesn't perform bound reduction of the resulting variable until the index -// variables are bound. -// Ownership of the 'values' callback is taken by the constraint. -template -class LightFunctionElement2Constraint : public Constraint { - public: - LightFunctionElement2Constraint(Solver* const solver, IntVar* const var, - IntVar* const index1, IntVar* const index2, - F values, - std::function deep_serialize) - : Constraint(solver), - var_(var), - index1_(index1), - index2_(index2), - values_(std::move(values)), - deep_serialize_(std::move(deep_serialize)) {} - ~LightFunctionElement2Constraint() override {} - void Post() override { - Demon* demon = MakeConstraintDemon0( - solver(), this, &LightFunctionElement2Constraint::IndexBound, - "IndexBound"); - index1_->WhenBound(demon); - index2_->WhenBound(demon); - } - void InitialPropagate() override { IndexBound(); } - - std::string DebugString() const override { - return "LightFunctionElement2Constraint"; - } - - void Accept(ModelVisitor* const visitor) const override { - visitor->BeginVisitConstraint(RoutingModelVisitor::kLightElement2, this); - visitor->VisitIntegerExpressionArgument(ModelVisitor::kTargetArgument, - var_); - visitor->VisitIntegerExpressionArgument(ModelVisitor::kIndexArgument, - index1_); - visitor->VisitIntegerExpressionArgument(ModelVisitor::kIndex2Argument, - index2_); - // Warning: This will expand all values into a vector. - const int64_t index1_min = index1_->Min(); - const int64_t index1_max = index1_->Max(); - visitor->VisitIntegerArgument(ModelVisitor::kMinArgument, index1_min); - visitor->VisitIntegerArgument(ModelVisitor::kMaxArgument, index1_max); - if (deep_serialize_()) { - for (int i = index1_min; i <= index1_max; ++i) { - visitor->VisitInt64ToInt64Extension( - [this, i](int64_t j) { return values_(i, j); }, index2_->Min(), - index2_->Max()); - } - } - visitor->EndVisitConstraint(RoutingModelVisitor::kLightElement2, this); - } - - private: - void IndexBound() { - if (index1_->Bound() && index2_->Bound()) { - var_->SetValue(values_(index1_->Min(), index2_->Min())); - } - } - - IntVar* const var_; - IntVar* const index1_; - IntVar* const index2_; - Solver::IndexEvaluator2 values_; - std::function deep_serialize_; -}; - -template -Constraint* MakeLightElement2(Solver* const solver, IntVar* const var, - IntVar* const index1, IntVar* const index2, - F values, std::function deep_serialize) { - return solver->RevAlloc(new LightFunctionElement2Constraint( - solver, var, index1, index2, std::move(values), - std::move(deep_serialize))); -} - // For each vehicle, computes information on the partially fixed start/end // chains (based on bound NextVar values): // - For every 'end_node', the last node of a start chain of a vehicle, @@ -970,9 +890,6 @@ class ResourceAssignmentConstraint : public Constraint { LocalDimensionCumulOptimizer* const mp_optimizer = model_.GetMutableLocalCumulMPOptimizer(dimension); DCHECK_NE(mp_optimizer, nullptr); - ResourceAssignmentOptimizer resource_assignment_optimizer( - &resource_group_, optimizer, mp_optimizer); - const auto transit = [&dimension](int64_t node, int64_t /*next*/) { // TODO(user): Get rid of this max() by only allowing resources on // dimensions with positive transits (model.AreVehicleTransitsPositive()). @@ -983,17 +900,22 @@ class ResourceAssignmentConstraint : public Constraint { }; std::vector> assignment_costs(model_.vehicles()); - for (int v : resource_group_.GetVehiclesRequiringAResource()) { - if (!resource_assignment_optimizer.ComputeAssignmentCostsForVehicle( - v, next, transit, /*optimize_vehicle_costs*/ false, + if (!ComputeVehicleToResourcesAssignmentCosts( + v, resource_group_, next, transit, + /*optimize_vehicle_costs*/ false, + model_.GetMutableLocalCumulLPOptimizer(dimension), + model_.GetMutableLocalCumulMPOptimizer(dimension), &assignment_costs[v], nullptr, nullptr)) { return false; } } - - return resource_assignment_optimizer.ComputeBestAssignmentCost( - assignment_costs, assignment_costs, [](int) { return true; }, + // TODO(user): Replace this call with a more efficient max-flow, instead + // of running the full min-cost flow. + return ComputeBestVehicleToResourceAssignment( + resource_group_.GetVehiclesRequiringAResource(), + resource_group_.Size(), + [&assignment_costs](int v) { return &assignment_costs[v]; }, nullptr) >= 0; } @@ -1030,8 +952,7 @@ class ResourceAssignmentConstraint : public Constraint { start_cumul ? vehicle_to_start_bound_vars_per_dimension_[v][d].lower_bound : vehicle_to_end_bound_vars_per_dimension_[v][d].lower_bound; - s->AddConstraint(MakeLightElement( - s, resource_lb_var, resource_var, + s->AddConstraint(s->MakeLightElement( [dim, start_cumul, &resource_group = resource_group_](int r) { if (r < 0) return std::numeric_limits::min(); return start_cumul ? resource_group.GetResources()[r] @@ -1043,6 +964,7 @@ class ResourceAssignmentConstraint : public Constraint { .end_domain() .Min(); }, + resource_lb_var, resource_var, [&model = model_]() { return model.enable_deep_serialization(); })); @@ -1052,8 +974,7 @@ class ResourceAssignmentConstraint : public Constraint { start_cumul ? vehicle_to_start_bound_vars_per_dimension_[v][d].upper_bound : vehicle_to_end_bound_vars_per_dimension_[v][d].upper_bound; - s->AddConstraint(MakeLightElement( - s, resource_ub_var, resource_var, + s->AddConstraint(s->MakeLightElement( [dim, start_cumul, &resource_group = resource_group_](int r) { if (r < 0) return std::numeric_limits::max(); return start_cumul ? resource_group.GetResources()[r] @@ -1065,6 +986,7 @@ class ResourceAssignmentConstraint : public Constraint { .end_domain() .Max(); }, + resource_ub_var, resource_var, [&model = model_]() { return model.enable_deep_serialization(); })); @@ -1152,8 +1074,7 @@ RoutingModel::RoutingModel(const RoutingIndexManager& index_manager, has_same_vehicle_type_requirements_(false), has_temporal_type_requirements_(false), num_visit_types_(0), - starts_(vehicles_), - ends_(vehicles_), + paths_metadata_(index_manager), manager_(index_manager) { // Initialize vehicle costs to the zero evaluator. vehicle_to_transit_cost_.assign( @@ -1177,14 +1098,6 @@ RoutingModel::RoutingModel(const RoutingIndexManager& index_manager, index_to_visit_type_.resize(index_manager.num_indices(), kUnassigned); index_to_type_policy_.resize(index_manager.num_indices()); - index_to_vehicle_.resize(index_manager.num_indices(), kUnassigned); - for (int v = 0; v < index_manager.num_vehicles(); ++v) { - starts_[v] = index_manager.GetStartIndex(v); - index_to_vehicle_[starts_[v]] = v; - ends_[v] = index_manager.GetEndIndex(v); - index_to_vehicle_[ends_[v]] = v; - } - const std::vector& index_to_node = index_manager.GetIndexToNodeMap(); index_to_equivalence_class_.resize(index_manager.num_indices()); @@ -1255,14 +1168,13 @@ int RegisterUnaryCallback(RoutingTransitCallback1 callback, bool is_positive, } // namespace int RoutingModel::RegisterUnaryTransitVector(std::vector values) { + bool is_positive = std::all_of(std::cbegin(values), std::cend(values), + [](int64_t transit) { return transit >= 0; }); return RegisterUnaryCallback( [this, values = std::move(values)](int64_t i) { return values[manager_.IndexToNode(i).value()]; }, - /*is_positive=*/ - std::all_of(std::cbegin(values), std::cend(values), - [](int64_t transit) { return transit >= 0; }), - this); + is_positive, this); } int RoutingModel::RegisterUnaryTransitCallback(TransitCallback1 callback) { @@ -1737,8 +1649,17 @@ const ResourceGroup::Attributes& ResourceGroup::Resource::GetDefaultAttributes() } int RoutingModel::AddResourceGroup() { + DCHECK_EQ(resource_groups_.size(), resource_vars_.size()); + // Create and add the resource group. resource_groups_.push_back(std::make_unique(this)); - return resource_groups_.size() - 1; + // Create and add the resource vars (the proper variable bounds and + // constraints are set up when closing the model). + const int rg_index = resource_groups_.size() - 1; + resource_vars_.push_back({}); + solver_->MakeIntVarArray(vehicles(), -1, std::numeric_limits::max(), + absl::StrCat("Resources[", rg_index, "]"), + &resource_vars_.back()); + return rg_index; } int RoutingModel::ResourceGroup::AddResource( @@ -2415,8 +2336,8 @@ void RoutingModel::AppendHomogeneousArcCosts( // the search in GLS in some cases (Solomon instances for instance). IntVar* const base_cost_var = solver_->MakeIntVar(0, std::numeric_limits::max()); - solver_->AddConstraint(MakeLightElement( - solver_.get(), base_cost_var, nexts_[node_index], arc_cost_evaluator, + solver_->AddConstraint(solver_->MakeLightElement( + arc_cost_evaluator, base_cost_var, nexts_[node_index], [this]() { return enable_deep_serialization_; })); IntVar* const var = solver_->MakeProd(base_cost_var, active_[node_index])->Var(); @@ -2440,12 +2361,11 @@ void RoutingModel::AppendArcCosts(const RoutingSearchParameters& parameters, // the search in GLS in some cases (Solomon instances for instance). IntVar* const base_cost_var = solver_->MakeIntVar(0, std::numeric_limits::max()); - solver_->AddConstraint(MakeLightElement2( - solver_.get(), base_cost_var, nexts_[node_index], - vehicle_vars_[node_index], + solver_->AddConstraint(solver_->MakeLightElement( [this, node_index](int64_t to, int64_t vehicle) { return GetArcCostForVehicle(node_index, to, vehicle); }, + base_cost_var, nexts_[node_index], vehicle_vars_[node_index], [this]() { return enable_deep_serialization_; })); IntVar* const var = solver_->MakeProd(base_cost_var, active_[node_index])->Var(); @@ -2470,7 +2390,7 @@ void RoutingModel::AppendArcCosts(const RoutingSearchParameters& parameters, } int RoutingModel::GetVehicleStartClass(int64_t start_index) const { - const int vehicle = index_to_vehicle_[start_index]; + const int vehicle = VehicleIndex(start_index); if (vehicle != kUnassigned) { return GetVehicleClassIndexOfVehicle(vehicle).value(); } @@ -2719,8 +2639,8 @@ void RoutingModel::CloseModelWithParameters( // Vehicle variable constraints for (int i = 0; i < vehicles_; ++i) { - const int64_t start = starts_[i]; - const int64_t end = ends_[i]; + const int64_t start = Start(i); + const int64_t end = End(i); solver_->AddConstraint( solver_->MakeEquality(vehicle_vars_[start], solver_->MakeIntConst(i))); solver_->AddConstraint( @@ -2787,8 +2707,8 @@ void RoutingModel::CloseModelWithParameters( // Reduce domain of next variables. for (int i = 0; i < size; ++i) { // No variable can point back to a start. - solver_->AddConstraint(solver_->RevAlloc( - new DifferentFromValues(solver_.get(), nexts_[i], starts_))); + solver_->AddConstraint(solver_->RevAlloc(new DifferentFromValues( + solver_.get(), nexts_[i], paths_metadata_.Starts()))); // Extra constraint to state an active node can't point to itself. solver_->AddConstraint( solver_->MakeIsDifferentCstCt(nexts_[i], i, active_[i])); @@ -2812,15 +2732,15 @@ void RoutingModel::CloseModelWithParameters( forbidden_ends.reserve(vehicles_ - 1); for (int j = 0; j < vehicles_; ++j) { if (i != j) { - forbidden_ends.push_back(ends_[j]); + forbidden_ends.push_back(End(j)); } } solver_->AddConstraint(solver_->RevAlloc(new DifferentFromValues( - solver_.get(), nexts_[starts_[i]], std::move(forbidden_ends)))); + solver_.get(), nexts_[Start(i)], std::move(forbidden_ends)))); } // Constraining is_bound_to_end_ variables. - for (const int64_t end : ends_) { + for (const int64_t end : paths_metadata_.Ends()) { is_bound_to_end_[end]->SetValue(1); } @@ -3077,13 +2997,14 @@ void RoutingModel::CloseModelWithParameters( } if (!resource_groups_.empty()) { - resource_vars_.resize(resource_groups_.size()); + DCHECK_EQ(resource_vars_.size(), resource_groups_.size()); for (int rg = 0; rg < resource_groups_.size(); ++rg) { const auto& resource_group = resource_groups_[rg]; + const int max_resource_index = resource_group->Size() - 1; std::vector& vehicle_res_vars = resource_vars_[rg]; - solver_->MakeIntVarArray(vehicles(), -1, resource_group->Size() - 1, - absl::StrCat("Resources[", rg, "]"), - &vehicle_res_vars); + for (IntVar* res_var : vehicle_res_vars) { + res_var->SetMax(max_resource_index); + } solver_->AddConstraint(MakeResourceConstraint(resource_group.get(), &vehicle_res_vars, this)); } @@ -3284,6 +3205,7 @@ class AtSolutionCallbackMonitor : public SearchMonitor { callback_(); return false; } + void Install() override { ListenToEvent(Solver::MonitorEvent::kAtSolution); } private: std::function callback_; @@ -3601,7 +3523,7 @@ int64_t RoutingModel::ComputeLowerBound() { // Therefore we are creating fake assignments for end nodes, forced to point // to the equivalent start node with a cost of 0. for (int tail = Size(); tail < num_nodes; ++tail) { - const ArcIndex arc = graph.AddArc(tail, num_nodes + starts_[tail - Size()]); + const ArcIndex arc = graph.AddArc(tail, num_nodes + Start(tail - Size())); linear_sum_assignment.SetArcCost(arc, 0); } if (linear_sum_assignment.ComputeAssignment()) { @@ -4083,11 +4005,11 @@ int64_t RoutingModel::GetArcCostForClassInternal( cost = CapAdd( evaluator(from_index, to_index), CapAdd(GetDimensionTransitCostSum(from_index, to_index, cost_class), - fixed_cost_of_vehicle_[index_to_vehicle_[from_index]])); + fixed_cost_of_vehicle_[VehicleIndex(from_index)])); } else { // If there's only the first and last nodes on the route, it is considered // as an empty route. - if (vehicle_used_when_empty_[index_to_vehicle_[from_index]]) { + if (vehicle_used_when_empty_[VehicleIndex(from_index)]) { cost = CapAdd(evaluator(from_index, to_index), GetDimensionTransitCostSum(from_index, to_index, cost_class)); @@ -4099,10 +4021,6 @@ int64_t RoutingModel::GetArcCostForClassInternal( return cost; } -bool RoutingModel::IsStart(int64_t index) const { - return !IsEnd(index) && index_to_vehicle_[index] != kUnassigned; -} - bool RoutingModel::IsVehicleUsed(const Assignment& assignment, int vehicle) const { CHECK_GE(vehicle, 0); @@ -4719,7 +4637,7 @@ void RoutingModel::CreateNeighborhoodOperators( const auto arc_cost_for_path_start = [this](int64_t before_node, int64_t after_node, int64_t start_index) { - const int vehicle = index_to_vehicle_[start_index]; + const int vehicle = VehicleIndex(start_index); const int64_t arc_cost = GetArcCostForVehicle(before_node, after_node, vehicle); return (before_node != start_index || IsEnd(after_node)) @@ -5034,8 +4952,8 @@ RoutingModel::CreateLocalSearchFilters( if (HasUnaryDimension(GetDimensions())) { std::vector path_starts; std::vector path_ends; - ConvertVectorInt64ToVectorInt(starts_, &path_starts); - ConvertVectorInt64ToVectorInt(ends_, &path_ends); + ConvertVectorInt64ToVectorInt(paths_metadata_.Starts(), &path_starts); + ConvertVectorInt64ToVectorInt(paths_metadata_.Ends(), &path_ends); auto path_state = std::make_unique( Size() + vehicles(), std::move(path_starts), std::move(path_ends)); @@ -5806,7 +5724,7 @@ void RoutingModel::SetupDecisionBuilders( void RoutingModel::SetupMetaheuristics( const RoutingSearchParameters& search_parameters) { - SearchMonitor* optimize; + SearchMonitor* optimize = nullptr; const LocalSearchMetaheuristic::Value metaheuristic = search_parameters.local_search_metaheuristic(); // Some metaheuristics will effectively never terminate; warn @@ -6508,11 +6426,6 @@ void RoutingDimension::InitializeCumuls() { } namespace { -template -int64_t IthElementOrValue(const std::vector& v, int64_t index) { - return index >= 0 ? v[index] : value; -} - void ComputeTransitClasses(const std::vector& evaluator_indices, std::vector* class_evaluators, std::vector* vehicle_to_class) { @@ -6917,8 +6830,8 @@ void RoutingDimension::CloseModel(bool use_light_propagation) { IntVar* const vehicle_var = model_->VehicleVar(i); IntVar* const capacity_var = capacity_vars_[i]; if (use_light_propagation) { - solver->AddConstraint(MakeLightElement( - solver, capacity_var, vehicle_var, capacity_lambda, + solver->AddConstraint(solver->MakeLightElement( + capacity_lambda, capacity_var, vehicle_var, [this]() { return model_->enable_deep_serialization_; })); } else { solver->AddConstraint(solver->MakeEquality( @@ -6939,19 +6852,19 @@ void RoutingDimension::CloseModel(bool use_light_propagation) { const auto& unary_callback = model_->UnaryTransitCallbackOrNull(class_evaluator_index); if (unary_callback == nullptr) { - solver->AddConstraint(MakeLightElement( - solver, fixed_transit, next_var, + solver->AddConstraint(solver->MakeLightElement( [this, i](int64_t to) { return model_->TransitCallback(class_evaluators_[0])(i, to); }, + fixed_transit, next_var, [this]() { return model_->enable_deep_serialization_; })); } else { fixed_transit->SetValue(unary_callback(i)); } } else { - solver->AddConstraint(MakeLightElement2( - solver, fixed_transit, next_var, model_->VehicleVar(i), - transit_vehicle_evaluator, + solver->AddConstraint(solver->MakeLightElement( + transit_vehicle_evaluator, fixed_transit, next_var, + model_->VehicleVar(i), [this]() { return model_->enable_deep_serialization_; })); } } else { @@ -7459,4 +7372,5 @@ void RoutingDimension::SetupSlackAndDependentTransitCosts() const { } } } + } // namespace operations_research diff --git a/ortools/constraint_solver/routing.h b/ortools/constraint_solver/routing.h index a2538e941a..a1f961f52d 100644 --- a/ortools/constraint_solver/routing.h +++ b/ortools/constraint_solver/routing.h @@ -158,6 +158,7 @@ #define OR_TOOLS_CONSTRAINT_SOLVER_ROUTING_H_ #include +#include #include #include #include @@ -169,6 +170,7 @@ #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" +#include "absl/container/inlined_vector.h" #include "absl/functional/bind_front.h" #include "absl/memory/memory.h" #include "absl/time/time.h" @@ -207,6 +209,44 @@ using util::ReverseArcListGraph; class SweepArranger; #endif +class PathsMetadata { + public: + explicit PathsMetadata(const RoutingIndexManager& manager) { + const int num_indices = manager.num_indices(); + const int num_paths = manager.num_vehicles(); + path_of_node_.resize(num_indices, -1); + is_start_.resize(num_indices, false); + is_end_.resize(num_indices, false); + start_of_path_.resize(num_paths); + end_of_path_.resize(num_paths); + for (int v = 0; v < num_paths; ++v) { + const int64_t start = manager.GetStartIndex(v); + start_of_path_[v] = start; + path_of_node_[start] = v; + is_start_[start] = true; + const int64_t end = manager.GetEndIndex(v); + end_of_path_[v] = end; + path_of_node_[end] = v; + is_end_[end] = true; + } + } + + bool IsStart(int64_t node) const { return is_start_[node]; } + bool IsEnd(int64_t node) const { return is_end_[node]; } + int GetPath(int64_t start_or_end_node) const { + return path_of_node_[start_or_end_node]; + } + const std::vector& Starts() const { return start_of_path_; } + const std::vector& Ends() const { return end_of_path_; } + + private: + std::vector is_start_; + std::vector is_end_; + std::vector start_of_path_; + std::vector end_of_path_; + std::vector path_of_node_; +}; + class RoutingModel { public: /// Status of the search. @@ -1300,6 +1340,65 @@ class RoutingModel { const Assignment* PackCumulsOfOptimizerDimensionsFromAssignment( const Assignment* original_assignment, absl::Duration duration_limit, bool* time_limit_was_reached = nullptr); + /// Contains the information needed by the solver to optimize a dimension's + /// cumuls with travel-start dependent transit values. + struct RouteDimensionTravelInfo { + /// Contains the information for a single transition on the route. + struct TransitionInfo { + /// The following struct defines a piecewise linear formulation, with + /// int64_t values for the "anchor" x and y values, and potential double + /// values for the slope of each linear function. + // TODO(user): Adjust the inlined vector sizes based on experiments. + struct PiecewiseLinearFormulation { + /// The set of *increasing* anchor cumul values for the interpolation. + absl::InlinedVector x_anchors; + /// The y values used for the interpolation: + /// For any x anchor value, let i be an index such that + /// x_anchors[i] ≤ x < x_anchors[i+1], then the y value for x is + /// y_anchors[i] * (1-λ) + y_anchors[i+1] * λ, with + /// λ = (x - x_anchors[i]) / (x_anchors[i+1] - x_anchors[i]). + absl::InlinedVector y_anchors; + + std::string DebugString(std::string line_prefix = "") const; + }; + + /// Models the (real) travel value Tᵣ, for this transition based on the + /// departure value of the travel. + PiecewiseLinearFormulation travel_start_dependent_travel; + + /// travel_compression_cost models the cost of the difference between the + /// (real) travel value Tᵣ given by travel_start_dependent_travel and the + /// compressed travel value considered in the scheduling. + PiecewiseLinearFormulation travel_compression_cost; + + /// The parts of the transit which occur pre/post travel between the + /// nodes. The total transit between the two nodes i and j is + /// = pre_travel_transit_value + travel(i, j) + post_travel_transit_value. + int64_t pre_travel_transit_value; + int64_t post_travel_transit_value; + + /// The hard lower bound of the compressed travel value that will be + /// enforced by the scheduling module. + int64_t compressed_travel_value_lower_bound; + + /// The hard upper bound of the (real) travel value Tᵣ (see + /// above). This value should be chosen so as to prevent + /// the overall cost of the model + /// (dimension costs + travel_compression_cost) to overflow. + int64_t travel_value_upper_bound; + + std::string DebugString(std::string line_prefix = "") const; + }; + + /// For each node #i on the route, transition_info[i] contains the relevant + /// information for the travel between nodes #i and #(i + 1) on the route. + std::vector transition_info; + /// The cost per unit of travel for this vehicle. + int64_t travel_cost_coefficient; + + std::string DebugString(std::string line_prefix = "") const; + }; + #ifndef SWIG // TODO(user): Revisit if coordinates are added to the RoutingModel class. void SetSweepArranger(SweepArranger* sweep_arranger); @@ -1322,16 +1421,18 @@ class RoutingModel { /// Model inspection. /// Returns the variable index of the starting node of a vehicle route. - int64_t Start(int vehicle) const { return starts_[vehicle]; } + int64_t Start(int vehicle) const { return paths_metadata_.Starts()[vehicle]; } /// Returns the variable index of the ending node of a vehicle route. - int64_t End(int vehicle) const { return ends_[vehicle]; } + int64_t End(int vehicle) const { return paths_metadata_.Ends()[vehicle]; } /// Returns true if 'index' represents the first node of a route. - bool IsStart(int64_t index) const; + bool IsStart(int64_t index) const { return paths_metadata_.IsStart(index); } /// Returns true if 'index' represents the last node of a route. - bool IsEnd(int64_t index) const { return index >= Size(); } + bool IsEnd(int64_t index) const { return paths_metadata_.IsEnd(index); } /// Returns the vehicle of the given start/end index, and -1 if the given /// index is not a vehicle start/end. - int VehicleIndex(int64_t index) const { return index_to_vehicle_[index]; } + int VehicleIndex(int64_t index) const { + return paths_metadata_.GetPath(index); + } /// Assignment inspection /// Returns the variable index of the node directly after the node /// corresponding to 'index' in 'assignment'. @@ -2083,9 +2184,7 @@ class RoutingModel { // Two indices are equivalent if they correspond to the same node (as given // to the constructors taking a RoutingIndexManager). std::vector index_to_equivalence_class_; - std::vector index_to_vehicle_; - std::vector starts_; - std::vector ends_; + const PathsMetadata paths_metadata_; // TODO(user): b/62478706 Once the port is done, this shouldn't be needed // anymore. RoutingIndexManager manager_; @@ -2302,7 +2401,6 @@ class GlobalVehicleBreaksConstraint : public Constraint { private: void PropagateNode(int node); void PropagateVehicle(int vehicle); - void PropagateMaxBreakDistance(int vehicle); const RoutingModel* model_; const RoutingDimension* const dimension_; @@ -2833,7 +2931,6 @@ class RoutingDimension { std::vector breaks, int vehicle, std::vector node_visit_transits, std::function delays); -#endif /// !defined(SWIGPYTHON) /// Returns the break intervals set by SetBreakIntervalsOfVehicle(). const std::vector& GetBreakIntervalsOfVehicle( @@ -2844,6 +2941,7 @@ class RoutingDimension { const std::vector >& GetBreakDistanceDurationOfVehicle(int vehicle) const; // clang-format on +#endif /// !defined(SWIGPYTHON) int GetPreTravelEvaluatorOfVehicle(int vehicle) const; int GetPostTravelEvaluatorOfVehicle(int vehicle) const; diff --git a/ortools/constraint_solver/routing_filters.cc b/ortools/constraint_solver/routing_filters.cc index caad525d96..6390df0426 100644 --- a/ortools/constraint_solver/routing_filters.cc +++ b/ortools/constraint_solver/routing_filters.cc @@ -1213,11 +1213,7 @@ class PathCumulFilter : public BasePathFilter { const bool can_use_lp_; const bool propagate_own_objective_value_; - // Used to do span lower bounding in presence of vehicle breaks. - DisjunctivePropagator disjunctive_propagator_; - DisjunctivePropagator::Tasks tasks_; - TravelBounds travel_bounds_; - std::vector current_path_; + std::vector min_path_cumuls_; }; PathCumulFilter::PathCumulFilter(const RoutingModel& routing_model, @@ -1428,9 +1424,8 @@ void PathCumulFilter::OnBeforeSynchronizePaths() { // Second pass: update cumul, transit and cost values. node = Start(r); int64_t cumul = cumuls_[node]->Min(); - std::vector min_path_cumuls; - min_path_cumuls.reserve(number_of_route_arcs + 1); - min_path_cumuls.push_back(cumul); + min_path_cumuls_.clear(); + min_path_cumuls_.push_back(cumul); int64_t current_cumul_cost_value = GetCumulSoftCost(node, cumul); current_cumul_cost_value = CapAdd( @@ -1447,7 +1442,7 @@ void PathCumulFilter::OnBeforeSynchronizePaths() { cumul = dimension_.GetFirstPossibleGreaterOrEqualValueForNode(next, cumul); cumul = std::max(cumuls_[next]->Min(), cumul); - min_path_cumuls.push_back(cumul); + min_path_cumuls_.push_back(cumul); node = next; current_cumul_cost_value = CapAdd(current_cumul_cost_value, GetCumulSoftCost(node, cumul)); @@ -1455,7 +1450,7 @@ void PathCumulFilter::OnBeforeSynchronizePaths() { current_cumul_cost_value, GetCumulPiecewiseLinearCost(node, cumul)); } if (FilterPrecedences()) { - StoreMinMaxCumulOfNodesOnPath(/*path=*/r, min_path_cumuls, + StoreMinMaxCumulOfNodesOnPath(/*path=*/r, min_path_cumuls_, /*is_delta=*/false); } if (number_of_route_arcs == 1 && @@ -1597,9 +1592,8 @@ bool PathCumulFilter::AcceptPath(int64_t path_start, int64_t /*chain_start*/, DCHECK_NE(node, kUnassigned); } delta_path_transits_.ReserveTransits(path, number_of_route_arcs); - std::vector min_path_cumuls; - min_path_cumuls.reserve(number_of_route_arcs + 1); - min_path_cumuls.push_back(cumul); + min_path_cumuls_.clear(); + min_path_cumuls_.push_back(cumul); // Check that the path is feasible with regards to cumul bounds, scanning // the paths from start to end (caching path node sequences and transits // for further span cost filtering). @@ -1616,7 +1610,7 @@ bool PathCumulFilter::AcceptPath(int64_t path_start, int64_t /*chain_start*/, return false; } cumul = std::max(cumuls_[next]->Min(), cumul); - min_path_cumuls.push_back(cumul); + min_path_cumuls_.push_back(cumul); node = next; if (filter_vehicle_costs) { cumul_cost_delta = @@ -1628,7 +1622,7 @@ bool PathCumulFilter::AcceptPath(int64_t path_start, int64_t /*chain_start*/, const int64_t min_end = cumul; if (!PickupToDeliveryLimitsRespected(delta_path_transits_, path, - min_path_cumuls)) { + min_path_cumuls_)) { return false; } if (FilterSlackCost() || FilterBreakCost(vehicle) || @@ -1712,7 +1706,7 @@ bool PathCumulFilter::AcceptPath(int64_t path_start, int64_t /*chain_start*/, GetPathCumulSoftLowerBoundCost(delta_path_transits_, path)); } if (FilterPrecedences()) { - StoreMinMaxCumulOfNodesOnPath(path, min_path_cumuls, /*is_delta=*/true); + StoreMinMaxCumulOfNodesOnPath(path, min_path_cumuls_, /*is_delta=*/true); } if (!filter_vehicle_costs) { // If this route's costs shouldn't be taken into account, reset the @@ -2093,15 +2087,15 @@ void AppendLightWeightDimensionFilters( const PathState* path_state, const std::vector& dimensions, std::vector* filters) { + using Interval = DimensionChecker::Interval; // For every dimension that fits, add a DimensionChecker. for (const RoutingDimension* dimension : dimensions) { // Skip dimension if not unary. if (dimension->GetUnaryTransitEvaluator(0) == nullptr) continue; - using Intervals = std::vector; // Fill path capacities and classes. const int num_vehicles = dimension->model()->vehicles(); - Intervals path_capacity(num_vehicles); + std::vector path_capacity(num_vehicles); std::vector path_class(num_vehicles); for (int v = 0; v < num_vehicles; ++v) { const auto& vehicle_capacities = dimension->vehicle_capacities(); @@ -2117,15 +2111,14 @@ void AppendLightWeightDimensionFilters( 1 + *std::max_element(path_class.begin(), path_class.end()); const int num_cumuls = dimension->cumuls().size(); const int num_slacks = dimension->slacks().size(); - using Interval = DimensionChecker::Interval; - std::vector> transits(num_vehicle_classes, - nullptr); + std::vector> transits( + num_vehicle_classes, nullptr); for (int vehicle = 0; vehicle < num_vehicles; ++vehicle) { const int vehicle_class = path_class[vehicle]; if (transits[vehicle_class] != nullptr) continue; const auto& evaluator = dimension->GetUnaryTransitEvaluator(vehicle); - transits[vehicle_class] = [&evaluator, dimension, - num_slacks](int64_t node) -> Interval { + transits[vehicle_class] = [&evaluator, dimension, num_slacks]( + int64_t node, int64_t next) -> Interval { if (node >= num_slacks) return {0, 0}; const int64_t min_transit = evaluator(node); const int64_t max_transit = @@ -2134,7 +2127,7 @@ void AppendLightWeightDimensionFilters( }; } // Fill node capacities. - Intervals node_capacity(num_cumuls); + std::vector node_capacity(num_cumuls); for (int node = 0; node < num_cumuls; ++node) { const IntVar* cumul = dimension->CumulVar(node); node_capacity[node] = {cumul->Min(), cumul->Max()}; @@ -2642,10 +2635,11 @@ bool LPCumulFilter::Accept(const Assignment* delta, // No need to compute the cost of the LP, only verify its feasibility. delta_cost_without_transit_ = 0; const DimensionSchedulingStatus status = - optimizer_.ComputeCumuls(next_accessor, nullptr, nullptr, nullptr); + optimizer_.ComputeCumuls(next_accessor, {}, nullptr, nullptr, nullptr); if (status == DimensionSchedulingStatus::OPTIMAL) return true; if (status == DimensionSchedulingStatus::RELAXED_OPTIMAL_ONLY && - mp_optimizer_.ComputeCumuls(next_accessor, nullptr, nullptr, nullptr) == + mp_optimizer_.ComputeCumuls(next_accessor, {}, nullptr, nullptr, + nullptr) == DimensionSchedulingStatus::OPTIMAL) { return true; } @@ -2727,7 +2721,7 @@ class ResourceGroupAssignmentFilter : public BasePathFilter { public: ResourceGroupAssignmentFilter(const std::vector& nexts, const ResourceGroup* resource_group, - LocalDimensionCumulOptimizer* optimizer, + LocalDimensionCumulOptimizer* lp_optimizer, LocalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost); bool InitializeAcceptPath() override; @@ -2746,15 +2740,15 @@ class ResourceGroupAssignmentFilter : public BasePathFilter { return synchronized_cost_without_transit_; } std::string DebugString() const override { - return "ResourceGroupAssignmentFilter(" + - resource_assignment_optimizer_.dimension()->name() + ")"; + return "ResourceGroupAssignmentFilter(" + dimension_.name() + ")"; } private: - ResourceAssignmentOptimizer resource_assignment_optimizer_; const RoutingModel& model_; const RoutingDimension& dimension_; const ResourceGroup& resource_group_; + LocalDimensionCumulOptimizer* lp_optimizer_; + LocalDimensionCumulOptimizer* mp_optimizer_; const bool filter_objective_cost_; bool current_synch_failed_; int64_t synchronized_cost_without_transit_; @@ -2765,13 +2759,14 @@ class ResourceGroupAssignmentFilter : public BasePathFilter { ResourceGroupAssignmentFilter::ResourceGroupAssignmentFilter( const std::vector& nexts, const ResourceGroup* resource_group, - LocalDimensionCumulOptimizer* optimizer, + LocalDimensionCumulOptimizer* lp_optimizer, LocalDimensionCumulOptimizer* mp_optimizer, bool filter_objective_cost) - : BasePathFilter(nexts, optimizer->dimension()->cumuls().size()), - resource_assignment_optimizer_(resource_group, optimizer, mp_optimizer), - model_(*optimizer->dimension()->model()), - dimension_(*optimizer->dimension()), + : BasePathFilter(nexts, lp_optimizer->dimension()->cumuls().size()), + model_(*lp_optimizer->dimension()->model()), + dimension_(*lp_optimizer->dimension()), resource_group_(*resource_group), + lp_optimizer_(lp_optimizer), + mp_optimizer_(mp_optimizer), filter_objective_cost_(filter_objective_cost), current_synch_failed_(false), synchronized_cost_without_transit_(-1), @@ -2801,19 +2796,25 @@ bool ResourceGroupAssignmentFilter::AcceptPath(int64_t path_start, int64_t /*chain_start*/, int64_t /*chain_end*/) { const int vehicle = model_.VehicleIndex(path_start); - return resource_assignment_optimizer_.ComputeAssignmentCostsForVehicle( - vehicle, [this](int64_t index) { return GetNext(index); }, + return ComputeVehicleToResourcesAssignmentCosts( + vehicle, resource_group_, + [this](int64_t index) { return GetNext(index); }, dimension_.transit_evaluator(vehicle), filter_objective_cost_, + lp_optimizer_, mp_optimizer_, &delta_vehicle_to_resource_assignment_costs_[vehicle], nullptr, nullptr); } bool ResourceGroupAssignmentFilter::FinalizeAcceptPath( int64_t /*objective_min*/, int64_t objective_max) { - delta_cost_without_transit_ = - resource_assignment_optimizer_.ComputeBestAssignmentCost( - delta_vehicle_to_resource_assignment_costs_, - vehicle_to_resource_assignment_costs_, - [this](int v) { return PathStartTouched(model_.Start(v)); }, nullptr); + delta_cost_without_transit_ = ComputeBestVehicleToResourceAssignment( + resource_group_.GetVehiclesRequiringAResource(), resource_group_.Size(), + /*vehicle_to_resource_assignment_costs=*/ + [this](int v) { + return PathStartTouched(model_.Start(v)) + ? &delta_vehicle_to_resource_assignment_costs_[v] + : &vehicle_to_resource_assignment_costs_[v]; + }, + nullptr); return delta_cost_without_transit_ >= 0 && delta_cost_without_transit_ <= objective_max; } @@ -2824,19 +2825,19 @@ void ResourceGroupAssignmentFilter::OnBeforeSynchronizePaths() { void ResourceGroupAssignmentFilter::OnSynchronizePathFromStart(int64_t start) { // NOTE(user): Even if filter_objective_cost_ is false, we still need to - // call ComputeAssignmentCostsForVehicle() for every vehicle to keep track - // of whether or not a given vehicle-to-resource assignment is possible by - // storing 0 or -1 in vehicle_to_resource_assignment_costs_. + // call ComputeVehicleToResourcesAssignmentCosts() for every vehicle to keep + // track of whether or not a given vehicle-to-resource assignment is possible + // by storing 0 or -1 in vehicle_to_resource_assignment_costs_. const auto& next_accessor = [this](int64_t index) { return IsVarSynced(index) ? Value(index) : model_.IsStart(index) ? model_.End(model_.VehicleIndex(index)) : index; }; const int v = model_.VehicleIndex(start); - if (!resource_assignment_optimizer_.ComputeAssignmentCostsForVehicle( - v, next_accessor, dimension_.transit_evaluator(v), - filter_objective_cost_, &vehicle_to_resource_assignment_costs_[v], - nullptr, nullptr)) { + if (!ComputeVehicleToResourcesAssignmentCosts( + v, resource_group_, next_accessor, dimension_.transit_evaluator(v), + filter_objective_cost_, lp_optimizer_, mp_optimizer_, + &vehicle_to_resource_assignment_costs_[v], nullptr, nullptr)) { vehicle_to_resource_assignment_costs_[v].assign(resource_group_.Size(), -1); current_synch_failed_ = true; } @@ -2846,9 +2847,12 @@ void ResourceGroupAssignmentFilter::OnAfterSynchronizePaths() { synchronized_cost_without_transit_ = (current_synch_failed_ || !filter_objective_cost_) ? 0 - : resource_assignment_optimizer_.ComputeBestAssignmentCost( - vehicle_to_resource_assignment_costs_, - vehicle_to_resource_assignment_costs_, [](int) { return true; }, + : ComputeBestVehicleToResourceAssignment( + resource_group_.GetVehiclesRequiringAResource(), + resource_group_.Size(), + [this](int v) { + return &vehicle_to_resource_assignment_costs_[v]; + }, nullptr); synchronized_cost_without_transit_ = std::max(synchronized_cost_without_transit_, 0); diff --git a/ortools/constraint_solver/routing_lp_scheduling.cc b/ortools/constraint_solver/routing_lp_scheduling.cc index 8deca00cbb..2021a211a7 100644 --- a/ortools/constraint_solver/routing_lp_scheduling.cc +++ b/ortools/constraint_solver/routing_lp_scheduling.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -32,6 +33,9 @@ #include "absl/strings/str_format.h" #include "absl/strings/str_join.h" #include "absl/time/time.h" +#include "ortools/base/check.h" +#include "ortools/base/dump_vars.h" +#include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/constraint_solver/constraint_solver.h" @@ -41,6 +45,7 @@ #include "ortools/graph/min_cost_flow.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/lp_utils.h" +#include "ortools/util/flat_matrix.h" #include "ortools/util/saturated_arithmetic.h" #include "ortools/util/sorted_interval_list.h" @@ -205,7 +210,7 @@ LocalDimensionCumulOptimizer::LocalDimensionCumulOptimizer( DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumulCost( int vehicle, const std::function& next_accessor, int64_t* optimal_cost) { - return optimizer_core_.OptimizeSingleRoute(vehicle, next_accessor, + return optimizer_core_.OptimizeSingleRoute(vehicle, next_accessor, {}, solver_[vehicle].get(), nullptr, nullptr, optimal_cost, nullptr); } @@ -217,8 +222,8 @@ LocalDimensionCumulOptimizer::ComputeRouteCumulCostWithoutFixedTransits( int64_t cost = 0; int64_t transit_cost = 0; const DimensionSchedulingStatus status = optimizer_core_.OptimizeSingleRoute( - vehicle, next_accessor, solver_[vehicle].get(), nullptr, nullptr, &cost, - &transit_cost); + vehicle, next_accessor, {}, solver_[vehicle].get(), nullptr, nullptr, + &cost, &transit_cost); if (status != DimensionSchedulingStatus::INFEASIBLE && optimal_cost_without_transits != nullptr) { *optimal_cost_without_transits = CapSub(cost, transit_cost); @@ -236,28 +241,57 @@ std::vector LocalDimensionCumulOptimizer:: std::vector>* optimal_cumuls, std::vector>* optimal_breaks) { return optimizer_core_.OptimizeSingleRouteWithResources( - vehicle, next_accessor, transit_accessor, resources, resource_indices, + vehicle, next_accessor, transit_accessor, {}, resources, resource_indices, optimize_vehicle_costs, solver_[vehicle].get(), optimal_costs_without_transits, optimal_cumuls, optimal_breaks); } DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputeRouteCumuls( int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, std::vector* optimal_cumuls, std::vector* optimal_breaks) { return optimizer_core_.OptimizeSingleRoute( - vehicle, next_accessor, solver_[vehicle].get(), optimal_cumuls, - optimal_breaks, nullptr, nullptr); + vehicle, next_accessor, dimension_travel_info, solver_[vehicle].get(), + optimal_cumuls, optimal_breaks, nullptr, nullptr); +} + +DimensionSchedulingStatus +LocalDimensionCumulOptimizer::ComputeRouteCumulsAndCost( + int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, + std::vector* optimal_cumuls, std::vector* optimal_breaks, + int64_t* optimal_cost) { + return optimizer_core_.OptimizeSingleRoute( + vehicle, next_accessor, dimension_travel_info, solver_[vehicle].get(), + optimal_cumuls, optimal_breaks, optimal_cost, nullptr); +} + +DimensionSchedulingStatus +LocalDimensionCumulOptimizer::ComputeRouteSolutionCost( + int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, + const std::vector& solution_cumul_values, + const std::vector& solution_break_values, int64_t* solution_cost, + int64_t* cost_offset, bool reuse_previous_model_if_possible, bool clear_lp, + absl::Duration* solve_duration) { + RoutingLinearSolverWrapper* solver = solver_[vehicle].get(); + return optimizer_core_.ComputeSingleRouteSolutionCost( + vehicle, next_accessor, dimension_travel_info, solver, + solution_cumul_values, solution_break_values, solution_cost, nullptr, + cost_offset, reuse_previous_model_if_possible, clear_lp, + /*clear_solution_constraints=*/true, solve_duration); } DimensionSchedulingStatus LocalDimensionCumulOptimizer::ComputePackedRouteCumuls( int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, const RoutingModel::ResourceGroup::Resource* resource, std::vector* packed_cumuls, std::vector* packed_breaks) { return optimizer_core_.OptimizeAndPackSingleRoute( - vehicle, next_accessor, resource, solver_[vehicle].get(), packed_cumuls, - packed_breaks); + vehicle, next_accessor, dimension_travel_info, resource, + solver_[vehicle].get(), packed_cumuls, packed_breaks); } const int CumulBoundsPropagator::kNoParent = -2; @@ -286,8 +320,9 @@ void CumulBoundsPropagator::AddArcs(int first_index, int second_index, } bool CumulBoundsPropagator::InitializeArcsAndBounds( - const std::function& next_accessor, - int64_t cumul_offset) { + const std::function& next_accessor, int64_t cumul_offset, + const std::vector* + dimension_travel_info_per_route) { propagated_bounds_.assign(num_nodes_, std::numeric_limits::min()); for (std::vector& arcs : outgoing_arcs_) { @@ -302,6 +337,7 @@ bool CumulBoundsPropagator::InitializeArcsAndBounds( dimension_.transit_evaluator(vehicle); int node = model->Start(vehicle); + int index_on_route = 0; while (true) { int64_t cumul_lb, cumul_ub; if (!GetCumulBoundsWithOffset(dimension_, node, cumul_offset, &cumul_lb, @@ -318,7 +354,17 @@ bool CumulBoundsPropagator::InitializeArcsAndBounds( } const int next = next_accessor(node); - const int64_t transit = transit_accessor(node, next); + int64_t transit = transit_accessor(node, next); + if (dimension_travel_info_per_route != nullptr && + !dimension_travel_info_per_route->empty()) { + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo& + transition_info = (*dimension_travel_info_per_route)[vehicle] + .transition_info[index_on_route]; + transit = transition_info.compressed_travel_value_lower_bound + + transition_info.pre_travel_transit_value + + transition_info.post_travel_transit_value; + ++index_on_route; + } const IntVar& slack_var = *dimension_.SlackVar(node); // node + transit + slack_var == next // Add arcs for node + transit + slack_min <= next @@ -427,14 +473,16 @@ bool CumulBoundsPropagator::DisassembleSubtree(int source, int target) { } bool CumulBoundsPropagator::PropagateCumulBounds( - const std::function& next_accessor, - int64_t cumul_offset) { + const std::function& next_accessor, int64_t cumul_offset, + const std::vector* + dimension_travel_info_per_route) { tree_parent_node_of_.assign(num_nodes_, kNoParent); DCHECK(std::none_of(node_in_queue_.begin(), node_in_queue_.end(), [](bool b) { return b; })); DCHECK(bf_queue_.empty()); - if (!InitializeArcsAndBounds(next_accessor, cumul_offset)) { + if (!InitializeArcsAndBounds(next_accessor, cumul_offset, + dimension_travel_info_per_route)) { return CleanupAndReturnFalse(); } @@ -510,35 +558,138 @@ DimensionCumulOptimizerCore::DimensionCumulOptimizerCore( } } -DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeSingleRoute( +bool DimensionCumulOptimizerCore::InitSingleRoute( int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, - std::vector* break_values, int64_t* cost, int64_t* transit_cost, - bool clear_lp) { + int64_t* cost, int64_t* transit_cost, int64_t* cumul_offset, + int64_t* const cost_offset) { InitOptimizer(solver); // Make sure SetRouteCumulConstraints will properly set the cumul bounds by // looking at this route only. DCHECK_EQ(propagator_.get(), nullptr); - RoutingModel* const model = dimension()->model(); + RoutingModel* model = dimension_->model(); const bool optimize_vehicle_costs = (cumul_values != nullptr || cost != nullptr) && (!model->IsEnd(next_accessor(model->Start(vehicle))) || model->IsVehicleUsedWhenEmpty(vehicle)); - const int64_t cumul_offset = - dimension_->GetLocalOptimizerOffsetForVehicle(vehicle); - int64_t cost_offset = 0; - if (!SetRouteCumulConstraints(vehicle, next_accessor, - dimension_->transit_evaluator(vehicle), - cumul_offset, optimize_vehicle_costs, solver, - transit_cost, &cost_offset)) { - return DimensionSchedulingStatus::INFEASIBLE; + *cumul_offset = dimension_->GetLocalOptimizerOffsetForVehicle(vehicle); + if (!SetRouteCumulConstraints( + vehicle, next_accessor, dimension_->transit_evaluator(vehicle), + dimension_travel_info, *cumul_offset, optimize_vehicle_costs, solver, + transit_cost, cost_offset)) { + return false; } if (model->CheckLimit()) { + return false; + } + return true; +} + +DimensionSchedulingStatus +DimensionCumulOptimizerCore::ComputeSingleRouteSolutionCost( + int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, + RoutingLinearSolverWrapper* solver, + const std::vector& solution_cumul_values, + const std::vector& solution_break_values, int64_t* cost, + int64_t* transit_cost, int64_t* cost_offset, + bool reuse_previous_model_if_possible, bool clear_lp, + bool clear_solution_constraints, absl::Duration* const solve_duration) { + absl::Duration solve_duration_value; + int64_t cost_offset_value; + if (!reuse_previous_model_if_possible || solver->ModelIsEmpty()) { + int64_t cumul_offset; + std::vector cumul_values; + if (!InitSingleRoute(vehicle, next_accessor, dimension_travel_info, solver, + &cumul_values, cost, transit_cost, &cumul_offset, + &cost_offset_value)) { + return DimensionSchedulingStatus::INFEASIBLE; + } + solve_duration_value = dimension_->model()->RemainingTime(); + if (solve_duration != nullptr) *solve_duration = solve_duration_value; + if (cost_offset != nullptr) *cost_offset = cost_offset_value; + } else { + CHECK(cost_offset != nullptr) + << "Cannot reuse model without the cost_offset"; + cost_offset_value = *cost_offset; + CHECK(solve_duration != nullptr) + << "Cannot reuse model without the solve_duration"; + solve_duration_value = *solve_duration; + } + + // Constrains the cumuls. + DCHECK_EQ(solution_cumul_values.size(), + current_route_cumul_variables_.size()); + for (int i = 0; i < current_route_cumul_variables_.size(); ++i) { + if (solution_cumul_values[i] < current_route_min_cumuls_[i] || + solution_cumul_values[i] > current_route_max_cumuls_[i]) { + return DimensionSchedulingStatus::INFEASIBLE; + } + solver->SetVariableBounds(current_route_cumul_variables_[i], + /*lower_bound=*/solution_cumul_values[i], + /*upper_bound=*/solution_cumul_values[i]); + } + + // Constrains the breaks. + DCHECK_EQ(solution_break_values.size(), + current_route_break_variables_.size()); + std::vector current_route_min_breaks( + current_route_break_variables_.size()); + std::vector current_route_max_breaks( + current_route_break_variables_.size()); + for (int i = 0; i < current_route_break_variables_.size(); ++i) { + current_route_min_breaks[i] = + solver->GetVariableLowerBound(current_route_break_variables_[i]); + current_route_max_breaks[i] = + solver->GetVariableUpperBound(current_route_break_variables_[i]); + solver->SetVariableBounds(current_route_break_variables_[i], + /*lower_bound=*/solution_break_values[i], + /*upper_bound=*/solution_break_values[i]); + } + + const DimensionSchedulingStatus status = solver->Solve(solve_duration_value); + if (status == DimensionSchedulingStatus::INFEASIBLE) { + solver->Clear(); + return status; + } + + if (cost != nullptr) { + *cost = CapAdd(cost_offset_value, solver->GetObjectiveValue()); + } + + if (clear_lp) { + solver->Clear(); + } else if (clear_solution_constraints) { + for (int i = 0; i < current_route_cumul_variables_.size(); ++i) { + solver->SetVariableBounds(current_route_cumul_variables_[i], + /*lower_bound=*/current_route_min_cumuls_[i], + /*upper_bound=*/current_route_max_cumuls_[i]); + } + for (int i = 0; i < current_route_break_variables_.size(); ++i) { + solver->SetVariableBounds(current_route_break_variables_[i], + /*lower_bound=*/current_route_min_breaks[i], + /*upper_bound=*/current_route_max_breaks[i]); + } + } + return status; +} + +DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeSingleRoute( + int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, + RoutingLinearSolverWrapper* solver, std::vector* cumul_values, + std::vector* break_values, int64_t* cost, int64_t* transit_cost, + bool clear_lp) { + int64_t cumul_offset, cost_offset; + if (!InitSingleRoute(vehicle, next_accessor, dimension_travel_info, solver, + cumul_values, cost, transit_cost, &cumul_offset, + &cost_offset)) { return DimensionSchedulingStatus::INFEASIBLE; } const DimensionSchedulingStatus status = - solver->Solve(model->RemainingTime()); + solver->Solve(dimension()->model()->RemainingTime()); if (status == DimensionSchedulingStatus::INFEASIBLE) { solver->Clear(); return status; @@ -624,6 +775,7 @@ std::vector DimensionCumulOptimizerCore::OptimizeSingleRouteWithResources( int vehicle, const std::function& next_accessor, const std::function& transit_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, const std::vector& resources, const std::vector& resource_indices, bool optimize_vehicle_costs, RoutingLinearSolverWrapper* solver, @@ -651,8 +803,9 @@ DimensionCumulOptimizerCore::OptimizeSingleRouteWithResources( int64_t cost_offset = 0; int64_t transit_cost = 0; if (!SetRouteCumulConstraints(vehicle, next_accessor, transit_accessor, - cumul_offset, optimize_vehicle_costs, solver, - &transit_cost, &cost_offset)) { + dimension_travel_info, cumul_offset, + optimize_vehicle_costs, solver, &transit_cost, + &cost_offset)) { return {DimensionSchedulingStatus::INFEASIBLE}; } @@ -719,6 +872,8 @@ DimensionCumulOptimizerCore::OptimizeSingleRouteWithResources( DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values, std::vector>* resource_indices_per_group, int64_t* cost, @@ -733,7 +888,8 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize( const int64_t cumul_offset = dimension_->GetGlobalOptimizerOffset(); if (propagator_ != nullptr && - !propagator_->PropagateCumulBounds(next_accessor, cumul_offset)) { + !propagator_->PropagateCumulBounds(next_accessor, cumul_offset, + &dimension_travel_info_per_route)) { return DimensionSchedulingStatus::INFEASIBLE; } @@ -747,10 +903,14 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize( !model->IsEnd(next_accessor(model->Start(vehicle))) || model->IsVehicleUsedWhenEmpty(vehicle); const bool optimize_vehicle_costs = optimize_costs && vehicle_is_used; - if (!SetRouteCumulConstraints(vehicle, next_accessor, - dimension_->transit_evaluator(vehicle), - cumul_offset, optimize_vehicle_costs, solver, - &route_transit_cost, &route_cost_offset)) { + const RouteDimensionTravelInfo& dimension_travel_info = + dimension_travel_info_per_route.empty() + ? RouteDimensionTravelInfo() + : dimension_travel_info_per_route[vehicle]; + if (!SetRouteCumulConstraints( + vehicle, next_accessor, dimension_->transit_evaluator(vehicle), + dimension_travel_info, cumul_offset, optimize_vehicle_costs, solver, + &route_transit_cost, &route_cost_offset)) { return DimensionSchedulingStatus::INFEASIBLE; } total_transit_cost = CapAdd(total_transit_cost, route_transit_cost); @@ -791,6 +951,8 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::Optimize( DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPack( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values, std::vector>* resource_indices_per_group) { @@ -806,8 +968,8 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPack( solver->SetParameters(packing_parameters.SerializeAsString()); } DimensionSchedulingStatus status = DimensionSchedulingStatus::OPTIMAL; - if (Optimize(next_accessor, solver, /*cumul_values=*/nullptr, - /*break_values=*/nullptr, + if (Optimize(next_accessor, dimension_travel_info_per_route, solver, + /*cumul_values=*/nullptr, /*break_values=*/nullptr, /*resource_indices_per_group=*/nullptr, &cost, /*transit_cost=*/nullptr, /*clear_lp=*/false) == DimensionSchedulingStatus::INFEASIBLE) { @@ -840,6 +1002,7 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPack( DimensionSchedulingStatus DimensionCumulOptimizerCore::OptimizeAndPackSingleRoute( int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, const RoutingModel::ResourceGroup::Resource* resource, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values) { @@ -856,8 +1019,8 @@ DimensionCumulOptimizerCore::OptimizeAndPackSingleRoute( // Note: We pass a non-nullptr cost to the OptimizeSingleRoute() method so // the costs are optimized by the LP. int64_t cost = 0; - if (OptimizeSingleRoute(vehicle, next_accessor, solver, - /*cumul_values=*/nullptr, + if (OptimizeSingleRoute(vehicle, next_accessor, dimension_travel_info, + solver, /*cumul_values=*/nullptr, /*break_values=*/nullptr, &cost, /*transit_cost=*/nullptr, /*clear_lp=*/false) == DimensionSchedulingStatus::INFEASIBLE) { @@ -868,8 +1031,9 @@ DimensionCumulOptimizerCore::OptimizeAndPackSingleRoute( const std::vector statuses = OptimizeSingleRouteWithResources( vehicle, next_accessor, dimension_->transit_evaluator(vehicle), - {*resource}, {0}, /*optimize_vehicle_costs=*/true, solver, - &costs_without_transits, /*cumul_values=*/nullptr, + dimension_travel_info, {*resource}, {0}, + /*optimize_vehicle_costs=*/true, solver, &costs_without_transits, + /*cumul_values=*/nullptr, /*break_values=*/nullptr, /*clear_lp=*/false); if (dimension_->model()->CheckLimit()) { status = DimensionSchedulingStatus::INFEASIBLE; @@ -967,33 +1131,44 @@ DimensionSchedulingStatus DimensionCumulOptimizerCore::PackRoutes( return status; } -#ifndef NDEBUG -#define SET_VARIABLE_NAME(solver, var, name) \ - do { \ - solver->SetVariableName(var, name); \ +#define SET_DEBUG_VARIABLE_NAME(solver, var, name) \ + do { \ + if (DEBUG_MODE) { \ + solver->SetVariableName(var, name); \ + } \ } while (false) -#else -#define SET_VARIABLE_NAME(solver, var, name) \ - do { \ - } while (false) -#endif void DimensionCumulOptimizerCore::InitOptimizer( RoutingLinearSolverWrapper* solver) { solver->Clear(); index_to_cumul_variable_.assign(dimension_->cumuls().size(), -1); max_end_cumul_ = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, max_end_cumul_, "max_end_cumul"); + SET_DEBUG_VARIABLE_NAME(solver, max_end_cumul_, "max_end_cumul"); min_start_cumul_ = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, min_start_cumul_, "min_start_cumul"); + SET_DEBUG_VARIABLE_NAME(solver, min_start_cumul_, "min_start_cumul"); } -bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds( - const std::vector& route, - const std::vector& fixed_transits, int64_t cumul_offset) { +bool DimensionCumulOptimizerCore::ExtractRouteCumulBounds( + const std::vector& route, int64_t cumul_offset) { const int route_size = route.size(); current_route_min_cumuls_.resize(route_size); current_route_max_cumuls_.resize(route_size); + + // Extract cumul min/max and fixed transits from CP. + for (int pos = 0; pos < route_size; ++pos) { + if (!GetCumulBoundsWithOffset(*dimension_, route[pos], cumul_offset, + ¤t_route_min_cumuls_[pos], + ¤t_route_max_cumuls_[pos])) { + return false; + } + } + return true; +} + +bool DimensionCumulOptimizerCore::TightenRouteCumulBounds( + const std::vector& route, const std::vector& min_transits, + int64_t cumul_offset) { + const int route_size = route.size(); if (propagator_ != nullptr) { for (int pos = 0; pos < route_size; pos++) { const int64_t node = route[pos]; @@ -1005,15 +1180,6 @@ bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds( return true; } - // Extract cumul min/max and fixed transits from CP. - for (int pos = 0; pos < route_size; ++pos) { - if (!GetCumulBoundsWithOffset(*dimension_, route[pos], cumul_offset, - ¤t_route_min_cumuls_[pos], - ¤t_route_max_cumuls_[pos])) { - return false; - } - } - // Refine cumul bounds using // cumul[i+1] >= cumul[i] + fixed_transit[i] + slack[i]. for (int pos = 1; pos < route_size; ++pos) { @@ -1021,7 +1187,7 @@ bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds( current_route_min_cumuls_[pos] = std::max( current_route_min_cumuls_[pos], CapAdd( - CapAdd(current_route_min_cumuls_[pos - 1], fixed_transits[pos - 1]), + CapAdd(current_route_min_cumuls_[pos - 1], min_transits[pos - 1]), slack_min)); current_route_min_cumuls_[pos] = GetFirstPossibleValueForCumulWithOffset( *dimension_, route[pos], current_route_min_cumuls_[pos], cumul_offset); @@ -1038,9 +1204,8 @@ bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds( const int64_t slack_min = dimension_->SlackVar(route[pos])->Min(); current_route_max_cumuls_[pos] = std::min( current_route_max_cumuls_[pos], - CapSub( - CapSub(current_route_max_cumuls_[pos + 1], fixed_transits[pos]), - slack_min)); + CapSub(CapSub(current_route_max_cumuls_[pos + 1], min_transits[pos]), + slack_min)); current_route_max_cumuls_[pos] = GetLastPossibleValueForCumulWithOffset( *dimension_, route[pos], current_route_max_cumuls_[pos], cumul_offset); @@ -1052,12 +1217,447 @@ bool DimensionCumulOptimizerCore::ComputeRouteCumulBounds( return true; } +std::vector SlopeAndYInterceptToConvexityRegions( + const std::vector& slope_and_y_intercept) { + CHECK(!slope_and_y_intercept.empty()); + std::vector convex(slope_and_y_intercept.size(), false); + double previous_slope = std::numeric_limits::max(); + for (int i = 0; i < slope_and_y_intercept.size(); ++i) { + const auto& pair = slope_and_y_intercept[i]; + if (pair.slope < previous_slope) { + convex[i] = true; + } + previous_slope = pair.slope; + } + return convex; +} + +std::vector PiecewiseLinearFormulationToSlopeAndYIntercept( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl_function, + int index_start, int index_end) { + if (index_end < 0) index_end = pwl_function.x_anchors.size() - 1; + const int num_segments = index_end - index_start; + DCHECK_GE(num_segments, 1); + std::vector slope_and_y_intercept(num_segments); + for (int seg = index_start; seg < index_end; ++seg) { + auto& [slope, y_intercept] = slope_and_y_intercept[seg - index_start]; + slope = (pwl_function.y_anchors[seg + 1] - pwl_function.y_anchors[seg]) / + static_cast(pwl_function.x_anchors[seg + 1] - + pwl_function.x_anchors[seg]); + y_intercept = + pwl_function.y_anchors[seg] - slope * pwl_function.x_anchors[seg]; + } + return slope_and_y_intercept; +} + +namespace { + +// Find a "good" scaling factor for constraints with non-integers coefficients. +// See sat::FindBestScalingAndComputeErrors() for more infos. +double FindBestScaling(const std::vector& coefficients, + const std::vector& lower_bounds, + const std::vector& upper_bounds, + int64_t max_absolute_activity, + double wanted_absolute_activity_precision) { + double unused_relative_coeff_error = 0; + double unused_scaled_sum_error = 0; + return sat::FindBestScalingAndComputeErrors( + coefficients, lower_bounds, upper_bounds, max_absolute_activity, + wanted_absolute_activity_precision, &unused_relative_coeff_error, + &unused_scaled_sum_error); +} + +// Returns the value of pwl(x) with pwl a PiecewiseLinearFormulation, knowing +// that x ∈ [pwl.x[upper_segment_index-1], pwl.x[upper_segment_index]]. +int64_t PieceWiseLinearFormulationValueKnownSegment( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl, + int64_t x, int upper_segment_index, double delta = 0) { + DCHECK_GE(upper_segment_index, 1); + DCHECK_LE(upper_segment_index, pwl.x_anchors.size() - 1); + const double alpha = + static_cast(pwl.y_anchors[upper_segment_index] - + pwl.y_anchors[upper_segment_index - 1]) / + (pwl.x_anchors[upper_segment_index] - + pwl.x_anchors[upper_segment_index - 1]); + const double beta = pwl.y_anchors[upper_segment_index] - + pwl.x_anchors[upper_segment_index] * alpha; + return std::ceil(alpha * x + beta + delta); +} + +} // namespace + +PiecewiseEvaluationStatus ComputePiecewiseLinearFormulationValue( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl, + int64_t x, int64_t* value, double delta) { + // Search for first element xi such that xi < x. + const auto upper_segment = + std::upper_bound(pwl.x_anchors.begin(), pwl.x_anchors.end(), x); + const int upper_segment_index = + std::distance(pwl.x_anchors.begin(), upper_segment); + + // Checking bounds + if (upper_segment_index == 0) { + return PiecewiseEvaluationStatus::SMALLER_THAN_LOWER_BOUND; + } else if (upper_segment == pwl.x_anchors.end()) { + if (x == pwl.x_anchors.back()) { + *value = std::ceil(pwl.y_anchors.back() + delta); + return PiecewiseEvaluationStatus::WITHIN_BOUNDS; + } + return PiecewiseEvaluationStatus::LARGER_THAN_UPPER_BOUND; + } + + *value = PieceWiseLinearFormulationValueKnownSegment( + pwl, x, upper_segment_index, delta); + return PiecewiseEvaluationStatus::WITHIN_BOUNDS; +} + +int64_t ComputeConvexPiecewiseLinearFormulationValue( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl, + int64_t x, double delta) { + int64_t y_value; + const PiecewiseEvaluationStatus status = + ComputePiecewiseLinearFormulationValue(pwl, x, &y_value, delta); + switch (status) { + case PiecewiseEvaluationStatus::UNSPECIFIED: + // The status should be specified. + LOG(FATAL) << "Unspecified PiecewiseEvaluationStatus."; + break; + case PiecewiseEvaluationStatus::WITHIN_BOUNDS: + // x is in the bounds, therefore, simply return the computed value. + return y_value; + case PiecewiseEvaluationStatus::SMALLER_THAN_LOWER_BOUND: + // In the convex case, if x <= lower_bound, the most restrictive + // constraint will be the first one. + return PieceWiseLinearFormulationValueKnownSegment(pwl, x, 1, delta); + case PiecewiseEvaluationStatus::LARGER_THAN_UPPER_BOUND: + // In the convex case, if x >= upper_bound, the most restrictive + // constraint will be the last one. + return PieceWiseLinearFormulationValueKnownSegment( + pwl, x, pwl.x_anchors.size() - 1, delta); + } +} + +bool DimensionCumulOptimizerCore::SetRouteTravelConstraints( + const RouteDimensionTravelInfo& dimension_travel_info, + const std::vector& lp_slacks, + const std::vector& fixed_transit, + RoutingLinearSolverWrapper* solver) { + const std::vector& lp_cumuls = current_route_cumul_variables_; + const int path_size = lp_cumuls.size(); + + if (dimension_travel_info.transition_info.empty()) { + // Travel is not travel-start dependent. + // Add all path constraints to LP: + // cumul[i] + fixed_transit[i] + slack[i] == cumul[i+1] + // <=> fixed_transit[i] == cumul[i+1] - cumul[i] - slack[i]. + for (int pos = 0; pos < path_size - 1; ++pos) { + const int ct = + solver->CreateNewConstraint(fixed_transit[pos], fixed_transit[pos]); + solver->SetCoefficient(ct, lp_cumuls[pos + 1], 1); + solver->SetCoefficient(ct, lp_cumuls[pos], -1); + solver->SetCoefficient(ct, lp_slacks[pos], -1); + } + return true; + } + + // See: + // https://docs.google.com/document/d/1niFfbCNjh70VepvVk4Ft8QgQAcrA5hrye-Vu_k0VgWw/edit?usp=sharing&resourcekey=0-rUJUHuWPjn1wWZaA0uDdtg + for (int pos = 0; pos < path_size - 1; ++pos) { + // Add a traffic-aware compression cost, for every path. + // compression_cost represents the cost of the ABSOLUTE compression of the + // travel. + const int compression_cost = solver->CreateNewPositiveVariable(); + SET_DEBUG_VARIABLE_NAME(solver, compression_cost, + absl::StrFormat("compression_cost(%ld)", pos)); + // relative_compression_cost represents the cost of the RELATIVE compression + // of the travel. This is the real cost used. In theory, + // relative_compression_cost = compression_cost / Tᵣ, where Tᵣ is the travel + // value (computed with the PWL). In practice, this requires a + // multiplication which is slow, so several approximations are implemented + // below. + const int relative_compression_cost = solver->CreateNewPositiveVariable(); + SET_DEBUG_VARIABLE_NAME( + solver, relative_compression_cost, + absl::StrFormat("relative_compression_cost(%ld)", pos)); + + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo& + transition_info = dimension_travel_info.transition_info[pos]; + const RouteDimensionTravelInfo::TransitionInfo::PiecewiseLinearFormulation& + travel_function = transition_info.travel_start_dependent_travel; + const int num_pwl_anchors = travel_function.x_anchors.size(); + DCHECK_GE(num_pwl_anchors, 2) + << "Travel value PWL must have at least 2 points"; + DCHECK_EQ(num_pwl_anchors, travel_function.y_anchors.size()) + << "Travel value PWL must have as many x anchors than y."; + + // 1. Create the travel value variable and set its constraints. + // 1.a. Create Variables for the start and value of a travel + const int64_t pre_travel_transit = transition_info.pre_travel_transit_value; + const int64_t post_travel_transit = + transition_info.post_travel_transit_value; + const int64_t compressed_travel_value_lower_bound = + transition_info.compressed_travel_value_lower_bound; + const int64_t travel_value_upper_bound = + dimension_travel_info.transition_info[pos].travel_value_upper_bound; + // The lower bound of travel_value is already implemented by constraints as + // travel_value >= compressed_travel_value (defined below) and + // compressed_travel_value has compressed_travel_value_lower_bound as a + // lower bound. The bound is added for the sake of clarity and being + // explicit. + const int travel_value = solver->AddVariable( + compressed_travel_value_lower_bound, travel_value_upper_bound); + SET_DEBUG_VARIABLE_NAME(solver, travel_value, + absl::StrFormat("travel_value(%ld)", pos)); + const int travel_start = solver->AddVariable( + current_route_min_cumuls_[pos] + pre_travel_transit, + current_route_max_cumuls_[pos + 1] - post_travel_transit - + compressed_travel_value_lower_bound); + SET_DEBUG_VARIABLE_NAME(solver, travel_start, + absl::StrFormat("travel_start(%ld)", pos)); + // travel_start = cumul[pos] + pre_travel[pos] + // This becomes: pre_travel[pos] = travel_start - cumul[pos] + solver->AddLinearConstraint(pre_travel_transit, pre_travel_transit, + {{travel_start, 1}, {lp_cumuls[pos], -1}}); + + // Find segments that are in bounds + // Only the segments in [index_anchor_start, index_anchor_end[ are in + // bounds, the others can therefore be discarded. + int index_anchor_start = 0; + while (index_anchor_start < num_pwl_anchors - 1 && + travel_function.x_anchors[index_anchor_start + 1] <= + current_route_min_cumuls_[pos] + pre_travel_transit) { + ++index_anchor_start; + } + int index_anchor_end = num_pwl_anchors - 1; + while (index_anchor_end > 0 && + travel_function.x_anchors[index_anchor_end - 1] >= + current_route_max_cumuls_[pos] + pre_travel_transit) { + --index_anchor_end; + } + // Check that there is at least one segment. + if (index_anchor_start >= index_anchor_end) return false; + + // Precompute the slopes and y-intercept as they will be used to detect + // convexities and in the constraints. + const std::vector slope_and_y_intercept = + PiecewiseLinearFormulationToSlopeAndYIntercept( + travel_function, index_anchor_start, index_anchor_end); + + // Optimize binary variables by detecting convexities + const std::vector convexities = + SlopeAndYInterceptToConvexityRegions(slope_and_y_intercept); + + int nb_bin_variables = 0; + for (const bool convexity : convexities) { + if (convexity) { + ++nb_bin_variables; + if (nb_bin_variables >= 2) break; + } + } + const bool need_bins = (nb_bin_variables > 1); + // 1.b. Create a constraint so that the start belongs to only one segment + const int travel_start_in_one_segment_ct = + need_bins ? solver->CreateNewConstraint(1, 1) + : -1; // -1 is a placeholder value here + + int belongs_to_this_segment_var; + for (int seg = 0; seg < convexities.size(); ++seg) { + if (need_bins && convexities[seg]) { + belongs_to_this_segment_var = solver->AddVariable(0, 1); + SET_DEBUG_VARIABLE_NAME( + solver, belongs_to_this_segment_var, + absl::StrFormat("travel_start(%ld)belongs_to_seg(%ld)", pos, seg)); + solver->SetCoefficient(travel_start_in_one_segment_ct, + belongs_to_this_segment_var, 1); + + // 1.c. Link the binary variable to the departure time + // If the travel_start value is outside the PWL, the closest segment + // will be used. This is why some bounds are infinite. + const int64_t lower_bound_interval = + seg > 0 ? travel_function.x_anchors[index_anchor_start + seg] + : current_route_min_cumuls_[pos] + pre_travel_transit; + int64_t end_of_seg = seg + 1; + while (end_of_seg < num_pwl_anchors - 1 && !convexities[end_of_seg]) { + ++end_of_seg; + } + const int64_t higher_bound_interval = + end_of_seg < num_pwl_anchors - 1 + ? travel_function.x_anchors[index_anchor_start + end_of_seg] + : current_route_max_cumuls_[pos] + pre_travel_transit; + const int travel_start_in_segment_ct = solver->AddLinearConstraint( + lower_bound_interval, higher_bound_interval, {{travel_start, 1}}); + solver->SetEnforcementLiteral(travel_start_in_segment_ct, + belongs_to_this_segment_var); + } + + // 1.d. Compute the slope and y_intercept, the + // coefficient used in the constraint, for each segment. + const auto [slope, y_intercept] = slope_and_y_intercept[seg]; + // Starting later should always mean arriving later + DCHECK_GE(slope, -1.0) << "Travel value PWL should have a slope >= -1"; + + // 1.e. Define the linearization of travel_value + // travel_value - slope * travel_start[pos] = y_intercept, for each + // segment. In order to have a softer constraint, we only impose: + // travel_value - slope * travel_start[pos] >= y_intercept + // and since the cost is increasing in the travel_value, it will + // Minimize it. In addition, since we are working with integers, we + // add a relaxation of 0.5 which gives: + // travel_value - slope * travel_start[pos] >= y_intercept - 0.5. + const double upper_bound = current_route_max_cumuls_[pos]; + const double factor = FindBestScaling( + {1.0, -slope, y_intercept - 0.5}, /*lower_bounds=*/ + {static_cast(compressed_travel_value_lower_bound), 0, 1}, + /*upper_bounds=*/ + {static_cast(travel_value_upper_bound), upper_bound, 1}, + /*max_absolute_activity=*/(int64_t{1} << 62), + /*wanted_absolute_activity_precision=*/1e-3); + // If no correct scaling is found, factor can be equal to 0. This will be + // translated as an unfeasible model as, it will not constraint the + // travel_value with a factor of 0. + if (factor <= 0) return false; + + const int linearization_ct = solver->AddLinearConstraint( + MathUtil::FastInt64Round(factor * (y_intercept - 0.5)), + std::numeric_limits::max(), + {{travel_value, MathUtil::FastInt64Round(factor)}, + {travel_start, MathUtil::FastInt64Round(-factor * slope)}}); + if (need_bins) { + solver->SetEnforcementLiteral(linearization_ct, + belongs_to_this_segment_var); + } + + // ====== UNCOMMENT TO USE ORDER0 ERROR APPROXIMATION ===== // + // Normally cost_scaled = C₂×(Tᵣ - T)²/Tᵣ + // but here we approximate it as cost_scaled = C₂×(Tᵣ - T)²/Tm + // with Tm the average travel value of this segment. + // Therefore, we compute cost_scaled = cost / Tm + // So the cost_function must be defined as cost = C₂×(Tᵣ - T)² + // const int64_t Tm = (transit_function.y_anchors[seg] + + // transit_function.y_anchors[seg + 1]) / 2; The constraint is + // implemented as: cost_scaled * Tm >= cost const int cost_ct = + // solver->AddLinearConstraint(0, std::numeric_limits::max(), + // {{cost_scaled, Tm}, {cost, -1}}); + // solver->SetEnforcementLiteral(cost_ct, belongs_to_this_segment_var); + } + + // 2. Create a variable for the compressed_travel_value. + // cumul[pos + 1] = cumul[pos] + slack[pos] + pre_travel_transit[pos] + + // compressed_travel_value[pos] + post_travel_transit[pos] This becomes: + // post_travel_transit[pos] + pre_travel_transit[pos] = cumul[pos + 1] - + // cumul[pos] - slack[pos] - compressed_travel_value[pos] The higher bound + // of compressed_travel_value is already implemented by constraints as + // travel_compression_absolute = travel_value - compressed_travel_value > 0 + // (see below) and travel_value has travel_value_upper_bound as an upper + // bound. The bound is added for the sake of clarity and being explicit. + const int compressed_travel_value = solver->AddVariable( + compressed_travel_value_lower_bound, travel_value_upper_bound); + SET_DEBUG_VARIABLE_NAME( + solver, compressed_travel_value, + absl::StrFormat("compressed_travel_value(%ld)", pos)); + solver->AddLinearConstraint(post_travel_transit + pre_travel_transit, + post_travel_transit + pre_travel_transit, + {{compressed_travel_value, -1}, + {lp_cumuls[pos + 1], 1}, + {lp_cumuls[pos], -1}, + {lp_slacks[pos], -1}}); + + // 2. Create the travel value compression variable + // travel_compression_absolute == travel_value - compressed_travel_value + // This becomes: 0 = travel_compression_absolute - travel_value + + // compressed_travel_value travel_compression_absolute must be positive or + // equal to 0. + const int travel_compression_absolute = solver->AddVariable( + 0, travel_value_upper_bound - compressed_travel_value_lower_bound); + SET_DEBUG_VARIABLE_NAME( + solver, travel_compression_absolute, + absl::StrFormat("travel_compression_absolute(%ld)", pos)); + + solver->AddLinearConstraint(0, 0, + {{travel_compression_absolute, 1}, + {travel_value, -1}, + {compressed_travel_value, 1}}); + + // 3. Add a cost per unit of travel + // The travel_cost_coefficient is set with the travel_value and not the + // compressed_travel_value to not give the incentive to compress a little + // bit in order to same some cost per travel. + solver->SetObjectiveCoefficient( + travel_value, dimension_travel_info.travel_cost_coefficient); + + // 4. Adds a convex cost in epsilon + // Here we DCHECK that the cost function is indeed convex + const RouteDimensionTravelInfo::TransitionInfo::PiecewiseLinearFormulation& + cost_function = + dimension_travel_info.transition_info[pos].travel_compression_cost; + const std::vector cost_slope_and_y_intercept = + PiecewiseLinearFormulationToSlopeAndYIntercept(cost_function); + const double cost_max = ComputeConvexPiecewiseLinearFormulationValue( + cost_function, + travel_value_upper_bound - compressed_travel_value_lower_bound); + double previous_slope = 0; + for (int seg = 0; seg < cost_function.x_anchors.size() - 1; ++seg) { + const auto [slope, y_intercept] = cost_slope_and_y_intercept[seg]; + // Check convexity + DCHECK_GE(slope, previous_slope) + << "Compression error is not convex. Segment " << (1 + seg) + << " out of " << (cost_function.x_anchors.size() - 1); + previous_slope = slope; + const double factor = FindBestScaling( + {1.0, -slope, y_intercept}, /*lower_bounds=*/ + {0, static_cast(compressed_travel_value_lower_bound), 1}, + /*upper_bounds=*/ + {cost_max, static_cast(travel_value_upper_bound), 1}, + /*max_absolute_activity=*/(static_cast(1) << 62), + /*wanted_absolute_activity_precision=*/1e-3); + // If no correct scaling is found, factor can be equal to 0. This will be + // translated as an unfeasible model as, it will not constraint the + // compression_cost with a factor of 0. + if (factor <= 0) return false; + + solver->AddLinearConstraint( + MathUtil::FastInt64Round(factor * y_intercept), + std::numeric_limits::max(), + {{compression_cost, std::round(factor)}, + {travel_compression_absolute, + MathUtil::FastInt64Round(-factor * slope)}}); + } + // ====== UNCOMMENT TO USE PRODUCT TO COMPUTE THE EXACT ERROR ===== // + // Normally cost_scaled = C₂×(Tᵣ - T)²/Tᵣ + // and here we compute cost_scaled = cost / Tᵣ + // So the cost_function must be defined as cost = C₂×(Tᵣ - T)² + // const int prod = solver->CreateNewPositiveVariable(); + // solver->AddProductConstraint(prod, {cost_scaled, travel_value}); + // The constraint is implemented as: cost_scaled * Tᵣ >= cost + // solver->AddLinearConstraint(0, std::numeric_limits::max(), + // {{prod, 1}, {cost, -1}}); + + // ====== UNCOMMENT TO USE AVERAGE ERROR APPROXIMATION ===== // + // Normally cost_scaled = C₂×(Tᵣ - T)²/Tᵣ + // but here we approximate it as cost_scaled = C₂×(Tᵣ - T)²/Tₐ + // with Tₐ the average travel value (on all the segments). + // Since we do not have access to Tₐ here, we define cost_scaled as + // cost_scaled = cost. So the cost_function must be defined as cost = + // C₂×(Tᵣ - T)²/Tₐ The constraint is implemented as: cost_scaled >= cost + solver->AddLinearConstraint( + 0, std::numeric_limits::max(), + {{relative_compression_cost, 1}, {compression_cost, -1}}); + + solver->SetObjectiveCoefficient(relative_compression_cost, 1.0); + } + return true; +} + bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( int vehicle, const std::function& next_accessor, const std::function& transit_accessor, - int64_t cumul_offset, bool optimize_costs, - RoutingLinearSolverWrapper* solver, int64_t* route_transit_cost, - int64_t* route_cost_offset) { + const RouteDimensionTravelInfo& dimension_travel_info, int64_t cumul_offset, + bool optimize_costs, RoutingLinearSolverWrapper* solver, + int64_t* route_transit_cost, int64_t* route_cost_offset) { RoutingModel* const model = dimension_->model(); // Extract the vehicle's path from next_accessor. std::vector path; @@ -1078,10 +1678,27 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( fixed_transit[pos - 1] = transit_accessor(path[pos - 1], path[pos]); } } - - if (!ComputeRouteCumulBounds(path, fixed_transit, cumul_offset)) { + if (!ExtractRouteCumulBounds(path, cumul_offset)) { return false; } + if (dimension_travel_info.transition_info.empty()) { + if (!TightenRouteCumulBounds(path, fixed_transit, cumul_offset)) { + return false; + } + } else { + // Tighten the bounds with the lower bound of the transit value + std::vector min_transit(path_size - 1); + for (int pos = 0; pos < path_size - 1; ++pos) { + const RouteDimensionTravelInfo::TransitionInfo& transition = + dimension_travel_info.transition_info[pos]; + min_transit[pos] = transition.pre_travel_transit_value + + transition.compressed_travel_value_lower_bound + + transition.post_travel_transit_value; + } + if (!TightenRouteCumulBounds(path, min_transit, cumul_offset)) { + return false; + } + } // LP Model variables, current_route_cumul_variables_ and lp_slacks. // Create LP variables for cumuls. @@ -1089,7 +1706,8 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( lp_cumuls.assign(path_size, -1); for (int pos = 0; pos < path_size; ++pos) { const int lp_cumul = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, lp_cumul, absl::StrFormat("lp_cumul(%ld)", pos)); + SET_DEBUG_VARIABLE_NAME(solver, lp_cumul, + absl::StrFormat("lp_cumul(%ld)", pos)); index_to_cumul_variable_[path[pos]] = lp_cumul; lp_cumuls[pos] = lp_cumul; if (!solver->SetVariableBounds(lp_cumul, current_route_min_cumuls_[pos], @@ -1116,25 +1734,19 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( for (int pos = 0; pos < path_size - 1; ++pos) { const IntVar* cp_slack = dimension_->SlackVar(path[pos]); lp_slacks[pos] = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, lp_slacks[pos], - absl::StrFormat("lp_slacks(%ld)", pos)); + SET_DEBUG_VARIABLE_NAME(solver, lp_slacks[pos], + absl::StrFormat("lp_slacks(%ld)", pos)); if (!solver->SetVariableBounds(lp_slacks[pos], cp_slack->Min(), cp_slack->Max())) { return false; } } - // LP Model constraints and costs. - // Add all path constraints to LP: - // cumul[i] + fixed_transit[i] + slack[i] == cumul[i+1] - // <=> fixed_transit[i] == cumul[i+1] - cumul[i] - slack[i]. - for (int pos = 0; pos < path_size - 1; ++pos) { - const int ct = - solver->CreateNewConstraint(fixed_transit[pos], fixed_transit[pos]); - solver->SetCoefficient(ct, lp_cumuls[pos + 1], 1); - solver->SetCoefficient(ct, lp_cumuls[pos], -1); - solver->SetCoefficient(ct, lp_slacks[pos], -1); + if (!SetRouteTravelConstraints(dimension_travel_info, lp_slacks, + fixed_transit, solver)) { + return false; } + if (route_cost_offset != nullptr) *route_cost_offset = 0; if (optimize_costs) { // Add soft upper bounds. @@ -1155,8 +1767,8 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( continue; } const int soft_ub_diff = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, soft_ub_diff, - absl::StrFormat("soft_ub_diff(%ld)", pos)); + SET_DEBUG_VARIABLE_NAME(solver, soft_ub_diff, + absl::StrFormat("soft_ub_diff(%ld)", pos)); solver->SetObjectiveCoefficient(soft_ub_diff, coef); // cumul - soft_ub_diff <= bound. const int ct = solver->CreateNewConstraint( @@ -1178,8 +1790,8 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( continue; } const int soft_lb_diff = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, soft_lb_diff, - absl::StrFormat("soft_lb_diff(%ld)", pos)); + SET_DEBUG_VARIABLE_NAME(solver, soft_lb_diff, + absl::StrFormat("soft_lb_diff(%ld)", pos)); solver->SetObjectiveCoefficient(soft_lb_diff, coef); // bound - cumul <= soft_lb_diff const int ct = solver->CreateNewConstraint( @@ -1241,7 +1853,7 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( if (bound_cost.bound < std::numeric_limits::max() && bound_cost.cost > 0) { const int span_violation = solver->CreateNewPositiveVariable(); - SET_VARIABLE_NAME(solver, span_violation, "span_violation"); + SET_DEBUG_VARIABLE_NAME(solver, span_violation, "span_violation"); // end - start <= bound + span_violation const int violation = solver->CreateNewConstraint( std::numeric_limits::min(), bound_cost.bound); @@ -1365,15 +1977,15 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( } lp_break_start[br] = solver->AddVariable(break_start_min, break_start_max); - SET_VARIABLE_NAME(solver, lp_break_start[br], - absl::StrFormat("lp_break_start(%ld)", br)); + SET_DEBUG_VARIABLE_NAME(solver, lp_break_start[br], + absl::StrFormat("lp_break_start(%ld)", br)); lp_break_end[br] = solver->AddVariable(break_end_min, break_end_max); - SET_VARIABLE_NAME(solver, lp_break_end[br], - absl::StrFormat("lp_break_end(%ld)", br)); + SET_DEBUG_VARIABLE_NAME(solver, lp_break_end[br], + absl::StrFormat("lp_break_end(%ld)", br)); lp_break_duration[br] = solver->AddVariable(break_duration_min, break_duration_max); - SET_VARIABLE_NAME(solver, lp_break_duration[br], - absl::StrFormat("lp_break_duration(%ld)", br)); + SET_DEBUG_VARIABLE_NAME(solver, lp_break_duration[br], + absl::StrFormat("lp_break_duration(%ld)", br)); // start + duration = end. solver->AddLinearConstraint(0, 0, {{lp_break_end[br], 1}, @@ -1409,8 +2021,9 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( 0, std::numeric_limits::max(), {{lp_cumuls.front(), 1}, {lp_break_end[br], -1}}); const int break_is_before_route = solver->AddVariable(0, 1); - SET_VARIABLE_NAME(solver, break_is_before_route, - absl::StrFormat("break_is_before_route(%ld)", br)); + SET_DEBUG_VARIABLE_NAME( + solver, break_is_before_route, + absl::StrFormat("break_is_before_route(%ld)", br)); solver->SetEnforcementLiteral(ct, break_is_before_route); solver->SetCoefficient(break_in_one_slack_ct, break_is_before_route, 1); } @@ -1420,8 +2033,9 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( 0, std::numeric_limits::max(), {{lp_break_start[br], 1}, {lp_cumuls.back(), -1}}); const int break_is_after_route = solver->AddVariable(0, 1); - SET_VARIABLE_NAME(solver, break_is_after_route, - absl::StrFormat("break_is_after_route(%ld)", br)); + SET_DEBUG_VARIABLE_NAME( + solver, break_is_after_route, + absl::StrFormat("break_is_after_route(%ld)", br)); solver->SetEnforcementLiteral(ct, break_is_after_route); solver->SetCoefficient(break_in_one_slack_ct, break_is_after_route, 1); } @@ -1450,15 +2064,26 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( // sum_br break_duration_min(br) * break_in_slack(br, pos) <= // lp_slacks(pos). const int break_in_slack = solver->AddVariable(0, 1); - SET_VARIABLE_NAME(solver, break_in_slack, - absl::StrFormat("break_in_slack(%ld, %ld)", br, pos)); - solver->SetCoefficient(break_in_one_slack_ct, break_in_slack, 1); + SET_DEBUG_VARIABLE_NAME( + solver, break_in_slack, + absl::StrFormat("break_in_slack(%ld, %ld)", br, pos)); if (slack_linear_lower_bound_ct[pos] == -1) { slack_linear_lower_bound_ct[pos] = solver->AddLinearConstraint( std::numeric_limits::min(), 0, {{lp_slacks[pos], -1}}); } - solver->SetCoefficient(slack_linear_lower_bound_ct[pos], break_in_slack, - break_duration_min); + // To keep the model clean + // (cf. glop::LinearProgram::NotifyThatColumnsAreClean), constraints on + // break_in_slack need to be in ascending order. + if (break_in_one_slack_ct < slack_linear_lower_bound_ct[pos]) { + solver->SetCoefficient(break_in_one_slack_ct, break_in_slack, 1); + solver->SetCoefficient(slack_linear_lower_bound_ct[pos], break_in_slack, + break_duration_min); + } else { + solver->SetCoefficient(slack_linear_lower_bound_ct[pos], break_in_slack, + break_duration_min); + solver->SetCoefficient(break_in_one_slack_ct, break_in_slack, 1); + } + if (solver->IsCPSATSolver()) { // Exact relation between breaks, slacks and cumul variables. // Make an exact slack lower bound (lazily), that represents @@ -1466,7 +2091,7 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( // lp_slacks(pos). const int break_duration_in_slack = solver->AddVariable(0, slack_duration_max); - SET_VARIABLE_NAME( + SET_DEBUG_VARIABLE_NAME( solver, break_duration_in_slack, absl::StrFormat("break_duration_in_slack(%ld, %ld)", br, pos)); solver->AddProductConstraint(break_duration_in_slack, @@ -1543,7 +2168,7 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( // the whole timeline (-infinity, +infinity). int previous_cover = solver->AddVariable(CapAdd(vehicle_start_min, limit), CapAdd(vehicle_start_max, limit)); - SET_VARIABLE_NAME(solver, previous_cover, "previous_cover"); + SET_DEBUG_VARIABLE_NAME(solver, previous_cover, "previous_cover"); solver->AddLinearConstraint(limit, limit, {{previous_cover, 1}, {lp_cumuls.front(), -1}}); for (int br = 0; br < num_breaks; ++br) { @@ -1553,11 +2178,12 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( // break_is_eligible <=> // break_end - break_start >= break_minimum_duration. const int break_is_eligible = solver->AddVariable(0, 1); - SET_VARIABLE_NAME(solver, break_is_eligible, - absl::StrFormat("break_is_eligible(%ld)", br)); + SET_DEBUG_VARIABLE_NAME(solver, break_is_eligible, + absl::StrFormat("break_is_eligible(%ld)", br)); const int break_is_not_eligible = solver->AddVariable(0, 1); - SET_VARIABLE_NAME(solver, break_is_not_eligible, - absl::StrFormat("break_is_not_eligible(%ld)", br)); + SET_DEBUG_VARIABLE_NAME( + solver, break_is_not_eligible, + absl::StrFormat("break_is_not_eligible(%ld)", br)); { solver->AddLinearConstraint( 1, 1, {{break_is_eligible, 1}, {break_is_not_eligible, 1}}); @@ -1578,8 +2204,8 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( const int break_cover = solver->AddVariable( CapAdd(std::min(vehicle_start_min, break_end_min), limit), CapAdd(std::max(vehicle_start_min, break_end_max), limit)); - SET_VARIABLE_NAME(solver, break_cover, - absl::StrFormat("break_cover(%ld)", br)); + SET_DEBUG_VARIABLE_NAME(solver, break_cover, + absl::StrFormat("break_cover(%ld)", br)); const int limit_cover_ct = solver->AddLinearConstraint( limit, limit, {{break_cover, 1}, {lp_break_end[br], -1}}); solver->SetEnforcementLiteral(limit_cover_ct, break_is_eligible); @@ -1591,7 +2217,7 @@ bool DimensionCumulOptimizerCore::SetRouteCumulConstraints( const int cover = solver->AddVariable(CapAdd(vehicle_start_min, limit), std::numeric_limits::max()); - SET_VARIABLE_NAME(solver, cover, absl::StrFormat("cover(%ld)", br)); + SET_DEBUG_VARIABLE_NAME(solver, cover, absl::StrFormat("cover(%ld)", br)); solver->AddMaximumConstraint(cover, {previous_cover, break_cover}); // Cover chaining. If route end is not covered, break start must be: // cover_{i-1} < route_end => s_i <= cover_{i-1} @@ -1725,12 +2351,21 @@ bool DimensionCumulOptimizerCore::SetGlobalConstraints( if (vehicle_constraints[v] == kNoConstraint) continue; const int assign_r_to_v = solver->AddVariable(0, 1); - SET_VARIABLE_NAME(solver, assign_r_to_v, - absl::StrFormat("assign_r_to_v(%ld, %ld)", r, v)); + SET_DEBUG_VARIABLE_NAME( + solver, assign_r_to_v, + absl::StrFormat("assign_r_to_v(%ld, %ld)", r, v)); resource_to_vehicle_assignment_variables[r * num_vehicles + v] = assign_r_to_v; - solver->SetCoefficient(vehicle_constraints[v], assign_r_to_v, 1); - solver->SetCoefficient(resource_constraints[r], assign_r_to_v, 1); + // To keep the model clean + // (cf. glop::LinearProgram::NotifyThatColumnsAreClean), constraints on + // assign_r_to_v need to be in ascending order. + if (vehicle_constraints[v] < resource_constraints[r]) { + solver->SetCoefficient(vehicle_constraints[v], assign_r_to_v, 1); + solver->SetCoefficient(resource_constraints[r], assign_r_to_v, 1); + } else { + solver->SetCoefficient(resource_constraints[r], assign_r_to_v, 1); + solver->SetCoefficient(vehicle_constraints[v], assign_r_to_v, 1); + } const auto& add_domain_constraint = [&solver, cumul_offset, assign_r_to_v](const Domain& domain, @@ -1758,7 +2393,7 @@ bool DimensionCumulOptimizerCore::SetGlobalConstraints( return true; } -#undef SET_VARIABLE_NAME +#undef SET_DEBUG_VARIABLE_NAME void DimensionCumulOptimizerCore::SetValuesFromLP( const std::vector& lp_variables, int64_t offset, @@ -1849,8 +2484,8 @@ GlobalDimensionCumulOptimizer::ComputeCumulCostWithoutFixedTransits( int64_t cost = 0; int64_t transit_cost = 0; DimensionSchedulingStatus status = - optimizer_core_.Optimize(next_accessor, solver_.get(), nullptr, nullptr, - nullptr, &cost, &transit_cost); + optimizer_core_.Optimize(next_accessor, {}, solver_.get(), nullptr, + nullptr, nullptr, &cost, &transit_cost); if (status != DimensionSchedulingStatus::INFEASIBLE && optimal_cost_without_transits != nullptr) { *optimal_cost_without_transits = CapSub(cost, transit_cost); @@ -1860,64 +2495,65 @@ GlobalDimensionCumulOptimizer::ComputeCumulCostWithoutFixedTransits( DimensionSchedulingStatus GlobalDimensionCumulOptimizer::ComputeCumuls( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, std::vector* optimal_cumuls, std::vector* optimal_breaks, std::vector>* optimal_resource_indices) { - return optimizer_core_.Optimize(next_accessor, solver_.get(), optimal_cumuls, - optimal_breaks, optimal_resource_indices, - nullptr, nullptr); + return optimizer_core_.Optimize(next_accessor, + dimension_travel_info_per_route, + solver_.get(), optimal_cumuls, optimal_breaks, + optimal_resource_indices, nullptr, nullptr); } DimensionSchedulingStatus GlobalDimensionCumulOptimizer::ComputePackedCumuls( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, std::vector* packed_cumuls, std::vector* packed_breaks, std::vector>* resource_indices) { - return optimizer_core_.OptimizeAndPack(next_accessor, solver_.get(), - packed_cumuls, packed_breaks, - resource_indices); + return optimizer_core_.OptimizeAndPack( + next_accessor, dimension_travel_info_per_route, solver_.get(), + packed_cumuls, packed_breaks, resource_indices); } -// ResourceAssignmentOptimizer - -ResourceAssignmentOptimizer::ResourceAssignmentOptimizer( - const RoutingModel::ResourceGroup* resource_group, - LocalDimensionCumulOptimizer* optimizer, - LocalDimensionCumulOptimizer* mp_optimizer) - : optimizer_(*optimizer), - mp_optimizer_(*mp_optimizer), - model_(*optimizer->dimension()->model()), - resource_group_(*resource_group) {} - -bool ResourceAssignmentOptimizer::ComputeAssignmentCostsForVehicle( - int v, const std::function& next_accessor, +bool ComputeVehicleToResourcesAssignmentCosts( + int v, const RoutingModel::ResourceGroup& resource_group, + const std::function& next_accessor, const std::function& transit_accessor, - bool optimize_vehicle_costs, std::vector* assignment_costs, + bool optimize_vehicle_costs, LocalDimensionCumulOptimizer* lp_optimizer, + LocalDimensionCumulOptimizer* mp_optimizer, + std::vector* assignment_costs, std::vector>* cumul_values, std::vector>* break_values) { + DCHECK(lp_optimizer != nullptr); + DCHECK(mp_optimizer != nullptr); + const RoutingDimension* dimension = lp_optimizer->dimension(); + DCHECK_EQ(dimension, mp_optimizer->dimension()); + RoutingModel* const model = dimension->model(); DCHECK_NE(assignment_costs, nullptr); - if (!resource_group_.VehicleRequiresAResource(v) || - (next_accessor(model_.Start(v)) == model_.End(v) && - !model_.IsVehicleUsedWhenEmpty(v))) { + if (!resource_group.VehicleRequiresAResource(v) || + (!model->IsVehicleUsedWhenEmpty(v) && + next_accessor(model->Start(v)) == model->End(v))) { assignment_costs->clear(); return true; } - const RoutingDimension& dimension = *optimizer_.dimension(); - if (dimension.model()->CheckLimit()) { + if (model->CheckLimit()) { // The model's time limit has been reached, stop everything. return false; } const std::vector& resources = - resource_group_.GetResources(); + resource_group.GetResources(); const int num_resources = resources.size(); std::vector all_resource_indices(num_resources); std::iota(all_resource_indices.begin(), all_resource_indices.end(), 0); const bool use_mp_optimizer = - dimension.HasBreakConstraints() && - !dimension.GetBreakIntervalsOfVehicle(v).empty(); - LocalDimensionCumulOptimizer& optimizer = - use_mp_optimizer ? mp_optimizer_ : optimizer_; + dimension->HasBreakConstraints() && + !dimension->GetBreakIntervalsOfVehicle(v).empty(); + LocalDimensionCumulOptimizer* optimizer = + use_mp_optimizer ? mp_optimizer : lp_optimizer; std::vector statuses = - optimizer.ComputeRouteCumulCostsForResourcesWithoutFixedTransits( + optimizer->ComputeRouteCumulCostsForResourcesWithoutFixedTransits( v, next_accessor, transit_accessor, resources, all_resource_indices, optimize_vehicle_costs, assignment_costs, cumul_values, break_values); @@ -1948,7 +2584,7 @@ bool ResourceAssignmentOptimizer::ComputeAssignmentCostsForVehicle( std::vector mp_assignment_costs; std::vector> mp_cumul_values; std::vector> mp_break_values; - mp_optimizer_.ComputeRouteCumulCostsForResourcesWithoutFixedTransits( + mp_optimizer->ComputeRouteCumulCostsForResourcesWithoutFixedTransits( v, next_accessor, transit_accessor, resources, mp_optimizer_resource_indices, optimize_vehicle_costs, &mp_assignment_costs, @@ -1979,19 +2615,82 @@ bool ResourceAssignmentOptimizer::ComputeAssignmentCostsForVehicle( [](int64_t cost) { return cost >= 0; }); } -int64_t ResourceAssignmentOptimizer::ComputeBestAssignmentCost( - const std::vector>& - primary_vehicle_to_resource_assignment_costs, - const std::vector>& - secondary_vehicle_to_resource_assignment_costs, - const std::function& use_primary_for_vehicle, - std::vector* resource_indices) const { - const int num_vehicles = model_.vehicles(); - DCHECK_EQ(primary_vehicle_to_resource_assignment_costs.size(), num_vehicles); - DCHECK_EQ(secondary_vehicle_to_resource_assignment_costs.size(), - num_vehicles); - const int num_resources = resource_group_.Size(); +int64_t ComputeBestVehicleToResourceAssignment( + std::vector vehicles, int num_resources, + std::function*(int)> + vehicle_to_resource_assignment_costs, + std::vector* resource_indices) { + DCHECK_GE(num_resources, 1); // Else this whole function doesn't make sense. + const int num_vehicles = vehicles.size(); + int num_total_vehicles = -1; + if (resource_indices != nullptr) { + num_total_vehicles = resource_indices->size(); + // When returning infeasible, 'resource_indices' must be cleared, so we do + // it here preemptively. + resource_indices->clear(); + DCHECK_GE(num_total_vehicles, num_vehicles); + for (int v : vehicles) { + DCHECK_GE(v, 0); + DCHECK_LT(v, num_total_vehicles); + } + } + // Collect vehicle_to_resource_assignment_costs(v) for all v ∈ vehicles. + // Then detect trivial infeasibility cases, before doing the min-cost-flow: + // - There are not enough resources overall. + // - There is no resource assignable to a vehicle that needs one. + std::vector*> vi_to_resource_cost(num_vehicles); + int num_vehicles_to_assign = 0; + for (int i = 0; i < num_vehicles; ++i) { + vi_to_resource_cost[i] = vehicle_to_resource_assignment_costs(vehicles[i]); + if (!vi_to_resource_cost[i]->empty()) { + DCHECK_EQ(vi_to_resource_cost[i]->size(), num_resources); + ++num_vehicles_to_assign; + } + } + if (num_vehicles_to_assign > num_resources) { + VLOG(3) << "Less resources (" << num_resources << ") than the vehicles" + << " requiring one (" << num_vehicles_to_assign << ")"; + return -1; // Infeasible. + } + // Catch infeasibility cases where ComputeVehicleToResourcesAssignmentCosts() + // hasn't "properly" initialized the vehicle to resource assignment costs + // (this can happen for instance in the ResourceGroupAssignmentFilter when + // routes are synchronized with an impossible first solution). + for (int i = 0; i < num_vehicles; ++i) { + if (!vi_to_resource_cost[i]->empty() && + *absl::c_max_element(*vi_to_resource_cost[i]) < 0) { + VLOG(3) << "Vehicle #" << vehicles[i] << " has no feasible resource"; + return -1; + } + } + + // We may need to apply some cost scaling when using SimpleMinCostFlow: + // 3 * max_arc_cost * num_nodes must be ≤ kint64max. + // To do that, we first find the maximum arc cost. + int64_t max_arc_cost = 0; + for (const std::vector* costs : vi_to_resource_cost) { + if (costs->empty()) continue; + max_arc_cost = std::max(max_arc_cost, *absl::c_max_element(*costs)); + } + // To avoid potential int64_t overflows, we tweak the above formula to: + // max_acceptable_arc_cost = kint64max / (3 * num_nodes) - 1. + // NOTE(user): SimpleMinCostFlow always adds a sink and source node (we + // probably shouldn't add a sink/source node ourselves in the graph). + const int real_num_nodes = 4 + num_vehicles + num_resources; + const int64_t max_acceptable_arc_cost = kint64max / (3 * real_num_nodes) - 1; + // We use a power of 2 for the cost scaling factor, to have clean (in)accuracy + // properties. Note also that we must round *down* the costs. + int cost_right_shift = 0; + while ((max_arc_cost >> cost_right_shift) > max_acceptable_arc_cost) { + ++cost_right_shift; + } + + // Then, we create the SimpleMinCostFlow and run the assignment algorithm. + // NOTE(user): We often don't create as many arcs as outlined below, + // especially when num_vehicles_to_assign < vehicles.size(). But since we + // want to eventually make this whole function incremental, we prefer sticking + // with the whole 'vehicles' set. SimpleMinCostFlow flow( /*reserve_num_nodes*/ 2 + num_vehicles + num_resources, /*reserve_num_arcs*/ num_vehicles + num_vehicles * num_resources + @@ -2002,51 +2701,28 @@ int64_t ResourceAssignmentOptimizer::ComputeBestAssignmentCost( return num_vehicles + r; }; - std::vector> vehicle_to_resource_arc_index; + // Used to store the arc indices, if we need to later recover the solution. + FlatMatrix vehicle_to_resource_arc_index; if (resource_indices != nullptr) { - vehicle_to_resource_arc_index.resize( - num_vehicles, std::vector(num_resources, -1)); + vehicle_to_resource_arc_index = + FlatMatrix(num_vehicles, num_resources, -1); } - int num_used_vehicles = 0; - for (int v : resource_group_.GetVehiclesRequiringAResource()) { - DCHECK(use_primary_for_vehicle(v) || - primary_vehicle_to_resource_assignment_costs[v].empty()); - const std::vector& assignment_costs = - use_primary_for_vehicle(v) - ? primary_vehicle_to_resource_assignment_costs[v] - : secondary_vehicle_to_resource_assignment_costs[v]; - if (assignment_costs.empty()) { - // We don't need a resource for this vehicle. - continue; - } - DCHECK_EQ(assignment_costs.size(), num_resources); - num_used_vehicles++; - if (num_used_vehicles > num_resources) { - // Not enough resources for all vehicles needing one. - return -1; - } + for (int vi = 0; vi < num_vehicles; ++vi) { + const std::vector& assignment_costs = *vi_to_resource_cost[vi]; + if (assignment_costs.empty()) continue; // Doesn't need resources. - // Add a source->vehicle arc to the flow. - flow.AddArcWithCapacityAndUnitCost(source_index, v, 1, 0); + // Add a source → vehicle arc to the min-cost-flow graph. + flow.AddArcWithCapacityAndUnitCost(source_index, vi, 1, 0); - // Add arcs to the min-cost-flow graph. - bool has_feasible_resource = false; + // Add vehicle → resource arcs to the min-cost-flow graph. for (int r = 0; r < num_resources; r++) { const int64_t assignment_cost = assignment_costs[r]; if (assignment_cost < 0) continue; - const ArcIndex arc_index = flow.AddArcWithCapacityAndUnitCost( - v, resource_index(r), 1, assignment_cost); - if (!vehicle_to_resource_arc_index.empty()) { - vehicle_to_resource_arc_index[v][r] = arc_index; + const ArcIndex arc = flow.AddArcWithCapacityAndUnitCost( + vi, resource_index(r), 1, assignment_cost >> cost_right_shift); + if (resource_indices != nullptr) { + vehicle_to_resource_arc_index[vi][r] = arc; } - has_feasible_resource = true; - } - if (!has_feasible_resource) { - // Catches cases of infeasibility where ComputeAssignmentCostsForVehicle() - // hasn't "properly" initialized the primary/secondary assignment costs - // (this can happen for instance in the ResourceGroupAssignmentFilter - // when routes are synchronized with an impossible first solution). - return -1; } } @@ -2056,32 +2732,374 @@ int64_t ResourceAssignmentOptimizer::ComputeBestAssignmentCost( } // Set the flow supply. - flow.SetNodeSupply(source_index, num_used_vehicles); - flow.SetNodeSupply(sink_index, -num_used_vehicles); + flow.SetNodeSupply(source_index, num_vehicles_to_assign); + flow.SetNodeSupply(sink_index, -num_vehicles_to_assign); // Solve the min-cost flow and return its cost. if (flow.Solve() != SimpleMinCostFlow::OPTIMAL) { - if (resource_indices != nullptr) resource_indices->clear(); + VLOG(3) << "Non-OPTIMAL flow result"; return -1; } - const int64_t cost = flow.OptimalCost(); - if (resource_indices == nullptr) { - return cost; - } - - // Fill the resource indices corresponding to the min-cost assignment. - resource_indices->assign(num_vehicles, -1); - for (int v : resource_group_.GetVehiclesRequiringAResource()) { - for (int r = 0; r < num_resources; r++) { - if (vehicle_to_resource_arc_index[v][r] >= 0 && - flow.Flow(vehicle_to_resource_arc_index[v][r]) > 0) { - resource_indices->at(v) = r; - break; + if (resource_indices != nullptr) { + // Fill the resource indices corresponding to the min-cost assignment. + resource_indices->assign(num_total_vehicles, -1); + for (int vi = 0; vi < num_vehicles; ++vi) { + for (int r = 0; r < num_resources; r++) { + const ArcIndex arc = vehicle_to_resource_arc_index[vi][r]; + if (arc >= 0 && flow.Flow(arc) > 0) { + resource_indices->at(vehicles[vi]) = r; + break; + } } } } - return cost; + + const int64_t cost = flow.OptimalCost(); + DCHECK_LE(cost, kint64max >> cost_right_shift); + return cost << cost_right_shift; +} + +std::string Int64ToStr(int64_t number) { + if (number == kint64min) return "-infty"; + if (number == kint64max) return "+infty"; + return std::to_string(number); +} + +std::string DomainToString( + const ::google::protobuf::RepeatedField* domain) { + if (domain->size() > 2 && domain->size() % 2 == 0) { + std::string s = "∈ "; + for (int i = 0; i < domain->size(); i += 2) { + s += absl::StrFormat("[%s, %s]", Int64ToStr(domain->Get(i)), + Int64ToStr(domain->Get(i + 1))); + if (i < domain->size() - 2) s += " ∪ "; + } + return s; + } else if (domain->size() == 2) { + if (domain->Get(0) == domain->Get(1)) { + return absl::StrFormat("= %s", Int64ToStr(domain->Get(0))); + } else if (domain->Get(0) == 0 && domain->Get(1) == 1) { + return "∈ Binary"; + } else if (domain->Get(0) == std::numeric_limits::min() && + domain->Get(1) == std::numeric_limits::max()) { + return "∈ ℝ"; + } else if (domain->Get(0) == std::numeric_limits::min()) { + return absl::StrFormat("≤ %s", Int64ToStr(domain->Get(1))); + } else if (domain->Get(1) == std::numeric_limits::max()) { + return absl::StrFormat("≥ %s", Int64ToStr(domain->Get(0))); + } + return absl::StrFormat("∈ [%ls, %s]", Int64ToStr(domain->Get(0)), + Int64ToStr(domain->Get(1))); + } else if (domain->size() == 1) { + return absl::StrFormat("= %s", Int64ToStr(domain->Get(0))); + } else { + return absl::StrFormat("∈ Unknown domain (size=%ld)", domain->size()); + } +} + +std::string VariableToString( + std::pair& variable_pair, + const sat::CpSolverResponse& response_) { + std::string s = ""; + sat::IntegerVariableProto& variable = variable_pair.first; + const int index = variable_pair.second; + if (response_.IsInitialized() && variable.IsInitialized() && + (response_.status() == sat::CpSolverStatus::OPTIMAL || + response_.status() == sat::CpSolverStatus::FEASIBLE)) { + const double lp_value_double = response_.solution(index); + const int64_t lp_value_int64 = + (lp_value_double >= std::numeric_limits::max()) + ? std::numeric_limits::max() + : MathUtil::FastInt64Round(lp_value_double); + s += Int64ToStr(lp_value_int64) + " "; + } else { + s += "? "; + } + s += DomainToString(variable.mutable_domain()); + return s; +} + +std::string ConstraintToString(const sat::ConstraintProto& constraint, + const sat::CpModelProto& model_, + bool show_enforcement = true) { + std::string s = ""; + if (constraint.has_linear()) { + const auto& linear = constraint.linear(); + for (int j = 0; j < linear.vars().size(); ++j) { + const std::string sign = linear.coeffs(j) > 0 ? "+" : "-"; + const std::string mult = + std::abs(linear.coeffs(j)) != 1 + ? std::to_string(std::abs(linear.coeffs(j))) + " * " + : ""; + if (j > 0 || sign != "+") s += sign + " "; + s += mult + model_.variables(linear.vars(j)).name() + " "; + } + s += DomainToString(&linear.domain()); + + // Enforcement literal. + if (show_enforcement) { + for (int j = 0; j < constraint.enforcement_literal_size(); ++j) { + s += (j == 0) ? "\t if " : " and "; + s += model_.variables(constraint.enforcement_literal(j)).name(); + } + } + } else { + s += constraint.ShortDebugString(); + } + return s; +} + +std::string VariablesToString( + absl::flat_hash_map>& + variables, + absl::flat_hash_map>& variable_instances, + absl::flat_hash_map>& + variable_childs, + const sat::CpSolverResponse& response_, const std::string& variable, + std::string prefix = "") { + if (variable.empty()) { + std::string s = ""; + const auto& childs = variable_childs[""]; + for (const std::string& child : childs) { + s += prefix + + VariablesToString(variables, variable_instances, variable_childs, + response_, child, prefix) + + prefix + "\n"; + } + return s; + } + + const auto& instances = variable_instances[variable]; + std::string variable_display = variable; + std::size_t bracket_pos = variable.find_last_of(')'); + if (bracket_pos != std::string::npos) { + variable_display = variable.substr(bracket_pos + 1); + } + std::string s = variable_display + " | "; + prefix += std::string(variable_display.length(), ' ') + " | "; + for (int i = 0; i < instances.size(); ++i) { + const std::string instance_name = + absl::StrFormat("%s(%ld)", variable, instances[i]); + if (i > 0) s += prefix; + s += absl::StrFormat("%ld: %s", instances[i], + VariableToString(variables[instance_name], response_)); + + // Childs + const auto& childs = variable_childs[instance_name]; + for (const std::string& child : childs) { + s += "\n" + prefix + "| "; + s += VariablesToString(variables, variable_instances, variable_childs, + response_, child, prefix + "| "); + } + if (childs.empty()) s += "\n"; + } + return s; +} + +std::string RoutingCPSatWrapper::PrintModel() const { + // Constraints you want to separate. + std::vector> constraints_apart; + constraints_apart.push_back( + {"compression_cost", "travel_compression_absolute"}); + + // variable_instances links the lemma of a variable to the different number of + // instantiation. For instance if you have in your model x(0), x(1) and x(4), + // the key "x" will be associated to {0,1,4}. + absl::flat_hash_map> variable_instances; + // variable_children links a variable to its children. That is, if you have in + // you model x(0), then typical childs would be {"x(0)in_segment(0)", + // "x(0)in_segment(1)", "x(0)scaled", ...} + absl::flat_hash_map> + variable_children; + // variables link the name of a variable to its Proto. + absl::flat_hash_map> + variables; + variable_children[""] = {}; + + const int num_constraints = model_.constraints_size(); + const int num_variables = model_.variables_size(); + int num_binary_variables = 0; + for (int i = 0; i < num_variables; ++i) { + const auto& variable = model_.variables(i); + const auto& name = variable.name(); + const int pos_bracket = name.find_last_of('('); + if (pos_bracket != std::string::npos) { + const std::string lemma = name.substr(0, pos_bracket); + const int pos_closing_bracket = name.find_last_of(')'); + CHECK_NE(pos_closing_bracket, std::string::npos); + const int index = + std::stoi(name.substr(pos_bracket + 1, pos_closing_bracket)); + std::vector* instances = gtl::FindOrNull(variable_instances, lemma); + if (instances != nullptr) { + instances->push_back(index); + } else { + variable_instances[lemma] = {index}; + } + variable_children[name] = {}; + + std::string parent = ""; + const int pos_parent_closing_bracket = lemma.find_last_of(')'); + if (pos_parent_closing_bracket != std::string::npos) { + parent = lemma.substr(0, pos_parent_closing_bracket + 1); + } + variable_children[parent].emplace(lemma); + variables[name] = std::make_pair(variable, i); + if (variable.domain(0) == 0 & variable.domain(1) == 1) { + ++num_binary_variables; + } + } + } + + // Preparing constraints + // the constraints hashmap associate enforcement to constraints. + // If the ket is "", then the constraint has no enforcement and if the key is + // "multiple", then the constraint has several enforcement. If the constraint + // has a single enforcement, then the key will be the variable name of the + // enforcement. + absl::flat_hash_map> + constraints; + absl::flat_hash_map, + std::vector> + constraint_groups; + for (int i = 0; i < num_constraints; ++i) { + const auto& constraint = model_.constraints(i); + std::string enforcement = ""; + if (constraint.enforcement_literal_size() == 1) { + enforcement = model_.variables(constraint.enforcement_literal(0)).name(); + } else if (constraint.enforcement_literal_size() > 1) { + enforcement = "multiple"; + } else { + if (constraint.has_linear()) { + const auto& linear = constraint.linear(); + std::vector key; + for (int j = 0; j < linear.vars().size(); ++j) { + std::string var_name = model_.variables(linear.vars(j)).name(); + std::string lemma = var_name.substr(0, var_name.find_last_of('(')); + key.push_back(lemma); + } + auto* constraint_group = gtl::FindOrNull(constraint_groups, key); + if (constraint_group != nullptr) { + constraint_group->push_back(constraint); + } else { + constraint_groups[key] = {constraint}; + } + } + } + auto* constraints_enforced = gtl::FindOrNull(constraints, enforcement); + if (constraints_enforced != nullptr) { + constraints[enforcement].push_back(constraint); + } else { + constraints[enforcement] = {constraint}; + } + } + + const std::string prefix_constraint = " • "; + std::string s = "Using RoutingCPSatWrapper.\n"; + s += absl::StrFormat("\nObjective = %f\n", this->GetObjectiveValue()); + + for (int i = 0; i < objective_coefficients_.size(); ++i) { + double coeff = objective_coefficients_[i]; + if (coeff != 0) { + s += absl::StrFormat(" | %f * %s\n", coeff, model_.variables(i).name()); + } + } + + s += absl::StrFormat("\nVariables %d (%d Binary - %d Non Binary)\n", + num_variables, num_binary_variables, + num_variables - num_binary_variables); + s += VariablesToString(variables, variable_instances, variable_children, + response_, "", " | "); + s += absl::StrFormat("\n\nConstraints (%d)\n", num_constraints); + + // Constraints NOT enforced + s += "\n- Not enforced\n"; + bool at_least_one_not_enforced = false; + for (const auto& pair : constraint_groups) { + if (!std::count(constraints_apart.begin(), constraints_apart.end(), + pair.first)) { + for (const auto& constraint : pair.second) { + s += prefix_constraint + ConstraintToString(constraint, model_, true) + + "\n"; + at_least_one_not_enforced = true; + } + } + } + if (!at_least_one_not_enforced) { + s += prefix_constraint + "None\n"; + } + + // Constraints with a SINGLE enforcement + s += "\n- Single enforcement\n"; + bool at_least_one_single_enforced = false; + for (const auto& pair : variable_instances) { + const std::string lemma = pair.first; + bool found_one_constraint = false; + std::string prefix = ""; + for (int instance : pair.second) { + const std::string enforcement = + absl::StrFormat("%s(%d)", lemma, instance); + auto* constraints_enforced = gtl::FindOrNull(constraints, enforcement); + std::string prefix_instance = ""; + if (constraints_enforced != nullptr) { + at_least_one_single_enforced = true; + if (!found_one_constraint) { + found_one_constraint = true; + s += prefix_constraint + "if " + lemma + " | "; + prefix = + std::string(prefix_constraint.size() + 1 + lemma.size(), ' ') + + " | "; + } else { + s += prefix; + } + s += absl::StrFormat("%d: | ", instance); + prefix_instance = prefix + " | "; + bool first = true; + for (const auto& constraint : *constraints_enforced) { + if (!first) + s += prefix_instance; + else + first = false; + s += ConstraintToString(constraint, model_, false) + "\n"; + } + } + } + } + if (!at_least_one_single_enforced) { + s += prefix_constraint + "None\n"; + } + + // Constraints with MULTIPLE enforcement + s += "\n- Multiple enforcement\n"; + auto* constraints_multiple_enforced = + gtl::FindOrNull(constraints, "multiple"); + if (constraints_multiple_enforced != nullptr) { + for (const auto& constraint : *constraints_multiple_enforced) { + s += prefix_constraint + ConstraintToString(constraint, model_, true) + + "\n"; + } + } else { + s += prefix_constraint + "None\n"; + } + + // Constraints apart + s += "\n- Set apart\n"; + bool at_least_one_apart = false; + for (const auto& pair : constraint_groups) { + if (std::count(constraints_apart.begin(), constraints_apart.end(), + pair.first)) { + for (const auto& constraint : pair.second) { + s += prefix_constraint + ConstraintToString(constraint, model_, true) + + "\n"; + at_least_one_apart = true; + } + } + } + if (!at_least_one_apart) { + s += prefix_constraint + "None\n"; + } + + return s; } } // namespace operations_research diff --git a/ortools/constraint_solver/routing_lp_scheduling.h b/ortools/constraint_solver/routing_lp_scheduling.h index eaa45342e3..9da95668ee 100644 --- a/ortools/constraint_solver/routing_lp_scheduling.h +++ b/ortools/constraint_solver/routing_lp_scheduling.h @@ -21,12 +21,14 @@ #include #include #include +#include #include #include #include #include "absl/container/flat_hash_map.h" #include "absl/time/time.h" +#include "ortools/base/dump_vars.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/constraint_solver/routing.h" @@ -60,7 +62,9 @@ class CumulBoundsPropagator { // bounds of an index. bool PropagateCumulBounds( const std::function& next_accessor, - int64_t cumul_offset); + int64_t cumul_offset, + const std::vector* + dimension_travel_info_per_route = nullptr); int64_t CumulMin(int index) const { return propagated_bounds_[PositiveNode(index)]; @@ -105,7 +109,9 @@ class CumulBoundsPropagator { bool InitializeArcsAndBounds( const std::function& next_accessor, - int64_t cumul_offset); + int64_t cumul_offset, + const std::vector* + dimension_travel_info_per_route = nullptr); bool UpdateCurrentLowerBoundOfNode(int node, int64_t new_lb, int64_t offset); @@ -186,6 +192,12 @@ class RoutingLinearSolverWrapper { // This function is meant to override the parameters of the solver. virtual void SetParameters(const std::string& parameters) = 0; + // Returns if the model is empty or not. + virtual bool ModelIsEmpty() const { return true; } + + // Prints an understandable view of the model. + virtual std::string PrintModel() const = 0; + // Adds a variable with bounds [lower_bound, upper_bound]. int AddVariable(int64_t lower_bound, int64_t upper_bound) { CHECK_LE(lower_bound, upper_bound); @@ -414,6 +426,10 @@ class RoutingGlopWrapper : public RoutingLinearSolverWrapper { lp_solver_.SetParameters(params); } + // Prints an understandable view of the model + // TODO(user): Improve output readability. + std::string PrintModel() const override { return linear_program_.Dump(); } + private: const bool is_relaxation_; glop::LinearProgram linear_program_; @@ -593,6 +609,11 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { DCHECK(false); } + bool ModelIsEmpty() const override { return model_.ByteSizeLong() == 0; } + + // Prints an understandable view of the model + std::string PrintModel() const override; + private: sat::CpModelProto model_; sat::CpSolverResponse response_; @@ -604,6 +625,8 @@ class RoutingCPSatWrapper : public RoutingLinearSolverWrapper { // Utility class used in Local/GlobalDimensionCumulOptimizer to set the linear // solver constraints and solve the problem. class DimensionCumulOptimizerCore { + using RouteDimensionTravelInfo = RoutingModel::RouteDimensionTravelInfo; + public: DimensionCumulOptimizerCore(const RoutingDimension* dimension, bool use_precedence_propagator); @@ -614,13 +637,29 @@ class DimensionCumulOptimizerCore { // feasibility is of interest). DimensionSchedulingStatus OptimizeSingleRoute( int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values, int64_t* cost, int64_t* transit_cost, bool clear_lp = true); + // Given some cumuls and breaks, computes the solution cost by solving the + // same model as in OptimizeSingleRoute() with the addition of constraints for + // cumuls and breaks. + DimensionSchedulingStatus ComputeSingleRouteSolutionCost( + int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, + RoutingLinearSolverWrapper* solver, + const std::vector& solution_cumul_values, + const std::vector& solution_break_values, int64_t* cost, + int64_t* transit_cost, int64_t* cost_offset = nullptr, + bool reuse_previous_model_if_possible = true, bool clear_lp = false, + bool clear_solution_constraints = true, + absl::Duration* const solve_duration = nullptr); + std::vector OptimizeSingleRouteWithResources( int vehicle, const std::function& next_accessor, const std::function& transit_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, const std::vector& resources, const std::vector& resource_indices, bool optimize_vehicle_costs, RoutingLinearSolverWrapper* solver, @@ -630,6 +669,8 @@ class DimensionCumulOptimizerCore { DimensionSchedulingStatus Optimize( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values, std::vector>* resource_indices_per_group, int64_t* cost, @@ -637,12 +678,15 @@ class DimensionCumulOptimizerCore { DimensionSchedulingStatus OptimizeAndPack( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values, std::vector>* resource_indices_per_group); DimensionSchedulingStatus OptimizeAndPackSingleRoute( int vehicle, const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, const RoutingModel::ResourceGroup::Resource* resource, RoutingLinearSolverWrapper* solver, std::vector* cumul_values, std::vector* break_values); @@ -654,14 +698,26 @@ class DimensionCumulOptimizerCore { // setting any constraints and solving. void InitOptimizer(RoutingLinearSolverWrapper* solver); + // Initializes the model for a single route and a given dimension and sets its + // constraints. + bool InitSingleRoute(int vehicle, + const std::function& next_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, + RoutingLinearSolverWrapper* solver, + std::vector* cumul_values, int64_t* cost, + int64_t* transit_cost, int64_t* cumul_offset, + int64_t* const cost_offset); + // Computes the minimum/maximum of cumuls for nodes on "route", and sets them // in current_route_[min|max]_cumuls_ respectively. + bool ExtractRouteCumulBounds(const std::vector& route, + int64_t cumul_offset); + + // Tighten the minimum/maximum of cumuls for nodes on "route" // If the propagator_ is not null, uses the bounds tightened by the - // propagator. - // Otherwise, the bounds are computed by going over the nodes on the route - // using the CP bounds, and the fixed transits are used to tighten them. - bool ComputeRouteCumulBounds(const std::vector& route, - const std::vector& fixed_transits, + // propagator. Otherwise, the minimum transits are used to tighten them. + bool TightenRouteCumulBounds(const std::vector& route, + const std::vector& min_transits, int64_t cumul_offset); // Sets the constraints for all nodes on "vehicle"'s route according to @@ -671,10 +727,20 @@ class DimensionCumulOptimizerCore { bool SetRouteCumulConstraints( int vehicle, const std::function& next_accessor, const std::function& transit_accessor, + const RouteDimensionTravelInfo& dimension_travel_info, int64_t cumul_offset, bool optimize_costs, RoutingLinearSolverWrapper* solver, int64_t* route_transit_cost, int64_t* route_cost_offset); + // Sets the constraints for all variables related to travel. Handles + // static or time-dependent travel values. + // Returns false if some infeasibility was detected, true otherwise. + bool SetRouteTravelConstraints( + const RouteDimensionTravelInfo& dimension_travel_info, + const std::vector& lp_slacks, + const std::vector& fixed_transit, + RoutingLinearSolverWrapper* solver); + // Sets the global constraints on the dimension, and adds global objective // cost coefficients if optimize_costs is true. // NOTE: When called, the call to this function MUST come after @@ -780,9 +846,28 @@ class LocalDimensionCumulOptimizer { // Returns false if the route is not feasible. DimensionSchedulingStatus ComputeRouteCumuls( int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, std::vector* optimal_cumuls, std::vector* optimal_breaks); + // Simple combination of ComputeRouteCumulCost() and ComputeRouteCumuls() + DimensionSchedulingStatus ComputeRouteCumulsAndCost( + int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, + std::vector* optimal_cumuls, + std::vector* optimal_breaks, int64_t* optimal_cost); + + // If feasible, computes the cost of a given route performed by a vehicle + // defined by its cumuls and breaks. + DimensionSchedulingStatus ComputeRouteSolutionCost( + int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, + const std::vector& solution_cumul_values, + const std::vector& solution_break_values, int64_t* solution_cost, + int64_t* cost_offset = nullptr, + bool reuse_previous_model_if_possible = false, bool clear_lp = true, + absl::Duration* solve_duration = nullptr); + // Similar to ComputeRouteCumuls, but also tries to pack the cumul values on // the route, such that the cost remains the same, the cumul of route end is // minimized, and then the cumul of the start of the route is maximized. @@ -790,6 +875,7 @@ class LocalDimensionCumulOptimizer { // time window. DimensionSchedulingStatus ComputePackedRouteCumuls( int vehicle, const std::function& next_accessor, + const RoutingModel::RouteDimensionTravelInfo& dimension_travel_info, const RoutingModel::ResourceGroup::Resource* resource, std::vector* packed_cumuls, std::vector* packed_breaks); @@ -822,6 +908,8 @@ class GlobalDimensionCumulOptimizer { // Returns false if the routes are not feasible. DimensionSchedulingStatus ComputeCumuls( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, std::vector* optimal_cumuls, std::vector* optimal_breaks, std::vector>* optimal_resource_indices_per_group); @@ -831,6 +919,8 @@ class GlobalDimensionCumulOptimizer { // minimized, and then the cumuls of the starts of the routes are maximized. DimensionSchedulingStatus ComputePackedCumuls( const std::function& next_accessor, + const std::vector& + dimension_travel_info_per_route, std::vector* packed_cumuls, std::vector* packed_breaks, std::vector>* resource_indices_per_group); @@ -843,56 +933,105 @@ class GlobalDimensionCumulOptimizer { DimensionCumulOptimizerCore optimizer_core_; }; -class ResourceAssignmentOptimizer { - public: - ResourceAssignmentOptimizer(const RoutingModel::ResourceGroup* resource_group, - LocalDimensionCumulOptimizer* optimizer, - LocalDimensionCumulOptimizer* mp_optimizer); +// Finds the approximate (*) min-cost (i.e. best) assignment of all vehicles +// v ∈ 'vehicles' to resources, i.e. indices in [0..num_resources), where the +// costs of assigning a vehicle v to a resource r is given by +// 'vehicle_to_resource_cost(v)[r]', unless 'vehicle_to_resource_cost(v)' is +// empty in which case vehicle v does not need a resource. +// +// Returns the cost of that optimal assignment, or -1 if it's infeasible. +// Moreover, if 'resource_indices' != nullptr, it assumes that its size is the +// global number of vehicles, and assigns its element #v with the resource r +// assigned to v, or -1 if none. +// +// (*) COST SCALING: When the costs are so large that they could possibly yield +// int64_t overflow, this method returns a *lower* bound of the actual optimal +// cost, and the assignment output in 'resource_indices' may be suboptimal if +// that lower bound isn't tight (but it should be very close). +// +// COMPLEXITY: in practice, should be roughly +// O(num_resources * vehicles.size() + resource_indices->size()). +int64_t ComputeBestVehicleToResourceAssignment( + std::vector vehicles, int num_resources, + std::function*(int)> + vehicle_to_resource_assignment_costs, + std::vector* resource_indices); - // Returns the cost resulting from the min-cost assignment of resources to - // (used) vehicles, or -1 if the assignment is impossible. - // For each vehicle v and resource r, the cost of assigning r to v is equal to - // - primary_vehicle_to_resource_assignment_costs[v][r] if - // primary_vehicle_to_resource_assignment_costs[v] is not empty, - // - secondary_vehicle_to_resource_assignment_costs[v][r] otherwise - // (secondary_vehicle_to_resource_assignment_costs[v] can never be empty). - // If non-null, 'resource_indices' contains the index of the resource assigned - // to each vehicle. - // TODO(user): Return std::optional to catch infeasibilities. - int64_t ComputeBestAssignmentCost( - const std::vector>& - primary_vehicle_to_resource_assignment_costs, - const std::vector>& - secondary_vehicle_to_resource_assignment_costs, - const std::function& use_primary_for_vehicle, - std::vector* resource_indices) const; +// Computes the vehicle-to-resource assignment costs for the given vehicle to +// all resources in the group, and sets these costs in 'assignment_costs' (if +// non-null). The latter is cleared and kept empty if the vehicle 'v' should not +// have a resource assigned to it. +// optimize_vehicle_costs indicates if the costs should be optimized or if +// we merely care about feasibility (cost of 0) and infeasibility (cost of -1) +// of the assignments. +// The cumul and break values corresponding to the assignment of each resource +// are also set in cumul_values and break_values, if non-null. +bool ComputeVehicleToResourcesAssignmentCosts( + int v, const RoutingModel::ResourceGroup& resource_group, + const std::function& next_accessor, + const std::function& transit_accessor, + bool optimize_vehicle_costs, LocalDimensionCumulOptimizer* lp_optimizer, + LocalDimensionCumulOptimizer* mp_optimizer, + std::vector* assignment_costs, + std::vector>* cumul_values, + std::vector>* break_values); - // Computes the vehicle-to-resource assignment costs for the given vehicle to - // all resources in the group, and sets these costs in assignment_costs (if - // non-null). - // optimize_vehicle_costs indicates if the costs should be optimized or if - // we merely care about feasibility (cost of 0) and infeasibility (cost of -1) - // of the assignments. - // The cumul and break values corresponding to the assignment of each resource - // are also set in cumul_values and break_values, if non-null. - bool ComputeAssignmentCostsForVehicle( - int v, const std::function& next_accessor, - const std::function& transit_accessor, - bool optimize_vehicle_costs, std::vector* assignment_costs, - std::vector>* cumul_values, - std::vector>* break_values); - - const RoutingDimension* const dimension() const { - return optimizer_.dimension(); - } - - private: - LocalDimensionCumulOptimizer& optimizer_; - LocalDimensionCumulOptimizer& mp_optimizer_; - const RoutingModel& model_; - const RoutingModel::ResourceGroup& resource_group_; +// Simple struct returned by ComputePiecewiseLinearFormulationValue() to +// indicate if the value could be computed and if not, on what side the value +// was from the definition interval. +enum class PiecewiseEvaluationStatus { + UNSPECIFIED = 0, + WITHIN_BOUNDS, + SMALLER_THAN_LOWER_BOUND, + LARGER_THAN_UPPER_BOUND }; +// Computes pwl(x) for pwl a PieceWiseLinearFormulation. +// Returns a PieceWiseEvaluationStatus to indicate if the value could be +// computed (filled in value) and if not, why. +PiecewiseEvaluationStatus ComputePiecewiseLinearFormulationValue( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl, + int64_t x, int64_t* value, double delta = 0); + +// Like ComputePiecewiseLinearFormulationValue(), computes pwl(x) for pwl a +// PiecewiseLinearFormulation. For convex PiecewiseLinearFormulations, if x is +// outside the bounds of the function, instead of returning an error like in +// PiecewiseLinearFormulation, the function will still be defined by its outer +// segments. +int64_t ComputeConvexPiecewiseLinearFormulationValue( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl, + int64_t x, double delta = 0); + +// Structure to store the slope and y_intercept of a segment. +struct SlopeAndYIntercept { + double slope; + double y_intercept; + + friend ::std::ostream& operator<<(::std::ostream& os, + const SlopeAndYIntercept& it) { + return os << "{" << it.slope << ", " << it.y_intercept << "}"; + } +}; + +// Converts a vector of SlopeAndYIntercept to a vector of convexity regions. +// Convexity regions are defined such that, all segment in a convexity region +// form a convex function. The boolean in the vector is set to true if the +// segment associated to it starts a new convexity region. Therefore, a convex +// function would yield {true, false, false, ...} and a concave function would +// yield {true, true, true, ...}. +std::vector SlopeAndYInterceptToConvexityRegions( + const std::vector& slope_and_y_intercept); + +// Given a PiecewiseLinearFormulation, returns a vector of slope and y-intercept +// corresponding to each segment. Only the segments in [index_start, index_end[ +// will be considered. +std::vector PiecewiseLinearFormulationToSlopeAndYIntercept( + const RoutingModel::RouteDimensionTravelInfo::TransitionInfo:: + PiecewiseLinearFormulation& pwl_function, + int index_start = 0, int index_end = -1); + } // namespace operations_research #endif // OR_TOOLS_CONSTRAINT_SOLVER_ROUTING_LP_SCHEDULING_H_ diff --git a/ortools/constraint_solver/routing_sat.cc b/ortools/constraint_solver/routing_sat.cc index 0cf822086a..e2993fde4a 100644 --- a/ortools/constraint_solver/routing_sat.cc +++ b/ortools/constraint_solver/routing_sat.cc @@ -527,7 +527,7 @@ std::vector CreateGeneralizedRanks(const RoutingModel& model, } for (const auto [arc, arc_var] : arc_vars) { const auto [cp_tail, cp_head] = arc; - if (model.IsEnd(cp_head - 1)) continue; + if (cp_head == 0 || model.IsEnd(cp_head - 1)) continue; if (cp_tail == cp_head || cp_head == depot) continue; // arc[tail][head] -> ranks[head] == ranks[tail] + 1. AddLinearConstraint(cp_model, 1, 1, diff --git a/ortools/constraint_solver/routing_search.cc b/ortools/constraint_solver/routing_search.cc index 2297f33782..16269e3d80 100644 --- a/ortools/constraint_solver/routing_search.cc +++ b/ortools/constraint_solver/routing_search.cc @@ -665,20 +665,6 @@ int64_t CheapestInsertionFilteredHeuristic::GetUnperformedValue( return std::numeric_limits::max(); } -namespace { -template -void SortAndExtractPairSeconds(std::vector>* pairs, - std::vector* sorted_seconds) { - CHECK(pairs != nullptr); - CHECK(sorted_seconds != nullptr); - std::sort(pairs->begin(), pairs->end()); - sorted_seconds->reserve(pairs->size()); - for (const std::pair& p : *pairs) { - sorted_seconds->push_back(p.second); - } -} -} // namespace - // Priority queue entries used by global cheapest insertion heuristic. // Entry in priority queue containing the insertion positions of a node pair. @@ -4680,7 +4666,7 @@ DecisionBuilder* RoutingModel::MakeSelfDependentDimensionFinalizer( solver_->MakeSolveOnce(guided_finalizer); std::vector start_cumuls(vehicles_, nullptr); for (int64_t vehicle_idx = 0; vehicle_idx < vehicles_; ++vehicle_idx) { - start_cumuls[vehicle_idx] = dimension->CumulVar(starts_[vehicle_idx]); + start_cumuls[vehicle_idx] = dimension->CumulVar(Start(vehicle_idx)); } LocalSearchOperator* const hill_climber = solver_->RevAlloc(new GreedyDescentLSOperator(start_cumuls)); diff --git a/ortools/constraint_solver/search.cc b/ortools/constraint_solver/search.cc index a43991a053..ed41ff6ffd 100644 --- a/ortools/constraint_solver/search.cc +++ b/ortools/constraint_solver/search.cc @@ -408,6 +408,7 @@ class AtSolutionCallback : public SearchMonitor { : SearchMonitor(solver), callback_(std::move(callback)) {} ~AtSolutionCallback() override {} bool AtSolution() override; + void Install() override; private: const std::function callback_; @@ -418,6 +419,10 @@ bool AtSolutionCallback::AtSolution() { return false; } +void AtSolutionCallback::Install() { + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + } // namespace SearchMonitor* Solver::MakeAtSolutionCallback(std::function callback) { @@ -431,6 +436,7 @@ class EnterSearchCallback : public SearchMonitor { : SearchMonitor(solver), callback_(std::move(callback)) {} ~EnterSearchCallback() override {} void EnterSearch() override; + void Install() override; private: const std::function callback_; @@ -438,6 +444,10 @@ class EnterSearchCallback : public SearchMonitor { void EnterSearchCallback::EnterSearch() { callback_(); } +void EnterSearchCallback::Install() { + ListenToEvent(Solver::MonitorEvent::kEnterSearch); +} + } // namespace SearchMonitor* Solver::MakeEnterSearchCallback(std::function callback) { @@ -451,6 +461,7 @@ class ExitSearchCallback : public SearchMonitor { : SearchMonitor(solver), callback_(std::move(callback)) {} ~ExitSearchCallback() override {} void ExitSearch() override; + void Install() override; private: const std::function callback_; @@ -458,6 +469,10 @@ class ExitSearchCallback : public SearchMonitor { void ExitSearchCallback::ExitSearch() { callback_(); } +void ExitSearchCallback::Install() { + ListenToEvent(Solver::MonitorEvent::kExitSearch); +} + } // namespace SearchMonitor* Solver::MakeExitSearchCallback(std::function callback) { @@ -1493,7 +1508,7 @@ int64_t StaticEvaluatorSelector::SelectValue(const IntVar* var, int64_t id) { int64_t StaticEvaluatorSelector::ChooseVariable() { if (first_ == -1) { // first call to select. update assignment costs - // Two phases: compute size then filland sort + // Two phases: compute size then fill and sort int64_t element_size = 0; for (int64_t i = 0; i < vars_.size(); ++i) { if (!vars_[i]->Bound()) { @@ -2279,6 +2294,10 @@ SolutionCollector::~SolutionCollector() { gtl::STLDeleteElements(&recycle_solutions_); } +void SolutionCollector::Install() { + ListenToEvent(Solver::MonitorEvent::kEnterSearch); +} + void SolutionCollector::Add(IntVar* const var) { if (prototype_ != nullptr) { prototype_->Add(var); @@ -2451,6 +2470,7 @@ class FirstSolutionCollector : public SolutionCollector { ~FirstSolutionCollector() override; void EnterSearch() override; bool AtSolution() override; + void Install() override; std::string DebugString() const override; private: @@ -2479,6 +2499,11 @@ bool FirstSolutionCollector::AtSolution() { return false; } +void FirstSolutionCollector::Install() { + SolutionCollector::Install(); + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + std::string FirstSolutionCollector::DebugString() const { if (prototype_ == nullptr) { return "FirstSolutionCollector()"; @@ -2507,6 +2532,7 @@ class LastSolutionCollector : public SolutionCollector { explicit LastSolutionCollector(Solver* const s); ~LastSolutionCollector() override; bool AtSolution() override; + void Install() override; std::string DebugString() const override; }; @@ -2525,6 +2551,11 @@ bool LastSolutionCollector::AtSolution() { return true; } +void LastSolutionCollector::Install() { + SolutionCollector::Install(); + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + std::string LastSolutionCollector::DebugString() const { if (prototype_ == nullptr) { return "LastSolutionCollector()"; @@ -2554,6 +2585,7 @@ class BestValueSolutionCollector : public SolutionCollector { ~BestValueSolutionCollector() override {} void EnterSearch() override; bool AtSolution() override; + void Install() override; std::string DebugString() const override; public: @@ -2600,6 +2632,11 @@ bool BestValueSolutionCollector::AtSolution() { return true; } +void BestValueSolutionCollector::Install() { + SolutionCollector::Install(); + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + std::string BestValueSolutionCollector::DebugString() const { if (prototype_ == nullptr) { return "BestValueSolutionCollector()"; @@ -2632,6 +2669,7 @@ class NBestValueSolutionCollector : public SolutionCollector { void EnterSearch() override; void ExitSearch() override; bool AtSolution() override; + void Install() override; std::string DebugString() const override; public: @@ -2696,6 +2734,12 @@ bool NBestValueSolutionCollector::AtSolution() { return true; } +void NBestValueSolutionCollector::Install() { + SolutionCollector::Install(); + ListenToEvent(Solver::MonitorEvent::kExitSearch); + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + std::string NBestValueSolutionCollector::DebugString() const { if (prototype_ == nullptr) { return "NBestValueSolutionCollector()"; @@ -2741,6 +2785,7 @@ class AllSolutionCollector : public SolutionCollector { explicit AllSolutionCollector(Solver* const s); ~AllSolutionCollector() override; bool AtSolution() override; + void Install() override; std::string DebugString() const override; }; @@ -2758,6 +2803,11 @@ bool AllSolutionCollector::AtSolution() { return true; } +void AllSolutionCollector::Install() { + SolutionCollector::Install(); + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + std::string AllSolutionCollector::DebugString() const { if (prototype_ == nullptr) { return "AllSolutionCollector()"; @@ -3099,11 +3149,9 @@ class TabuSearch : public Metaheuristic { protected: struct VarValue { - VarValue(IntVar* const var, int64_t value, int64_t stamp) - : var_(var), value_(value), stamp_(stamp) {} - IntVar* const var_; - const int64_t value_; - const int64_t stamp_; + IntVar* const var; + const int64_t value; + const int64_t stamp; }; typedef std::list TabuList; @@ -3210,11 +3258,11 @@ std::vector TabuSearch::CreateTabuVars() { // the tabu criterion which is tolerated; a factor of 1 means no violations // allowed, a factor of 0 means all violations allowed. std::vector tabu_vars; - for (const VarValue& vv : keep_tabu_list_) { - tabu_vars.push_back(s->MakeIsEqualCstVar(vv.var_, vv.value_)); + for (const auto [var, value, unused_stamp] : keep_tabu_list_) { + tabu_vars.push_back(s->MakeIsEqualCstVar(var, value)); } - for (const VarValue& vv : forbid_tabu_list_) { - tabu_vars.push_back(s->MakeIsDifferentCstVar(vv.var_, vv.value_)); + for (const auto [var, value, unused_stamp] : forbid_tabu_list_) { + tabu_vars.push_back(s->MakeIsDifferentCstVar(var, value)); } return tabu_vars; } @@ -3235,12 +3283,10 @@ bool TabuSearch::AtSolution() { const int64_t new_value = var->Value(); if (old_value != new_value) { if (keep_tenure_ > 0) { - VarValue keep_value(var, new_value, stamp_); - keep_tabu_list_.push_front(keep_value); + keep_tabu_list_.push_front({var, new_value, stamp_}); } if (forbid_tenure_ > 0) { - VarValue forbid_value(var, old_value, stamp_); - forbid_tabu_list_.push_front(forbid_value); + forbid_tabu_list_.push_front({var, old_value, stamp_}); } } } @@ -3267,7 +3313,7 @@ void TabuSearch::AcceptNeighbor() { } void TabuSearch::AgeList(int64_t tenure, TabuList* list) { - while (!list->empty() && list->back().stamp_ < stamp_ - tenure) { + while (!list->empty() && list->back().stamp < stamp_ - tenure) { list->pop_back(); } } @@ -3296,8 +3342,8 @@ std::vector GenericTabuSearch::CreateTabuVars() { // Tabu criterion // At least one element of the forbid_tabu_list must change value. std::vector forbid_values; - for (const VarValue& vv : forbid_tabu_list()) { - forbid_values.push_back(s->MakeIsDifferentCstVar(vv.var_, vv.value_)); + for (const auto [var, value, unused_stamp] : forbid_tabu_list()) { + forbid_values.push_back(s->MakeIsDifferentCstVar(var, value)); } std::vector tabu_vars; if (!forbid_values.empty()) { @@ -3432,45 +3478,39 @@ SearchMonitor* Solver::MakeSimulatedAnnealing(bool maximize, IntVar* const v, // ---------- Guided Local Search ---------- -typedef std::pair Arc; - namespace { -// Base GLS penalties abstract class. Maintains the penalty frequency for each -// (variable, value) arc. -class GuidedLocalSearchPenalties { - public: - virtual ~GuidedLocalSearchPenalties() {} - virtual bool HasValues() const = 0; - virtual void Increment(const Arc& arc) = 0; - virtual int64_t Value(const Arc& arc) const = 0; - virtual void Reset() = 0; -}; +// GLS penalty management classes. Maintains the penalty frequency for each +// (variable, value) pair. // Dense GLS penalties implementation using a matrix to store penalties. -class GuidedLocalSearchPenaltiesTable : public GuidedLocalSearchPenalties { +class GuidedLocalSearchPenaltiesTable { public: - explicit GuidedLocalSearchPenaltiesTable(int size); - ~GuidedLocalSearchPenaltiesTable() override {} - bool HasValues() const override { return has_values_; } - void Increment(const Arc& arc) override; - int64_t Value(const Arc& arc) const override; - void Reset() override; + struct VarValue { + int64_t var; + int64_t value; + }; + explicit GuidedLocalSearchPenaltiesTable(int num_vars); + bool HasPenalties() const { return has_values_; } + void IncrementPenalty(const VarValue& var_value); + int64_t GetPenalty(const VarValue& var_value) const; + void Reset(); private: std::vector> penalties_; bool has_values_; }; -GuidedLocalSearchPenaltiesTable::GuidedLocalSearchPenaltiesTable(int size) - : penalties_(size), has_values_(false) {} +GuidedLocalSearchPenaltiesTable::GuidedLocalSearchPenaltiesTable(int num_vars) + : penalties_(num_vars), has_values_(false) {} -void GuidedLocalSearchPenaltiesTable::Increment(const Arc& arc) { - std::vector& first_penalties = penalties_[arc.first]; - const int64_t second = arc.second; - if (second >= first_penalties.size()) { - first_penalties.resize(second + 1, 0); +void GuidedLocalSearchPenaltiesTable::IncrementPenalty( + const VarValue& var_value) { + std::vector& var_penalties = penalties_[var_value.var]; + const int64_t value = var_value.value; + if (value >= var_penalties.size()) { + var_penalties.resize(value + 1, 0); } - ++first_penalties[second]; + ++var_penalties[value]; has_values_ = true; } @@ -3481,37 +3521,46 @@ void GuidedLocalSearchPenaltiesTable::Reset() { } } -int64_t GuidedLocalSearchPenaltiesTable::Value(const Arc& arc) const { - const std::vector& first_penalties = penalties_[arc.first]; - const int64_t second = arc.second; - if (second >= first_penalties.size()) { - return 0; - } else { - return first_penalties[second]; - } +int64_t GuidedLocalSearchPenaltiesTable::GetPenalty( + const VarValue& var_value) const { + const std::vector& var_penalties = penalties_[var_value.var]; + const int64_t value = var_value.value; + return (value >= var_penalties.size()) ? 0 : var_penalties[value]; } // Sparse GLS penalties implementation using hash_map to store penalties. -class GuidedLocalSearchPenaltiesMap : public GuidedLocalSearchPenalties { +class GuidedLocalSearchPenaltiesMap { public: - explicit GuidedLocalSearchPenaltiesMap(int size); - ~GuidedLocalSearchPenaltiesMap() override {} - bool HasValues() const override { return (!penalties_.empty()); } - void Increment(const Arc& arc) override; - int64_t Value(const Arc& arc) const override; - void Reset() override; + struct VarValue { + int64_t var; + int64_t value; + + friend bool operator==(const VarValue& lhs, const VarValue& rhs) { + return lhs.var == rhs.var && lhs.value == rhs.value; + } + template + friend H AbslHashValue(H h, const VarValue& var_value) { + return H::combine(std::move(h), var_value.var, var_value.value); + } + }; + explicit GuidedLocalSearchPenaltiesMap(int num_vars); + bool HasPenalties() const { return (!penalties_.empty()); } + void IncrementPenalty(const VarValue& var_value); + int64_t GetPenalty(const VarValue& var_value) const; + void Reset(); private: Bitmap penalized_; - absl::flat_hash_map penalties_; + absl::flat_hash_map penalties_; }; -GuidedLocalSearchPenaltiesMap::GuidedLocalSearchPenaltiesMap(int size) - : penalized_(size, false) {} +GuidedLocalSearchPenaltiesMap::GuidedLocalSearchPenaltiesMap(int num_vars) + : penalized_(num_vars, false) {} -void GuidedLocalSearchPenaltiesMap::Increment(const Arc& arc) { - ++penalties_[arc]; - penalized_.Set(arc.first, true); +void GuidedLocalSearchPenaltiesMap::IncrementPenalty( + const VarValue& var_value) { + ++penalties_[var_value]; + penalized_.Set(var_value.var, true); } void GuidedLocalSearchPenaltiesMap::Reset() { @@ -3519,13 +3568,14 @@ void GuidedLocalSearchPenaltiesMap::Reset() { penalized_.Clear(); } -int64_t GuidedLocalSearchPenaltiesMap::Value(const Arc& arc) const { - if (penalized_.Get(arc.first)) { - return gtl::FindWithDefault(penalties_, arc); - } - return 0; +int64_t GuidedLocalSearchPenaltiesMap::GetPenalty( + const VarValue& var_value) const { + return (penalized_.Get(var_value.var)) + ? gtl::FindWithDefault(penalties_, var_value) + : 0; } +template class GuidedLocalSearch : public Metaheuristic { public: GuidedLocalSearch(Solver* const s, IntVar* objective, bool maximize, @@ -3537,68 +3587,124 @@ class GuidedLocalSearch : public Metaheuristic { bool AtSolution() override; void EnterSearch() override; bool LocalOptimum() override; - virtual int64_t AssignmentElementPenalty(const Assignment& assignment, - int index) = 0; - virtual int64_t AssignmentPenalty(const Assignment& assignment, int index, - int64_t next) = 0; - virtual bool EvaluateElementValue(const Assignment::IntContainer& container, - int64_t index, int* container_index, - int64_t* penalty) = 0; + virtual int64_t AssignmentElementPenalty(int index) const = 0; + virtual int64_t AssignmentPenalty(int64_t var, int64_t value) const = 0; + virtual int64_t Evaluate(const Assignment* delta, int64_t current_penalty, + bool incremental) = 0; virtual IntExpr* MakeElementPenalty(int index) = 0; std::string DebugString() const override { return "Guided Local Search"; } protected: - struct Comparator { - bool operator()(const std::pair& i, - const std::pair& j) { - return i.second > j.second; + // Array which keeps track of modifications done. This allows to effectively + // revert or commit modifications. + // TODO(user): Expose this in a utility file. + template + class DirtyArray { + public: + explicit DirtyArray(IndexType size) + : base_data_(size), modified_data_(size), touched_(size) {} + // Sets a value in the array. This value will be reverted if Revert() is + // called. + void Set(IndexType i, const T& value) { + modified_data_[i] = value; + touched_.Set(i); } + // Same as Set() but modifies all values of the array. + void SetAll(const T& value) { + for (IndexType i = 0; i < modified_data_.size(); ++i) { + Set(i, value); + } + } + // Returns the modified value in the array. + T Get(IndexType i) const { return modified_data_[i]; } + // Commits all modifications done to the array, effectively copying all + // modifications to the base values. + void Commit() { + for (const IndexType index : touched_.PositionsSetAtLeastOnce()) { + base_data_[index] = modified_data_[index]; + } + touched_.SparseClearAll(); + } + // Reverts all modified values in the array. + void Revert() { + for (const IndexType index : touched_.PositionsSetAtLeastOnce()) { + modified_data_[index] = base_data_[index]; + } + touched_.SparseClearAll(); + } + // Returns the number of values modified since the last call to Commit or + // Revert. + int NumSetValues() const { + return touched_.NumberOfSetCallsWithDifferentArguments(); + } + + private: + std::vector base_data_; + std::vector modified_data_; + SparseBitset touched_; }; - int64_t Evaluate(const Assignment* delta, int64_t current_penalty, - const int64_t* const out_values, bool cache_delta_values); + int64_t GetValue(int64_t index) const { + return assignment_.Element(index).Value(); + } + IntVar* GetVar(int64_t index) const { + return assignment_.Element(index).Var(); + } + void AddVars(const std::vector& vars); + int NumPrimaryVars() const { return num_vars_; } + int GetLocalIndexFromVar(IntVar* var) const { + const int var_index = var->index(); + return (var_index < var_index_to_local_index_.size()) + ? var_index_to_local_index_[var_index] + : -1; + } IntVar* penalized_objective_; - Assignment assignment_; + Assignment::IntContainer assignment_; int64_t assignment_penalized_value_; int64_t old_penalized_value_; - const std::vector vars_; - absl::flat_hash_map indices_; + const int num_vars_; + std::vector var_index_to_local_index_; const double penalty_factor_; - std::unique_ptr penalties_; - std::unique_ptr current_penalized_values_; - std::unique_ptr delta_cache_; + P penalties_; + DirtyArray penalized_values_; bool incremental_; }; -GuidedLocalSearch::GuidedLocalSearch(Solver* const s, IntVar* objective, - bool maximize, int64_t step, - const std::vector& vars, - double penalty_factor) +template +GuidedLocalSearch

::GuidedLocalSearch(Solver* const s, IntVar* objective, + bool maximize, int64_t step, + const std::vector& vars, + double penalty_factor) : Metaheuristic(s, maximize, objective, step), penalized_objective_(nullptr), - assignment_(s), assignment_penalized_value_(0), old_penalized_value_(0), - vars_(vars), + num_vars_(vars.size()), penalty_factor_(penalty_factor), + penalties_(vars.size()), + penalized_values_(vars.size()), incremental_(false) { - if (!vars.empty()) { - // TODO(user): Remove scoped_array. - assignment_.Add(vars_); - current_penalized_values_ = std::make_unique(vars_.size()); - delta_cache_ = std::make_unique(vars_.size()); - memset(current_penalized_values_.get(), 0, - vars_.size() * sizeof(*current_penalized_values_.get())); + AddVars(vars); +} + +template +void GuidedLocalSearch

::AddVars(const std::vector& vars) { + const int offset = assignment_.Size(); + if (vars.empty()) return; + assignment_.Resize(offset + vars.size()); + for (int i = 0; i < vars.size(); ++i) { + assignment_.AddAtPosition(vars[i], offset + i); } - for (int i = 0; i < vars_.size(); ++i) { - indices_[vars_[i]] = i; + const int max_var_index = + (*std::max_element(vars.begin(), vars.end(), [](IntVar* a, IntVar* b) { + return a->index() < b->index(); + }))->index(); + if (max_var_index >= var_index_to_local_index_.size()) { + var_index_to_local_index_.resize(max_var_index + 1, -1); } - if (absl::GetFlag(FLAGS_cp_use_sparse_gls_penalties)) { - penalties_ = std::make_unique(vars_.size()); - } else { - penalties_ = - std::make_unique(vars_.size()); + for (int i = 0; i < vars.size(); ++i) { + var_index_to_local_index_[vars[i]->index()] = offset + i; } } @@ -3609,26 +3715,27 @@ GuidedLocalSearch::GuidedLocalSearch(Solver* const s, IntVar* objective, // if maximizing, // objective >= Min(current penalized cost - penalized_objective + step, // best solution cost + step) -void GuidedLocalSearch::ApplyDecision(Decision* const d) { +template +void GuidedLocalSearch

::ApplyDecision(Decision* const d) { if (d == solver()->balancing_decision()) { return; } assignment_penalized_value_ = 0; - if (penalties_->HasValues()) { + if (penalties_.HasPenalties()) { // Computing sum of penalties expression. // Scope needed to avoid potential leak of elements. { std::vector elements; - for (int i = 0; i < vars_.size(); ++i) { + for (int i = 0; i < num_vars_; ++i) { elements.push_back(MakeElementPenalty(i)->Var()); - const int64_t penalty = AssignmentElementPenalty(assignment_, i); - current_penalized_values_[i] = penalty; - delta_cache_[i] = penalty; + const int64_t penalty = AssignmentElementPenalty(i); + penalized_values_.Set(i, penalty); assignment_penalized_value_ = CapAdd(assignment_penalized_value_, penalty); } penalized_objective_ = solver()->MakeSum(elements)->Var(); } + penalized_values_.Commit(); old_penalized_value_ = assignment_penalized_value_; incremental_ = false; if (maximize_) { @@ -3659,124 +3766,101 @@ void GuidedLocalSearch::ApplyDecision(Decision* const d) { } } -bool GuidedLocalSearch::AtSolution() { +template +bool GuidedLocalSearch

::AtSolution() { if (!Metaheuristic::AtSolution()) { return false; } if (penalized_objective_ != nullptr) { // In case no move has been found - current_ += penalized_objective_->Value(); + current_ = CapAdd(current_, penalized_objective_->Value()); } assignment_.Store(); return true; } -void GuidedLocalSearch::EnterSearch() { +template +void GuidedLocalSearch

::EnterSearch() { Metaheuristic::EnterSearch(); penalized_objective_ = nullptr; assignment_penalized_value_ = 0; old_penalized_value_ = 0; - memset(current_penalized_values_.get(), 0, - vars_.size() * sizeof(*current_penalized_values_.get())); - penalties_->Reset(); + penalized_values_.SetAll(0); + penalized_values_.Commit(); + penalties_.Reset(); } // GLS filtering; compute the penalized value corresponding to the delta and // modify objective bound accordingly. -bool GuidedLocalSearch::AcceptDelta(Assignment* delta, Assignment* deltadelta) { - if (delta != nullptr || deltadelta != nullptr) { - if (!penalties_->HasValues()) { - return Metaheuristic::AcceptDelta(delta, deltadelta); - } - int64_t penalty = 0; - if (!deltadelta->Empty()) { - if (!incremental_) { - penalty = Evaluate(delta, assignment_penalized_value_, - current_penalized_values_.get(), true); - } else { - penalty = Evaluate(deltadelta, old_penalized_value_, delta_cache_.get(), - true); - } - incremental_ = true; +template +bool GuidedLocalSearch

::AcceptDelta(Assignment* delta, + Assignment* deltadelta) { + if (delta == nullptr && deltadelta == nullptr) return true; + if (!penalties_.HasPenalties()) { + return Metaheuristic::AcceptDelta(delta, deltadelta); + } + int64_t penalty = 0; + if (!deltadelta->Empty()) { + if (!incremental_) { + DCHECK_EQ(penalized_values_.NumSetValues(), 0); + penalty = Evaluate(delta, assignment_penalized_value_, true); } else { - if (incremental_) { - for (int i = 0; i < vars_.size(); ++i) { - delta_cache_[i] = current_penalized_values_[i]; - } - old_penalized_value_ = assignment_penalized_value_; - } - incremental_ = false; - penalty = Evaluate(delta, assignment_penalized_value_, - current_penalized_values_.get(), false); + penalty = Evaluate(deltadelta, old_penalized_value_, true); } - old_penalized_value_ = penalty; - if (!delta->HasObjective()) { - delta->AddObjective(objective_); + incremental_ = true; + } else { + if (incremental_) { + penalized_values_.Revert(); } - if (delta->Objective() == objective_) { - if (maximize_) { - delta->SetObjectiveMin( - std::max(std::min(CapSub(CapAdd(current_, step_), penalty), - CapAdd(best_, step_)), - delta->ObjectiveMin())); - } else { - delta->SetObjectiveMax( - std::min(std::max(CapSub(CapSub(current_, step_), penalty), - CapSub(best_, step_)), - delta->ObjectiveMax())); - } + incremental_ = false; + DCHECK_EQ(penalized_values_.NumSetValues(), 0); + penalty = Evaluate(delta, assignment_penalized_value_, false); + } + old_penalized_value_ = penalty; + if (!delta->HasObjective()) { + delta->AddObjective(objective_); + } + if (delta->Objective() == objective_) { + if (maximize_) { + delta->SetObjectiveMin( + std::max(std::min(CapSub(CapAdd(current_, step_), penalty), + CapAdd(best_, step_)), + delta->ObjectiveMin())); + } else { + delta->SetObjectiveMax( + std::min(std::max(CapSub(CapSub(current_, step_), penalty), + CapSub(best_, step_)), + delta->ObjectiveMax())); } } return true; } -int64_t GuidedLocalSearch::Evaluate(const Assignment* delta, - int64_t current_penalty, - const int64_t* const out_values, - bool cache_delta_values) { - int64_t penalty = current_penalty; - const Assignment::IntContainer& container = delta->IntVarContainer(); - const int size = container.Size(); - for (int i = 0; i < size; ++i) { - const IntVarElement& new_element = container.Element(i); - IntVar* var = new_element.Var(); - int64_t index = -1; - if (gtl::FindCopy(indices_, var, &index)) { - penalty = CapSub(penalty, out_values[index]); - int64_t new_penalty = 0; - if (EvaluateElementValue(container, index, &i, &new_penalty)) { - penalty = CapAdd(penalty, new_penalty); - if (cache_delta_values) { - delta_cache_[index] = new_penalty; - } - } - } - } - return penalty; -} - -// Penalize all the most expensive arcs (var, value) according to their utility: -// utility(i, j) = cost(i, j) / (1 + penalty(i, j)) -bool GuidedLocalSearch::LocalOptimum() { - std::vector> utility(vars_.size()); - for (int i = 0; i < vars_.size(); ++i) { - if (!assignment_.Bound(vars_[i])) { +// Penalize (var, value) pairs of maximum utility, with +// utility(var, value) = cost(var, value) / (1 + penalty(var, value)) +template +bool GuidedLocalSearch

::LocalOptimum() { + std::vector utilities(num_vars_); + double max_utility = -std::numeric_limits::infinity(); + for (int var = 0; var < num_vars_; ++var) { + const IntVarElement& element = assignment_.Element(var); + if (!element.Bound()) { // Never synced with a solution, problem infeasible. return false; } - const int64_t var_value = assignment_.Value(vars_[i]); - const int64_t value = - (var_value != i) ? AssignmentPenalty(assignment_, i, var_value) : 0; - const Arc arc(i, var_value); - const int64_t penalty = penalties_->Value(arc); - utility[i] = std::pair(arc, value / (penalty + 1.0)); + const int64_t value = element.Value(); + // The fact that we do not penalize loops is influenced by vehicle routing. + // Assuming a cost of 0 in that case. + const int64_t cost = (value != var) ? AssignmentPenalty(var, value) : 0; + const double utility = cost / (penalties_.GetPenalty({var, value}) + 1.0); + utilities[var] = utility; + if (utility > max_utility) max_utility = utility; } - Comparator comparator; - std::sort(utility.begin(), utility.end(), comparator); - int64_t utility_value = utility[0].second; - penalties_->Increment(utility[0].first); - for (int i = 1; i < utility.size() && utility_value == utility[i].second; - ++i) { - penalties_->Increment(utility[i].first); + for (int var = 0; var < num_vars_; ++var) { + if (utilities[var] == max_utility) { + const IntVarElement& element = assignment_.Element(var); + DCHECK(element.Bound()); + penalties_.IncrementPenalty({var, element.Value()}); + } } if (maximize_) { current_ = std::numeric_limits::min(); @@ -3786,7 +3870,8 @@ bool GuidedLocalSearch::LocalOptimum() { return true; } -class BinaryGuidedLocalSearch : public GuidedLocalSearch { +template +class BinaryGuidedLocalSearch : public GuidedLocalSearch

{ public: BinaryGuidedLocalSearch( Solver* const solver, IntVar* const objective, @@ -3795,76 +3880,81 @@ class BinaryGuidedLocalSearch : public GuidedLocalSearch { double penalty_factor); ~BinaryGuidedLocalSearch() override {} IntExpr* MakeElementPenalty(int index) override; - int64_t AssignmentElementPenalty(const Assignment& assignment, - int index) override; - int64_t AssignmentPenalty(const Assignment& assignment, int index, - int64_t next) override; - bool EvaluateElementValue(const Assignment::IntContainer& container, - int64_t index, int* container_index, - int64_t* penalty) override; + int64_t AssignmentElementPenalty(int index) const override; + int64_t AssignmentPenalty(int64_t var, int64_t value) const override; + int64_t Evaluate(const Assignment* delta, int64_t current_penalty, + bool incremental) override; private: - int64_t PenalizedValue(int64_t i, int64_t j); + int64_t PenalizedValue(int64_t i, int64_t j) const; std::function objective_function_; }; -BinaryGuidedLocalSearch::BinaryGuidedLocalSearch( +template +BinaryGuidedLocalSearch

::BinaryGuidedLocalSearch( Solver* const solver, IntVar* const objective, std::function objective_function, bool maximize, int64_t step, const std::vector& vars, double penalty_factor) - : GuidedLocalSearch(solver, objective, maximize, step, vars, - penalty_factor), + : GuidedLocalSearch

(solver, objective, maximize, step, vars, + penalty_factor), objective_function_(std::move(objective_function)) {} -IntExpr* BinaryGuidedLocalSearch::MakeElementPenalty(int index) { - return solver()->MakeElement( +template +IntExpr* BinaryGuidedLocalSearch

::MakeElementPenalty(int index) { + return this->solver()->MakeElement( [this, index](int64_t i) { return PenalizedValue(index, i); }, - vars_[index]); + this->GetVar(index)); } -int64_t BinaryGuidedLocalSearch::AssignmentElementPenalty( - const Assignment& assignment, int index) { - return PenalizedValue(index, assignment.Value(vars_[index])); +template +int64_t BinaryGuidedLocalSearch

::AssignmentElementPenalty(int index) const { + return PenalizedValue(index, this->GetValue(index)); } -int64_t BinaryGuidedLocalSearch::AssignmentPenalty(const Assignment& assignment, - int index, int64_t next) { - return objective_function_(index, next); +template +int64_t BinaryGuidedLocalSearch

::AssignmentPenalty(int64_t var, + int64_t value) const { + return objective_function_(var, value); } -bool BinaryGuidedLocalSearch::EvaluateElementValue( - const Assignment::IntContainer& container, int64_t index, - int* container_index, int64_t* penalty) { - const IntVarElement& element = container.Element(*container_index); - if (element.Activated()) { - *penalty = PenalizedValue(index, element.Value()); - return true; +template +int64_t BinaryGuidedLocalSearch

::Evaluate(const Assignment* delta, + int64_t current_penalty, + bool incremental) { + int64_t penalty = current_penalty; + const Assignment::IntContainer& container = delta->IntVarContainer(); + for (const IntVarElement& new_element : container.elements()) { + const int index = this->GetLocalIndexFromVar(new_element.Var()); + if (index == -1) continue; + penalty = CapSub(penalty, this->penalized_values_.Get(index)); + if (new_element.Activated()) { + const int64_t new_penalty = PenalizedValue(index, new_element.Value()); + penalty = CapAdd(penalty, new_penalty); + if (incremental) { + this->penalized_values_.Set(index, new_penalty); + } + } } - return false; + return penalty; } // Penalized value for (i, j) = penalty_factor_ * penalty(i, j) * cost (i, j) -int64_t BinaryGuidedLocalSearch::PenalizedValue(int64_t i, int64_t j) { - const Arc arc(i, j); - const int64_t penalty = penalties_->Value(arc); - if (penalty != 0) { // objective_function_->Run(i, j) can be costly - const double penalized_value_fp = - penalty_factor_ * penalty * objective_function_(i, j); - const int64_t penalized_value = - (penalized_value_fp <= std::numeric_limits::max()) - ? static_cast(penalized_value_fp) - : std::numeric_limits::max(); - if (maximize_) { - return -penalized_value; - } else { - return penalized_value; - } - } else { - return 0; - } +template +int64_t BinaryGuidedLocalSearch

::PenalizedValue(int64_t i, int64_t j) const { + const int64_t penalty = this->penalties_.GetPenalty({i, j}); + // Calls to objective_function_(i, j) can be costly. + if (penalty == 0) return 0; + const double penalized_value_fp = + this->penalty_factor_ * penalty * objective_function_(i, j); + const int64_t penalized_value = + (penalized_value_fp <= std::numeric_limits::max()) + ? static_cast(penalized_value_fp) + : std::numeric_limits::max(); + return this->maximize_ ? -penalized_value : penalized_value; } -class TernaryGuidedLocalSearch : public GuidedLocalSearch { +template +class TernaryGuidedLocalSearch : public GuidedLocalSearch

{ public: TernaryGuidedLocalSearch( Solver* const solver, IntVar* const objective, @@ -3873,104 +3963,103 @@ class TernaryGuidedLocalSearch : public GuidedLocalSearch { const std::vector& secondary_vars, double penalty_factor); ~TernaryGuidedLocalSearch() override {} IntExpr* MakeElementPenalty(int index) override; - int64_t AssignmentElementPenalty(const Assignment& assignment, - int index) override; - int64_t AssignmentPenalty(const Assignment& assignment, int index, - int64_t next) override; - bool EvaluateElementValue(const Assignment::IntContainer& container, - int64_t index, int* container_index, - int64_t* penalty) override; + int64_t AssignmentElementPenalty(int index) const override; + int64_t AssignmentPenalty(int64_t var, int64_t value) const override; + int64_t Evaluate(const Assignment* delta, int64_t current_penalty, + bool incremental) override; private: - int64_t PenalizedValue(int64_t i, int64_t j, int64_t k); - int64_t GetAssignmentSecondaryValue(const Assignment::IntContainer& container, - int index, int* container_index) const; + int64_t PenalizedValue(int64_t i, int64_t j, int64_t k) const; - const std::vector secondary_vars_; std::function objective_function_; }; -TernaryGuidedLocalSearch::TernaryGuidedLocalSearch( +template +TernaryGuidedLocalSearch

::TernaryGuidedLocalSearch( Solver* const solver, IntVar* const objective, std::function objective_function, bool maximize, int64_t step, const std::vector& vars, const std::vector& secondary_vars, double penalty_factor) - : GuidedLocalSearch(solver, objective, maximize, step, vars, - penalty_factor), - secondary_vars_(secondary_vars), + : GuidedLocalSearch

(solver, objective, maximize, step, vars, + penalty_factor), objective_function_(std::move(objective_function)) { - if (!secondary_vars.empty()) { - assignment_.Add(secondary_vars); - } + this->AddVars(secondary_vars); } -IntExpr* TernaryGuidedLocalSearch::MakeElementPenalty(int index) { - return solver()->MakeElement( - [this, index](int64_t i, int64_t j) { - return PenalizedValue(index, i, j); +template +IntExpr* TernaryGuidedLocalSearch

::MakeElementPenalty(int index) { + return this->solver()->MakeElement( + [this, index](int64_t j, int64_t k) { + return PenalizedValue(index, j, k); }, - vars_[index], secondary_vars_[index]); + this->GetVar(index), this->GetVar(this->NumPrimaryVars() + index)); } -int64_t TernaryGuidedLocalSearch::AssignmentElementPenalty( - const Assignment& assignment, int index) { - return PenalizedValue(index, assignment.Value(vars_[index]), - assignment.Value(secondary_vars_[index])); +template +int64_t TernaryGuidedLocalSearch

::AssignmentElementPenalty(int index) const { + return PenalizedValue(index, this->GetValue(index), + this->GetValue(this->NumPrimaryVars() + index)); } -int64_t TernaryGuidedLocalSearch::AssignmentPenalty( - const Assignment& assignment, int index, int64_t next) { - return objective_function_(index, next, - assignment.Value(secondary_vars_[index])); +template +int64_t TernaryGuidedLocalSearch

::AssignmentPenalty(int64_t var, + int64_t value) const { + return objective_function_(var, value, + this->GetValue(this->NumPrimaryVars() + var)); } -bool TernaryGuidedLocalSearch::EvaluateElementValue( - const Assignment::IntContainer& container, int64_t index, - int* container_index, int64_t* penalty) { - const IntVarElement& element = container.Element(*container_index); - if (element.Activated()) { - *penalty = PenalizedValue( - index, element.Value(), - GetAssignmentSecondaryValue(container, index, container_index)); - return true; - } - return false; -} - -// Penalized value for (i, j) = penalty_factor_ * penalty(i, j) * cost (i, j) -int64_t TernaryGuidedLocalSearch::PenalizedValue(int64_t i, int64_t j, - int64_t k) { - const Arc arc(i, j); - const int64_t penalty = penalties_->Value(arc); - if (penalty != 0) { // objective_function_(i, j, k) can be costly - const double penalized_value_fp = - penalty_factor_ * penalty * objective_function_(i, j, k); - const int64_t penalized_value = - (penalized_value_fp < std::numeric_limits::max()) - ? static_cast(penalized_value_fp) - : std::numeric_limits::max(); - if (maximize_) { - return -penalized_value; - } else { - return penalized_value; +template +int64_t TernaryGuidedLocalSearch

::Evaluate(const Assignment* delta, + int64_t current_penalty, + bool incremental) { + int64_t penalty = current_penalty; + const Assignment::IntContainer& container = delta->IntVarContainer(); + // Collect values for each secondary variable, matching them with their + // corresponding primary variable. + std::vector secondary_values(this->NumPrimaryVars(), -1); + for (int i = 0; i < container.Size(); ++i) { + const IntVarElement& new_element = container.Element(i); + if (!new_element.Activated()) continue; + const int index = this->GetLocalIndexFromVar(new_element.Var()); + if (index != -1 && index >= this->NumPrimaryVars()) { // secondary variable + secondary_values[index - this->NumPrimaryVars()] = new_element.Value(); } - } else { - return 0; } + for (int i = 0; i < container.Size(); ++i) { + const IntVarElement& new_element = container.Element(i); + const int index = this->GetLocalIndexFromVar(new_element.Var()); + // Only process primary variables. + if (index == -1 || index >= this->NumPrimaryVars()) { + continue; + } + penalty = CapSub(penalty, this->penalized_values_.Get(index)); + // Performed and active. + if (new_element.Activated() && secondary_values[index] != -1) { + const int64_t new_penalty = + PenalizedValue(index, new_element.Value(), secondary_values[index]); + penalty = CapAdd(penalty, new_penalty); + if (incremental) { + this->penalized_values_.Set(index, new_penalty); + } + } + } + return penalty; } -int64_t TernaryGuidedLocalSearch::GetAssignmentSecondaryValue( - const Assignment::IntContainer& container, int index, - int* container_index) const { - const IntVar* secondary_var = secondary_vars_[index]; - int hint_index = *container_index + 1; - if (hint_index > 0 && hint_index < container.Size() && - secondary_var == container.Element(hint_index).Var()) { - *container_index = hint_index; - return container.Element(hint_index).Value(); - } else { - return container.Element(secondary_var).Value(); - } +// Penalized value for (i, j) = penalty_factor_ * penalty(i, j) * cost (i, j, k) +template +int64_t TernaryGuidedLocalSearch

::PenalizedValue(int64_t i, int64_t j, + int64_t k) const { + const int64_t penalty = this->penalties_.GetPenalty({i, j}); + // Calls to objective_function_(i, j, k) can be costly. + if (penalty == 0) return 0; + const double penalized_value_fp = + this->penalty_factor_ * penalty * objective_function_(i, j, k); + const int64_t penalized_value = + (penalized_value_fp < std::numeric_limits::max()) + ? static_cast(penalized_value_fp) + : std::numeric_limits::max(); + return this->maximize_ ? -penalized_value : penalized_value; } } // namespace @@ -3978,9 +4067,16 @@ SearchMonitor* Solver::MakeGuidedLocalSearch( bool maximize, IntVar* const objective, Solver::IndexEvaluator2 objective_function, int64_t step, const std::vector& vars, double penalty_factor) { - return RevAlloc(new BinaryGuidedLocalSearch( - this, objective, std::move(objective_function), maximize, step, vars, - penalty_factor)); + if (absl::GetFlag(FLAGS_cp_use_sparse_gls_penalties)) { + return RevAlloc(new BinaryGuidedLocalSearch( + this, objective, std::move(objective_function), maximize, step, vars, + penalty_factor)); + } else { + return RevAlloc( + new BinaryGuidedLocalSearch( + this, objective, std::move(objective_function), maximize, step, + vars, penalty_factor)); + } } SearchMonitor* Solver::MakeGuidedLocalSearch( @@ -3988,9 +4084,16 @@ SearchMonitor* Solver::MakeGuidedLocalSearch( Solver::IndexEvaluator3 objective_function, int64_t step, const std::vector& vars, const std::vector& secondary_vars, double penalty_factor) { - return RevAlloc(new TernaryGuidedLocalSearch( - this, objective, std::move(objective_function), maximize, step, vars, - secondary_vars, penalty_factor)); + if (absl::GetFlag(FLAGS_cp_use_sparse_gls_penalties)) { + return RevAlloc(new TernaryGuidedLocalSearch( + this, objective, std::move(objective_function), maximize, step, vars, + secondary_vars, penalty_factor)); + } else { + return RevAlloc( + new TernaryGuidedLocalSearch( + this, objective, std::move(objective_function), maximize, step, + vars, secondary_vars, penalty_factor)); + } } // ---------- Search Limits ---------- @@ -3999,6 +4102,13 @@ SearchMonitor* Solver::MakeGuidedLocalSearch( SearchLimit::~SearchLimit() {} +void SearchLimit::Install() { + ListenToEvent(Solver::MonitorEvent::kEnterSearch); + ListenToEvent(Solver::MonitorEvent::kBeginNextDecision); + ListenToEvent(Solver::MonitorEvent::kPeriodicCheck); + ListenToEvent(Solver::MonitorEvent::kRefuteDecision); +} + void SearchLimit::EnterSearch() { crossed_ = false; Init(); @@ -4050,6 +4160,14 @@ RegularLimit::RegularLimit(Solver* const s, absl::Duration time, RegularLimit::~RegularLimit() {} +void RegularLimit::Install() { + SearchLimit::Install(); + ListenToEvent(Solver::MonitorEvent::kExitSearch); + ListenToEvent(Solver::MonitorEvent::kIsUncheckedSolutionLimitReached); + ListenToEvent(Solver::MonitorEvent::kProgressPercent); + ListenToEvent(Solver::MonitorEvent::kAccept); +} + void RegularLimit::Copy(const SearchLimit* const limit) { const RegularLimit* const regular = reinterpret_cast(limit); @@ -4252,6 +4370,11 @@ ImprovementSearchLimit::ImprovementSearchLimit( ImprovementSearchLimit::~ImprovementSearchLimit() {} +void ImprovementSearchLimit::Install() { + SearchLimit::Install(); + ListenToEvent(Solver::MonitorEvent::kAtSolution); +} + void ImprovementSearchLimit::Init() { best_objective_ = maximize_ ? -std::numeric_limits::infinity() : std::numeric_limits::infinity(); @@ -4721,6 +4844,8 @@ class LubyRestart : public SearchMonitor { } } + void Install() override { ListenToEvent(Solver::MonitorEvent::kBeginFail); } + std::string DebugString() const override { return absl::StrFormat("LubyRestart(%i)", scale_factor_); } @@ -4756,6 +4881,8 @@ class ConstantRestart : public SearchMonitor { } } + void Install() override { ListenToEvent(Solver::MonitorEvent::kBeginFail); } + std::string DebugString() const override { return absl::StrFormat("ConstantRestart(%i)", frequency_); }