[CP-SAT] improve memory usage for LNS

lperron · lperron · commit 113c4fa8a1d0 · 2025-02-10T14:48:48.000+01:00
diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc
@@ -19,6 +19,7 @@
 #include <deque>
 #include <functional>
 #include <limits>
+#include <memory>
 #include <random>
 #include <string>
 #include <tuple>
@@ -77,16 +78,19 @@ NeighborhoodGeneratorHelper::NeighborhoodGeneratorHelper(
       shared_bounds_(shared_bounds),
       shared_response_(shared_response) {
   // Initialize proto memory.
+  local_arena_storage_.assign(Neighborhood::kDefaultArenaSizePerVariable *
+                                  model_proto_.variables_size(),
+                              0);
+  local_arena_ = std::make_unique<google::protobuf::Arena>(
+      local_arena_storage_.data(), local_arena_storage_.size());
   simplified_model_proto_ =
-      google::protobuf::Arena::Create<CpModelProto>(&local_arena_);
-  model_proto_with_only_variables_ =
-      google::protobuf::Arena::Create<CpModelProto>(&local_arena_);
+      google::protobuf::Arena::Create<CpModelProto>(local_arena_.get());
 
   CHECK(shared_response_ != nullptr);
   if (shared_bounds_ != nullptr) {
     shared_bounds_id_ = shared_bounds_->RegisterNewId();
   }
-  *model_proto_with_only_variables_->mutable_variables() =
+  *model_proto_with_only_variables_.mutable_variables() =
       model_proto_.variables();
   InitializeHelperData();
   RecomputeHelperData();
@@ -112,15 +116,15 @@ void NeighborhoodGeneratorHelper::Synchronize() {
         const int64_t new_ub = new_upper_bounds[i];
         if (VLOG_IS_ON(3)) {
           const auto& domain =
-              model_proto_with_only_variables_->variables(var).domain();
+              model_proto_with_only_variables_.variables(var).domain();
           const int64_t old_lb = domain.Get(0);
           const int64_t old_ub = domain.Get(domain.size() - 1);
           VLOG(3) << "Variable: " << var << " old domain: [" << old_lb << ", "
                   << old_ub << "] new domain: [" << new_lb << ", " << new_ub
                   << "]";
         }
         const Domain old_domain = ReadDomainFromProto(
-            model_proto_with_only_variables_->variables(var));
+            model_proto_with_only_variables_.variables(var));
         const Domain new_domain =
             old_domain.IntersectionWith(Domain(new_lb, new_ub));
         if (new_domain.IsEmpty()) {
@@ -141,7 +145,7 @@ void NeighborhoodGeneratorHelper::Synchronize() {
         }
         FillDomainInProto(
             new_domain,
-            model_proto_with_only_variables_->mutable_variables(var));
+            model_proto_with_only_variables_.mutable_variables(var));
         new_variables_have_been_fixed |= new_domain.IsFixed();
       }
     }
