diff --git a/ortools/sat/all_different.cc b/ortools/sat/all_different.cc index 3f861e1b73..6b47c5d96f 100644 --- a/ortools/sat/all_different.cc +++ b/ortools/sat/all_different.cc @@ -74,8 +74,9 @@ std::function AllDifferentBinary( } std::function AllDifferentOnBounds( - const std::vector& expressions) { - return [=](Model* model) { + absl::Span expressions) { + return [=, expressions = std::vector( + expressions.begin(), expressions.end())](Model* model) { if (expressions.empty()) return; auto* constraint = new AllDifferentBoundsPropagator( expressions, model->GetOrCreate()); diff --git a/ortools/sat/all_different.h b/ortools/sat/all_different.h index 0d833e7f1c..4b202720f9 100644 --- a/ortools/sat/all_different.h +++ b/ortools/sat/all_different.h @@ -49,7 +49,7 @@ std::function AllDifferentBinary( std::function AllDifferentOnBounds( absl::Span vars); std::function AllDifferentOnBounds( - const std::vector& expressions); + absl::Span expressions); // This constraint forces all variables to take different values. This is meant // to be used as a complement to an alldifferent decomposition like diff --git a/ortools/sat/cp_model_solver_helpers.cc b/ortools/sat/cp_model_solver_helpers.cc index 9e013c5d33..73e6933cdd 100644 --- a/ortools/sat/cp_model_solver_helpers.cc +++ b/ortools/sat/cp_model_solver_helpers.cc @@ -430,7 +430,7 @@ void InitializeLinearConstraintSymmetrizerIfRequested( orbit_is_ok.resize(orbit_index + 1, true); } - // In linerization level >=2, all variables should have a view. + // In linearization level >=2, all variables should have a view. // Otherwise revisit and skip orbit without a full view. const IntegerVariable var = mapping->Integer(proto_var); CHECK_NE(var, kNoIntegerVariable); @@ -1191,9 +1191,13 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { LoadBooleanSymmetries(model_proto, model); } + TimeLimit* time_limit = model->GetOrCreate(); + if (time_limit->LimitReached()) return; + ExtractEncoding(model_proto, model); PropagateEncodingFromEquivalenceRelations(model_proto, model); + if (time_limit->LimitReached()) return; // Check the model is still feasible before continuing. if (sat_solver->ModelIsUnsat()) return unsat(); @@ -1207,8 +1211,12 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { FillBinaryRelationRepository(model_proto, model); + if (time_limit->LimitReached()) return; + // Load the constraints. int num_ignored_constraints = 0; + int time_limit_check_count = 0; + static const int kTimeLimitCheckFrequency = 1000; absl::flat_hash_set unsupported_types; for (const ConstraintProto& ct : model_proto.constraints()) { if (mapping->ConstraintIsAlreadyLoaded(&ct)) { @@ -1221,6 +1229,12 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { continue; } + ++time_limit_check_count; + if (time_limit_check_count >= kTimeLimitCheckFrequency) { + if (time_limit->LimitReached()) return; + time_limit_check_count = 0; + } + // We propagate after each new Boolean constraint but not the integer // ones. So we call FinishPropagation() manually here. // @@ -1277,6 +1291,8 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { void LoadFeasibilityPump(const CpModelProto& model_proto, Model* model) { LoadBaseModel(model_proto, model); + if (model->GetOrCreate()->LimitReached()) return; + auto* mapping = model->GetOrCreate(); const SatParameters& parameters = *(model->GetOrCreate()); if (parameters.linearization_level() == 0) return; @@ -1309,6 +1325,8 @@ void LoadFeasibilityPump(const CpModelProto& model_proto, Model* model) { void LoadCpModel(const CpModelProto& model_proto, Model* model) { LoadBaseModel(model_proto, model); + if (model->GetOrCreate()->LimitReached()) return; + // We want to load the debug solution before the initial propag. // But at this point the objective is not loaded yet, so we will not have // a value for the objective integer variable, so we do it again later. @@ -1571,6 +1589,7 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, Model* model) { auto* shared_response_manager = model->GetOrCreate(); if (shared_response_manager->ProblemIsSolved()) return; + if (model->GetOrCreate()->LimitReached()) return; const SatParameters& parameters = *model->GetOrCreate(); if (parameters.stop_after_root_propagation()) return; @@ -1702,6 +1721,8 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, Model* model) { void QuickSolveWithHint(const CpModelProto& model_proto, Model* model) { if (!model_proto.has_solution_hint()) return; + if (model->GetOrCreate()->LimitReached()) return; + auto* shared_response_manager = model->GetOrCreate(); if (shared_response_manager->ProblemIsSolved()) return; diff --git a/ortools/sat/python/linear_expr.cc b/ortools/sat/python/linear_expr.cc index 66d2214ae2..e9c58c3aae 100644 --- a/ortools/sat/python/linear_expr.cc +++ b/ortools/sat/python/linear_expr.cc @@ -25,6 +25,7 @@ #include "absl/log/check.h" #include "absl/strings/str_cat.h" #include "absl/strings/str_join.h" +#include "absl/types/span.h" #include "ortools/util/fp_roundtrip_conv.h" #include "ortools/util/sorted_interval_list.h" @@ -403,8 +404,8 @@ std::string SumArray::DebugString() const { } FloatWeightedSum::FloatWeightedSum( - const std::vector>& exprs, - const std::vector& coeffs, double offset) + absl::Span> exprs, + absl::Span coeffs, double offset) : exprs_(exprs.begin(), exprs.end()), coeffs_(coeffs.begin(), coeffs.end()), offset_(offset) { @@ -480,8 +481,8 @@ std::string FloatWeightedSum::DebugString() const { } IntWeightedSum::IntWeightedSum( - const std::vector>& exprs, - const std::vector& coeffs, int64_t offset) + absl::Span> exprs, + absl::Span coeffs, int64_t offset) : exprs_(exprs.begin(), exprs.end()), coeffs_(coeffs.begin(), coeffs.end()), offset_(offset) {} diff --git a/ortools/sat/python/linear_expr.h b/ortools/sat/python/linear_expr.h index 1f3b07eed0..2f021d9114 100644 --- a/ortools/sat/python/linear_expr.h +++ b/ortools/sat/python/linear_expr.h @@ -24,6 +24,7 @@ #include "absl/container/fixed_array.h" #include "absl/log/check.h" #include "absl/strings/str_cat.h" +#include "absl/types/span.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/util/sorted_interval_list.h" @@ -299,8 +300,8 @@ class SumArray : public LinearExpr { /// A class to hold a weighted sum of floating point linear expressions. class FloatWeightedSum : public LinearExpr { public: - FloatWeightedSum(const std::vector>& exprs, - const std::vector& coeffs, double offset); + FloatWeightedSum(absl::Span> exprs, + absl::Span coeffs, double offset); ~FloatWeightedSum() override = default; void VisitAsFloat(FloatExprVisitor& lin, double c) override; @@ -320,8 +321,8 @@ class FloatWeightedSum : public LinearExpr { /// A class to hold a weighted sum of integer linear expressions. class IntWeightedSum : public LinearExpr { public: - IntWeightedSum(const std::vector>& exprs, - const std::vector& coeffs, int64_t offset); + IntWeightedSum(absl::Span> exprs, + absl::Span coeffs, int64_t offset); ~IntWeightedSum() override = default; void VisitAsFloat(FloatExprVisitor& lin, double c) override; diff --git a/ortools/sat/routing_cuts.cc b/ortools/sat/routing_cuts.cc index dc98c6928a..e98c62d7fc 100644 --- a/ortools/sat/routing_cuts.cc +++ b/ortools/sat/routing_cuts.cc @@ -146,6 +146,8 @@ MinOutgoingFlowHelper::~MinOutgoingFlowHelper() { stats.push_back( {"RoutingDp/num_full_dp_work_abort", num_full_dp_work_abort_}); stats.push_back({"RoutingDp/num_full_dp_rc_skip", num_full_dp_rc_skip_}); + stats.push_back( + {"RoutingDp/num_full_dp_unique_arc", num_full_dp_unique_arc_}); for (const auto& [name, count] : num_by_type_) { stats.push_back({absl::StrCat("RoutingDp/num_bounds_", name), count}); } @@ -170,6 +172,8 @@ void MinOutgoingFlowHelper::PrecomputeDataForSubset( has_incoming_arcs_from_outside_.assign(num_nodes, false); has_outgoing_arcs_to_outside_.assign(num_nodes, false); + incoming_arcs_from_outside_.assign(num_nodes, {}); + outgoing_arcs_to_outside_.assign(num_nodes, {}); for (auto& v : incoming_arc_indices_) v.clear(); for (auto& v : outgoing_arc_indices_) v.clear(); @@ -186,8 +190,10 @@ void MinOutgoingFlowHelper::PrecomputeDataForSubset( outgoing_arc_indices_[tail].push_back(i); incoming_arc_indices_[head].push_back(i); } else if (tail_in && !head_in) { + outgoing_arcs_to_outside_[tail].Add(i); has_outgoing_arcs_to_outside_[tail] = true; } else if (!tail_in && head_in) { + incoming_arcs_from_outside_[head].Add(i); has_incoming_arcs_from_outside_[head] = true; } } @@ -261,6 +267,42 @@ MinOutgoingFlowHelper::RelaxIntoSpecialBinPackingProblem( std::vector max_upper_bound_of_others(num_nodes, kMinIntegerValue); + // Lower bound on the node expression when it is first to serve the subset. + const auto lb_when_first = [&helper, this](int n, int dimension, + bool negate_expressions) { + AffineExpression expr = helper.GetNodeExpression(n, dimension); + if (negate_expressions) expr = expr.Negated(); + IntegerValue lb = integer_trail_.LevelZeroLowerBound(expr); + if (incoming_arcs_from_outside_[n].IsUnique()) { + const int a = incoming_arcs_from_outside_[n].Get(); + AffineExpression tail_expr = + helper.GetNodeExpression(tails_[a], dimension); + if (negate_expressions) tail_expr = tail_expr.Negated(); + const IntegerValue offset = + helper.GetArcOffsetLowerBound(a, dimension, negate_expressions); + lb = std::max(lb, integer_trail_.LevelZeroLowerBound(tail_expr) + offset); + } + return lb; + }; + + // Upper bound on the node expression when it is last to serve the subset. + const auto ub_when_last = [&helper, this](int n, int dimension, + bool negate_expressions) { + AffineExpression expr = helper.GetNodeExpression(n, dimension); + if (negate_expressions) expr = expr.Negated(); + IntegerValue ub = integer_trail_.LevelZeroUpperBound(expr); + if (outgoing_arcs_to_outside_[n].IsUnique()) { + const int a = outgoing_arcs_to_outside_[n].Get(); + AffineExpression head_expr = + helper.GetNodeExpression(heads_[a], dimension); + if (negate_expressions) head_expr = head_expr.Negated(); + const IntegerValue offset = + helper.GetArcOffsetLowerBound(a, dimension, negate_expressions); + ub = std::min(ub, integer_trail_.LevelZeroUpperBound(head_expr) - offset); + } + return ub; + }; + // We do a forward and a backward pass in order to compute the min/max // of all the other nodes for each node. for (const bool forward_pass : {true, false}) { @@ -269,8 +311,6 @@ MinOutgoingFlowHelper::RelaxIntoSpecialBinPackingProblem( const int size = subset.size(); for (int i = 0; i < size; ++i) { const int n = forward_pass ? subset[i] : subset[size - 1 - i]; - AffineExpression expr = helper.GetNodeExpression(n, dimension); - if (negate_expressions) expr = expr.Negated(); // The local_min/max contains the min/max of all nodes strictly before 'n' // in the forward pass (resp. striclty after). So after two passes, @@ -283,20 +323,17 @@ MinOutgoingFlowHelper::RelaxIntoSpecialBinPackingProblem( // Update local_min/max. if (has_incoming_arcs_from_outside_[n]) { - local_min = - std::min(local_min, integer_trail_.LevelZeroLowerBound(expr)); + local_min = std::min(local_min, + lb_when_first(n, dimension, negate_expressions)); } if (has_outgoing_arcs_to_outside_[n]) { local_max = - std::max(local_max, integer_trail_.LevelZeroUpperBound(expr)); + std::max(local_max, ub_when_last(n, dimension, negate_expressions)); } } } for (const int n : subset) { - AffineExpression expr = helper.GetNodeExpression(n, dimension); - if (negate_expressions) expr = expr.Negated(); - const absl::Span arcs = use_incoming ? incoming_arc_indices_[n] : outgoing_arc_indices_[n]; IntegerValue demand = kMaxIntegerValue; @@ -313,16 +350,14 @@ MinOutgoingFlowHelper::RelaxIntoSpecialBinPackingProblem( // a non-negative capacity (see above). obj.capacity = 0; if (use_incoming) { - if (max_upper_bound_of_others[n] > - integer_trail_.LevelZeroLowerBound(expr)) { - obj.capacity = max_upper_bound_of_others[n] - - integer_trail_.LevelZeroLowerBound(expr); + const IntegerValue lb = lb_when_first(n, dimension, negate_expressions); + if (max_upper_bound_of_others[n] > lb) { + obj.capacity = max_upper_bound_of_others[n] - lb; } } else { - if (integer_trail_.LevelZeroUpperBound(expr) > - min_lower_bound_of_others[n]) { - obj.capacity = integer_trail_.LevelZeroUpperBound(expr) - - min_lower_bound_of_others[n]; + const IntegerValue ub = ub_when_last(n, dimension, negate_expressions); + if (ub > min_lower_bound_of_others[n]) { + obj.capacity = ub - min_lower_bound_of_others[n]; } } @@ -753,7 +788,8 @@ int MinOutgoingFlowHelper::ComputeTightMinOutgoingFlow( bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( int k, absl::Span subset, RouteRelationsHelper* helper, - LinearConstraintManager* manager, int special_node, bool use_outgoing) { + LinearConstraintManager* manager, int special_node, + bool use_forward_direction) { cannot_be_first_.clear(); cannot_be_last_.clear(); if (k >= subset.size()) return true; @@ -789,7 +825,7 @@ bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( adjacency.Add({}); const int node = subset[i]; if (helper != nullptr) { - CHECK(use_outgoing); + CHECK(use_forward_direction); for (int j = 0; j < subset.size(); ++j) { if (i == j) continue; if (helper->PathExists(node, subset[j])) { @@ -797,7 +833,7 @@ bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( } } } else { - if (use_outgoing) { + if (use_forward_direction) { for (const int arc : outgoing_arc_indices_[node]) { adjacency.AppendToLastVector({heads_[arc], literals_[arc]}); } @@ -842,10 +878,12 @@ bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( const int size = subset.size(); const uint32_t final_mask = (1 << size) - 1; - const auto& can_be_last_set = use_outgoing ? has_outgoing_arcs_to_outside_ - : has_incoming_arcs_from_outside_; - const auto& can_be_first_set = use_outgoing ? has_incoming_arcs_from_outside_ - : has_outgoing_arcs_to_outside_; + const auto& can_be_last_set = use_forward_direction + ? has_outgoing_arcs_to_outside_ + : has_incoming_arcs_from_outside_; + const auto& can_be_first_set = use_forward_direction + ? has_incoming_arcs_from_outside_ + : has_outgoing_arcs_to_outside_; uint32_t can_be_last_mask = 0; for (int i = 0; i < subset.size(); ++i) { @@ -862,10 +900,38 @@ bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( std::vector states; states.push_back(State()); + // If a state has an unique arc to enter/leave it, we can update it further. + const auto update_using_unique_arc = [manager, this](UniqueArc unique_arc, + State& state) { + if (!unique_arc.IsUnique()) return true; + if (manager == nullptr) return true; + + ++num_full_dp_unique_arc_; + const Literal unique_lit = literals_[unique_arc.Get()]; + state.sum_of_reduced_costs += manager->GetLiteralReducedCost(unique_lit); + if (state.sum_of_reduced_costs > manager->ReducedCostsGap()) { + ++num_full_dp_rc_skip_; + return false; + } + + absl::flat_hash_map copy = state.lbs; + return binary_relation_repository_.PropagateLocalBounds( + integer_trail_, unique_lit, copy, &state.lbs); + }; + // We always start with the first node in this case. if (special_node >= 0) { states.back().node_set = 1; states.back().last_nodes_set = 1; + + // If there is a unique way to enter (resp. leave) that subset, we can start + // with tighter bounds! + const UniqueArc unique_arc = use_forward_direction + ? incoming_arcs_from_outside_[subset[0]] + : outgoing_arcs_to_outside_[subset[0]]; + if (!update_using_unique_arc(unique_arc, states.back())) { + return false; // We can't even start! + } } while (!states.empty()) { @@ -893,11 +959,20 @@ bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( const uint32_t head_mask = (1 << i); to_state.node_set = from_state.node_set | head_mask; to_state.last_nodes_set = from_state.last_nodes_set | head_mask; + to_state.lbs = from_state.lbs; + to_state.sum_of_reduced_costs = from_state.sum_of_reduced_costs; + + const UniqueArc unique_arc = + use_forward_direction ? incoming_arcs_from_outside_[subset[i]] + : outgoing_arcs_to_outside_[subset[i]]; + if (!update_using_unique_arc(unique_arc, to_state)) { + continue; + } + if (to_state.node_set == final_mask) { ++num_full_dp_early_abort_; return true; // All served! } - states.push_back(std::move(to_state)); } continue; @@ -949,6 +1024,25 @@ bool MinOutgoingFlowHelper::SubsetMightBeServedWithKRoutes( // One of the last node has no arc to outside, this is not possible. if (to_state.last_nodes_set & ~can_be_last_mask) continue; + // Add the constraints implied by each last node that has an unique + // way to leave the subset. + bool infeasible = false; + int last_mask = to_state.last_nodes_set; + for (int j = 0; last_mask; j++, last_mask /= 2) { + if ((last_mask & 1) == 0) continue; + + const UniqueArc unique_arc = + use_forward_direction ? outgoing_arcs_to_outside_[subset[j]] + : incoming_arcs_from_outside_[subset[j]]; + if (!update_using_unique_arc(unique_arc, to_state)) { + infeasible = true; + break; + } + } + if (infeasible) { + continue; + } + ++num_full_dp_early_abort_; return true; } @@ -2535,14 +2629,14 @@ bool RoutingCutHelper::TrySubsetCut(int known_bound, std::string name, total_incoming_weight - nodes_incoming_weight_[n] < threshold && !min_outgoing_flow_helper_.SubsetMightBeServedWithKRoutes( best_bound.bound(), subset, nullptr, manager, /*special_node=*/n, - /*use_outgoing=*/true)) { + /*use_forward_direction=*/true)) { best_bound.Update(best_bound.bound(), "", {n}, {}); } if (nodes_outgoing_weight_[n] > 0.1 && total_outgoing_weight - nodes_outgoing_weight_[n] < threshold && !min_outgoing_flow_helper_.SubsetMightBeServedWithKRoutes( best_bound.bound(), subset, nullptr, manager, /*special_node=*/n, - /*use_outgoing=*/false)) { + /*use_forward_direction=*/false)) { best_bound.Update(best_bound.bound(), "", {}, {n}); } } diff --git a/ortools/sat/routing_cuts.h b/ortools/sat/routing_cuts.h index 39e85bd307..31b07dd781 100644 --- a/ortools/sat/routing_cuts.h +++ b/ortools/sat/routing_cuts.h @@ -449,8 +449,8 @@ class MinOutgoingFlowHelper { // can stop as soon as one solution is found. // // If special_node is non-negative, we will only look for routes that - // - Start at this special_node if use_outgoing = true - // - End at this special_node if use_outgoing = false + // - Start at this special_node if use_forward_direction = true + // - End at this special_node if use_forward_direction = false // // If the RouteRelationsHelper is non null, then we will use "shortest path" // bounds instead of recomputing them from the binary relation of the model. @@ -461,7 +461,7 @@ class MinOutgoingFlowHelper { int k, absl::Span subset, RouteRelationsHelper* helper = nullptr, LinearConstraintManager* manager = nullptr, int special_node = -1, - bool use_outgoing = true); + bool use_forward_direction = true); // Advanced. If non-empty, and one of the functions above proved that a subset // needs at least k vehicles to serve it, then these vector list the nodes @@ -579,6 +579,25 @@ class MinOutgoingFlowHelper { std::vector has_incoming_arcs_from_outside_; std::vector has_outgoing_arcs_to_outside_; + // If a subset has an unique arc arriving or leaving at a given node, we can + // derive tighter bounds. + class UniqueArc { + public: + bool IsUnique() const { return num_seen_ == 1; } + int Get() const { return arc_; } + + void Add(int arc) { + ++num_seen_; + arc_ = arc; + } + + private: + int num_seen_ = 0; + int arc_ = 0; + }; + std::vector incoming_arcs_from_outside_; + std::vector outgoing_arcs_to_outside_; + // For each node n, whether it can appear at the current and next position in // a feasible path. std::vector reachable_; @@ -607,6 +626,7 @@ class MinOutgoingFlowHelper { int64_t num_full_dp_early_abort_ = 0; int64_t num_full_dp_work_abort_ = 0; int64_t num_full_dp_rc_skip_ = 0; + int64_t num_full_dp_unique_arc_ = 0; absl::flat_hash_map num_by_type_; };