From 8d4b32967c678f4bb0df169a3e44c50c7e498401 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Wed, 31 Jan 2024 14:44:02 +0100 Subject: [PATCH] [CP-SAT] cleanup diffn propagation; remove dead code --- ortools/sat/2d_orthogonal_packing.cc | 283 +++++++++++++++------- ortools/sat/2d_orthogonal_packing.h | 102 ++++++-- ortools/sat/BUILD.bazel | 5 +- ortools/sat/cp_model_lns.cc | 1 - ortools/sat/cp_model_presolve.cc | 1 + ortools/sat/cp_model_search.cc | 14 ++ ortools/sat/diffn.cc | 289 +++++++---------------- ortools/sat/diffn.h | 17 +- ortools/sat/integer_expr.h | 30 ++- ortools/sat/linear_constraint_manager.cc | 10 +- ortools/sat/linear_propagation.cc | 157 +++++++++--- ortools/sat/linear_propagation.h | 24 +- ortools/util/bitset.h | 50 +++- 13 files changed, 609 insertions(+), 374 deletions(-) diff --git a/ortools/sat/2d_orthogonal_packing.cc b/ortools/sat/2d_orthogonal_packing.cc index 627cfdca40..4f50402ea3 100644 --- a/ortools/sat/2d_orthogonal_packing.cc +++ b/ortools/sat/2d_orthogonal_packing.cc @@ -23,9 +23,9 @@ #include "absl/log/check.h" #include "absl/numeric/bits.h" +#include "absl/random/distributions.h" #include "absl/types/span.h" #include "ortools/base/logging.h" -#include "ortools/base/vlog_is_on.h" #include "ortools/sat/integer.h" #include "ortools/sat/util.h" #include "ortools/util/bitset.h" @@ -220,47 +220,19 @@ std::vector GetDffConflict( return result; } -// Check for conflict all combinations of the two Dual Feasible Functions f_0 -// (see documentation for GetDffConflict()) and f_2 (see documentation for -// RoundingDualFeasibleFunction). More precisely, check whether there exist `l` -// and `k` so that -// -// sum_i f_2^k(f_0^l(sizes_x[i])) * sizes_y[i] > f_2^k(f_0^l(x_bb_size)) * -// y_bb_size -// -// The function returns the smallest subset of items enough to make the -// inequality above true or an empty vector if impossible. -std::vector CheckFeasibilityWithDualFunction2( - absl::Span sizes_x, - absl::Span sizes_y, - absl::Span index_by_decreasing_x_size, IntegerValue x_bb_size, - IntegerValue y_bb_size) { - if (x_bb_size == 1) { - return {}; - } - std::vector sizes_x_rescaled; - if (x_bb_size >= std::numeric_limits::max()) { - // To do fast division we want our sizes to fit in a uint16_t. The simplest - // way of doing that is to just first apply this DFF with the right - // power-of-two value of the parameter. - const int log2_k = - absl::bit_width(static_cast(x_bb_size.value() + 1)) - 16 + 1; - const RoundingDualFeasibleFunctionPowerOfTwo dff(x_bb_size, log2_k); - sizes_x_rescaled.resize(sizes_x.size()); - for (int i = 0; i < sizes_x.size(); i++) { - sizes_x_rescaled[i] = dff(sizes_x[i]); - } - x_bb_size = dff(x_bb_size); - CHECK_LT(x_bb_size, std::numeric_limits::max()); - sizes_x = sizes_x_rescaled; - } +} // namespace - // We now want to find the minimum set of values of `k` that would always - // find a conflict if there is a `k` for which it exists. In the literature - // it is often implied (but not stated) that it is sufficient to test the - // values of `k` that correspond to the size of an item. This is not true. To - // find the minimum set of values of `k` we look for all values of `k` that - // are "extreme": ie., the rounding on the division truncates the most (or the +void OrthogonalPackingInfeasibilityDetector::GetAllCandidatesForKForDff2( + absl::Span sizes_x, + absl::Span sizes_y, IntegerValue x_bb_size, + IntegerValue sqrt_bb_size, IntegerValue y_bb_size, + Bitset64& candidates) { + // We want to find the minimum set of values of `k` that would always find a + // conflict if there is a `k` for which it exists. In the literature it is + // often implied (but not stated) that it is sufficient to test the values of + // `k` that correspond to the size of an item. This is not true. To find the + // minimum set of values of `k` we look for all values of `k` that are + // "extreme": ie., the rounding on the division truncates the most (or the // least) amount, depending on the sign it appears in the formula. IntegerValue sum_widths = 0; for (int i = 0; i < sizes_x.size(); i++) { @@ -273,10 +245,9 @@ std::vector CheckFeasibilityWithDualFunction2( } const IntegerValue round_up = sum_widths > 2 * y_bb_size ? 0 : 1; // x_bb_size is less than 65536, so this fits in only 4kib. - Bitset64 candidates(x_bb_size / 2 + 2); + candidates.ClearAndResize(x_bb_size / 2 + 2); // `sqrt_bb_size` is lower than 256. - const IntegerValue sqrt_bb_size = FloorSquareRoot(x_bb_size.value()); for (IntegerValue i = 2; i <= sqrt_bb_size; i++) { candidates.Set(i); } @@ -296,9 +267,6 @@ std::vector CheckFeasibilityWithDualFunction2( } } - std::vector result; - int num_items = sizes_x.size(); - // Remove some bogus candidates added by the logic above. candidates.Clear(0); candidates.Clear(1); @@ -318,6 +286,89 @@ std::vector CheckFeasibilityWithDualFunction2( if (x_bb_size > 3) { candidates.Set(x_bb_size / 3 + 1); } +} + +// Check for conflict all combinations of the two Dual Feasible Functions f_0 +// (see documentation for GetDffConflict()) and f_2 (see documentation for +// RoundingDualFeasibleFunction). More precisely, check whether there exist `l` +// and `k` so that +// +// sum_i f_2^k(f_0^l(sizes_x[i])) * sizes_y[i] > f_2^k(f_0^l(x_bb_size)) * +// y_bb_size +// +// The function returns the smallest subset of items enough to make the +// inequality above true or an empty vector if impossible. +std::vector +OrthogonalPackingInfeasibilityDetector::CheckFeasibilityWithDualFunction2( + absl::Span sizes_x, + absl::Span sizes_y, + absl::Span index_by_decreasing_x_size, IntegerValue x_bb_size, + IntegerValue y_bb_size, int max_number_of_parameters_to_check) { + if (x_bb_size == 1) { + return {}; + } + std::vector sizes_x_rescaled; + if (x_bb_size >= std::numeric_limits::max()) { + // To do fast division we want our sizes to fit in a uint16_t. The simplest + // way of doing that is to just first apply this DFF with the right + // power-of-two value of the parameter. + const int log2_k = + absl::bit_width(static_cast(x_bb_size.value() + 1)) - 16 + 1; + const RoundingDualFeasibleFunctionPowerOfTwo dff(x_bb_size, log2_k); + sizes_x_rescaled.resize(sizes_x.size()); + for (int i = 0; i < sizes_x.size(); i++) { + sizes_x_rescaled[i] = dff(sizes_x[i]); + } + x_bb_size = dff(x_bb_size); + CHECK_LT(x_bb_size, std::numeric_limits::max()); + sizes_x = sizes_x_rescaled; + } + + Bitset64 candidates; + const IntegerValue sqrt_bb_size = FloorSquareRoot(x_bb_size.value()); + int num_items = sizes_x.size(); + const IntegerValue max_possible_number_of_parameters = + std::min(x_bb_size / 4 + 1, sqrt_bb_size * num_items); + if (5ull * max_number_of_parameters_to_check < + max_possible_number_of_parameters) { + // There are many more possible values than what we want to sample. It is + // not worth to pay the price of computing all optimal values to drop most + // of them, so let's just pick it randomly. + candidates.Resize(x_bb_size / 4 + 1); + int num_candidates = 0; + while (num_candidates < max_number_of_parameters_to_check) { + const IntegerValue pick = + absl::Uniform(random_, 1, x_bb_size.value() / 4); + if (!candidates.IsSet(pick)) { + candidates.Set(pick); + num_candidates++; + } + } + } else { + GetAllCandidatesForKForDff2(sizes_x, sizes_y, x_bb_size, sqrt_bb_size, + y_bb_size, candidates); + + if (max_number_of_parameters_to_check < max_possible_number_of_parameters) { + // We might have produced too many candidates. Let's count them and if it + // is the case, sample them. + int count = 0; + for (auto it = candidates.begin(); it != candidates.end(); ++it) { + count++; + } + if (count > max_number_of_parameters_to_check) { + std::vector sampled_candidates( + max_number_of_parameters_to_check); + std::sample(candidates.begin(), candidates.end(), + sampled_candidates.begin(), + max_number_of_parameters_to_check, random_); + candidates.ClearAll(); + for (const IntegerValue k : sampled_candidates) { + candidates.Set(k); + } + } + } + } + std::vector result; // Finally run our small loop to look for the conflict! std::vector modified_sizes(num_items); @@ -341,21 +392,13 @@ std::vector CheckFeasibilityWithDualFunction2( return result; } -} // namespace - -OrthogonalPackingInfeasibilityDetector::Result -OrthogonalPackingInfeasibilityDetector::TestFeasibility( +OrthogonalPackingResult +OrthogonalPackingInfeasibilityDetector::TestFeasibilityImpl( absl::Span sizes_x, absl::Span sizes_y, - std::pair bounding_box_size) { - enum class ConflictType { - NO_CONFLICT, - TRIVIAL, - DFF_F0, - DFF_F2, - }; - ConflictType conflict_type = ConflictType::NO_CONFLICT; - num_calls_++; + std::pair bounding_box_size, + const OrthogonalPackingOptions& options) { + using ConflictType = OrthogonalPackingResult::ConflictType; const int num_items = sizes_x.size(); DCHECK_EQ(num_items, sizes_y.size()); @@ -363,6 +406,11 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility( bounding_box_size.first * bounding_box_size.second; IntegerValue total_energy = 0; + auto make_item = [sizes_x, sizes_y](int i) { + return OrthogonalPackingResult::Item{ + .index = i, .size_x = sizes_x[i], .size_y = sizes_y[i]}; + }; + index_by_decreasing_x_size_.resize(num_items); index_by_decreasing_y_size_.resize(num_items); for (int i = 0; i < num_items; i++) { @@ -371,14 +419,16 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility( index_by_decreasing_y_size_[i] = i; if (sizes_x[i] > bounding_box_size.first || sizes_y[i] > bounding_box_size.second) { - num_conflicts_++; - return Result{.result = Status::INFEASIBLE, - .minimum_unfeasible_subproblem = {i}}; + OrthogonalPackingResult result( + OrthogonalPackingResult::Status::INFEASIBLE); + result.conflict_type_ = ConflictType::TRIVIAL; + result.items_participating_on_conflict_ = {make_item(i)}; + return result; } } if (num_items <= 1) { - return Result{.result = Status::FEASIBLE}; + return OrthogonalPackingResult(OrthogonalPackingResult::Status::FEASIBLE); } std::sort(index_by_decreasing_x_size_.begin(), @@ -389,25 +439,27 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility( [&sizes_y](int a, int b) { return sizes_y[a] > sizes_y[b]; }); // First look for pairwise incompatible pairs. - if (options_.use_pairwise) { + if (options.use_pairwise) { if (auto pair = FindPairwiseConflict(sizes_x, sizes_y, bounding_box_size, index_by_decreasing_x_size_, index_by_decreasing_y_size_); pair.has_value()) { - num_conflicts_++; - num_conflicts_two_items_++; - return Result{.result = Status::INFEASIBLE, - .minimum_unfeasible_subproblem = {pair.value().first, - pair.value().second}}; + OrthogonalPackingResult result( + OrthogonalPackingResult::Status::INFEASIBLE); + result.conflict_type_ = ConflictType::PAIRWISE; + result.items_participating_on_conflict_ = { + make_item(pair.value().first), make_item(pair.value().second)}; + return result; } if (num_items == 2) { - return Result{.result = Status::FEASIBLE}; + return OrthogonalPackingResult(OrthogonalPackingResult::Status::FEASIBLE); } } - std::vector result; + OrthogonalPackingResult result(OrthogonalPackingResult::Status::UNKNOWN); if (total_energy > bb_area) { - conflict_type = ConflictType::TRIVIAL; + result.conflict_type_ = ConflictType::TRIVIAL; + result.result_ = OrthogonalPackingResult::Status::INFEASIBLE; std::vector> index_to_energy; index_to_energy.reserve(num_items); for (int i = 0; i < num_items; i++) { @@ -422,26 +474,39 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility( for (int i = 0; i < index_to_energy.size(); i++) { recomputed_energy += index_to_energy[i].second; if (recomputed_energy > bb_area) { - result.resize(i + 1); + result.items_participating_on_conflict_.resize(i + 1); for (int j = 0; j <= i; j++) { - result[j] = index_to_energy[j].first; + result.items_participating_on_conflict_[j] = + make_item(index_to_energy[j].first); } + result.slack_ = recomputed_energy - bb_area - 1; break; } } } - auto set_conflict_if_better = [&result, &conflict_type]( + const int minimum_conflict_size = options.use_pairwise ? 3 : 2; + if (result.items_participating_on_conflict_.size() == minimum_conflict_size) { + return result; + } + + auto set_conflict_if_better = [&result, &make_item]( const std::vector& conflict, ConflictType type) { if (!conflict.empty() && - (result.empty() || conflict.size() < result.size())) { - conflict_type = type; - result = conflict; + (result.result_ != OrthogonalPackingResult::Status::INFEASIBLE || + conflict.size() < result.items_participating_on_conflict_.size())) { + result.result_ = OrthogonalPackingResult::Status::INFEASIBLE; + result.conflict_type_ = type; + result.items_participating_on_conflict_.clear(); + for (int i : conflict) { + result.items_participating_on_conflict_.push_back(make_item(i)); + } + result.slack_ = 0; // Only supported for trivial for now. } }; - if (options_.use_dff_f0) { + if (options.use_dff_f0) { // If there is no pairwise incompatible pairs, this DFF cannot find a // conflict by enlarging a item on both x and y directions: this would // create an item as long as the whole box and another item as high as the @@ -460,29 +525,56 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility( set_conflict_if_better(conflict, ConflictType::DFF_F0); } - if (options_.use_dff_f2) { + if (result.items_participating_on_conflict_.size() == minimum_conflict_size) { + return result; + } + + if (options.use_dff_f2) { // We only check for conflicts applying this DFF on heights and widths, but // not on both, which would be too expensive if done naively. auto conflict = CheckFeasibilityWithDualFunction2( sizes_x, sizes_y, index_by_decreasing_x_size_, bounding_box_size.first, - bounding_box_size.second); + bounding_box_size.second, + options.dff2_max_number_of_parameters_to_check); set_conflict_if_better(conflict, ConflictType::DFF_F2); + if (result.items_participating_on_conflict_.size() == + minimum_conflict_size) { + return result; + } conflict = CheckFeasibilityWithDualFunction2( sizes_y, sizes_x, index_by_decreasing_y_size_, bounding_box_size.second, - bounding_box_size.first); + bounding_box_size.first, + options.dff2_max_number_of_parameters_to_check); set_conflict_if_better(conflict, ConflictType::DFF_F2); } - if (!result.empty()) { + return result; +} + +OrthogonalPackingResult OrthogonalPackingInfeasibilityDetector::TestFeasibility( + absl::Span sizes_x, + absl::Span sizes_y, + std::pair bounding_box_size, + const OrthogonalPackingOptions& options) { + using ConflictType = OrthogonalPackingResult::ConflictType; + + num_calls_++; + OrthogonalPackingResult result = + TestFeasibilityImpl(sizes_x, sizes_y, bounding_box_size, options); + + if (result.result_ == OrthogonalPackingResult::Status::INFEASIBLE) { num_conflicts_++; - switch (conflict_type) { + switch (result.conflict_type_) { case ConflictType::DFF_F0: num_conflicts_dff0_++; break; case ConflictType::DFF_F2: num_conflicts_dff2_++; break; + case ConflictType::PAIRWISE: + num_conflicts_two_items_++; + break; case ConflictType::TRIVIAL: // The total area of the items was larger than the area of the box. num_trivial_conflicts_++; @@ -491,10 +583,27 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility( LOG(FATAL) << "Should never happen"; break; } - return Result{.result = Status::INFEASIBLE, - .minimum_unfeasible_subproblem = result}; } - return Result{.result = Status::UNKNOWN}; + return result; +} + +bool OrthogonalPackingResult::TryUseSlackToReduceItemSize( + int i, Coord coord, IntegerValue lower_bound) { + Item& item = items_participating_on_conflict_[i]; + IntegerValue& size = coord == Coord::kCoordX ? item.size_x : item.size_y; + const IntegerValue orthogonal_size = + coord == Coord::kCoordX ? item.size_y : item.size_x; + + if (size <= lower_bound || orthogonal_size > slack_) { + return false; + } + const IntegerValue new_size = + std::max(lower_bound, size - slack_ / orthogonal_size); + slack_ -= (size - new_size) * orthogonal_size; + DCHECK_NE(size, new_size); + DCHECK_GE(slack_, 0); + size = new_size; + return true; } } // namespace sat diff --git a/ortools/sat/2d_orthogonal_packing.h b/ortools/sat/2d_orthogonal_packing.h index 6168d53f5a..84274e5d7b 100644 --- a/ortools/sat/2d_orthogonal_packing.h +++ b/ortools/sat/2d_orthogonal_packing.h @@ -20,6 +20,7 @@ #include #include "absl/log/check.h" +#include "absl/random/bit_gen_ref.h" #include "absl/types/span.h" #include "ortools/sat/integer.h" #include "ortools/sat/synchronization.h" @@ -31,6 +32,60 @@ struct OrthogonalPackingOptions { bool use_pairwise = true; bool use_dff_f0 = true; bool use_dff_f2 = true; + int dff2_max_number_of_parameters_to_check = std::numeric_limits::max(); +}; + +class OrthogonalPackingResult { + public: + explicit OrthogonalPackingResult() : result_(Status::UNKNOWN) {} + enum class Status { + INFEASIBLE, + FEASIBLE, + UNKNOWN, + }; + + Status GetResult() const { return result_; } + + struct Item { + // Index of the item on the original sizes_x/sizes_y input. + int index; + // New size for item of index `i` which is smaller or equal to the initial + // size. The subproblem remain infeasible if every item is shrinked to its + // new size. + IntegerValue size_x; + IntegerValue size_y; + }; + const std::vector& GetItemsParticipatingOnConflict() const { + return items_participating_on_conflict_; + } + + bool HasSlack() const { return slack_ > IntegerValue(0); } + + enum class Coord { kCoordX, kCoordY }; + // Use an eventual slack to reduce the size of item corresponding to the + // `i`-th element on GetItemsParticipatingOnConflict(). It will not use any + // slack to reduce it beyond lower_bound. This is a no-op if HasSlack() is + // false. + bool TryUseSlackToReduceItemSize(int i, Coord coord, + IntegerValue lower_bound = 0); + + private: + friend class OrthogonalPackingInfeasibilityDetector; + + explicit OrthogonalPackingResult(Status result) : result_(result) {} + + enum class ConflictType { + NO_CONFLICT, + TRIVIAL, + PAIRWISE, + DFF_F0, + DFF_F2, + }; + + Status result_; + ConflictType conflict_type_ = ConflictType::NO_CONFLICT; + IntegerValue slack_ = 0; + std::vector items_participating_on_conflict_; }; // Class for solving the orthogonal packing problem when it can be done @@ -38,29 +93,43 @@ struct OrthogonalPackingOptions { class OrthogonalPackingInfeasibilityDetector { public: explicit OrthogonalPackingInfeasibilityDetector( - SharedStatistics* shared_stats, - const OrthogonalPackingOptions& options = OrthogonalPackingOptions()) - : options_(options), shared_stats_(shared_stats) {} + absl::BitGenRef random, SharedStatistics* shared_stats) + : random_(random), shared_stats_(shared_stats) {} ~OrthogonalPackingInfeasibilityDetector(); - enum class Status { - INFEASIBLE, - FEASIBLE, - UNKNOWN, - }; - struct Result { - Status result; - std::vector minimum_unfeasible_subproblem; - }; - - Result TestFeasibility( + OrthogonalPackingResult TestFeasibility( absl::Span sizes_x, absl::Span sizes_y, - std::pair bounding_box_size); + std::pair bounding_box_size, + const OrthogonalPackingOptions& options = OrthogonalPackingOptions()); private: - const OrthogonalPackingOptions options_; + OrthogonalPackingResult TestFeasibilityImpl( + absl::Span sizes_x, + absl::Span sizes_y, + std::pair bounding_box_size, + const OrthogonalPackingOptions& options); + + // Returns a minimum set of values of the parameter `k` of f_2^k that is + // sufficient to find a conflict. This function runs in + // O(num_items * sqrt(bb_size)) operations. + // All sizes must be positive values less than UINT16_MAX. + // The returned bitset will contain less elements than + // min(sqrt_bb_size * num_items, x_bb_size/4+1). + void GetAllCandidatesForKForDff2(absl::Span sizes_x, + absl::Span sizes_y, + IntegerValue x_bb_size, + IntegerValue sqrt_bb_size, + IntegerValue y_bb_size, + Bitset64& candidates); + + std::vector CheckFeasibilityWithDualFunction2( + absl::Span sizes_x, + absl::Span sizes_y, + absl::Span index_by_decreasing_x_size, IntegerValue x_bb_size, + IntegerValue y_bb_size, int max_number_of_parameters_to_check); + std::vector index_by_decreasing_x_size_; std::vector index_by_decreasing_y_size_; @@ -71,6 +140,7 @@ class OrthogonalPackingInfeasibilityDetector { int64_t num_conflicts_dff2_ = 0; int64_t num_conflicts_dff0_ = 0; + absl::BitGenRef random_; SharedStatistics* shared_stats_; }; diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index 3d2bb71ca6..9484ada2d6 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -266,6 +266,7 @@ cc_library( ":cp_model_utils", ":integer", ":integer_search", + ":linear_propagation", ":model", ":sat_base", ":sat_parameters_cc_proto", @@ -1823,6 +1824,8 @@ cc_library( "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/numeric:bits", + "@com_google_absl//absl/random:bit_gen_ref", + "@com_google_absl//absl/random:distributions", "@com_google_absl//absl/types:span", ], ) @@ -1930,7 +1933,6 @@ cc_library( ":subsolver", ":synchronization", ":util", - "//ortools/base", "//ortools/base:stl_util", "//ortools/graph:connected_components", "//ortools/util:adaptative_parameter_value", @@ -1943,6 +1945,7 @@ cc_library( "@com_google_absl//absl/base:core_headers", "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/container:flat_hash_set", + "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/meta:type_traits", "@com_google_absl//absl/random:bit_gen_ref", diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index db2cd889d9..5a78cf33f3 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -38,7 +38,6 @@ #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" -#include "ortools/base/vlog_is_on.h" #include "ortools/graph/connected_components.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_mapping.h" diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 35abe4782c..e8b576cdfa 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -5438,6 +5438,7 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) { return RemoveConstraint(ct); } + // TODO(user): handle this case. See issue #4068. if (!has_zero_sized_interval && (x_constant || y_constant)) { context_->UpdateRuleStats( "no_overlap_2d: a dimension is constant, splitting into many " diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index c1c8588277..8059c372b6 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -34,6 +34,7 @@ #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_search.h" +#include "ortools/sat/linear_propagation.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" #include "ortools/sat/sat_parameters.pb.h" @@ -327,12 +328,25 @@ std::function ConstructHeuristicSearchStrategy( if (ModelHasSchedulingConstraints(cp_model_proto)) { std::vector> heuristics; const auto& params = *model->GetOrCreate(); + bool possible_new_constraints = false; if (params.use_dynamic_precedence_in_disjunctive()) { + possible_new_constraints = true; heuristics.push_back(DisjunctivePrecedenceSearchHeuristic(model)); } if (params.use_dynamic_precedence_in_cumulative()) { + possible_new_constraints = true; heuristics.push_back(CumulativePrecedenceSearchHeuristic(model)); } + + // Tricky: we need to create this at level zero in case there are no linear + // constraint in the model at the beginning. + // + // TODO(user): Alternatively, support creation of SatPropagator at positive + // level. + if (possible_new_constraints && params.new_linear_propagation()) { + model->GetOrCreate(); + } + heuristics.push_back(SchedulingSearchHeuristic(model)); return SequentialSearch(std::move(heuristics)); } diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index a1c3a6ae7c..0a752ddd1d 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -246,17 +246,14 @@ NonOverlappingRectanglesEnergyPropagator:: {"NonOverlappingRectanglesEnergyPropagator/called", num_calls_}); stats.push_back( {"NonOverlappingRectanglesEnergyPropagator/conflicts", num_conflicts_}); - stats.push_back( - {"NonOverlappingRectanglesEnergyPropagator/multiple_conflicts", - num_multiple_conflicts_}); stats.push_back( {"NonOverlappingRectanglesEnergyPropagator/conflicts_two_boxes", num_conflicts_two_boxes_}); stats.push_back({"NonOverlappingRectanglesEnergyPropagator/refined", num_refined_conflicts_}); stats.push_back( - {"NonOverlappingRectanglesEnergyPropagator/orthogonal_packing_conflict", - num_orthogonal_packing_conflict_}); + {"NonOverlappingRectanglesEnergyPropagator/conflicts_with_slack", + num_conflicts_with_slack_}); shared_stats_->AddStats(stats); } @@ -267,7 +264,6 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() { if (!x_.SynchronizeAndSetTimeDirection(true)) return false; if (!y_.SynchronizeAndSetTimeDirection(true)) return false; - num_calls_++; Rectangle bounding_box = {.x_min = std::numeric_limits::max(), .x_max = std::numeric_limits::min(), .y_min = std::numeric_limits::max(), @@ -310,6 +306,8 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() { return true; } + num_calls_++; + std::optional best_conflict = FindConflict(std::move(active_box_ranges)); if (!best_conflict.has_value()) { @@ -321,18 +319,30 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() { // search for a best explanation. This is specially the case since we only // want to re-run it over the items that participate in the conflict, so it is // a much smaller problem. - IntegerValue best_explanation_size = best_conflict->items.size(); + IntegerValue best_explanation_size = + best_conflict->opp_result.GetItemsParticipatingOnConflict().size(); bool refined = false; while (true) { - const std::optional conflict = FindConflict(best_conflict->items); + std::vector items_participating_in_conflict; + items_participating_in_conflict.reserve( + best_conflict->items_for_opp.size()); + for (const auto& item : + best_conflict->opp_result.GetItemsParticipatingOnConflict()) { + items_participating_in_conflict.push_back( + best_conflict->items_for_opp[item.index]); + } + std::optional conflict = + FindConflict(items_participating_in_conflict); if (!conflict.has_value()) break; // We prefer an explanation with the least number of boxes. - if (conflict->items.size() >= best_explanation_size) { + if (conflict->opp_result.GetItemsParticipatingOnConflict().size() >= + best_explanation_size) { // The new explanation isn't better than the old one. Stop trying. break; } - best_explanation_size = conflict->items.size(); - best_conflict = conflict; + best_explanation_size = + conflict->opp_result.GetItemsParticipatingOnConflict().size(); + best_conflict = std::move(conflict); refined = true; } @@ -342,9 +352,6 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() { if (best_explanation_size == 2) { num_conflicts_two_boxes_++; } - if (best_conflict->opp_problem_size > 0) { - num_orthogonal_packing_conflict_++; - } BuildAndReportEnergyTooLarge(generalized_explanation); return false; } @@ -359,28 +366,13 @@ NonOverlappingRectanglesEnergyPropagator::FindConflict( rectangles_with_too_much_energy.candidates.empty()) { return std::nullopt; } - num_multiple_conflicts_ += - rectangles_with_too_much_energy.conflicts.size() > 1; - std::vector best_explanation; - Rectangle best_rectangle; - for (const Rectangle& r : rectangles_with_too_much_energy.conflicts) { - const std::vector range_for_explanation = - GetEnergyConflictForRectangle(r, active_box_ranges); - CheckPropagationIsValid(range_for_explanation, r); - if (best_explanation.empty() || - best_explanation.size() > range_for_explanation.size()) { - best_explanation = range_for_explanation; - best_rectangle = r; - } - } + Conflict best_conflict; - // Now try with a DFF on the candidates. - int opp_problem_size = 0; - // Also try the DFF on known conflicts, hopefully we will get a smaller - // explanation. - // Sample 10 rectangles for looking for a conflict with DFF. This avoids - // making this heuristic too costly. + // Sample 10 rectangles (at least five among the ones for which we already + // detected an energy overflow), extract an orthogonal packing subproblem for + // each and look for conflict. Sampling avoids making this heuristic too + // costly. constexpr int kSampleSize = 10; absl::InlinedVector sampled_rectangles; std::sample(rectangles_with_too_much_energy.conflicts.begin(), @@ -430,158 +422,87 @@ NonOverlappingRectanglesEnergyPropagator::FindConflict( // our rectangle is the bounding box, and we need to fit inside it a set of // items corresponding to the minimum intersection of the original items // with this bounding box. + const bool use_expensive_heuristics = + (sizes_x.size() < 10) || + // Propagating a large conflict is expensive, so it's worth running an + // expensive heuristic to make it smaller. + best_conflict.opp_result.GetResult() == + OrthogonalPackingResult::Status::INFEASIBLE || + !rectangles_with_too_much_energy.conflicts.empty(); const auto opp_result = orthogonal_packing_checker_.TestFeasibility( - sizes_x, sizes_y, {r.SizeX(), r.SizeY()}); - if (opp_result.result == - OrthogonalPackingInfeasibilityDetector::Status::INFEASIBLE && - (best_explanation.empty() || - best_explanation.size() > - opp_result.minimum_unfeasible_subproblem.size())) { - best_explanation.clear(); - for (int idx : opp_result.minimum_unfeasible_subproblem) { - best_explanation.push_back(filtered_items[idx]); - } - best_rectangle = r; - opp_problem_size = static_cast(active_box_ranges.size()); + sizes_x, sizes_y, {r.SizeX(), r.SizeY()}, + OrthogonalPackingOptions{.use_pairwise = true, + .use_dff_f0 = true, + .use_dff_f2 = true, + .dff2_max_number_of_parameters_to_check = + use_expensive_heuristics ? 25 : 3}); + if (opp_result.GetResult() == OrthogonalPackingResult::Status::INFEASIBLE && + (best_conflict.opp_result.GetResult() != + OrthogonalPackingResult::Status::INFEASIBLE || + best_conflict.opp_result.GetItemsParticipatingOnConflict().size() > + opp_result.GetItemsParticipatingOnConflict().size())) { + best_conflict.items_for_opp = filtered_items; + best_conflict.opp_result = opp_result; + best_conflict.rectangle_with_too_much_energy = r; } // Use the fact that our rectangles are ordered in shrinking order to remove // all items that will never contribute again. filtered_items.swap(active_box_ranges); } - if (!best_explanation.empty()) { - return Conflict{.items = std::move(best_explanation), - .rectangle_with_too_much_energy = best_rectangle, - .opp_problem_size = opp_problem_size}; + if (best_conflict.opp_result.GetResult() == + OrthogonalPackingResult::Status::INFEASIBLE) { + return best_conflict; + } else { + return std::nullopt; } - return std::nullopt; } std::vector NonOverlappingRectanglesEnergyPropagator::GeneralizeExplanation( const Conflict& conflict) { const Rectangle& rectangle = conflict.rectangle_with_too_much_energy; - IntegerValue available_energy = rectangle.Area(); - IntegerValue used_energy = 0; - for (const RectangleInRange& box : conflict.items) { - IntegerValue energy = box.GetMinimumIntersectionArea(rectangle); - used_energy += energy; - } - std::vector ranges_for_explanation(conflict.items); + OrthogonalPackingResult relaxed_result = conflict.opp_result; - IntegerValue slack = used_energy - available_energy - 1; - - if (conflict.opp_problem_size > 0) { - VLOG_EVERY_N_SEC(2, 3) - << "Rectangle with energy overflow after solving OPP: " << rectangle - << " using " << ranges_for_explanation.size() << "/" - << conflict.opp_problem_size << " items."; - } else { - VLOG_EVERY_N_SEC(2, 3) << "Rectangle with energy overflow: " << rectangle - << " (" << used_energy << "/" << available_energy - << "), with " << ranges_for_explanation.size() - << " boxes inside."; - } - - for (auto& range : ranges_for_explanation) { - Rectangle overlap = range.GetMinimumIntersection(rectangle); - DCHECK_GT(overlap.SizeX(), 0); - DCHECK_GT(overlap.SizeY(), 0); - range = RectangleInRange::BiggestWithMinIntersection( - rectangle, range, overlap.SizeX(), overlap.SizeY()); - - // Now use the energy slack to make the intervals even bigger. - if (overlap.SizeX() < slack) { - IntegerValue y_slack = slack / overlap.SizeX(); - // We know we can remove an area as big as y_slack * x_size. That means - // we can reduce the start of an interval of y_slack or increase the end - // by the same amount. But there is no point wasting the slack on a - // interval larger than the zero-level interval, so we check for that. - const IntegerValue y_min_slack = std::max( - IntegerValue(0), - std::min(y_slack, range.bounding_area.y_min - - y_.LevelZeroStartMin(range.box_index))); - if (y_min_slack > 0) { - range.bounding_area.y_min -= y_min_slack; - slack -= y_min_slack * overlap.SizeX(); - y_slack -= y_min_slack; - } - const IntegerValue y_max_slack = - std::min(y_slack, y_.LevelZeroStartMax(range.box_index) - - range.bounding_area.y_max); - if (y_max_slack > 0) { - range.bounding_area.y_max += y_max_slack; - slack -= y_max_slack * overlap.SizeX(); - } - // Recompute the overlap, since we changed the item dimensions. - overlap = range.GetMinimumIntersection(rectangle); - } - if (overlap.SizeY() < slack) { - IntegerValue x_slack = slack / overlap.SizeY(); - const IntegerValue x_min_slack = std::max( - IntegerValue(0), - std::min(x_slack, range.bounding_area.x_min - - x_.LevelZeroStartMax(range.box_index))); - if (x_min_slack > 0) { - range.bounding_area.x_min -= x_min_slack; - slack -= x_min_slack * overlap.SizeY(); - x_slack -= x_min_slack; - } - const IntegerValue x_max_slack = - std::min(x_slack, x_.LevelZeroStartMax(range.box_index) - - range.bounding_area.x_max); - if (x_max_slack > 0) { - range.bounding_area.x_max += x_max_slack; - slack -= x_max_slack * overlap.SizeY(); - } - } - } - if (conflict.opp_problem_size == 0) { - CHECK_GT(used_energy, available_energy); - CHECK_GT(available_energy, 0); - } - - return ranges_for_explanation; -} - -std::vector -NonOverlappingRectanglesEnergyPropagator::GetEnergyConflictForRectangle( - const Rectangle& rectangle, - const std::vector& active_box_ranges) { - struct OverlapPerBox { - IntegerValue energy; - IntegerValue overlap_x_size; - IntegerValue overlap_y_size; - int index; - }; - - // First select the minimum amount of elements that would be enough to exceed - // the energy. - std::vector energy_per_box(active_box_ranges.size()); - for (int i = 0; i < active_box_ranges.size(); ++i) { - const Rectangle overlap = - active_box_ranges[i].GetMinimumIntersection(rectangle); - OverlapPerBox& box = energy_per_box[i]; - box.overlap_x_size = overlap.SizeX(); - box.overlap_y_size = overlap.SizeY(); - box.index = i; - box.energy = box.overlap_x_size * box.overlap_y_size; - } - std::sort(energy_per_box.begin(), energy_per_box.end(), - [](const OverlapPerBox& a, const OverlapPerBox& b) { - return a.energy > b.energy; - }); - const IntegerValue available_energy = rectangle.Area(); - IntegerValue used_energy = 0; - std::vector ranges_for_explanation; - ranges_for_explanation.reserve(energy_per_box.size()); - for (int i = 0; i < energy_per_box.size(); ++i) { - const OverlapPerBox& box = energy_per_box[i]; - used_energy += box.energy; - CHECK_GT(box.energy, 0); - ranges_for_explanation.push_back(active_box_ranges[box.index]); - if (used_energy > available_energy) { + // Use the potential slack to have a stronger reason. + int used_slack_count = 0; + const auto& items = relaxed_result.GetItemsParticipatingOnConflict(); + for (int i = 0; i < items.size(); ++i) { + if (!relaxed_result.HasSlack()) { break; } + const RectangleInRange& range = conflict.items_for_opp[items[i].index]; + const RectangleInRange item_in_zero_level_range = { + .bounding_area = {.x_min = x_.LevelZeroStartMin(range.box_index), + .x_max = x_.LevelZeroStartMax(range.box_index) + + range.x_size, + .y_min = y_.LevelZeroStartMin(range.box_index), + .y_max = y_.LevelZeroStartMax(range.box_index) + + range.y_size}, + .x_size = range.x_size, + .y_size = range.y_size}; + // There is no point trying to intersect less the item with the rectangle + // than it would on zero-level. So do not waste the slack with it. + const Rectangle max_overlap = + item_in_zero_level_range.GetMinimumIntersection(rectangle); + used_slack_count += relaxed_result.TryUseSlackToReduceItemSize( + i, OrthogonalPackingResult::Coord::kCoordX, max_overlap.SizeX()); + used_slack_count += relaxed_result.TryUseSlackToReduceItemSize( + i, OrthogonalPackingResult::Coord::kCoordY, max_overlap.SizeY()); + } + num_conflicts_with_slack_ += (used_slack_count > 0); + VLOG_EVERY_N_SEC(2, 3) + << "Found a conflict on the OPP sub-problem of rectangle: " << rectangle + << " using " + << conflict.opp_result.GetItemsParticipatingOnConflict().size() << "/" + << conflict.items_for_opp.size() << " items."; + + std::vector ranges_for_explanation; + ranges_for_explanation.reserve(conflict.items_for_opp.size()); + for (const auto& item : relaxed_result.GetItemsParticipatingOnConflict()) { + const RectangleInRange& range = conflict.items_for_opp[item.index]; + ranges_for_explanation.push_back( + RectangleInRange::BiggestWithMinIntersection(rectangle, range, + item.size_x, item.size_y)); } return ranges_for_explanation; } @@ -596,34 +517,6 @@ int NonOverlappingRectanglesEnergyPropagator::RegisterWith( return id; } -void NonOverlappingRectanglesEnergyPropagator::CheckPropagationIsValid( - const std::vector& ranges, - const Rectangle& rectangle_with_too_much_energy) { - const IntegerValue available_energy = rectangle_with_too_much_energy.Area(); - IntegerValue used_energy = 0; - for (const auto& range : ranges) { - const int b = range.box_index; - // Check that we are explaining intervals at least as large as the current. - CHECK_GE(x_.StartMin(b), range.bounding_area.x_min); - CHECK_GE(range.bounding_area.x_max - range.x_size, x_.StartMax(b)); - CHECK_GE(y_.StartMin(b), range.bounding_area.y_min); - CHECK_GE(range.bounding_area.y_max - range.y_size, y_.StartMax(b)); - - // Each one of the boxes-in-range that we found on the cut does intersect - // the rectangle we found. - const auto intersection = - range.GetMinimumIntersection(rectangle_with_too_much_energy); - CHECK_GT(intersection.Area(), 0); - // It cannot intersect more than the size of the object. - CHECK_GE(x_.SizeMin(b), intersection.SizeX()); - CHECK_GE(y_.SizeMin(b), intersection.SizeY()); - used_energy += intersection.Area(); - } - - // We have more area inside the rectangle than the area of the rectangle. - CHECK_GT(used_energy, available_energy); -} - bool NonOverlappingRectanglesEnergyPropagator::BuildAndReportEnergyTooLarge( const std::vector& ranges) { if (ranges.size() == 2) { diff --git a/ortools/sat/diffn.h b/ortools/sat/diffn.h index 3792d405b7..8fd0353975 100644 --- a/ortools/sat/diffn.h +++ b/ortools/sat/diffn.h @@ -45,7 +45,7 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { y_(*y), random_(model->GetOrCreate()), shared_stats_(model->GetOrCreate()), - orthogonal_packing_checker_(shared_stats_) {} + orthogonal_packing_checker_(*random_, shared_stats_) {} ~NonOverlappingRectanglesEnergyPropagator() override; @@ -54,11 +54,10 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { private: struct Conflict { - std::vector items; + // The Orthogonal Packing subproblem we used. + std::vector items_for_opp; Rectangle rectangle_with_too_much_energy; - // Size of the Orthogonal Packing Problem we used. Can be smaller than - // `items.size()` if we found a OPP conflict that didn't need all items. - int opp_problem_size; + OrthogonalPackingResult opp_result; }; std::optional FindConflict( std::vector active_box_ranges); @@ -67,11 +66,6 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { bool BuildAndReportEnergyTooLarge( const std::vector& ranges); - void CheckPropagationIsValid(const std::vector& ranges, - const Rectangle& rectangle_with_too_much_energy); - std::vector GetEnergyConflictForRectangle( - const Rectangle& rectangle, - const std::vector& active_box_ranges); SchedulingConstraintHelper& x_; SchedulingConstraintHelper& y_; @@ -81,10 +75,9 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface { int64_t num_calls_ = 0; int64_t num_conflicts_ = 0; - int64_t num_multiple_conflicts_ = 0; int64_t num_conflicts_two_boxes_ = 0; int64_t num_refined_conflicts_ = 0; - int64_t num_orthogonal_packing_conflict_ = 0; + int64_t num_conflicts_with_slack_ = 0; NonOverlappingRectanglesEnergyPropagator( const NonOverlappingRectanglesEnergyPropagator&) = delete; diff --git a/ortools/sat/integer_expr.h b/ortools/sat/integer_expr.h index 764b543ecb..2c9e099dea 100644 --- a/ortools/sat/integer_expr.h +++ b/ortools/sat/integer_expr.h @@ -455,10 +455,19 @@ inline std::function WeightedSumLowerOrEqual( } if (params.new_linear_propagation()) { - model->GetOrCreate()->AddConstraint( + const bool ok = model->GetOrCreate()->AddConstraint( {}, vars, std::vector(coefficients.begin(), coefficients.end()), IntegerValue(upper_bound)); + if (!ok) { + auto* sat_solver = model->GetOrCreate(); + if (sat_solver->CurrentDecisionLevel() == 0) { + sat_solver->NotifyThatModelIsUnsat(); + } else { + LOG(FATAL) << "We currently do not support adding conflicting " + "constraint at positive level."; + } + } } else { IntegerSumLE* constraint = new IntegerSumLE( {}, vars, @@ -544,8 +553,8 @@ inline std::function ConditionalWeightedSumLowerOrEqual( for (int i = 0; i < vars.size(); ++i) { expression_min += coefficients[i] * (coefficients[i] >= 0 - ? integer_trail->LowerBound(vars[i]) - : integer_trail->UpperBound(vars[i])); + ? integer_trail->LevelZeroLowerBound(vars[i]) + : integer_trail->LevelZeroUpperBound(vars[i])); } if (expression_min == upper_bound) { // Tricky: as we create integer literal, we might propagate stuff and @@ -554,12 +563,12 @@ inline std::function ConditionalWeightedSumLowerOrEqual( IntegerValue non_cached_min; for (int i = 0; i < vars.size(); ++i) { if (coefficients[i] > 0) { - const IntegerValue lb = integer_trail->LowerBound(vars[i]); + const IntegerValue lb = integer_trail->LevelZeroLowerBound(vars[i]); non_cached_min += coefficients[i] * lb; model->Add(Implication(enforcement_literals, IntegerLiteral::LowerOrEqual(vars[i], lb))); } else if (coefficients[i] < 0) { - const IntegerValue ub = integer_trail->UpperBound(vars[i]); + const IntegerValue ub = integer_trail->LevelZeroUpperBound(vars[i]); non_cached_min += coefficients[i] * ub; model->Add(Implication(enforcement_literals, IntegerLiteral::GreaterOrEqual(vars[i], ub))); @@ -574,10 +583,19 @@ inline std::function ConditionalWeightedSumLowerOrEqual( } } else { if (params.new_linear_propagation()) { - model->GetOrCreate()->AddConstraint( + const bool ok = model->GetOrCreate()->AddConstraint( enforcement_literals, vars, std::vector(coefficients.begin(), coefficients.end()), IntegerValue(upper_bound)); + if (!ok) { + auto* sat_solver = model->GetOrCreate(); + if (sat_solver->CurrentDecisionLevel() == 0) { + sat_solver->NotifyThatModelIsUnsat(); + } else { + LOG(FATAL) << "We currently do not support adding conflicting " + "constraint at positive level."; + } + } } else { IntegerSumLE* constraint = new IntegerSumLE( enforcement_literals, vars, diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index d60f0f1dcb..ffc597ec59 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -447,8 +447,6 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { {CeilRatio(threshold, IntegerValue(2)), threshold - min_magnitude, std::min(threshold_lb, threshold_ub)}); if (max_magnitude > second_threshold) { - term_changed = true; - ++num_coeff_strenghtening_; const int num_terms = ct->num_terms; for (int i = 0; i < num_terms; ++i) { // In all cases, we reason on a transformed constraint where the term @@ -460,18 +458,26 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); if (coeff > threshold) { + term_changed = true; + ++num_coeff_strenghtening_; ct->coeffs[i] = threshold; ct->ub -= (coeff - threshold) * ub; ct->lb -= (coeff - threshold) * lb; } else if (coeff > second_threshold && coeff < threshold) { + term_changed = true; + ++num_coeff_strenghtening_; ct->coeffs[i] = second_threshold; ct->ub -= (coeff - second_threshold) * ub; ct->lb -= (coeff - second_threshold) * lb; } else if (coeff < -threshold) { + term_changed = true; + ++num_coeff_strenghtening_; ct->coeffs[i] = -threshold; ct->ub -= (coeff + threshold) * lb; ct->lb -= (coeff + threshold) * ub; } else if (coeff < -second_threshold && coeff > -threshold) { + term_changed = true; + ++num_coeff_strenghtening_; ct->coeffs[i] = -second_threshold; ct->ub -= (coeff + second_threshold) * lb; ct->lb -= (coeff + second_threshold) * ub; diff --git a/ortools/sat/linear_propagation.cc b/ortools/sat/linear_propagation.cc index 7b47b4418a..87c5fc7f33 100644 --- a/ortools/sat/linear_propagation.cc +++ b/ortools/sat/linear_propagation.cc @@ -46,12 +46,24 @@ namespace sat { void CustomFifoQueue::IncreaseSize(int n) { CHECK_GE(n, pos_.size()); - CHECK_LE(left_, right_); pos_.resize(n, -1); tmp_positions_.resize(n, 0); // We need 1 more space since we can add at most n element. + if (n + 1 < queue_.size()) return; + const int old_size = queue_.size(); queue_.resize(n + 1, 0); + queue_.resize(queue_.capacity(), 0); + + // We need to reconstruct the queue in this case. + if (left_ > right_) { + const int diff = queue_.size() - old_size; + for (int i = old_size - 1; i >= left_; --i) { + pos_[queue_[i]] = i + diff; + queue_[i + diff] = queue_[i]; + } + left_ += diff; + } } int CustomFifoQueue::Pop() { @@ -230,10 +242,11 @@ void EnforcementPropagator::Untrail(const Trail& /*trail*/, int trail_index) { EnforcementId EnforcementPropagator::Register( absl::Span enforcement, std::function callback) { - CHECK_EQ(trail_.CurrentDecisionLevel(), 0); int num_true = 0; int num_false = 0; + bool is_always_false = false; temp_literals_.clear(); + const int level = trail_.CurrentDecisionLevel(); for (const Literal l : enforcement) { // Make sure we always have enough room for the literal and its negation. const int size = std::max(l.Index().value(), l.NegatedIndex().value()) + 1; @@ -241,23 +254,25 @@ EnforcementId EnforcementPropagator::Register( watcher_.resize(size); } if (assignment_.LiteralIsTrue(l)) { + if (level == 0 || trail_.Info(l.Variable()).level == 0) continue; ++num_true; - continue; - } - if (assignment_.LiteralIsFalse(l)) { + } else if (assignment_.LiteralIsFalse(l)) { + if (level == 0 || trail_.Info(l.Variable()).level == 0) { + is_always_false = true; + break; + } ++num_false; - continue; } temp_literals_.push_back(l); } gtl::STLSortAndRemoveDuplicates(&temp_literals_); // Return special indices if never/always enforced. - if (num_false > 0) { + if (is_always_false) { if (callback != nullptr) callback(EnforcementStatus::IS_FALSE); return EnforcementId(-1); } - if (num_true == enforcement.size()) { + if (temp_literals_.empty()) { if (callback != nullptr) callback(EnforcementStatus::IS_ENFORCED); return EnforcementId(-1); } @@ -267,15 +282,56 @@ EnforcementId EnforcementPropagator::Register( CHECK(!temp_literals_.empty()); buffer_.insert(buffer_.end(), temp_literals_.begin(), temp_literals_.end()); - starts_.push_back(buffer_.size()); // Sentinel. - statuses_.push_back(EnforcementStatus::CANNOT_PROPAGATE); + starts_.push_back(buffer_.size()); // Sentinel/next-start. + + // The default status at level zero. + statuses_.push_back(temp_literals_.size() == 1 + ? EnforcementStatus::CAN_PROPAGATE + : EnforcementStatus::CANNOT_PROPAGATE); if (temp_literals_.size() == 1) { watcher_[temp_literals_[0].Index()].push_back(id); - ChangeStatus(id, EnforcementStatus::CAN_PROPAGATE); } else { - watcher_[temp_literals_[0].Index()].push_back(id); - watcher_[temp_literals_[1].Index()].push_back(id); + // Make sure we watch correct literals. + const auto span = GetSpan(id); + int num_not_true = 0; + for (int i = 0; i < span.size(); ++i) { + if (assignment_.LiteralIsTrue(span[i])) continue; + std::swap(span[num_not_true], span[i]); + ++num_not_true; + if (num_not_true == 2) break; + } + + // We need to watch one of the literals at highest level. + if (num_not_true == 1) { + int max_level = trail_.Info(span[1].Variable()).level; + for (int i = 2; i < span.size(); ++i) { + const int level = trail_.Info(span[i].Variable()).level; + if (level > max_level) { + max_level = level; + std::swap(span[1], span[i]); + } + } + } + + watcher_[span[0].Index()].push_back(id); + watcher_[span[1].Index()].push_back(id); + } + + // Change status, call callback and set up untrail if the status is different + // from EnforcementStatus::CANNOT_PROPAGATE. + if (num_false > 0) { + ChangeStatus(id, EnforcementStatus::IS_FALSE); + } else if (num_true == temp_literals_.size()) { + ChangeStatus(id, EnforcementStatus::IS_ENFORCED); + } else if (num_true + 1 == temp_literals_.size()) { + ChangeStatus(id, EnforcementStatus::CAN_PROPAGATE); + // Because this is the default status, we still need to call the callback. + if (temp_literals_.size() == 1) { + if (callbacks_[id] != nullptr) { + callbacks_[id](EnforcementStatus::CAN_PROPAGATE); + } + } } return id; } @@ -390,6 +446,22 @@ void EnforcementPropagator::ChangeStatus(EnforcementId id, if (callbacks_[id] != nullptr) callbacks_[id](new_status); } +EnforcementStatus EnforcementPropagator::DebugStatus(EnforcementId id) { + if (id < 0) return EnforcementStatus::IS_ENFORCED; + + int num_true = 0; + for (const Literal l : GetSpan(id)) { + if (assignment_.LiteralIsFalse(l)) { + return EnforcementStatus::IS_FALSE; + } + if (assignment_.LiteralIsTrue(l)) ++num_true; + } + const int size = GetSpan(id).size(); + if (num_true == size) return EnforcementStatus::IS_ENFORCED; + if (num_true + 1 == size) return EnforcementStatus::CAN_PROPAGATE; + return EnforcementStatus::CANNOT_PROPAGATE; +} + LinearPropagator::LinearPropagator(Model* model) : trail_(model->GetOrCreate()), integer_trail_(model->GetOrCreate()), @@ -438,13 +510,14 @@ void LinearPropagator::SetLevel(int level) { while (!propagation_queue_.empty()) { in_queue_[propagation_queue_.Pop()] = false; } - for (int i = rev_at_false_size_; i < in_queue_and_at_false_.size(); ++i) { - in_queue_[in_queue_and_at_false_[i]] = false; + for (int i = rev_unenforced_size_; i < unenforced_constraints_.size(); + ++i) { + in_queue_[unenforced_constraints_[i]] = false; } - in_queue_and_at_false_.resize(rev_at_false_size_); + unenforced_constraints_.resize(rev_unenforced_size_); } else if (level > previous_level_) { - rev_at_false_size_ = in_queue_and_at_false_.size(); - rev_int_repository_->SaveState(&rev_at_false_size_); + rev_unenforced_size_ = unenforced_constraints_.size(); + rev_int_repository_->SaveState(&rev_unenforced_size_); } previous_level_ = level; @@ -468,12 +541,12 @@ bool LinearPropagator::Propagate() { AddWatchedToQueue(var); } - // TODO(user): Abort this propagator as soon as a Boolean is propagated ? so - // that we always finish the Boolean propagation first. This can happen when - // we push a bound that has associated Booleans. The idea is to resume from - // our current state when we are called again. Note however that we have to - // clear the propagated_by_ info has other propagator might have pushed the - // same variable further. + // We abort this propagator as soon as a Boolean is propagated, so that we + // always finish the Boolean propagation first. This can happen when we push a + // bound that has associated Booleans or push enforcement to false. The idea + // is to resume from our current state when we are called again. Note however + // that we have to clear the propagated_by_ info (as done above) has other + // propagator might have pushed the same variable further. // // Empty FIFO queue. const int saved_index = trail_->Index(); @@ -497,13 +570,15 @@ bool LinearPropagator::Propagate() { } // Adds a new constraint to the propagator. -void LinearPropagator::AddConstraint( +bool LinearPropagator::AddConstraint( absl::Span enforcement_literals, absl::Span vars, absl::Span coeffs, IntegerValue upper_bound) { - if (vars.empty()) return; - for (const Literal l : enforcement_literals) { - if (trail_->Assignment().LiteralIsFalse(l)) return; + if (vars.empty()) return true; + if (trail_->CurrentDecisionLevel() == 0) { + for (const Literal l : enforcement_literals) { + if (trail_->Assignment().LiteralIsFalse(l)) return true; + } } // Make sure max_variations_ is of correct size. @@ -551,12 +626,13 @@ void LinearPropagator::AddConstraint( propagation_queue_.IncreaseSize(in_queue_.size()); if (!enforcement_literals.empty()) { - infos_.back().enf_status = EnforcementStatus::CANNOT_PROPAGATE; + infos_.back().enf_status = + static_cast(EnforcementStatus::CANNOT_PROPAGATE); infos_.back().enf_id = enforcement_propagator_->Register( enforcement_literals, [this, id](EnforcementStatus status) { - infos_[id].enf_status = status; + infos_[id].enf_status = static_cast(status); // TODO(user): With some care, when we cannot propagate or the - // constraint is not enforced, we could live in_queue_[] at true but + // constraint is not enforced, we could leave in_queue_[] at true but // not put the constraint in the queue. if (status == EnforcementStatus::CAN_PROPAGATE || status == EnforcementStatus::IS_ENFORCED) { @@ -567,7 +643,7 @@ void LinearPropagator::AddConstraint( } else { AddToQueueIfNeeded(id); infos_.back().enf_id = -1; - infos_.back().enf_status = EnforcementStatus::IS_ENFORCED; + infos_.back().enf_status = static_cast(EnforcementStatus::IS_ENFORCED); } for (const IntegerVariable var : GetVariables(infos_[id])) { @@ -593,6 +669,9 @@ void LinearPropagator::AddConstraint( watcher_->WatchLowerBound(var, watcher_id_); } } + + // Propagate this new constraint. + return PropagateOneConstraint(id); } absl::Span LinearPropagator::GetCoeffs( @@ -635,14 +714,16 @@ bool LinearPropagator::PropagateOneConstraint(int id) { // Skip constraint not enforced or that cannot propagate if false. ConstraintInfo& info = infos_[id]; - if (info.enf_status == EnforcementStatus::IS_FALSE || - info.enf_status == EnforcementStatus::CANNOT_PROPAGATE) { + const EnforcementStatus enf_status = EnforcementStatus(info.enf_status); + DCHECK_EQ(enf_status, enforcement_propagator_->DebugStatus(info.enf_id)); + if (enf_status == EnforcementStatus::IS_FALSE || + enf_status == EnforcementStatus::CANNOT_PROPAGATE) { DCHECK(!in_queue_[id]); - if (info.enf_status == EnforcementStatus::IS_FALSE) { + if (enf_status == EnforcementStatus::IS_FALSE) { // We mark this constraint as in the queue but will never inspect it // again until we backtrack over this time. in_queue_[id] = true; - in_queue_and_at_false_.push_back(id); + unenforced_constraints_.push_back(id); } ++num_ignored_; return true; @@ -706,7 +787,9 @@ bool LinearPropagator::PropagateOneConstraint(int id) { } // We can only propagate more if all the enforcement literals are true. - if (info.enf_status != EnforcementStatus::IS_ENFORCED) return true; + if (info.enf_status != static_cast(EnforcementStatus::IS_ENFORCED)) { + return true; + } // The lower bound of all the variables except one can be used to update the // upper bound of the last one. diff --git a/ortools/sat/linear_propagation.h b/ortools/sat/linear_propagation.h index 64de1c755b..17d7aa6e08 100644 --- a/ortools/sat/linear_propagation.h +++ b/ortools/sat/linear_propagation.h @@ -48,8 +48,6 @@ class CustomFifoQueue { public: CustomFifoQueue() = default; - // Note that this requires the queue to be empty or to never have been poped - // before. void IncreaseSize(int n); int Pop(); @@ -83,7 +81,7 @@ class CustomFifoQueue { // An enforced constraint can be in one of these 4 states. // Note that we rely on the integer encoding to take 2 bits for optimization. -enum EnforcementStatus { +enum class EnforcementStatus { // One enforcement literal is false. IS_FALSE = 0, // More than two literals are unassigned. @@ -128,6 +126,10 @@ class EnforcementPropagator : SatPropagator { EnforcementStatus Status(EnforcementId id) const { return statuses_[id]; } + // Recompute the status from the current assignment. + // This should only used in DCHECK(). + EnforcementStatus DebugStatus(EnforcementId id); + private: absl::Span GetSpan(EnforcementId id); absl::Span GetSpan(EnforcementId id) const; @@ -181,7 +183,10 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface { void SetLevel(int level) final; // Adds a new constraint to the propagator. - void AddConstraint(absl::Span enforcement_literals, + // We support adding constraint at a positive level: + // - This will push new propagation right away. + // - This will returns false if the constraint is currently a conflict. + bool AddConstraint(absl::Span enforcement_literals, absl::Span vars, absl::Span coeffs, IntegerValue upper_bound); @@ -244,8 +249,9 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface { // IntegerVariable that changed since the last clear. SparseBitset modified_vars_; - // Per constraint info used during propagation. - std::vector infos_; + // Per constraint info used during propagation. Note that we keep pointer for + // the rev_size/rhs there, so we do need a deque. + std::deque infos_; // Buffer of the constraints data. // @@ -267,8 +273,10 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface { std::vector in_queue_; CustomFifoQueue propagation_queue_; - int rev_at_false_size_ = 0; - std::vector in_queue_and_at_false_; + // Unenforced constraints are marked as "in_queue_" but not actually added + // to the propagation_queue_. + int rev_unenforced_size_ = 0; + std::vector unenforced_constraints_; // Watchers. absl::StrongVector is_watched_; diff --git a/ortools/util/bitset.h b/ortools/util/bitset.h index 45be84d84e..5a8191fff6 100644 --- a/ortools/util/bitset.h +++ b/ortools/util/bitset.h @@ -19,10 +19,14 @@ #include #include +#include #include +#include #include +#include #include +#include "absl/log/check.h" #include "ortools/base/logging.h" #include "ortools/base/types.h" @@ -412,6 +416,8 @@ inline uint64_t TwoBitsFromPos64(uint64_t pos) { template class Bitset64 { public: + using value_type = IndexType; + // When speed matter, caching the base pointer like this to access this class // in a read only mode help. class ConstView { @@ -592,9 +598,21 @@ class Bitset64 { // // IMPORTANT: Because the iterator "caches" the current uint64_t bucket, this // will probably not do what you want if Bitset64 is modified while iterating. - class EndIterator {}; class Iterator { public: + // Make this iterator a std::forward_iterator, so it works with std::sample, + // std::max_element, etc. + Iterator() : data_(nullptr), size_(0) {} + Iterator(Iterator&& other) = default; + Iterator(const Iterator& other) = default; + Iterator& operator=(const Iterator& other) = default; + using difference_type = std::ptrdiff_t; + using iterator_category = std::forward_iterator_tag; + using value_type = IndexType; + using size_type = std::size_t; + using reference = value_type&; + using pointer = value_type*; + explicit Iterator(const Bitset64& bitset) : data_(bitset.data_.data()), size_(bitset.data_.size()) { if (!bitset.data_.empty()) { @@ -603,17 +621,34 @@ class Bitset64 { } } - bool operator!=(const EndIterator&) const { return size_ != 0; } + static Iterator EndIterator(const Bitset64& bitset) { + return Iterator(bitset.data_.data()); + } + + bool operator==(const Iterator& other) const { return !(*this == other); } + bool operator!=(const Iterator& other) const { + if (other.size_ == 0) { + return size_ != 0; + } + return std::tie(index_, current_) != + std::tie(other.index_, other.current_); + } IndexType operator*() const { return IndexType(index_); } - void operator++() { + Iterator operator++(int) { + Iterator other = *this; + (*this)++; + return other; + } + + Iterator& operator++() { int bucket = BitOffset64(index_); while (current_ == 0) { bucket++; if (bucket == size_) { size_ = 0; - return; + return *this; } current_ = data_[bucket]; } @@ -621,10 +656,13 @@ class Bitset64 { // Computes the index and clear the least significant bit of current_. index_ = BitShift64(bucket) | LeastSignificantBitPosition64(current_); current_ &= current_ - 1; + return *this; } private: - const uint64_t* const data_; + explicit Iterator(const uint64_t* data) : data_(data), size_(0) {} + + const uint64_t* data_; int size_; int index_ = 0; uint64_t current_ = 0; @@ -633,7 +671,7 @@ class Bitset64 { // Allows range-based "for" loop on the non-zero positions: // for (const IndexType index : bitset) {} Iterator begin() const { return Iterator(*this); } - EndIterator end() const { return EndIterator(); } + Iterator end() const { return Iterator::EndIterator(*this); } // Cryptic function! This is just an optimized version of a given piece of // code and has probably little general use.