@@ -164,7 +168,7 @@ bool NeighborhoodGeneratorHelper::ObjectiveDomainIsConstraining() const {
     const int var = PositiveRef(model_proto_.objective().vars(i));
     const int64_t coeff = model_proto_.objective().coeffs(i);
     const auto& var_domain =
-        model_proto_with_only_variables_->variables(var).domain();
+        model_proto_with_only_variables_.variables(var).domain();
     const int64_t v1 = coeff * var_domain[0];
     const int64_t v2 = coeff * var_domain[var_domain.size() - 1];
     min_activity += std::min(v1, v2);
@@ -218,9 +222,20 @@ void NeighborhoodGeneratorHelper::RecomputeHelperData() {
   {
     Model local_model;
     CpModelProto mapping_proto;
+    // We want to replace the simplified_model_proto_ by a new one. Since
+    // deleting an object in the arena doesn't free the memory, we also delete
+    // and recreate the arena, but reusing the same storage.
+    int64_t new_size = local_arena_->SpaceUsed();
+    new_size += new_size / 2;
     simplified_model_proto_->Clear();
+    local_arena_.reset();
+    local_arena_storage_.resize(new_size);
+    local_arena_ = std::make_unique<google::protobuf::Arena>(
+        local_arena_storage_.data(), local_arena_storage_.size());
+    simplified_model_proto_ =
+        google::protobuf::Arena::Create<CpModelProto>(local_arena_.get());
     *simplified_model_proto_->mutable_variables() =
-        model_proto_with_only_variables_->variables();
+        model_proto_with_only_variables_.variables();
     PresolveContext context(&local_model, simplified_model_proto_,
                             &mapping_proto);
     ModelCopy copier(&context);
@@ -383,7 +398,7 @@ bool NeighborhoodGeneratorHelper::IsActive(int var) const {
 }
 
 bool NeighborhoodGeneratorHelper::IsConstant(int var) const {
-  const auto& var_proto = model_proto_with_only_variables_->variables(var);
+  const auto& var_proto = model_proto_with_only_variables_.variables(var);
   return var_proto.domain_size() == 2 &&
          var_proto.domain(0) == var_proto.domain(1);
 }
@@ -395,7 +410,7 @@ Neighborhood NeighborhoodGeneratorHelper::FullNeighborhood() const {
   {
     absl::ReaderMutexLock lock(&domain_mutex_);
     *neighborhood.delta.mutable_variables() =
-        model_proto_with_only_variables_->variables();
+        model_proto_with_only_variables_.variables();
   }
   return neighborhood;
 }
@@ -1057,7 +1072,7 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables(
     absl::ReaderMutexLock domain_lock(&domain_mutex_);
     for (int i = 0; i < num_variables; ++i) {
       const IntegerVariableProto& current_var =
-          model_proto_with_only_variables_->variables(i);
+          model_proto_with_only_variables_.variables(i);
       IntegerVariableProto* new_var = neighborhood.delta.add_variables();
 
       // We only copy the name in debug mode.
@@ -1203,7 +1218,7 @@ CpModelProto NeighborhoodGeneratorHelper::UpdatedModelProtoCopy() const {
   {
     absl::MutexLock domain_lock(&domain_mutex_);
     *updated_model.mutable_variables() =
-        model_proto_with_only_variables_->variables();
+        model_proto_with_only_variables_.variables();
   }
   return updated_model;
 }
diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h
@@ -319,14 +319,15 @@ class NeighborhoodGeneratorHelper : public SubSolver {
   // Arena holding the memory of the CpModelProto* of this class. This saves the
   // destruction cost that can take time on problem with millions of
   // variables/constraints.
-  google::protobuf::Arena local_arena_;
+  std::vector<char> local_arena_storage_;
+  std::unique_ptr<google::protobuf::Arena> local_arena_;
 
   // This proto will only contain the field variables() with an updated version
   // of the domains compared to model_proto_.variables(). We do it like this to
   // reduce the memory footprint of the helper when the model is large.
   //
   // TODO(user): Use custom domain repository rather than a proto?
-  CpModelProto* model_proto_with_only_variables_ ABSL_GUARDED_BY(domain_mutex_);
+  CpModelProto model_proto_with_only_variables_ ABSL_GUARDED_BY(domain_mutex_);
 
   // Constraints by types. This never changes.
   std::vector<std::vector<int>> type_to_constraints_;
diff --git a/ortools/sat/diffn_cuts.cc b/ortools/sat/diffn_cuts.cc
@@ -391,6 +391,28 @@ std::string DiffnCtEvent::DebugString() const {
 // The original cut is:
 //    sum(end_min_i * duration_min_i) >=
 //        (sum(duration_min_i^2) + sum(duration_min_i)^2) / 2
+//
+// Let's build a figure where each horizontal rectangle represent a task. It
+// ends at the end of the task, and its height is the duration of the task.
+// For a given order, we pack each rectangle to the left while not overlapping,
+// that is one rectangle starts when the previous one ends.
+//
+//     e1
+// -----
+// :\  | s1
+// :  \|       e2
+// -------------
+//     :\      |
+//     :   \   | s2
+//     :      \|  e3
+// ----------------
+//             : \| s3
+// ----------------
+//
+// We can notice that the total area is independent of the order of tasks.
+// The first term of the rhs is the area above the diagonal.
+// The second term of the rhs is the area below the diagonal.
+//
 // We apply the following changes (see the code for cumulative constraints):
 //   - we strengthen this cuts by noticing that if all tasks starts after S,
 //     then replacing end_min_i by (end_min_i - S) is still valid.
@@ -461,12 +483,11 @@ void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
     // Best cut so far for this loop.
     int best_end = -1;
     double best_efficacy = 0.01;
-    IntegerValue best_min_contrib = 0;
-    IntegerValue best_capacity = 0;
-    IntegerValue best_y_range = 0;
+    IntegerValue best_min_rhs = 0;
+    bool best_use_subset_sum = false;
 
     // Used in the first term of the rhs of the equation.
-    IntegerValue sum_event_contributions = 0;
+    IntegerValue sum_event_areas = 0;
     // Used in the second term of the rhs of the equation.
     IntegerValue sum_energy = 0;
     // For normalization.
@@ -486,8 +507,7 @@ void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
       DCHECK_GE(event.x_start_min, sequence_start_min);
       // Make sure we do not overflow.
       if (!AddTo(event.energy_min, &sum_energy)) break;
-      if (!AddProductTo(event.energy_min, event.x_size_min,
-                        &sum_event_contributions)) {
+      if (!AddProductTo(event.energy_min, event.x_size_min, &sum_event_areas)) {
         break;
       }
       if (!AddSquareTo(event.energy_min, &sum_square_energy)) break;
@@ -505,7 +525,8 @@ void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
         if (i == 0) {
           dp.Reset((y_max_of_subset - y_min_of_subset).value());
         } else {
-          if (y_max_of_subset - y_min_of_subset != dp.Bound()) {
+          // TODO(user): Can we increase the bound dynamically ?
+          if (y_max_of_subset - y_min_of_subset > dp.Bound()) {
             use_dp = false;
           }
         }
@@ -522,23 +543,21 @@ void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
       if (sum_of_y_size_min <= reachable_capacity) continue;
 
       // Do we have a violated cut ?
-      const IntegerValue large_rectangle_contrib =
-          CapProdI(sum_energy, sum_energy);
-      if (AtMinOrMaxInt64I(large_rectangle_contrib)) break;
-      const IntegerValue mean_large_rectangle_contrib =
-          CeilRatio(large_rectangle_contrib, reachable_capacity);
+      const IntegerValue square_sum_energy = CapProdI(sum_energy, sum_energy);
+      if (AtMinOrMaxInt64I(square_sum_energy)) break;
+      const IntegerValue rhs_second_term =
+          CeilRatio(square_sum_energy, reachable_capacity);
 
-      IntegerValue min_contrib =
-          CapAddI(sum_event_contributions, mean_large_rectangle_contrib);
-      if (AtMinOrMaxInt64I(min_contrib)) break;
-      min_contrib = CeilRatio(min_contrib, 2);
+      IntegerValue min_rhs = CapAddI(sum_event_areas, rhs_second_term);
+      if (AtMinOrMaxInt64I(min_rhs)) break;
+      min_rhs = CeilRatio(min_rhs, 2);
 
       // shift contribution by current_start_min.
-      if (!AddProductTo(sum_energy, current_start_min, &min_contrib)) break;
+      if (!AddProductTo(sum_energy, current_start_min, &min_rhs)) break;
 
       // The efficacy of the cut is the normalized violation of the above
       // equation. We will normalize by the sqrt of the sum of squared energies.
-      const double efficacy = (ToDouble(min_contrib) - lp_contrib) /
+      const double efficacy = (ToDouble(min_rhs) - lp_contrib) /
                               std::sqrt(ToDouble(sum_square_energy));
 
       // For a given start time, we only keep the best cut.
@@ -550,13 +569,13 @@ void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
       if (efficacy > best_efficacy) {
         best_efficacy = efficacy;
         best_end = i;
-        best_min_contrib = min_contrib;
-        best_capacity = reachable_capacity;
-        best_y_range = y_max_of_subset - y_min_of_subset;
+        best_min_rhs = min_rhs;
+        best_use_subset_sum =
+            reachable_capacity < y_max_of_subset - y_min_of_subset;
       }
     }
     if (best_end != -1) {
-      LinearConstraintBuilder cut(model, best_min_contrib, kMaxIntegerValue);
+      LinearConstraintBuilder cut(model, best_min_rhs, kMaxIntegerValue);
       bool is_lifted = false;
       bool add_energy_to_name = false;
       for (int i = 0; i <= best_end; ++i) {
@@ -568,9 +587,7 @@ void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
       std::string full_name(cut_name);
       if (is_lifted) full_name.append("_lifted");
       if (add_energy_to_name) full_name.append("_energy");
-      if (best_capacity < best_y_range) {
-        full_name.append("_subsetsum");
-      }
+      if (best_use_subset_sum) full_name.append("_subsetsum");
       top_n_cuts.AddCut(cut.Build(), full_name, manager->LpValues());
     }
   }
diff --git a/ortools/util/fixed_shape_binary_tree.h b/ortools/util/fixed_shape_binary_tree.h
@@ -218,6 +218,7 @@ class FixedShapeBinaryTree {
   template <typename TypeWithPushBack>
   void PartitionIntervalIntoNodes(LeafIndex first_leaf, LeafIndex last_leaf,
                                   TypeWithPushBack* result) const {
+    DCHECK_LE(first_leaf, last_leaf);
     TreeNodeIndex prev(0);
     TreeNodeIndex current = GetNodeStartOfRange(first_leaf, last_leaf);
     if (current == Root()) {
