From 823ce1188c87c8522994e2cda351f5016722762d Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 4 Jul 2023 14:11:25 +0200 Subject: [PATCH] [CP-SAT] better management of cuts; better overflow detection; do not run feasibility_jump if search is globally finished; fix #3842 --- ortools/sat/BUILD.bazel | 1 + ortools/sat/cuts.cc | 30 +- ortools/sat/feasibility_jump.cc | 3 + ortools/sat/integer.h | 44 ++- ortools/sat/integer_expr.cc | 37 +-- ortools/sat/linear_constraint_manager.cc | 28 ++ ortools/sat/linear_constraint_manager.h | 5 + ortools/sat/linear_programming_constraint.cc | 172 ---------- ortools/sat/linear_programming_constraint.h | 13 - ortools/sat/routing_cuts.cc | 315 ++++++++++++++----- ortools/sat/routing_cuts.h | 33 +- ortools/sat/scheduling_cuts.cc | 21 +- 12 files changed, 356 insertions(+), 346 deletions(-) diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index 61eebee3d4..0f10cb1673 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -1229,6 +1229,7 @@ cc_library( ":linear_constraint", ":model", "//ortools/base", + "//ortools/base:cleanup", "//ortools/base:mathutil", "//ortools/graph", "//ortools/graph:max_flow", diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 7fa9a313bf..1e395aa1c8 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -274,15 +274,14 @@ void CutDataBuilder::AddOrMergeTerm(const CutTerm& term, IntegerValue t, // // If we cannot merge the term, we will keep them separate. The produced cut // will be less strong, but can still be used. - const int64_t new_coeff = - CapAdd(cut->terms[entry_index].coeff.value(), term.coeff.value()); - const int64_t overflow_check = CapProd(t.value(), new_coeff); - if (AtMinOrMaxInt64(new_coeff) || AtMinOrMaxInt64(overflow_check)) { + const IntegerValue new_coeff = + CapAddI(cut->terms[entry_index].coeff, term.coeff); + if (AtMinOrMaxInt64I(new_coeff) || ProdOverflow(t, new_coeff)) { // If we cannot merge the term, we keep them separate. cut->terms.push_back(term); } else { ++num_merges_; - cut->terms[entry_index].coeff = IntegerValue(new_coeff); + cut->terms[entry_index].coeff = new_coeff; } } } @@ -328,22 +327,6 @@ namespace { // is needed to avoid numerical issues and adding cuts with minor effect. const double kMinCutViolation = 1e-4; -IntegerValue CapProdI(IntegerValue a, IntegerValue b) { - return IntegerValue(CapProd(a.value(), b.value())); -} - -IntegerValue CapSubI(IntegerValue a, IntegerValue b) { - return IntegerValue(CapSub(a.value(), b.value())); -} - -IntegerValue CapAddI(IntegerValue a, IntegerValue b) { - return IntegerValue(CapAdd(a.value(), b.value())); -} - -bool ProdOverflow(IntegerValue t, IntegerValue value) { - return AtMinOrMaxInt64(CapProd(t.value(), value.value())); -} - IntegerValue PositiveRemainder(absl::int128 dividend, IntegerValue positive_divisor) { DCHECK_GT(positive_divisor, 0); @@ -933,7 +916,7 @@ bool IntegerRoundingCutHelper::ComputeCut( } // Use implied bounds to "lift" Booleans into the cut. - // This should lead to stronger cuts even if the norms migth be worse. + // This should lead to stronger cuts even if the norms might be worse. num_ib_used_ = 0; if (ib_processor != nullptr) { const auto [num_lb, num_ub] = ib_processor->PostprocessWithImpliedBound( @@ -2643,8 +2626,7 @@ bool BuildMaxAffineUpConstraint( // -delta_y * var + delta_x * target <= delta_x * y_at_min - delta_y * x_min // // Checks the rhs for overflows. - if (AtMinOrMaxInt64(CapProd(delta_x.value(), y_at_min.value())) || - AtMinOrMaxInt64(CapProd(delta_y.value(), x_min.value()))) { + if (ProdOverflow(delta_x, y_at_min) || ProdOverflow(delta_y, x_min)) { return false; } diff --git a/ortools/sat/feasibility_jump.cc b/ortools/sat/feasibility_jump.cc index fe124032fc..353efbfe8b 100644 --- a/ortools/sat/feasibility_jump.cc +++ b/ortools/sat/feasibility_jump.cc @@ -348,6 +348,9 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { if (shared_bounds_ != nullptr) { shared_bounds_->UpdateDomains(&var_domains_); for (int var = 0; var < var_domains_.size(); ++var) { + // We abort if the problem is trivially UNSAT. This might happen while + // we are cleaning up all workers at the end of a search. + if (var_domains_[var].IsEmpty()) return; var_has_two_values_[var] = var_domains_[var].HasTwoValues(); } } diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index 82e0ef2a49..a1acdec320 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -105,6 +105,28 @@ inline IntegerValue FloorRatio(IntegerValue dividend, return result - adjust; } +// Overflows and saturated arithmetic. + +inline IntegerValue CapProdI(IntegerValue a, IntegerValue b) { + return IntegerValue(CapProd(a.value(), b.value())); +} + +inline IntegerValue CapSubI(IntegerValue a, IntegerValue b) { + return IntegerValue(CapSub(a.value(), b.value())); +} + +inline IntegerValue CapAddI(IntegerValue a, IntegerValue b) { + return IntegerValue(CapAdd(a.value(), b.value())); +} + +inline bool ProdOverflow(IntegerValue t, IntegerValue value) { + return AtMinOrMaxInt64(CapProd(t.value(), value.value())); +} + +inline bool AtMinOrMaxInt64I(IntegerValue t) { + return AtMinOrMaxInt64(t.value()); +} + // Returns dividend - FloorRatio(dividend, divisor) * divisor; // // This function is around the same speed than the computation above, but it @@ -118,17 +140,21 @@ inline IntegerValue PositiveRemainder(IntegerValue dividend, return m < 0 ? m + positive_divisor : m; } +inline bool AddTo(IntegerValue a, IntegerValue* result) { + if (AtMinOrMaxInt64I(a)) return false; + const IntegerValue add = CapAddI(a, *result); + if (AtMinOrMaxInt64I(add)) return false; + *result = add; + return true; +} + // Computes result += a * b, and return false iff there is an overflow. inline bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue* result) { - const int64_t prod = CapProd(a.value(), b.value()); - if (prod == std::numeric_limits::min() || - prod == std::numeric_limits::max()) - return false; - const int64_t add = CapAdd(prod, result->value()); - if (add == std::numeric_limits::min() || - add == std::numeric_limits::max()) - return false; - *result = IntegerValue(add); + const IntegerValue prod = CapProdI(a, b); + if (AtMinOrMaxInt64I(prod)) return false; + const IntegerValue add = CapAddI(prod, *result); + if (AtMinOrMaxInt64I(add)) return false; + *result = add; return true; } diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index 9aa1e3234a..0ac02b7028 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -851,7 +851,7 @@ bool ProductPropagator::PropagateWhenAllNonNegative() { { const IntegerValue max_a = integer_trail_->UpperBound(a_); const IntegerValue max_b = integer_trail_->UpperBound(b_); - const IntegerValue new_max(CapProd(max_a.value(), max_b.value())); + const IntegerValue new_max = CapProdI(max_a, max_b); if (new_max < integer_trail_->UpperBound(p_)) { if (!integer_trail_->SafeEnqueue( p_.LowerOrEqual(new_max), @@ -866,7 +866,7 @@ bool ProductPropagator::PropagateWhenAllNonNegative() { { const IntegerValue min_a = integer_trail_->LowerBound(a_); const IntegerValue min_b = integer_trail_->LowerBound(b_); - const IntegerValue new_min(CapProd(min_a.value(), min_b.value())); + const IntegerValue new_min = CapProdI(min_a, min_b); // The conflict test is needed because when new_min is large, we could // have an overflow in p_.GreaterOrEqual(new_min); @@ -893,7 +893,7 @@ bool ProductPropagator::PropagateWhenAllNonNegative() { const IntegerValue min_b = integer_trail_->LowerBound(b); const IntegerValue min_p = integer_trail_->LowerBound(p_); const IntegerValue max_p = integer_trail_->UpperBound(p_); - const IntegerValue prod(CapProd(max_a.value(), min_b.value())); + const IntegerValue prod = CapProdI(max_a, min_b); if (prod > max_p) { if (!integer_trail_->SafeEnqueue(a.LowerOrEqual(FloorRatio(max_p, min_b)), {integer_trail_->LowerBoundAsLiteral(b), @@ -980,12 +980,12 @@ bool ProductPropagator::Propagate() { // // TODO(user): In the reasons, including all 4 bounds is always correct, but // we might be able to relax some of them. - const int64_t max_a = integer_trail_->UpperBound(a_).value(); - const int64_t max_b = integer_trail_->UpperBound(b_).value(); - const IntegerValue p1(CapProd(max_a, max_b)); - const IntegerValue p2(CapProd(max_a, min_b)); - const IntegerValue p3(CapProd(min_a, max_b)); - const IntegerValue p4(CapProd(min_a, min_b)); + const IntegerValue max_a = integer_trail_->UpperBound(a_); + const IntegerValue max_b = integer_trail_->UpperBound(b_); + const IntegerValue p1 = CapProdI(max_a, max_b); + const IntegerValue p2 = CapProdI(max_a, min_b); + const IntegerValue p3 = CapProdI(min_a, max_b); + const IntegerValue p4 = CapProdI(min_a, min_b); const IntegerValue new_max_p = std::max({p1, p2, p3, p4}); if (new_max_p < integer_trail_->UpperBound(p_)) { if (!integer_trail_->SafeEnqueue( @@ -1124,7 +1124,7 @@ SquarePropagator::SquarePropagator(AffineExpression x, AffineExpression s, bool SquarePropagator::Propagate() { const IntegerValue min_x = integer_trail_->LowerBound(x_); const IntegerValue min_s = integer_trail_->LowerBound(s_); - const IntegerValue min_x_square(CapProd(min_x.value(), min_x.value())); + const IntegerValue min_x_square = CapProdI(min_x, min_x); if (min_x_square > min_s) { if (!integer_trail_->SafeEnqueue(s_.GreaterOrEqual(min_x_square), {x_.GreaterOrEqual(min_x)})) { @@ -1141,7 +1141,7 @@ bool SquarePropagator::Propagate() { const IntegerValue max_x = integer_trail_->UpperBound(x_); const IntegerValue max_s = integer_trail_->UpperBound(s_); - const IntegerValue max_x_square(CapProd(max_x.value(), max_x.value())); + const IntegerValue max_x_square = CapProdI(max_x, max_x); if (max_x_square < max_s) { if (!integer_trail_->SafeEnqueue(s_.LowerOrEqual(max_x_square), {x_.LowerOrEqual(max_x)})) { @@ -1151,9 +1151,7 @@ bool SquarePropagator::Propagate() { const IntegerValue new_max(FloorSquareRoot(max_s.value())); if (!integer_trail_->SafeEnqueue( x_.LowerOrEqual(new_max), - {s_.LowerOrEqual(IntegerValue(CapProd(new_max.value() + 1, - new_max.value() + 1)) - - 1)})) { + {s_.LowerOrEqual(CapProdI(new_max + 1, new_max + 1) - 1)})) { return false; } } @@ -1273,7 +1271,7 @@ bool DivisionPropagator::PropagateUpperBounds(AffineExpression num, // num < (max_div + 1) * denom // num + 1 <= (max_div + 1) * max_denom. const IntegerValue new_max_num = - IntegerValue(CapAdd(CapProd(max_div.value() + 1, max_denom.value()), -1)); + CapAddI(CapProdI(max_div + 1, max_denom), -1); if (max_num > new_max_num) { if (!integer_trail_->SafeEnqueue( num.LowerOrEqual(new_max_num), @@ -1309,8 +1307,7 @@ bool DivisionPropagator::PropagatePositiveDomains(AffineExpression num, // We start from num / denom >= min_div. // num >= min_div * denom. // num >= min_div * min_denom. - const IntegerValue new_min_num = - IntegerValue(CapProd(min_denom.value(), min_div.value())); + const IntegerValue new_min_num = CapProdI(min_denom, min_div); if (min_num < new_min_num) { if (!integer_trail_->SafeEnqueue( num.GreaterOrEqual(new_min_num), @@ -1382,8 +1379,7 @@ bool FixedDivisionPropagator::Propagate() { } } else if (max_a / b_ > max_c) { const IntegerValue new_max_a = - max_c >= 0 ? max_c * b_ + b_ - 1 - : IntegerValue(CapProd(max_c.value(), b_.value())); + max_c >= 0 ? max_c * b_ + b_ - 1 : CapProdI(max_c, b_); CHECK_LT(new_max_a, max_a); if (!integer_trail_->SafeEnqueue( a_.LowerOrEqual(new_max_a), @@ -1401,8 +1397,7 @@ bool FixedDivisionPropagator::Propagate() { } } else if (min_a / b_ < min_c) { const IntegerValue new_min_a = - min_c > 0 ? IntegerValue(CapProd(min_c.value(), b_.value())) - : min_c * b_ - b_ + 1; + min_c > 0 ? CapProdI(min_c, b_) : min_c * b_ - b_ + 1; CHECK_GT(new_min_a, min_a); if (!integer_trail_->SafeEnqueue( a_.GreaterOrEqual(new_min_a), diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 83479a9e59..7f8bce7408 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -46,6 +46,28 @@ namespace operations_research { namespace sat { +bool PossibleOverflow(const IntegerTrail& integer_trail, + const LinearConstraint& constraint) { + IntegerValue min_activity(0); + IntegerValue max_activity(0); + const int size = constraint.vars.size(); + for (int i = 0; i < size; ++i) { + const IntegerVariable var = constraint.vars[i]; + const IntegerValue coeff = constraint.coeffs[i]; + CHECK_NE(coeff, 0); + const IntegerValue lb = integer_trail.LevelZeroLowerBound(var); + const IntegerValue ub = integer_trail.LevelZeroUpperBound(var); + if (coeff > 0) { + if (!AddProductTo(lb, coeff, &min_activity)) return true; + if (!AddProductTo(ub, coeff, &max_activity)) return true; + } else { + if (!AddProductTo(ub, coeff, &min_activity)) return true; + if (!AddProductTo(lb, coeff, &max_activity)) return true; + } + } + return AtMinOrMaxInt64(CapSub(max_activity.value(), min_activity.value())); +} + namespace { const LinearConstraintManager::ConstraintIndex kInvalidConstraintIndex(-1); @@ -290,6 +312,12 @@ bool LinearConstraintManager::AddCut(const LinearConstraint& ct, return false; } + // TODO(user): We could prevent overflow by dividing more. Note that mainly + // happen with super large variable domain since we usually restrict the size + // of the generated coefficients in our cuts. So it shouldn't be that + // important. + if (PossibleOverflow(integer_trail_, ct)) return false; + bool added = false; const ConstraintIndex ct_index = Add(ct, &added); diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h index d50260731c..819a04c528 100644 --- a/ortools/sat/linear_constraint_manager.h +++ b/ortools/sat/linear_constraint_manager.h @@ -41,6 +41,11 @@ namespace operations_research { namespace sat { +// Tests for possible overflow in the simplification of the given linear +// constraint. +bool PossibleOverflow(const IntegerTrail& integer_trail, + const LinearConstraint& constraint); + // Stores for each IntegerVariable its temporary LP solution. // // This is shared between all LinearProgrammingConstraint because in the corner diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 7f4fca44e1..f58672813c 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -1888,178 +1888,6 @@ absl::int128 LinearProgrammingConstraint::GetImpliedLowerBound( return lower_bound; } -bool PossibleOverflow(const IntegerTrail& integer_trail, - const LinearConstraint& constraint) { - IntegerValue lower_bound(0); - const int size = constraint.vars.size(); - for (int i = 0; i < size; ++i) { - const IntegerVariable var = constraint.vars[i]; - const IntegerValue coeff = constraint.coeffs[i]; - CHECK_NE(coeff, 0); - const IntegerValue bound = coeff > 0 - ? integer_trail.LevelZeroLowerBound(var) - : integer_trail.LevelZeroUpperBound(var); - if (!AddProductTo(bound, coeff, &lower_bound)) { - return true; - } - } - const int64_t slack = CapSub(constraint.ub.value(), lower_bound.value()); - return slack == std::numeric_limits::min() || - slack == std::numeric_limits::max(); -} - -namespace { - -absl::int128 FloorRatio128(absl::int128 x, IntegerValue positive_div) { - absl::int128 div128(positive_div.value()); - absl::int128 result = x / div128; - if (result * div128 > x) return result - 1; - return result; -} - -absl::int128 CeilRatio128(absl::int128 x, absl::int128 div128) { - absl::int128 result = x / div128; - if (result * div128 < x) return result + 1; - return result; -} - -// TODO(user): This code is tricky and similar to the one to generate cuts. -// Maybe reduce the duplication? note however that here we use int128 to deal -// with potential overflow. -void DivideConstraint(const IntegerTrail& integer_trail, IntegerValue divisor, - LinearConstraint* constraint) { - // To be correct, we need to shift all variable so that they are positive. - // - // Important: One might be tempted to think that using the current variable - // bounds is okay here since we only use this to derive cut/constraint that - // only needs to be locally valid. However, in some corner cases (like when - // one term become zero), we might loose the fact that we used one of the - // variable bound to derive the new constraint, so we will miss it in the - // explanation !! - int new_size = 0; - absl::int128 adjust = 0; - const int size = constraint->vars.size(); - for (int i = 0; i < size; ++i) { - const IntegerValue old_coeff = constraint->coeffs[i]; - const IntegerValue new_coeff = FloorRatio(old_coeff, divisor); - - // Compute the rhs adjustement. - const absl::int128 remainder = - absl::int128(old_coeff.value()) - - absl::int128(new_coeff.value()) * absl::int128(divisor.value()); - adjust += - remainder * - absl::int128( - integer_trail.LevelZeroLowerBound(constraint->vars[i]).value()); - - if (new_coeff == 0) continue; - constraint->vars[new_size] = constraint->vars[i]; - constraint->coeffs[new_size] = new_coeff; - ++new_size; - } - constraint->vars.resize(new_size); - constraint->coeffs.resize(new_size); - - // TODO(user): I am not 100% sure this cannot overflow. If it does it means - // our reduced constraint is trivial though, and we can cap it. - constraint->ub = IntegerValue(static_cast( - FloorRatio128(absl::int128(constraint->ub.value()) - adjust, divisor))); -} - -} // namespace - -// The goal here is to prevent overflow in the IntegerSumLE propagation code. -// We want to be as tight as possible. We do it in two steps, which is not ideal -// but easier. -bool PreventOverflow(const IntegerTrail& integer_trail, - LinearConstraint* constraint) { - bool prevented = false; - - // We use kint64max - 1 so that PossibleOverflow() can distinguish overflow - // for a sum exactly equal to kint64max. - const absl::int128 threshold(std::numeric_limits::max() - 1); - - // First, make all coefficient positive. - MakeAllCoefficientsPositive(constraint); - - // First step is to make sure coeff * (ub - lb) and coeff * lb will not - // overflow. Note that we already know (ub - lb) cannot overflow. - { - absl::int128 max_delta = 0; - const int size = constraint->vars.size(); - for (int i = 0; i < size; ++i) { - const IntegerVariable var = constraint->vars[i]; - const IntegerValue lb = integer_trail.LevelZeroLowerBound(var); - const IntegerValue ub = integer_trail.LevelZeroUpperBound(var); - const absl::int128 coeff(constraint->coeffs[i].value()); - const absl::int128 diff( - std::max({IntTypeAbs(lb), IntTypeAbs(ub), ub - lb}).value()); - max_delta = std::max(max_delta, coeff * diff); - } - if (max_delta > threshold) { - prevented = true; - const IntegerValue divisor( - static_cast(CeilRatio128(max_delta, threshold))); - DivideConstraint(integer_trail, divisor, constraint); - } - } - - // Second step is to make sure computing the lower bound will not overflow - // whatever the order and at whatever level. And also that computing the slack - // will not overflow. - // - // Note that because each term fit on an int64_t per first step, we will not - // have int128 overflow. - // - // TODO(user): We could change the propag to detect a conflict without - // computing the full activity, and thus avoid some overflow. Like precompute - // a base lb and then compute the activity from there? Or we could have - // a custom code here, actually we only propagate this with the current lb - // instead of the LevelZeroUpperBound(). Or we could just propagate using - // int128 arithmetic. - { - absl::int128 sum_min_neg = 0; - absl::int128 sum_min_pos = 0; - absl::int128 sum_max_neg = 0; - absl::int128 sum_max_pos = 0; - const int size = constraint->vars.size(); - for (int i = 0; i < size; ++i) { - const IntegerVariable var = constraint->vars[i]; - const absl::int128 coeff(constraint->coeffs[i].value()); - const absl::int128 lb(integer_trail.LevelZeroLowerBound(var).value()); - if (lb > 0) { - sum_min_pos += coeff * lb; - } else { - sum_min_neg += coeff * lb; - } - const absl::int128 ub(integer_trail.LevelZeroUpperBound(var).value()); - if (ub > 0) { - sum_max_pos += coeff * ub; - } else { - sum_max_neg += coeff * ub; - } - } - const absl::int128 min_slack = - static_cast(constraint->ub.value()) - - (sum_min_pos + sum_min_neg); - const absl::int128 max_slack = - static_cast(constraint->ub.value()) - - (sum_max_pos + sum_max_neg); - const absl::int128 max_value = - std::max({-sum_min_neg, sum_min_pos, sum_min_pos + sum_min_neg, - -sum_max_neg, sum_max_pos, sum_max_pos + sum_max_neg, - min_slack, -min_slack, max_slack, -max_slack}); - if (max_value > threshold) { - prevented = true; - const IntegerValue divisor( - static_cast(CeilRatio128(max_value, threshold))); - DivideConstraint(integer_trail, divisor, constraint); - } - } - - return prevented; -} - // TODO(user): combine this with RelaxLinearReason() to avoid the extra // magnitude vector and the weird precondition of RelaxLinearReason(). void LinearProgrammingConstraint::SetImpliedLowerBoundReason( diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 258c36042b..6f69844fe2 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -610,19 +610,6 @@ class LinearProgrammingConstraintCollection } }; -// Tests for possible overflow in the propagation of the given linear -// constraint. -bool PossibleOverflow(const IntegerTrail& integer_trail, - const LinearConstraint& constraint); - -// Reduces the coefficient of the constraint so that we cannot have overflow -// in the propagation of the given linear constraint. Note that we may loose -// some strength by doing so. -// -// Returns true if we changed the constraint. -bool PreventOverflow(const IntegerTrail& integer_trail, - LinearConstraint* constraint); - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/routing_cuts.cc b/ortools/sat/routing_cuts.cc index 9b7e3c34d2..aec4c48b64 100644 --- a/ortools/sat/routing_cuts.cc +++ b/ortools/sat/routing_cuts.cc @@ -19,12 +19,14 @@ #include #include #include +#include #include #include #include "absl/container/inlined_vector.h" #include "absl/log/check.h" #include "absl/types/span.h" +#include "ortools/base/cleanup.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/base/strong_vector.h" @@ -74,6 +76,21 @@ class OutgoingCutHelper { // Try to add an outgoing cut from the given subset. bool TrySubsetCut(std::string name, absl::Span subset); + // If we look at the symmetrized version (tail <-> head = tail->head + + // head->tail) and we split all the edges between a subset of nodes S and the + // outside into a set A and the other d(S)\A, and |A| is odd, we have a + // constraint of the form: + // "all edge of A at 1" => sum other edges >= 1. + // This is because a cycle or multiple-cycle must go in/out an even number + // of time. This enforced constraint simply linearize to: + // sum_d(S)\A x_e + sum_A (1 - x_e) >= 1. + // + // Given a subset of nodes, it is easy to identify the best subset A of edge + // to consider. + bool TryBlossomSubsetCut(std::string name, + const std::vector& symmetrized_edges, + absl::Span subset); + private: // Add a cut of the form Sum_{outgoing arcs from S} lp >= rhs_lower_bound. // @@ -252,16 +269,109 @@ bool OutgoingCutHelper::TrySubsetCut(std::string name, return result; } +bool OutgoingCutHelper::TryBlossomSubsetCut( + std::string name, const std::vector& symmetrized_edges, + absl::Span subset) { + DCHECK_GE(subset.size(), 1); + DCHECK_LT(subset.size(), num_nodes_); + + // Initialize "in_subset" and the subset demands. + for (const int n : subset) in_subset_[n] = true; + auto cleanup = ::absl::MakeCleanup([subset, this]() { + for (const int n : subset) in_subset_[n] = false; + }); + + // The heuristic assumes non-duplicates arcs, otherwise they are all bundled + // together in the same symmetric edge, and the result is probably wrong. + absl::flat_hash_set> special_edges; + int num_inverted = 0; + double sum_inverted = 0.0; + double sum = 0.0; + double best_change = 1.0; + ArcWithLpValue const* best_swap = nullptr; + for (const ArcWithLpValue& arc : symmetrized_edges) { + if (in_subset_[arc.tail] == in_subset_[arc.head]) continue; + + if (arc.lp_value > 0.5) { + ++num_inverted; + special_edges.insert({arc.tail, arc.head}); + sum_inverted += 1.0 - arc.lp_value; + } else { + sum += arc.lp_value; + } + + const double change = std::abs(2 * arc.lp_value - 1.0); + if (change < best_change) { + best_change = change; + best_swap = &arc; + } + } + + // If the we don't have an odd number, we move the best edge from one set to + // the other. + if (num_inverted % 2 == 0) { + if (best_swap == nullptr) return false; + if (special_edges.contains({best_swap->tail, best_swap->head})) { + --num_inverted; + special_edges.erase({best_swap->tail, best_swap->head}); + sum_inverted -= (1.0 - best_swap->lp_value); + sum += best_swap->lp_value; + } else { + ++num_inverted; + special_edges.insert({best_swap->tail, best_swap->head}); + sum_inverted += (1.0 - best_swap->lp_value); + sum -= best_swap->lp_value; + } + } + if (sum + sum_inverted > 0.99) return false; + + // For the route constraint, it is actually allowed to have circuit of size + // 2, so the reasoning is wrong if one of the edges touches the depot. + if (!demands_.empty()) { + for (const auto [tail, head] : special_edges) { + if (tail == 0) return false; + } + } + + // Try to generate the cut. + // + // We deal with the corner case with duplicate arcs, or just one side of a + // "symmetric" edge present. + int num_actual_inverted = 0; + absl::flat_hash_set> processed_arcs; + LinearConstraintBuilder builder(encoder_, IntegerValue(1), kMaxIntegerValue); + for (int i = 0; i < tails_.size(); ++i) { + if (tails_[i] == heads_[i]) continue; + if (in_subset_[tails_[i]] == in_subset_[heads_[i]]) continue; + + const std::pair key = {tails_[i], heads_[i]}; + const std::pair r_key = {heads_[i], tails_[i]}; + const std::pair s_key = std::min(key, r_key); + if (special_edges.contains(s_key) && !processed_arcs.contains(key)) { + processed_arcs.insert(key); + CHECK(builder.AddLiteralTerm(literals_[i], IntegerValue(-1))); + if (!processed_arcs.contains(r_key)) { + ++num_actual_inverted; + } + continue; + } + + // Normal edge. + CHECK(builder.AddLiteralTerm(literals_[i], IntegerValue(1))); + } + builder.AddConstant(IntegerValue(num_actual_inverted)); + if (num_actual_inverted % 2 == 0) return false; + + return manager_->AddCut(builder.Build(), name); +} + } // namespace void GenerateInterestingSubsets(int num_nodes, const std::vector>& arcs, - int min_subset_size, int stop_at_num_components, + int stop_at_num_components, std::vector* subset_data, std::vector>* subsets) { - subset_data->resize(num_nodes); - subsets->clear(); - // We will do a union-find by adding one by one the arc of the lp solution // in the order above. Every intermediate set during this construction will // be a candidate for a cut. @@ -313,71 +423,28 @@ void GenerateInterestingSubsets(int num_nodes, // a consecutive span in the subset_data vector. This vector just lists the // nodes in the "pre-order" graph traversal order. The Spans will point inside // the subset_data vector, it is why we initialize it once and for all. - int new_size = 0; - { - std::vector> graph(parent.size()); - for (int i = 0; i < parent.size(); ++i) { - if (parent[i] != i) graph[parent[i]].push_back(i); - } - std::vector queue; - std::vector seen(graph.size(), false); - std::vector start_index(parent.size()); - for (int i = 0; i < parent.size(); ++i) { - // Note that because of the way we constructed 'parent', the graph is a - // binary tree. This is not required for the correctness of the algorithm - // here though. - CHECK(graph[i].empty() || graph[i].size() == 2); - if (parent[i] != i) continue; - - // Explore the subtree rooted at node i. - CHECK(!seen[i]); - queue.push_back(i); - while (!queue.empty()) { - const int node = queue.back(); - if (seen[node]) { - queue.pop_back(); - // All the children of node are in the span [start, end) of the - // subset_data vector. - const int start = start_index[node]; - if (new_size - start >= min_subset_size) { - subsets->emplace_back(&(*subset_data)[start], new_size - start); - } - continue; - } - seen[node] = true; - start_index[node] = new_size; - if (node < num_nodes) (*subset_data)[new_size++] = node; - for (const int child : graph[node]) { - if (!seen[child]) queue.push_back(child); - } - } - } - } - - DCHECK_EQ(new_size, num_nodes); + ExtractAllSubsetsFromForest(parent, subset_data, subsets, + /*node_limit=*/num_nodes); } -void ExtractAllSubsetsFromTree(const std::vector& parent, - std::vector* subset_data, - std::vector>* subsets) { +void ExtractAllSubsetsFromForest(const std::vector& parent, + std::vector* subset_data, + std::vector>* subsets, + int node_limit) { // To not reallocate memory since we need the span to point inside this // vector, we resize subset_data right away. int out_index = 0; const int num_nodes = parent.size(); - subset_data->resize(num_nodes); + subset_data->resize(std::min(num_nodes, node_limit)); subsets->clear(); // Starts by creating the corresponding graph and find the root. - int root = -1; util::StaticGraph graph(num_nodes, num_nodes - 1); for (int i = 0; i < num_nodes; ++i) { - if (parent[i] == i) { - root = i; - } else { + if (parent[i] != i) { graph.AddArc(parent[i], i); } } - if (root == -1) return; graph.Build(); // Perform a dfs on the rooted tree. @@ -385,24 +452,30 @@ void ExtractAllSubsetsFromTree(const std::vector& parent, std::vector subtree_starts(num_nodes, -1); std::vector stack; stack.reserve(num_nodes); - stack.push_back(root); - while (!stack.empty()) { - const int node = stack.back(); + for (int i = 0; i < parent.size(); ++i) { + if (parent[i] != i) continue; - // The node was already explored, output its subtree and pop it. - if (subtree_starts[node] >= 0) { - stack.pop_back(); - (*subset_data)[out_index++] = node; - const int start = subtree_starts[node]; - const int size = out_index - start; - subsets->push_back(absl::MakeSpan(&(*subset_data)[start], size)); - continue; - } + stack.push_back(i); // root. + while (!stack.empty()) { + const int node = stack.back(); - // Explore. - subtree_starts[node] = out_index; - for (const int child : graph[node]) { - stack.push_back(child); + // The node was already explored, output its subtree and pop it. + if (subtree_starts[node] >= 0) { + stack.pop_back(); + if (node < node_limit) { + (*subset_data)[out_index++] = node; + } + const int start = subtree_starts[node]; + const int size = out_index - start; + subsets->push_back(absl::MakeSpan(&(*subset_data)[start], size)); + continue; + } + + // Explore. + subtree_starts[node] = out_index; + for (const int child : graph[node]) { + stack.push_back(child); + } } } } @@ -439,6 +512,31 @@ std::vector ComputeGomoryHuTree( return parent; } +void SymmetrizeArcs(std::vector* arcs) { + for (ArcWithLpValue& arc : *arcs) { + if (arc.tail <= arc.head) continue; + std::swap(arc.tail, arc.head); + } + std::sort(arcs->begin(), arcs->end(), + [](const ArcWithLpValue& a, const ArcWithLpValue& b) { + return std::tie(a.tail, a.head) < std::tie(b.tail, b.head); + }); + + int new_size = 0; + int last_tail = -1; + int last_head = -1; + for (const ArcWithLpValue& arc : *arcs) { + if (arc.tail == last_tail && arc.head == last_head) { + (*arcs)[new_size - 1].lp_value += arc.lp_value; + continue; + } + (*arcs)[new_size++] = arc; + last_tail = arc.tail; + last_head = arc.head; + } + arcs->resize(new_size); +} + // We roughly follow the algorithm described in section 6 of "The Traveling // Salesman Problem, A computational Study", David L. Applegate, Robert E. // Bixby, Vasek Chvatal, William J. Cook. @@ -486,7 +584,6 @@ void SeparateSubtourInequalities( std::vector subset_data; std::vector> subsets; GenerateInterestingSubsets(num_nodes, ordered_arcs, - /*min_subset_size=*/2, /*stop_at_num_components=*/2, &subset_data, &subsets); @@ -500,10 +597,23 @@ void SeparateSubtourInequalities( OutgoingCutHelper helper(num_nodes, capacity, demands, tails, heads, literals, literal_lp_values, relevant_arcs, manager, model); + // Hack/optim: we exploit the tree structure of the subsets to not add a cut + // for a larger subset if we added a cut from one included in it. + // + // TODO(user): Currently if we add too many not so relevant cuts, our generic + // MIP cut heuritic are way too slow on TSP/VRP problems. + int last_added_start = -1; + // Process each subsets and add any violated cut. int num_added = 0; for (const absl::Span subset : subsets) { - if (helper.TrySubsetCut("Circuit", subset)) ++num_added; + if (subset.size() <= 1) continue; + const int start = static_cast(subset.data() - subset_data.data()); + if (start <= last_added_start) continue; + if (helper.TrySubsetCut("Circuit", subset)) { + ++num_added; + last_added_start = start; + } } // If there were no cut added by the heuristic above, we try exact separation. @@ -516,19 +626,53 @@ void SeparateSubtourInequalities( // Note(user): Compared to any min-cut, these cut have some nice properties // since they are "included" in each other. This might help with combining // them within our generic IP cuts framework. - if (num_added == 0) { - // TODO(user): I had an older version that tried the n-cuts generated during - // the course of the algorithm. This could also be interesting. But it is - // hard to tell with our current benchmark setup. - const std::vector parent = - ComputeGomoryHuTree(num_nodes, relevant_arcs); + // + // TODO(user): I had an older version that tried the n-cuts generated during + // the course of the algorithm. This could also be interesting. But it is + // hard to tell with our current benchmark setup. + if (num_added != 0) return; + SymmetrizeArcs(&relevant_arcs); + std::vector parent = ComputeGomoryHuTree(num_nodes, relevant_arcs); - // Try all interesting subset from the Gomory-Hu tree. - ExtractAllSubsetsFromTree(parent, &subset_data, &subsets); - for (const absl::Span subset : subsets) { - if (subset.size() == 1) continue; - if (subset.size() == num_nodes) continue; - helper.TrySubsetCut("CircuitExact", subset); + // Try all interesting subset from the Gomory-Hu tree. + ExtractAllSubsetsFromForest(parent, &subset_data, &subsets); + last_added_start = -1; + for (const absl::Span subset : subsets) { + if (subset.size() <= 1) continue; + if (subset.size() == num_nodes) continue; + const int start = static_cast(subset.data() - subset_data.data()); + if (start <= last_added_start) continue; + if (helper.TrySubsetCut("CircuitExact", subset)) { + ++num_added; + last_added_start = start; + } + } + + // Exact separation of symmetric Blossom cut. We use the algorithm in the + // paper: "A Faster Exact Separation Algorithm for Blossom Inequalities", Adam + // N. Letchford, Gerhard Reinelt, Dirk Oliver Theis, 2004. + // + // Note that the "relevant_arcs" were symmetrized above. + if (num_added != 0) return; + if (num_nodes <= 2) return; + std::vector for_blossom; + for_blossom.reserve(relevant_arcs.size()); + for (ArcWithLpValue arc : relevant_arcs) { + if (arc.lp_value > 0.5) arc.lp_value = 1.0 - arc.lp_value; + if (arc.lp_value < 1e-6) continue; + for_blossom.push_back(arc); + } + parent = ComputeGomoryHuTree(num_nodes, for_blossom); + ExtractAllSubsetsFromForest(parent, &subset_data, &subsets); + last_added_start = -1; + for (const absl::Span subset : subsets) { + if (subset.size() <= 1) continue; + if (subset.size() == num_nodes) continue; + const int start = static_cast(subset.data() - subset_data.data()); + if (start <= last_added_start) continue; + if (helper.TryBlossomSubsetCut("CircuitBlossom", relevant_arcs, subset)) { + ++num_added; + last_added_start = start; } } } @@ -662,7 +806,6 @@ void SeparateFlowInequalities( std::vector subset_data; std::vector> subsets; GenerateInterestingSubsets(num_nodes, ordered_arcs, - /*min_subset_size=*/1, /*stop_at_num_components=*/1, &subset_data, &subsets); diff --git a/ortools/sat/routing_cuts.h b/ortools/sat/routing_cuts.h index af603650e5..be724332c4 100644 --- a/ortools/sat/routing_cuts.h +++ b/ortools/sat/routing_cuts.h @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -41,9 +42,6 @@ namespace sat { // subsets spans will point in the subset_data vector (which will be of size // exactly num_nodes). // -// Only subsets of size >= min_subset_size will be returned. This is mainly here -// to exclude subsets of size 1. -// // This is an heuristic to generate interesting cuts for TSP or other graph // based constraints. We roughly follow the algorithm described in section 6 of // "The Traveling Salesman Problem, A computational Study", David L. Applegate, @@ -53,25 +51,27 @@ namespace sat { // the asymmetric case. void GenerateInterestingSubsets(int num_nodes, const std::vector>& arcs, - int min_subset_size, int stop_at_num_components, + int stop_at_num_components, std::vector* subset_data, std::vector>* subsets); -// Given a rooted tree on n nodes represented by the parent vector, returns the -// n sets of nodes corresponding to all the possible subtree. Note that the -// output memory is just n as all subset will point into the same vector. +// Given a set of rooted tree on n nodes represented by the parent vector, +// returns the n sets of nodes corresponding to all the possible subtree. Note +// that the output memory is just n as all subset will point into the same +// vector. // -// This assumes a single rooted tree, otherwise it will not crash but the result -// will not be correct. +// This assumes no cycles, otherwise it will not crash but the result will not +// be correct. // // In the TSP context, if the tree is a Gomory-Hu cut tree, this will returns // a set of "min-cut" that contains a min-cut for all node pairs. // // TODO(user): This also allocate O(n) memory internally, we could reuse it from // call to call if needed. -void ExtractAllSubsetsFromTree(const std::vector& parent, - std::vector* subset_data, - std::vector>* subsets); +void ExtractAllSubsetsFromForest( + const std::vector& parent, std::vector* subset_data, + std::vector>* subsets, + int node_limit = std::numeric_limits::max()); // In the routing context, we usually always have lp_value in [0, 1] and only // looks at arcs with a lp_value that is not too close to zero. @@ -79,8 +79,17 @@ struct ArcWithLpValue { int tail; int head; double lp_value; + + bool operator==(const ArcWithLpValue& o) const { + return tail == o.tail && head == o.head && lp_value == o.lp_value; + } }; +// Regroups and sum the lp values on duplicate arcs or reversed arcs +// (tail->head) and (head->tail). As a side effect, we will always +// have tail <= head. +void SymmetrizeArcs(std::vector* arcs); + // Given a set of arcs on a directed graph with n nodes (in [0, num_nodes)), // returns a "parent" vector of size n encoding a rooted Gomory-Hu tree. // diff --git a/ortools/sat/scheduling_cuts.cc b/ortools/sat/scheduling_cuts.cc index 6b3169d704..d73a29f2c8 100644 --- a/ortools/sat/scheduling_cuts.cc +++ b/ortools/sat/scheduling_cuts.cc @@ -1631,7 +1631,8 @@ void GenerateCompletionTimeCutsWithEnergy(const std::string& cut_name, DCHECK_GE(event.x_start_min, sequence_start_min); const IntegerValue energy = event.energy_min; sum_duration += energy; - sum_square_duration += energy * energy; + if (!AddProductTo(energy, energy, &sum_square_duration)) break; + unscaled_lp_contrib += event.x_lp_end * ToDouble(energy); current_start_min = std::min(current_start_min, event.x_start_min); @@ -1668,14 +1669,16 @@ void GenerateCompletionTimeCutsWithEnergy(const std::string& cut_name, // duration actually equal to energy / capacity. But to keep the // computation in the integer domain, we multiply by capacity // everywhere instead. - if (AtMinOrMaxInt64( - CapAdd(CapProd(sum_duration.value(), sum_duration.value()), - sum_square_duration.value()))) { - break; // Overflow, we exit the loop. - } - const IntegerValue min_contrib = - (sum_duration * sum_duration + sum_square_duration) / 2 + - current_start_min * sum_duration * capacity; + IntegerValue min_contrib = 0; + if (!AddProductTo(sum_duration, sum_duration, &min_contrib)) break; + if (!AddTo(sum_square_duration, &min_contrib)) break; + min_contrib = min_contrib / 2; // The above is the double of the area. + + const IntegerValue intermediate = CapProdI(sum_duration, capacity); + if (AtMinOrMaxInt64I(intermediate)) break; + const IntegerValue offset = CapProdI(current_start_min, intermediate); + if (AtMinOrMaxInt64I(offset)) break; + if (!AddTo(offset, &min_contrib)) break; // We compute the efficacity in the unscaled domain where the l2 norm of // the cuts is exactly the sqrt of the sum of squared duration.