[CP-SAT] fix cumulative cuts

Author: lperron

ortools/sat/diffn_cuts.cc

@@ -314,8 +314,12 @@ CutGenerator CreateNoOverlap2dEnergyCutGenerator(
     NoOverlap2DConstraintHelper* helper, Model* model) {
   CutGenerator result;
   result.only_run_at_level_zero = true;
-  AddIntegerVariableFromIntervals(&helper->x_helper(), model, &result.vars);
-  AddIntegerVariableFromIntervals(&helper->y_helper(), model, &result.vars);
+  AddIntegerVariableFromIntervals(
+      &helper->x_helper(), model, &result.vars,
+      IntegerVariablesToAddMask::kSize | IntegerVariablesToAddMask::kPresence);
+  AddIntegerVariableFromIntervals(
+      &helper->y_helper(), model, &result.vars,
+      IntegerVariablesToAddMask::kSize | IntegerVariablesToAddMask::kPresence);
   gtl::STLSortAndRemoveDuplicates(&result.vars);
   ProductDecomposer* product_decomposer =
       model->GetOrCreate<ProductDecomposer>();
@@ -558,8 +562,12 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator(
     NoOverlap2DConstraintHelper* helper, Model* model) {
   CutGenerator result;
   result.only_run_at_level_zero = true;
-  AddIntegerVariableFromIntervals(&helper->x_helper(), model, &result.vars);
-  AddIntegerVariableFromIntervals(&helper->y_helper(), model, &result.vars);
+  AddIntegerVariableFromIntervals(
+      &helper->x_helper(), model, &result.vars,
+      IntegerVariablesToAddMask::kEnd | IntegerVariablesToAddMask::kSize);
+  AddIntegerVariableFromIntervals(
+      &helper->y_helper(), model, &result.vars,
+      IntegerVariablesToAddMask::kEnd | IntegerVariablesToAddMask::kSize);
   gtl::STLSortAndRemoveDuplicates(&result.vars);
 
   auto* product_decomposer = model->GetOrCreate<ProductDecomposer>();

ortools/sat/integer_base.h

@@ -140,6 +140,11 @@ inline bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue* result) {
   return true;
 }
 
+// Computes result += a * a, and return false iff there is an overflow.
+inline bool AddSquareTo(IntegerValue a, IntegerValue* result) {
+  return AddProductTo(a, a, result);
+}
+
 // Index of an IntegerVariable.
 //
 // Each time we create an IntegerVariable we also create its negation. This is

ortools/sat/integer_expr.h

@@ -456,9 +456,10 @@ inline std::function<void(Model*)> WeightedSumGreaterOrEqual(
 // Weighted sum == constant.
 template <typename VectorInt>
 inline std::function<void(Model*)> FixedWeightedSum(
-    const std::vector<IntegerVariable>& vars, const VectorInt& coefficients,
+    absl::Span<const IntegerVariable> vars, const VectorInt& coefficients,
     int64_t value) {
-  return [=](Model* model) {
+  return [=, vars = std::vector<IntegerVariable>(vars.begin(), vars.end())](
+             Model* model) {
     model->Add(WeightedSumGreaterOrEqual(vars, coefficients, value));
     model->Add(WeightedSumLowerOrEqual(vars, coefficients, value));
   };

ortools/sat/linear_relaxation.cc

@@ -1070,10 +1070,13 @@ void AppendNoOverlap2dRelaxation(const ConstraintProto& ct, Model* model,
 
   CutGenerator& result = relaxation->cut_generators.emplace_back();
   result.only_run_at_level_zero = true;
-  AddIntegerVariableFromIntervals(&no_overlap_helper->x_helper(), model,
-                                  &result.vars);
-  AddIntegerVariableFromIntervals(&no_overlap_helper->y_helper(), model,
-                                  &result.vars);
+
+  AddIntegerVariableFromIntervals(
+      &no_overlap_helper->x_helper(), model, &result.vars,
+      IntegerVariablesToAddMask::kSize | IntegerVariablesToAddMask::kPresence);
+  AddIntegerVariableFromIntervals(
+      &no_overlap_helper->y_helper(), model, &result.vars,
+      IntegerVariablesToAddMask::kSize | IntegerVariablesToAddMask::kPresence);
   result.generate_cuts = [no_overlap_helper, model, product_decomposer,
                           last_level_zero_bound_change_idx = int64_t{-1}](
                              LinearConstraintManager* manager) mutable {

ortools/sat/routing_cuts.cc

@@ -47,6 +47,8 @@
 namespace operations_research {
 namespace sat {
 
+namespace {
+
 // Helper to compute lower bounds based on binary relations enforced by arc
 // literals in a RoutesConstraint.
 class LowerBoundsHelper {
@@ -171,6 +173,8 @@ class LowerBoundsHelper {
       lower_bound_by_var_and_arc_index_;
 };
 
+}  // namespace
+
 MinOutgoingFlowHelper::MinOutgoingFlowHelper(
     int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
     const std::vector<Literal>& literals, Model* model)
@@ -310,6 +314,7 @@ class OutgoingCutHelper {
     // Compute the total demands in order to know the minimum incoming/outgoing
     // flow.
     for (const int64_t demand : demands) total_demand_ += demand;
+    complement_of_subset_.reserve(num_nodes_);
   }
 
   // Try to add an outgoing cut from the given subset.
@@ -356,6 +361,7 @@ class OutgoingCutHelper {
 
   int64_t total_demand_ = 0;
   std::vector<bool> in_subset_;
+  std::vector<int> complement_of_subset_;
   MinOutgoingFlowHelper min_outgoing_flow_helper_;
 };
 
@@ -469,15 +475,23 @@ bool OutgoingCutHelper::TrySubsetCut(std::string name,
   // Bounds inferred automatically from the enforced binary relation of the
   // model.
   //
-  // TODO(user): use the complement of subset if contain_depot?
-  //
   // TODO(user): This is still not as good as the "capacity" bounds below in
   // some cases. Fix! we should be able to use the same relation to infer the
   // capacity bounds somehow.
-  if (subset.size() <= routing_cut_subset_size_for_binary_relation_bound_ &&
-      !contain_depot) {
-    const int64_t automatic_bound =
-        min_outgoing_flow_helper_.ComputeMinOutgoingFlow(subset);
+  if ((contain_depot ? num_nodes_ - subset.size() : subset.size()) <=
+      routing_cut_subset_size_for_binary_relation_bound_) {
+    int automatic_bound;
+    if (contain_depot) {
+      complement_of_subset_.clear();
+      for (int i = 0; i < num_nodes_; ++i) {
+        if (!in_subset_[i]) complement_of_subset_.push_back(i);
+      }
+      automatic_bound = min_outgoing_flow_helper_.ComputeMinOutgoingFlow(
+          complement_of_subset_);
+    } else {
+      automatic_bound =
+          min_outgoing_flow_helper_.ComputeMinOutgoingFlow(subset);
+    }
     if (automatic_bound > min_outgoing_flow) {
       absl::StrAppend(&name, "Automatic");
       min_outgoing_flow = automatic_bound;

ortools/sat/sat_parameters.proto

@@ -914,7 +914,9 @@ message SatParameters {
   // If the size of a subset of nodes of a RoutesConstraint is less than this
   // value, use linear constraints of size 1 and 2 (such as capacity and time
   // window constraints) enforced by the arc literals to compute cuts for this
-  // subset.
+  // subset. The algorithm for these cuts has a O(n^3) complexity, where n is
+  // the subset size. Hence the value of this parameter should not be too large
+  // (e.g. 10 or 20).
   optional int32 routing_cut_subset_size_for_binary_relation_bound = 312
       [default = 0];