[CP-SAT] more fuzzer bugs fixed; reorganize scheduling and packing code; speed up presolve on knapsack constraints

Author: lperron

ortools/sat/2d_orthogonal_packing.cc

@@ -26,6 +26,7 @@
 #include "absl/numeric/bits.h"
 #include "absl/random/distributions.h"
 #include "absl/types/span.h"
+#include "ortools/base/constant_divisor.h"
 #include "ortools/base/logging.h"
 #include "ortools/sat/2d_packing_brute_force.h"
 #include "ortools/sat/integer_base.h"
@@ -404,16 +405,16 @@ void OrthogonalPackingInfeasibilityDetector::GetAllCandidatesForKForDff2(
     candidates.Set(i);
   }
   for (int i = 1; i <= sqrt_bb_size; i++) {
-    const QuickSmallDivision div(i);
+    const ::util::math::ConstantDivisor<uint16_t> div(i);
     if (i > 1) {
-      candidates.Set(div.DivideByDivisor(bb_size.value()));
+      candidates.Set(bb_size.value() / div);
     }
     for (int k = 0; k < sizes.size(); k++) {
       IntegerValue size = sizes[k];
       if (2 * size > bb_size && size < bb_size) {
-        candidates.Set(div.DivideByDivisor(bb_size.value() - size.value() + 1));
+        candidates.Set((bb_size.value() - size.value() + 1) / div);
       } else if (2 * size < bb_size) {
-        candidates.Set(div.DivideByDivisor(size.value()));
+        candidates.Set(size.value() / div);
       }
     }
   }

ortools/sat/2d_orthogonal_packing.h

@@ -22,6 +22,7 @@
 #include "absl/log/check.h"
 #include "absl/random/bit_gen_ref.h"
 #include "absl/types/span.h"
+#include "ortools/base/constant_divisor.h"
 #include "ortools/sat/integer_base.h"
 #include "ortools/sat/synchronization.h"
 #include "ortools/util/bitset.h"
