From 64230c64234427cab691690ace578480162d0d90 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 16 Apr 2024 16:17:44 +0200 Subject: [PATCH] [CP-SAT] more work on 2d packing; fix bugs --- ortools/sat/2d_packing_brute_force.cc | 321 ++++++++++++++++++-------- ortools/sat/2d_packing_brute_force.h | 15 ++ ortools/sat/cp_model_expand.cc | 19 +- ortools/sat/cp_model_presolve.cc | 109 +-------- ortools/sat/cp_model_presolve.h | 6 +- ortools/sat/integer_expr.h | 27 +-- ortools/sat/presolve_context.cc | 159 +++++++++++-- ortools/sat/presolve_context.h | 12 + 8 files changed, 420 insertions(+), 248 deletions(-) diff --git a/ortools/sat/2d_packing_brute_force.cc b/ortools/sat/2d_packing_brute_force.cc index d2aa9f0196e..8f258a603e5 100644 --- a/ortools/sat/2d_packing_brute_force.cc +++ b/ortools/sat/2d_packing_brute_force.cc @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -352,48 +353,46 @@ bool BruteForceOrthogonalPackingImpl( } bool BruteForceOrthogonalPackingNoPreprocessing( - absl::Span sizes_x, - absl::Span sizes_y, - const std::pair bounding_box_size, - absl::Span result) { + const absl::Span items, + const std::pair bounding_box_size) { IntegerValue smallest_x = std::numeric_limits::max(); IntegerValue smallest_y = std::numeric_limits::max(); - const int num_items = sizes_x.size(); + const int num_items = items.size(); CHECK_LE(num_items, kMaxProblemSize); - std::vector item_index_sorted_by_area_desc(num_items); - std::array, kMaxProblemSize> - potential_item_positions_storage; - absl::Span> - potential_item_positions(potential_item_positions_storage.data(), - num_items); + IntegerValue slack = bounding_box_size.first * bounding_box_size.second; - for (int i = 0; i < num_items; ++i) { - smallest_x = std::min(smallest_x, sizes_x[i]); - smallest_y = std::min(smallest_y, sizes_y[i]); - item_index_sorted_by_area_desc[i] = i; - potential_item_positions[i].push_back({0, 0, false}); - if (sizes_x[i] > bounding_box_size.first || - sizes_y[i] > bounding_box_size.second) { + for (const PermutableItem& item : items) { + smallest_x = std::min(smallest_x, item.size_x); + smallest_y = std::min(smallest_y, item.size_y); + slack -= item.size_x * item.size_y; + if (item.size_x > bounding_box_size.first || + item.size_y > bounding_box_size.second) { return false; } } - std::sort(item_index_sorted_by_area_desc.begin(), - item_index_sorted_by_area_desc.end(), - [sizes_x, sizes_y](int a, int b) { - return sizes_x[a] * sizes_y[a] > sizes_x[b] * sizes_y[b]; + if (slack < 0) { + return false; + } + + std::sort(items.begin(), items.end(), + [](const PermutableItem& a, const PermutableItem& b) { + return std::make_tuple(a.size_x * a.size_y, a.index) > + std::make_tuple(b.size_x * b.size_y, b.index); }); + std::array new_sizes_x_storage, new_sizes_y_storage; absl::Span new_sizes_x(new_sizes_x_storage.data(), num_items); absl::Span new_sizes_y(new_sizes_y_storage.data(), num_items); - IntegerValue slack = bounding_box_size.first * bounding_box_size.second; + std::array, kMaxProblemSize> + potential_item_positions_storage; + absl::Span> + potential_item_positions(potential_item_positions_storage.data(), + num_items); for (int i = 0; i < num_items; ++i) { - new_sizes_x[i] = sizes_x[item_index_sorted_by_area_desc[i]]; - new_sizes_y[i] = sizes_y[item_index_sorted_by_area_desc[i]]; - slack -= sizes_x[i] * sizes_y[i]; - } - if (slack < 0) { - return {}; + new_sizes_x[i] = items[i].size_x; + new_sizes_y[i] = items[i].size_y; + potential_item_positions[i].push_back({0, 0, false}); } std::array item_positions_storage; absl::Span item_positions(item_positions_storage.data(), @@ -406,23 +405,169 @@ bool BruteForceOrthogonalPackingNoPreprocessing( return false; } for (int i = 0; i < num_items; ++i) { - result[item_index_sorted_by_area_desc[i]] = item_positions[i]; + items[i].position = item_positions[i]; } return true; } +// This function tries to pack a set of "tall" items with all the "shallow" +// items that can fit in the area around them. See discussion about +// the identically-feasible function v_1 in [1]. For example, for the packing: +// +// +----------------------------+ +// |------+ | +// |888888+----| | +// |888888|&&&&| | +// |------+&&&&| | +// |####| |&&&&| | +// |####| +----+ | +// |####| +------+ | +// |####| |......| | +// |####| |......|@@@ | +// |####+---------+......|@@@ | +// |####|OOOOOOOOO|......|@@@ | +// |####|OOOOOOOOO|......|@@@ | +// |####|OOOOOOOOO|......|@@@ | +// |####|OOOOOOOOO|......|@@@ | +// |####|OOOOOOOOO|......|@@@ | +// +----+---------+------+------+ +// +// We can move all "tall" and "shallow" items to the left and pack all the +// shallow items on the space remaining on the top of the tall items: +// +// +----------------------------+ +// |------+ | +// |888888+----| | +// |888888|&&&&| | +// |------+&&&&| | +// |####| |&&&&| | +// |####| +----+ | +// |####+------+ | +// |####|......| | +// |####|......| @@@ | +// |####|......+---------+@@@ | +// |####|......|OOOOOOOOO|@@@ | +// |####|......|OOOOOOOOO|@@@ | +// |####|......|OOOOOOOOO|@@@ | +// |####|......|OOOOOOOOO|@@@ | +// |####|......|OOOOOOOOO|@@@ | +// +----+------+---------+------+ +// +// If the packing is successful, we can remove both set from the problem: +// +----------------+ +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// | @@@ | +// |---------+@@@ | +// |OOOOOOOOO|@@@ | +// |OOOOOOOOO|@@@ | +// |OOOOOOOOO|@@@ | +// |OOOOOOOOO|@@@ | +// |OOOOOOOOO|@@@ | +// +---------+------+ +// +// [1] Carlier, Jacques, François Clautiaux, and Aziz Moukrim. "New reduction +// procedures and lower bounds for the two-dimensional bin packing problem with +// fixed orientation." Computers & Operations Research 34.8 (2007): 2223-2250. +// +// See Preprocess() for the API documentation. +bool PreprocessLargeWithSmallX( + absl::Span& items, + std::pair& bounding_box_size, + int max_complexity) { + // The simplest way to implement this is to sort the shallow items alongside + // their corresponding tall items. More precisely, the smallest and largest + // items are at the end of the list. The reason we put the most interesting + // values in the end even if this means we want to iterate on the list + // backward is that if the heuristic is successful we will trim them from the + // back of the list of unfixed items. + std::sort( + items.begin(), items.end(), + [&bounding_box_size](const PermutableItem& a, const PermutableItem& b) { + const bool a_is_small = 2 * a.size_x <= bounding_box_size.first; + const bool b_is_small = 2 * b.size_x <= bounding_box_size.first; + const IntegerValue a_size_for_comp = + a_is_small ? a.size_x : bounding_box_size.first - a.size_x; + const IntegerValue b_size_for_comp = + b_is_small ? b.size_x : bounding_box_size.first - b.size_x; + return std::make_tuple(a_size_for_comp, !a_is_small, a.index) > + std::make_tuple(b_size_for_comp, !b_is_small, b.index); + }); + IntegerValue total_large_item_width = 0; + IntegerValue largest_small_item_height = 0; + for (int i = items.size() - 1; i >= 1; --i) { + if (2 * items[i].size_x <= bounding_box_size.first) { + largest_small_item_height = items[i].size_x; + continue; + } + DCHECK_LE(items[i].size_x + largest_small_item_height, + bounding_box_size.first); + // We found a big item. So all values that we visited before are either + // taller than it or shallow enough to fit on top of it. Try to fit all that + // together! + if (items.subspan(i).size() > max_complexity) { + return false; + } + total_large_item_width += items[i].size_y; + if (BruteForceOrthogonalPackingNoPreprocessing( + items.subspan(i), + {bounding_box_size.first, total_large_item_width})) { + bounding_box_size.second -= total_large_item_width; + for (auto& item : items.subspan(i)) { + item.position.y_min += bounding_box_size.second; + item.position.y_max += bounding_box_size.second; + } + items = items.subspan(0, i); + return true; + } + } + + return false; +} + +// Same API as Preprocess(). +bool PreprocessLargeWithSmallY( + absl::Span& items, + std::pair& bounding_box_size, + int max_complexity) { + absl::Span orig_items = items; + for (PermutableItem& item : orig_items) { + std::swap(item.size_x, item.size_y); + std::swap(item.position.x_min, item.position.y_min); + std::swap(item.position.x_max, item.position.y_max); + } + std::swap(bounding_box_size.first, bounding_box_size.second); + const bool ret = + PreprocessLargeWithSmallX(items, bounding_box_size, max_complexity); + std::swap(bounding_box_size.first, bounding_box_size.second); + for (PermutableItem& item : orig_items) { + std::swap(item.size_x, item.size_y); + std::swap(item.position.x_min, item.position.y_min); + std::swap(item.position.x_max, item.position.y_max); + } + return ret; +} + +} // namespace + // Try to find an equivalent smaller OPP problem by fixing large items. // The API is a bit unusual: it takes a reference to a mutable Span of sizes and -// rectangles. When this function finds an item that can be fixed, it first adds -// it fixed position to `result` then reorders `sizes_x`, `sizes_y`, -// `result_index_map` and `result` to put that item in the end of the span and -// then resizes the span so it contain only non-fixed items. -bool Preprocess(absl::Span& sizes_x, - absl::Span& sizes_y, +// rectangles. When this function finds an item that can be fixed, it sets +// the position of the PermutableItem, reorders `items` to put that item in the +// end of the span and then resizes the span so it contain only non-fixed items. +// +// Note that the position of input items is not used and the position of +// non-fixed items will not be modified by this function. +bool Preprocess(absl::Span& items, std::pair& bounding_box_size, - absl::Span& result_index_map, - absl::Span& result) { - const int num_items = sizes_x.size(); + int max_complexity) { + const int num_items = items.size(); if (num_items == 1) { return false; } @@ -433,16 +578,16 @@ bool Preprocess(absl::Span& sizes_x, int largest_x_idx = -1; int largest_y_idx = -1; for (int i = 0; i < num_items; ++i) { - if (sizes_x[i] > largest_x) { - largest_x = sizes_x[i]; + if (items[i].size_x > largest_x) { + largest_x = items[i].size_x; largest_x_idx = i; } - if (sizes_y[i] > largest_y) { - largest_y = sizes_y[i]; + if (items[i].size_y > largest_y) { + largest_y = items[i].size_y; largest_y_idx = i; } - smallest_x = std::min(smallest_x, sizes_x[i]); - smallest_y = std::min(smallest_y, sizes_y[i]); + smallest_x = std::min(smallest_x, items[i].size_x); + smallest_y = std::min(smallest_y, items[i].size_y); } if (largest_x > bounding_box_size.first || @@ -450,45 +595,48 @@ bool Preprocess(absl::Span& sizes_x, // No point in optimizing obviously infeasible instance. return false; } - const auto remove_item = [&sizes_x, &sizes_y, - &result_index_map](int index_to_remove) { - std::swap(sizes_x[index_to_remove], sizes_x.back()); - sizes_x.remove_suffix(1); - std::swap(sizes_y[index_to_remove], sizes_y.back()); - sizes_y.remove_suffix(1); - std::swap(result_index_map[index_to_remove], result_index_map.back()); - result_index_map.remove_suffix(1); + const auto remove_item = [&items](int index_to_remove) { + std::swap(items[index_to_remove], items.back()); + items.remove_suffix(1); }; if (largest_x + smallest_x > bounding_box_size.first) { // No item (not even the narrowest one) fit alongside the widest item. So we // care only about fitting the remaining items in the remaining space. - bounding_box_size.second -= sizes_y[largest_x_idx]; - result.back() = { + bounding_box_size.second -= items[largest_x_idx].size_y; + items[largest_x_idx].position = { .x_min = 0, .x_max = largest_x, .y_min = bounding_box_size.second, - .y_max = bounding_box_size.second + sizes_y[largest_x_idx]}; - result.remove_suffix(1); + .y_max = bounding_box_size.second + items[largest_x_idx].size_y}; remove_item(largest_x_idx); - Preprocess(sizes_x, sizes_y, bounding_box_size, result_index_map, result); + Preprocess(items, bounding_box_size, max_complexity); return true; } if (largest_y + smallest_y > bounding_box_size.second) { - bounding_box_size.first -= sizes_x[largest_y_idx]; - result.back() = {.x_min = bounding_box_size.first, - .x_max = bounding_box_size.first + sizes_x[largest_y_idx], - .y_min = 0, - .y_max = largest_y}; - result.remove_suffix(1); + bounding_box_size.first -= items[largest_y_idx].size_x; + items[largest_y_idx].position = { + .x_min = bounding_box_size.first, + .x_max = bounding_box_size.first + items[largest_y_idx].size_x, + .y_min = 0, + .y_max = largest_y}; remove_item(largest_y_idx); - Preprocess(sizes_x, sizes_y, bounding_box_size, result_index_map, result); + Preprocess(items, bounding_box_size, max_complexity); + return true; + } + + if (PreprocessLargeWithSmallX(items, bounding_box_size, max_complexity)) { + Preprocess(items, bounding_box_size, max_complexity); return true; } + + if (PreprocessLargeWithSmallY(items, bounding_box_size, max_complexity)) { + Preprocess(items, bounding_box_size, max_complexity); + return true; + } + return false; } -} // namespace - BruteForceResult BruteForceOrthogonalPacking( absl::Span sizes_x, absl::Span sizes_y, @@ -502,40 +650,27 @@ BruteForceResult BruteForceOrthogonalPacking( return {.status = BruteForceResult::Status::kTooBig}; } CHECK_LE(num_items, kMaxProblemSize); - std::vector preprocessing_order(num_items); - std::array new_sizes_x_storage, - new_sizes_y_storage; - // We need a mutable array of sizes for preprocessing. - absl::Span new_sizes_x(new_sizes_x_storage.data(), num_items); - absl::Span new_sizes_y(new_sizes_y_storage.data(), num_items); + + std::array items_storage; + absl::Span items(items_storage.data(), num_items); for (int i = 0; i < num_items; ++i) { - preprocessing_order[i] = i; - new_sizes_x[i] = sizes_x[i]; - new_sizes_y[i] = sizes_y[i]; + items[i] = { + .size_x = sizes_x[i], .size_y = sizes_y[i], .index = i, .position = {}}; } - std::array unordered_result_storage; - absl::Span unordered_result(unordered_result_storage.data(), - num_items); - absl::Span post_processed_sizes_x = new_sizes_x; - absl::Span post_processed_sizes_y = new_sizes_y; - absl::Span result_after_preprocessing = unordered_result; + absl::Span post_processed_items = items; std::pair post_processed_bounding_box_size = bounding_box_size; - absl::Span permutation = absl::MakeSpan(preprocessing_order); const bool post_processed = - Preprocess(post_processed_sizes_x, post_processed_sizes_y, - post_processed_bounding_box_size, permutation, - result_after_preprocessing); - DCHECK_EQ(result_after_preprocessing.size(), post_processed_sizes_x.size()); - if (result_after_preprocessing.size() > max_complexity) { + Preprocess(post_processed_items, post_processed_bounding_box_size, + max_complexity - 1); + if (post_processed_items.size() > max_complexity) { return {.status = BruteForceResult::Status::kTooBig}; } const bool is_feasible = BruteForceOrthogonalPackingNoPreprocessing( - post_processed_sizes_x, post_processed_sizes_y, - post_processed_bounding_box_size, result_after_preprocessing); + post_processed_items, post_processed_bounding_box_size); VLOG_EVERY_N_SEC(2, 3) << "Solved by brute force a problem of " << num_items << " items" - << (post_processed ? absl::StrCat(" (", post_processed_sizes_x.size(), + << (post_processed ? absl::StrCat(" (", post_processed_items.size(), " after preprocessing)") : "") << ": solution " << (is_feasible ? "found" : "not found") << "."; @@ -543,8 +678,8 @@ BruteForceResult BruteForceOrthogonalPacking( return {.status = BruteForceResult::Status::kNoSolutionExists}; } std::vector result(num_items); - for (int i = 0; i < num_items; ++i) { - result[preprocessing_order[i]] = unordered_result[i]; + for (const PermutableItem& item : items) { + result[item.index] = item.position; } VLOG_EVERY_N_SEC(3, 3) << "Found a feasible packing by brute force. Dot:\n " << RenderDot(bounding_box_size, result); diff --git a/ortools/sat/2d_packing_brute_force.h b/ortools/sat/2d_packing_brute_force.h index b3c9d681480..3a8129ca476 100644 --- a/ortools/sat/2d_packing_brute_force.h +++ b/ortools/sat/2d_packing_brute_force.h @@ -48,6 +48,21 @@ BruteForceResult BruteForceOrthogonalPacking( std::pair bounding_box_size, int max_complexity); +// Note that functions taking a Span are free to permute them +// as they see fit unless documented otherwise. +struct PermutableItem { + IntegerValue size_x; + IntegerValue size_y; + // Position of the item in the original input. + int index; + Rectangle position; +}; + +// Exposed for testing +bool Preprocess(absl::Span& items, + std::pair& bounding_box_size, + int max_complexity); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index 627a29a5a10..ed013b63402 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -191,18 +191,20 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) { CapAdd(CapSub(reservoir.min_level(), demand_i), offset)); level->mutable_linear()->add_domain( CapAdd(CapSub(reservoir.max_level(), demand_i), offset)); + context->CanonicalizeLinearConstraint(level); } } else { // If all level_changes have the same sign, we do not care about the order, // just the sum. - auto* const sum = - context->working_model->add_constraints()->mutable_linear(); + ConstraintProto* new_ct = context->working_model->add_constraints(); + auto* const sum = new_ct->mutable_linear(); for (int i = 0; i < num_events; ++i) { sum->add_vars(is_active_literal(i)); sum->add_coeffs(context->FixedValue(reservoir.level_changes(i))); } sum->add_domain(reservoir.min_level()); sum->add_domain(reservoir.max_level()); + context->CanonicalizeLinearConstraint(new_ct); } ct->Clear(); @@ -464,17 +466,22 @@ void ExpandLinMax(ConstraintProto* ct, PresolveContext* context) { const int num_exprs = ct->lin_max().exprs().size(); if (num_exprs < 2) return; + // We have a special treatment for Abs, Earlyness, Tardiness, and all + // affine_max where there is only one variable present in all the expressions. + if (ExpressionsContainsOnlyOneVar(ct->lin_max().exprs())) return; + // We will create 2 * num_exprs constraints for target = max(a1, .., an). // First. // - target >= ai for (const LinearExpressionProto& expr : ct->lin_max().exprs()) { - LinearConstraintProto* lin = - context->working_model->add_constraints()->mutable_linear(); + ConstraintProto* new_ct = context->working_model->add_constraints(); + LinearConstraintProto* lin = ct->mutable_linear(); lin->add_domain(0); lin->add_domain(std::numeric_limits::max()); AddLinearExpressionToLinearConstraint(ct->lin_max().target(), 1, lin); AddLinearExpressionToLinearConstraint(expr, -1, lin); + context->CanonicalizeLinearConstraint(new_ct); } // Second, for each expr, create a new boolean bi, and add bi => target >= ai @@ -497,16 +504,16 @@ void ExpandLinMax(ConstraintProto* ct, PresolveContext* context) { for (int i = 0; i < num_exprs; ++i) { ConstraintProto* new_ct = context->working_model->add_constraints(); new_ct->add_enforcement_literal(enforcement_literals[i]); - LinearConstraintProto* lin = new_ct->mutable_linear(); lin->add_domain(std::numeric_limits::min()); lin->add_domain(0); AddLinearExpressionToLinearConstraint(ct->lin_max().target(), 1, lin); AddLinearExpressionToLinearConstraint(ct->lin_max().exprs(i), -1, lin); + context->CanonicalizeLinearConstraint(new_ct); } - ct->Clear(); context->UpdateRuleStats("lin_max: expanded lin_max"); + ct->Clear(); } // A[V] == V means for all i, V == i => A_i == i diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 03d73663309..990ddd04ae8 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2045,107 +2045,9 @@ bool CpModelPresolver::DivideLinearByGcd(ConstraintProto* ct) { return false; } -template -bool CpModelPresolver::CanonicalizeLinearExpressionInternal( - const ConstraintProto& ct, ProtoWithVarsAndCoeffs* proto, int64_t* offset) { - // First regroup the terms on the same variables and sum the fixed ones. - // - // TODO(user): Add a quick pass to skip most of the work below if the - // constraint is already in canonical form? - tmp_terms_.clear(); - int64_t sum_of_fixed_terms = 0; - bool remapped = false; - const int old_size = proto->vars().size(); - DCHECK_EQ(old_size, proto->coeffs().size()); - for (int i = 0; i < old_size; ++i) { - // Remove fixed variable and take affine representative. - // - // Note that we need to do that before we test for equality with an - // enforcement (they should already have been mapped). - int new_var; - int64_t new_coeff; - { - const int ref = proto->vars(i); - const int var = PositiveRef(ref); - const int64_t coeff = - RefIsPositive(ref) ? proto->coeffs(i) : -proto->coeffs(i); - if (coeff == 0) continue; - - if (context_->IsFixed(var)) { - sum_of_fixed_terms += coeff * context_->FixedValue(var); - continue; - } - - const AffineRelation::Relation r = context_->GetAffineRelation(var); - if (r.representative != var) { - remapped = true; - sum_of_fixed_terms += coeff * r.offset; - } - - new_var = r.representative; - new_coeff = coeff * r.coeff; - } - - // TODO(user): Avoid the quadratic loop for the corner case of many - // enforcement literal (this should be pretty rare though). - bool removed = false; - for (const int enf : ct.enforcement_literal()) { - if (new_var == PositiveRef(enf)) { - if (RefIsPositive(enf)) { - // If the constraint is enforced, we can assume the variable is at 1. - sum_of_fixed_terms += new_coeff; - } else { - // We can assume the variable is at zero. - } - removed = true; - break; - } - } - if (removed) { - context_->UpdateRuleStats("linear: enforcement literal in expression"); - continue; - } - - tmp_terms_.push_back({new_var, new_coeff}); - } - proto->clear_vars(); - proto->clear_coeffs(); - std::sort(tmp_terms_.begin(), tmp_terms_.end()); - int current_var = 0; - int64_t current_coeff = 0; - for (const auto& entry : tmp_terms_) { - CHECK(RefIsPositive(entry.first)); - if (entry.first == current_var) { - current_coeff += entry.second; - } else { - if (current_coeff != 0) { - proto->add_vars(current_var); - proto->add_coeffs(current_coeff); - } - current_var = entry.first; - current_coeff = entry.second; - } - } - if (current_coeff != 0) { - proto->add_vars(current_var); - proto->add_coeffs(current_coeff); - } - if (remapped) { - context_->UpdateRuleStats("linear: remapped using affine relations"); - } - if (proto->vars().size() < old_size) { - context_->UpdateRuleStats("linear: fixed or dup variables"); - } - *offset = sum_of_fixed_terms; - return remapped || proto->vars().size() < old_size; -} - bool CpModelPresolver::CanonicalizeLinearExpression( const ConstraintProto& ct, LinearExpressionProto* exp) { - int64_t offset = 0; - const bool result = CanonicalizeLinearExpressionInternal(ct, exp, &offset); - exp->set_offset(exp->offset() + offset); - return result; + return context_->CanonicalizeLinearExpression(ct.enforcement_literal(), exp); } bool CpModelPresolver::CanonicalizeLinear(ConstraintProto* ct) { @@ -2157,14 +2059,7 @@ bool CpModelPresolver::CanonicalizeLinear(ConstraintProto* ct) { return MarkConstraintAsFalse(ct); } - int64_t offset = 0; - bool changed = - CanonicalizeLinearExpressionInternal(*ct, ct->mutable_linear(), &offset); - if (offset != 0) { - FillDomainInProto( - ReadDomainFromProto(ct->linear()).AdditionWith(Domain(-offset)), - ct->mutable_linear()); - } + bool changed = context_->CanonicalizeLinearConstraint(ct); changed |= DivideLinearByGcd(ct); // For duplicate detection, we always make the first coeff positive. diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h index b5b5f817eaa..39c3377a6b8 100644 --- a/ortools/sat/cp_model_presolve.h +++ b/ortools/sat/cp_model_presolve.h @@ -146,12 +146,8 @@ class CpModelPresolver { // Regroups terms and substitute affine relations. // Returns true if the set of variables in the expression changed. - template - bool CanonicalizeLinearExpressionInternal(const ConstraintProto& ct, - ProtoWithVarsAndCoeffs* proto, - int64_t* offset); bool CanonicalizeLinearExpression(const ConstraintProto& ct, - LinearExpressionProto* exp); + LinearExpressionProto* proto); bool CanonicalizeLinearArgument(const ConstraintProto& ct, LinearArgumentProto* proto); diff --git a/ortools/sat/integer_expr.h b/ortools/sat/integer_expr.h index 1c5a7496bda..7f17b88e6ce 100644 --- a/ortools/sat/integer_expr.h +++ b/ortools/sat/integer_expr.h @@ -720,27 +720,16 @@ inline std::function IsEqualToMinOf( min_var = integer_trail->AddIntegerVariable(min_lb, min_ub); // min_var = min_expr - std::vector min_sum_vars = min_expr.vars; - std::vector min_sum_coeffs; - for (IntegerValue coeff : min_expr.coeffs) { - min_sum_coeffs.push_back(coeff.value()); - } - min_sum_vars.push_back(min_var); - min_sum_coeffs.push_back(-1); - - model->Add(FixedWeightedSum(min_sum_vars, min_sum_coeffs, - -min_expr.offset.value())); + LinearConstraintBuilder builder(0, 0); + builder.AddLinearExpression(min_expr, 1); + builder.AddTerm(min_var, -1); + LoadLinearConstraint(builder.Build(), model); } for (const LinearExpression& expr : exprs) { - // min_var <= expr - std::vector vars = expr.vars; - std::vector coeffs; - for (IntegerValue coeff : expr.coeffs) { - coeffs.push_back(coeff.value()); - } - vars.push_back(min_var); - coeffs.push_back(-1); - model->Add(WeightedSumGreaterOrEqual(vars, coeffs, -expr.offset.value())); + LinearConstraintBuilder builder(0, kMaxIntegerValue); + builder.AddLinearExpression(expr, 1); + builder.AddTerm(min_var, -1); + LoadLinearConstraint(builder.Build(), model); } LinMinPropagator* constraint = new LinMinPropagator(exprs, min_var, model); diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 28d6c0df0c8..e6a7ea9920a 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -2149,6 +2149,8 @@ int PresolveContext::GetOrCreateReifiedPrecedenceLiteral( (IsFixed(time_j) ? FixedValue(time_j) : time_j.offset()); lesseq->mutable_linear()->add_domain(offset); lesseq->mutable_linear()->add_domain(std::numeric_limits::max()); + CanonicalizeLinearConstraint(lesseq); + if (!LiteralIsTrue(active_i)) { AddImplication(result, active_i); } @@ -2157,25 +2159,27 @@ int PresolveContext::GetOrCreateReifiedPrecedenceLiteral( } // Not(result) && active_i && active_j => (time_i > time_j) - ConstraintProto* const greater = working_model->add_constraints(); - if (!IsFixed(time_i)) { - greater->mutable_linear()->add_vars(time_i.vars(0)); - greater->mutable_linear()->add_coeffs(-time_i.coeffs(0)); - } - if (!IsFixed(time_j)) { - greater->mutable_linear()->add_vars(time_j.vars(0)); - greater->mutable_linear()->add_coeffs(time_j.coeffs(0)); - } - greater->mutable_linear()->add_domain(std::numeric_limits::min()); - greater->mutable_linear()->add_domain(offset - 1); + { + ConstraintProto* const greater = working_model->add_constraints(); + if (!IsFixed(time_i)) { + greater->mutable_linear()->add_vars(time_i.vars(0)); + greater->mutable_linear()->add_coeffs(-time_i.coeffs(0)); + } + if (!IsFixed(time_j)) { + greater->mutable_linear()->add_vars(time_j.vars(0)); + greater->mutable_linear()->add_coeffs(time_j.coeffs(0)); + } + greater->mutable_linear()->add_domain(std::numeric_limits::min()); + greater->mutable_linear()->add_domain(offset - 1); - // Manages enforcement literal. - greater->add_enforcement_literal(NegatedRef(result)); - if (!LiteralIsTrue(active_i)) { - greater->add_enforcement_literal(active_i); - } - if (!LiteralIsTrue(active_j)) { - greater->add_enforcement_literal(active_j); + greater->add_enforcement_literal(NegatedRef(result)); + if (!LiteralIsTrue(active_i)) { + greater->add_enforcement_literal(active_i); + } + if (!LiteralIsTrue(active_j)) { + greater->add_enforcement_literal(active_j); + } + CanonicalizeLinearConstraint(greater); } // This is redundant but should improves performance. @@ -2314,5 +2318,124 @@ int PresolveContext::GetIntervalRepresentative(int index) { return index; } +template +bool CanonicalizeLinearExpressionInternal( + absl::Span enforcements, ProtoWithVarsAndCoeffs* proto, + int64_t* offset, std::vector>* tmp_terms, + PresolveContext* context) { + // First regroup the terms on the same variables and sum the fixed ones. + // + // TODO(user): Add a quick pass to skip most of the work below if the + // constraint is already in canonical form? + tmp_terms->clear(); + int64_t sum_of_fixed_terms = 0; + bool remapped = false; + const int old_size = proto->vars().size(); + DCHECK_EQ(old_size, proto->coeffs().size()); + for (int i = 0; i < old_size; ++i) { + // Remove fixed variable and take affine representative. + // + // Note that we need to do that before we test for equality with an + // enforcement (they should already have been mapped). + int new_var; + int64_t new_coeff; + { + const int ref = proto->vars(i); + const int var = PositiveRef(ref); + const int64_t coeff = + RefIsPositive(ref) ? proto->coeffs(i) : -proto->coeffs(i); + if (coeff == 0) continue; + + if (context->IsFixed(var)) { + sum_of_fixed_terms += coeff * context->FixedValue(var); + continue; + } + + const AffineRelation::Relation r = context->GetAffineRelation(var); + if (r.representative != var) { + remapped = true; + sum_of_fixed_terms += coeff * r.offset; + } + + new_var = r.representative; + new_coeff = coeff * r.coeff; + } + + // TODO(user): Avoid the quadratic loop for the corner case of many + // enforcement literal (this should be pretty rare though). + bool removed = false; + for (const int enf : enforcements) { + if (new_var == PositiveRef(enf)) { + if (RefIsPositive(enf)) { + // If the constraint is enforced, we can assume the variable is at 1. + sum_of_fixed_terms += new_coeff; + } else { + // We can assume the variable is at zero. + } + removed = true; + break; + } + } + if (removed) { + context->UpdateRuleStats("linear: enforcement literal in expression"); + continue; + } + + tmp_terms->push_back({new_var, new_coeff}); + } + proto->clear_vars(); + proto->clear_coeffs(); + std::sort(tmp_terms->begin(), tmp_terms->end()); + int current_var = 0; + int64_t current_coeff = 0; + for (const auto& entry : *tmp_terms) { + CHECK(RefIsPositive(entry.first)); + if (entry.first == current_var) { + current_coeff += entry.second; + } else { + if (current_coeff != 0) { + proto->add_vars(current_var); + proto->add_coeffs(current_coeff); + } + current_var = entry.first; + current_coeff = entry.second; + } + } + if (current_coeff != 0) { + proto->add_vars(current_var); + proto->add_coeffs(current_coeff); + } + if (remapped) { + context->UpdateRuleStats("linear: remapped using affine relations"); + } + if (proto->vars().size() < old_size) { + context->UpdateRuleStats("linear: fixed or dup variables"); + } + *offset = sum_of_fixed_terms; + return remapped || proto->vars().size() < old_size; +} + +bool PresolveContext::CanonicalizeLinearConstraint(ConstraintProto* ct) { + int64_t offset = 0; + const bool result = CanonicalizeLinearExpressionInternal( + ct->enforcement_literal(), ct->mutable_linear(), &offset, &tmp_terms_, + this); + if (offset != 0) { + FillDomainInProto( + ReadDomainFromProto(ct->linear()).AdditionWith(Domain(-offset)), + ct->mutable_linear()); + } + return result; +} + +bool PresolveContext::CanonicalizeLinearExpression( + absl::Span enforcements, LinearExpressionProto* expr) { + int64_t offset = 0; + const bool result = CanonicalizeLinearExpressionInternal( + enforcements, expr, &offset, &tmp_terms_, this); + expr->set_offset(expr->offset() + offset); + return result; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index ceb97f1e3a4..1fa50cea078 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -185,6 +185,15 @@ class PresolveContext { } } + // Canonicalization of linear constraint. This might also be needed when + // creating new constraint to make sure there are no duplicate variables. + // Returns true if the set of variables in the expression changed. + // + // This uses affine relation and regroup duplicate/fixed terms. + bool CanonicalizeLinearConstraint(ConstraintProto* ct); + bool CanonicalizeLinearExpression(absl::Span enforcements, + LinearExpressionProto* expr); + // This methods only works for affine expressions (checked). bool DomainContains(const LinearExpressionProto& expr, int64_t value) const; @@ -731,6 +740,9 @@ class PresolveContext { // Serialized proto (should be small) to index. absl::flat_hash_map interval_representative_; + // Used by CanonicalizeLinearExpressionInternal(). + std::vector> tmp_terms_; + bool model_is_expanded_ = false; };