diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 432e673304..726cef5184 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -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( + local_arena_storage_.data(), local_arena_storage_.size()); simplified_model_proto_ = - google::protobuf::Arena::Create(&local_arena_); - model_proto_with_only_variables_ = - google::protobuf::Arena::Create(&local_arena_); + google::protobuf::Arena::Create(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,7 +116,7 @@ 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 << ", " @@ -120,7 +124,7 @@ void NeighborhoodGeneratorHelper::Synchronize() { << "]"; } 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( + local_arena_storage_.data(), local_arena_storage_.size()); + simplified_model_proto_ = + google::protobuf::Arena::Create(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 index a92f63291d..ca0d365e5f 100644 --- 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 local_arena_storage_; + std::unique_ptr 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> type_to_constraints_; diff --git a/ortools/sat/diffn_cuts.cc b/ortools/sat/diffn_cuts.cc index f372ae3ef4..a89675c518 100644 --- 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 index 86a5587e5a..17008245c5 100644 --- a/ortools/util/fixed_shape_binary_tree.h +++ b/ortools/util/fixed_shape_binary_tree.h @@ -218,6 +218,7 @@ class FixedShapeBinaryTree { template 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()) {