@@ -281,10 +282,7 @@ class RoundingDualFeasibleFunction {
  public:
   // `max_x` must fit in a uint16_t and `k` in [0, max_x/2].
   RoundingDualFeasibleFunction(IntegerValue max_x, IntegerValue k)
-      : div_(k.value()),
-        max_x_(max_x),
-        c_k_(div_.DivideByDivisor(max_x_.value())),
-        k_(k) {
+      : div_(k.value()), max_x_(max_x), c_k_(max_x_.value() / div_), k_(k) {
     DCHECK_GT(k, 0);
     DCHECK_LE(2 * k, max_x_);
     DCHECK_LE(max_x_, std::numeric_limits<uint16_t>::max());
@@ -296,11 +294,11 @@ class RoundingDualFeasibleFunction {
     DCHECK_LE(x, max_x_);
 
     if (2 * x > max_x_) {
-      return 2 * (c_k_ - div_.DivideByDivisor(max_x_.value() - x.value()));
+      return 2 * (c_k_ - (max_x_.value() - x.value()) / div_);
     } else if (2 * x == max_x_) {
       return c_k_;
     } else {
-      return 2 * div_.DivideByDivisor(x.value());
+      return 2 * (x.value() / div_);
     }
   }
 
@@ -309,7 +307,7 @@ class RoundingDualFeasibleFunction {
   IntegerValue LowestInverse(IntegerValue y) const;
 
  private:
-  const QuickSmallDivision div_;
+  const ::util::math::ConstantDivisor<uint16_t> div_;
   const IntegerValue max_x_;
   const IntegerValue c_k_;
   const IntegerValue k_;

ortools/sat/2d_orthogonal_packing_testing.cc

@@ -178,13 +178,12 @@ std::vector<RectangleInRange> MakeItemsFromRectangles(
   return ranges;
 }
 
-std::vector<ItemForPairwiseRestriction>
-GenerateItemsRectanglesWithNoPairwiseConflict(
+std::vector<ItemWithVariableSize> GenerateItemsRectanglesWithNoPairwiseConflict(
     absl::Span<const Rectangle> rectangles, double slack_factor,
     absl::BitGenRef random) {
   const std::vector<RectangleInRange> range_items =
       MakeItemsFromRectangles(rectangles, slack_factor, random);
-  std::vector<ItemForPairwiseRestriction> items;
+  std::vector<ItemWithVariableSize> items;
   items.reserve(rectangles.size());
   for (int i = 0; i < range_items.size(); ++i) {
     const RectangleInRange& rec = range_items[i];
@@ -201,13 +200,13 @@ GenerateItemsRectanglesWithNoPairwiseConflict(
   return items;
 }
 
-std::vector<ItemForPairwiseRestriction>
+std::vector<ItemWithVariableSize>
 GenerateItemsRectanglesWithNoPairwisePropagation(int num_rectangles,
                                                  double slack_factor,
                                                  absl::BitGenRef random) {
   const std::vector<Rectangle> rectangles =
       GenerateNonConflictingRectangles(num_rectangles, random);
-  std::vector<ItemForPairwiseRestriction> items =
+  std::vector<ItemWithVariableSize> items =
       GenerateItemsRectanglesWithNoPairwiseConflict(rectangles, slack_factor,
                                                     random);
   bool done = false;

ortools/sat/2d_orthogonal_packing_testing.h

@@ -38,12 +38,11 @@ std::vector<RectangleInRange> MakeItemsFromRectangles(
     absl::Span<const Rectangle> rectangles, double slack_factor,
     absl::BitGenRef random);
 
-std::vector<ItemForPairwiseRestriction>
-GenerateItemsRectanglesWithNoPairwiseConflict(
+std::vector<ItemWithVariableSize> GenerateItemsRectanglesWithNoPairwiseConflict(
     absl::Span<const Rectangle> rectangles, double slack_factor,
     absl::BitGenRef random);
 
-std::vector<ItemForPairwiseRestriction>
+std::vector<ItemWithVariableSize>
 GenerateItemsRectanglesWithNoPairwisePropagation(int num_rectangles,
                                                  double slack_factor,
                                                  absl::BitGenRef random);

ortools/sat/BUILD.bazel

@@ -1444,6 +1444,7 @@ cc_library(
         "//ortools/util:saturated_arithmetic",
         "//ortools/util:sorted_interval_list",
         "//ortools/util:strong_integers",
+        "@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/log:check",
@@ -2957,6 +2958,7 @@ cc_library(
         ":integer_base",
         ":synchronization",
         ":util",
+        "//ortools/base:constant_divisor",
         "//ortools/util:bitset",
         "@com_google_absl//absl/log",
         "@com_google_absl//absl/log:check",

ortools/sat/all_different.cc

@@ -83,8 +83,9 @@ std::function<void(Model*)> AllDifferentOnBounds(
 }
 
 std::function<void(Model*)> AllDifferentOnBounds(
-    const std::vector<IntegerVariable>& vars) {
-  return [=](Model* model) {
+    absl::Span<const IntegerVariable> vars) {
+  return [=, vars = std::vector<IntegerVariable>(vars.begin(), vars.end())](
+             Model* model) {
     if (vars.empty()) return;
     std::vector<AffineExpression> expressions;
     expressions.reserve(vars.size());

ortools/sat/all_different.h

@@ -45,7 +45,7 @@ std::function<void(Model*)> AllDifferentBinary(
 // this will not remove already taken values from inside a domain, but it will
 // propagates more the domain bounds.
 std::function<void(Model*)> AllDifferentOnBounds(
-    const std::vector<IntegerVariable>& vars);
+    absl::Span<const IntegerVariable> vars);
 std::function<void(Model*)> AllDifferentOnBounds(
     const std::vector<AffineExpression>& expressions);
 

ortools/sat/cp_model_presolve.cc

@@ -42,7 +42,9 @@
 #include "absl/status/statusor.h"
 #include "absl/strings/str_cat.h"
 #include "absl/types/span.h"
+#include "google/protobuf/arena.h"
 #include "google/protobuf/repeated_field.h"
+#include "google/protobuf/repeated_ptr_field.h"
 #include "google/protobuf/text_format.h"
 #include "ortools/base/logging.h"
 #include "ortools/base/mathutil.h"
@@ -5836,9 +5838,8 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
 
   // Filter absent boxes.
   int new_size = 0;
-  std::vector<Rectangle> bounding_boxes, fixed_boxes;
+  std::vector<Rectangle> bounding_boxes, fixed_boxes, non_fixed_bounding_boxes;
   std::vector<RectangleInRange> non_fixed_boxes;
-  std::vector<int> active_boxes;
   absl::flat_hash_set<int> fixed_item_indexes;
   for (int i = 0; i < proto.x_intervals_size(); ++i) {
     const int x_interval_index = proto.x_intervals(i);
@@ -5877,14 +5878,14 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
          IntegerValue(context_->EndMax(x_interval_index)),
          IntegerValue(context_->StartMin(y_interval_index)),
          IntegerValue(context_->EndMax(y_interval_index))});
-    active_boxes.push_back(new_size);
     if (context_->IntervalIsConstant(x_interval_index) &&
         context_->IntervalIsConstant(y_interval_index) &&
         context_->SizeMax(x_interval_index) > 0 &&
         context_->SizeMax(y_interval_index) > 0) {
       fixed_boxes.push_back(bounding_boxes.back());
       fixed_item_indexes.insert(new_size);
     } else {
+      non_fixed_bounding_boxes.push_back(bounding_boxes.back());
       non_fixed_boxes.push_back(
           {.box_index = new_size,
            .bounding_area = bounding_boxes.back(),
@@ -5909,15 +5910,25 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
     }
   }
 
-  const CompactVectorVector<int> components = GetOverlappingRectangleComponents(
-      bounding_boxes, absl::MakeSpan(active_boxes));
-  // The result of GetOverlappingRectangleComponents() omit singleton components
-  // thus to check whether a graph is fully connected we must check also the
-  // size of the unique component.
-  const bool is_fully_connected =
-      (components.size() == 1 && components[0].size() == active_boxes.size()) ||
-      (active_boxes.size() <= 1);
-  if (!is_fully_connected) {
+  if (new_size < initial_num_boxes) {
+    context_->UpdateRuleStats("no_overlap_2d: removed inactive boxes");
+    ct->mutable_no_overlap_2d()->mutable_x_intervals()->Truncate(new_size);
+    ct->mutable_no_overlap_2d()->mutable_y_intervals()->Truncate(new_size);
+  }
+
+  if (new_size == 0) {
+    context_->UpdateRuleStats("no_overlap_2d: no boxes");
+    return RemoveConstraint(ct);
+  }
+
+  if (new_size == 1) {
+    context_->UpdateRuleStats("no_overlap_2d: only one box");
+    return RemoveConstraint(ct);
+  }
+
+  const CompactVectorVector<int> components =
+      GetOverlappingRectangleComponents(bounding_boxes);
+  if (components.size() > 1) {
     for (int i = 0; i < components.size(); ++i) {
       absl::Span<const int> boxes = components[i];
       if (boxes.size() <= 1) continue;
@@ -5962,30 +5973,14 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
     return RemoveConstraint(ct);
   }
 
-  if (new_size < initial_num_boxes) {
-    context_->UpdateRuleStats("no_overlap_2d: removed inactive boxes");
-    ct->mutable_no_overlap_2d()->mutable_x_intervals()->Truncate(new_size);
-    ct->mutable_no_overlap_2d()->mutable_y_intervals()->Truncate(new_size);
-  }
-
-  if (new_size == 0) {
-    context_->UpdateRuleStats("no_overlap_2d: no boxes");
-    return RemoveConstraint(ct);
-  }
-
-  if (new_size == 1) {
-    context_->UpdateRuleStats("no_overlap_2d: only one box");
-    return RemoveConstraint(ct);
-  }
-
   // We check if the fixed boxes are not overlapping so downstream code can
   // assume it to be true.
   if (!FindPartialRectangleIntersections(fixed_boxes).empty()) {
     return context_->NotifyThatModelIsUnsat(
         "Two fixed boxes in no_overlap_2d overlap");
   }
 
-  if (fixed_boxes.size() == active_boxes.size()) {
+  if (non_fixed_bounding_boxes.empty()) {
     context_->UpdateRuleStats("no_overlap_2d: all boxes are fixed");
     return RemoveConstraint(ct);
   }
@@ -6035,6 +6030,36 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
       return RemoveConstraint(ct);
     }
   }
+  // If the non-fixed boxes are disjoint but connected by fixed boxes, we can
+  // split the constraint and duplicate the fixed boxes. To avoid duplicating
+  // too many fixed boxes, we do this after we we applied the presolve reducing
+  // their number to as few as possible.
+  const CompactVectorVector<int> non_fixed_components =
+      GetOverlappingRectangleComponents(non_fixed_bounding_boxes);
+  if (non_fixed_components.size() > 1) {
+    for (int i = 0; i < non_fixed_components.size(); ++i) {
+      // Note: we care about components of size 1 because they might be
+      // overlapping with the fixed boxes.
+      absl::Span<const int> indexes = non_fixed_components[i];
+
+      NoOverlap2DConstraintProto* new_no_overlap_2d =
+          context_->working_model->add_constraints()->mutable_no_overlap_2d();
+      for (const int idx : indexes) {
+        const int b = non_fixed_boxes[idx].box_index;
+        new_no_overlap_2d->add_x_intervals(proto.x_intervals(b));
+        new_no_overlap_2d->add_y_intervals(proto.y_intervals(b));
+      }
+      for (const int b : fixed_item_indexes) {
+        new_no_overlap_2d->add_x_intervals(proto.x_intervals(b));
+        new_no_overlap_2d->add_y_intervals(proto.y_intervals(b));
+      }
+    }
+    context_->UpdateNewConstraintsVariableUsage();
+    context_->UpdateRuleStats(
+        "no_overlap_2d: split into disjoint components duplicating fixed "
+        "boxes");
+    return RemoveConstraint(ct);
+  }
   RunPropagatorsForConstraint(*ct);
   return new_size < initial_num_boxes;
 }
@@ -9295,8 +9320,9 @@ void CpModelPresolver::DetectDuplicateColumns() {
     if (rep_to_dups[var].empty()) continue;
 
     // Since columns are the same, we can introduce a new variable = sum all
-    // columns. Note that we shouldn't have any overflow here by the
-    // precondition on our variable domains.
+    // columns. Note that the linear expression will not overflow, but the
+    // overflow check also requires that max_sum < int_max/2, which might
+    // happen.
     //
     // In the corner case where there is a lot of holes in the domain, and the
     // sum domain is too complex, we skip. Hopefully this should be rare.
@@ -9318,7 +9344,10 @@ void CpModelPresolver::DetectDuplicateColumns() {
     }
     const int new_var = context_->NewIntVarWithDefinition(
         domain, definition, /*append_constraint_to_mapping_model=*/true);
-    CHECK_NE(new_var, -1);
+    if (new_var == -1) {
+      context_->UpdateRuleStats("TODO duplicate: possible overflow");
+      continue;
+    }
 
     var_to_remove.push_back(var);
     CHECK_EQ(var_to_rep[var], -1);
@@ -14357,19 +14386,25 @@ void ApplyVariableMapping(absl::Span<int> mapping,
   }
 
   // Move the variable definitions.
-  std::vector<IntegerVariableProto> new_variables;
+  google::protobuf::RepeatedPtrField<IntegerVariableProto>
+      new_variables_storage;
+  google::protobuf::RepeatedPtrField<IntegerVariableProto>* new_variables;
+  if (proto->GetArena() == nullptr) {
+    new_variables = &new_variables_storage;
+  } else {
+    new_variables = google::protobuf::Arena::Create<
+        google::protobuf::RepeatedPtrField<IntegerVariableProto>>(
+        proto->GetArena());
+  }
   for (int i = 0; i < mapping.size(); ++i) {
     const int image = mapping[i];
     if (image < 0) continue;
-    if (image >= new_variables.size()) {
-      new_variables.resize(image + 1, IntegerVariableProto());
+    while (image >= new_variables->size()) {
+      new_variables->Add();
     }
-    new_variables[image].Swap(proto->mutable_variables(i));
-  }
-  proto->clear_variables();
-  for (IntegerVariableProto& proto_ref : new_variables) {
-    proto->add_variables()->Swap(&proto_ref);
+    (*new_variables)[image].Swap(proto->mutable_variables(i));
   }
+  proto->mutable_variables()->Swap(new_variables);
 
   // Check that all variables have a non-empty domain.
   for (const IntegerVariableProto& v : proto->variables()) {

ortools/sat/cuts.cc

@@ -2382,7 +2382,7 @@ IntegerValue SumOfAllDiffLowerBounder::SumOfMinDomainValues() {
   int count = 0;
   IntegerValue sum = 0;
   for (const IntegerValue value : min_values_) {
-    sum += value;
+    sum = CapAddI(sum, value);
     if (++count >= expr_mins_.size()) return sum;
   }
   return sum;
@@ -2439,6 +2439,8 @@ void TryToGenerateAllDiffCut(
     std::string max_suffix;
     const IntegerValue required_max_sum =
         -negated_diff_maxes.GetBestLowerBound(max_suffix);
+    if (required_max_sum == std::numeric_limits<IntegerValue>::max()) continue;
+    DCHECK_LE(required_min_sum, required_max_sum);
     if (sum < ToDouble(required_min_sum) - kMinCutViolation ||
         sum > ToDouble(required_max_sum) + kMinCutViolation) {
       LinearConstraintBuilder cut(model, required_min_sum, required_max_sum);

ortools/sat/cuts.h

@@ -727,6 +727,7 @@ class SumOfAllDiffLowerBounder {
   void Add(const AffineExpression& expr, int num_expr,
            const IntegerTrail& integer_trail);
 
+  // Return int_max if the sum overflows.
   IntegerValue SumOfMinDomainValues();
   IntegerValue SumOfDifferentMins();
   IntegerValue GetBestLowerBound(std::string& suffix);