From ebb88d9c50b1b572b071eb6de46d848ea2fe8a7e Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 30 May 2023 23:25:31 +0200 Subject: [PATCH] [CP-SAT] fix corner cases in presolve; rewrite part of the cut management; automatic fixes --- ortools/sat/cp_model_presolve.cc | 71 +++++++++++++------ ortools/sat/cp_model_search.cc | 4 -- ortools/sat/cp_model_search.h | 4 ++ ortools/sat/cp_model_solver.cc | 4 +- ortools/sat/cp_model_solver.h | 2 +- ortools/sat/csharp/CpModel.cs | 2 +- ortools/sat/cuts.cc | 15 ++-- ortools/sat/diffn_util.cc | 1 + ortools/sat/integer.cc | 1 + ortools/sat/lp_utils.h | 2 +- ortools/sat/precedences.cc | 2 +- ortools/sat/presolve_context.cc | 12 ++++ ortools/sat/presolve_context.h | 5 ++ ortools/sat/probing.cc | 2 +- .../earliness_tardiness_cost_sample_sat.py | 2 +- ortools/sat/util.cc | 4 +- ortools/sat/util.h | 11 +-- ortools/sat/var_domination.cc | 4 +- 18 files changed, 98 insertions(+), 50 deletions(-) diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 767c2bb54e..3e79641ed9 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1192,21 +1192,18 @@ bool CpModelPresolver::PresolveIntProd(ConstraintProto* ct) { LinearArgumentProto* proto = ct->mutable_int_prod(); for (int i = 0; i < ct->int_prod().exprs().size(); ++i) { LinearExpressionProto expr = ct->int_prod().exprs(i); + const int64_t local_constant_factor = context_->ExpressionDivisor(expr); if (context_->IsFixed(expr)) { context_->UpdateRuleStats("int_prod: removed constant expressions."); changed = true; - constant_factor = CapProd(constant_factor, context_->FixedValue(expr)); + constant_factor = CapProd(constant_factor, local_constant_factor); continue; } else { - const int64_t coeff = expr.coeffs(0); - const int64_t offset = expr.offset(); - const int64_t gcd = - MathUtil::GCD64(static_cast(std::abs(coeff)), - static_cast(std::abs(offset))); - if (gcd != 1) { - constant_factor = CapProd(constant_factor, gcd); - expr.set_coeffs(0, coeff / gcd); - expr.set_offset(offset / gcd); + // Transfer the local_constant_factor to the constant_factor. + if (local_constant_factor != 1) { + constant_factor = CapProd(constant_factor, local_constant_factor); + expr.set_coeffs(0, expr.coeffs(0) / local_constant_factor); + expr.set_offset(expr.offset() / local_constant_factor); } } *proto->mutable_exprs(new_size++) = expr; @@ -1233,8 +1230,7 @@ bool CpModelPresolver::PresolveIntProd(ConstraintProto* ct) { // In this case, the only possible value that fit in the domains is zero. // We will check for UNSAT if zero is not achievable by the rhs below. - if (constant_factor == std::numeric_limits::min() || - constant_factor == std::numeric_limits::max()) { + if (AtMinOrMaxInt64(constant_factor)) { context_->UpdateRuleStats("int_prod: overflow if non zero"); if (!context_->IntersectDomainWith(ct->int_prod().target(), Domain(0))) { return false; @@ -1242,18 +1238,47 @@ bool CpModelPresolver::PresolveIntProd(ConstraintProto* ct) { constant_factor = 1; } - // Replace by linear! + // Replace by linear if it cannot overflow. if (ct->int_prod().exprs().size() == 1) { - LinearConstraintProto* const lin = - context_->working_model->add_constraints()->mutable_linear(); - lin->add_domain(0); - lin->add_domain(0); - AddLinearExpressionToLinearConstraint(ct->int_prod().target(), 1, lin); - AddLinearExpressionToLinearConstraint(ct->int_prod().exprs(0), - -constant_factor, lin); - context_->UpdateNewConstraintsVariableUsage(); - context_->UpdateRuleStats("int_prod: linearize product by constant."); - return RemoveConstraint(ct); + const LinearExpressionProto& target = ct->int_prod().target(); + const int64_t target_constant_factor = context_->ExpressionDivisor(target); + const int64_t gcd = MathUtil::GCD64( + static_cast(std::abs(constant_factor)), + static_cast(std::abs(target_constant_factor))); + + if (gcd != 1) { + constant_factor /= gcd; + ct->mutable_int_prod()->mutable_target()->set_offset(target.offset() / + gcd); + if (target.coeffs_size() == 1) { + ct->mutable_int_prod()->mutable_target()->set_coeffs( + 0, target.coeffs(0) / gcd); + } + } + + // Overflow ? + const Domain expr_domain = + context_->DomainSuperSetOf(ct->int_prod().exprs(0)); + const Domain expanded = + expr_domain.ContinuousMultiplicationBy(constant_factor) + .AdditionWith(context_->DomainSuperSetOf(target)); + const bool possible_overflow = + AtMinOrMaxInt64(expanded.Min()) || AtMinOrMaxInt64(expanded.Max()); + if (possible_overflow) { + // Re-add a new term with the constant factor. + ct->mutable_int_prod()->add_exprs()->set_offset(constant_factor); + } else { // Replace with a linear equation. + LinearConstraintProto* const lin = + context_->working_model->add_constraints()->mutable_linear(); + lin->add_domain(0); + lin->add_domain(0); + AddLinearExpressionToLinearConstraint(ct->int_prod().target(), 1, lin); + AddLinearExpressionToLinearConstraint(ct->int_prod().exprs(0), + -constant_factor, lin); + context_->UpdateNewConstraintsVariableUsage(); + context_->UpdateRuleStats("int_prod: linearize product by constant."); + return RemoveConstraint(ct); + } } if (constant_factor != 1) { diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index e5f9036aba..5526a91d85 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -417,8 +417,6 @@ std::function InstrumentSearchStrategy( }; } -namespace { - // This generates a valid random seed (base_seed + delta) without overflow. // We assume |delta| is small. int ValidSumSeed(int base_seed, int delta) { @@ -431,8 +429,6 @@ int ValidSumSeed(int base_seed, int delta) { return static_cast(result); } -} // namespace - // Note: in flatzinc setting, we know we always have a fixed search defined. // // Things to try: diff --git a/ortools/sat/cp_model_search.h b/ortools/sat/cp_model_search.h index 4fc38a47fe..b652a721ce 100644 --- a/ortools/sat/cp_model_search.h +++ b/ortools/sat/cp_model_search.h @@ -112,6 +112,10 @@ std::vector GetWorkSharingParams( const SatParameters& base_params, const CpModelProto& cp_model, int num_params_to_generate); +// This generates a valid random seed (base_seed + delta) without overflow. +// We assume |delta| is small. +int ValidSumSeed(int base_seed, int delta); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index bd33c7a968..e33bd1ea49 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -3219,7 +3219,7 @@ void SolveCpModelParallel(const CpModelProto& model_proto, if (params.test_feasibility_jump()) { for (int i = 0; i < params.num_workers(); ++i) { SatParameters local_params = params; - local_params.set_random_seed(params.random_seed() + i); + local_params.set_random_seed(ValidSumSeed(params.random_seed(), i)); local_params.set_feasibility_jump_decay(((i / 2) % 2) == 1 ? 0.95 : 1.0); incomplete_subsolvers.push_back(std::make_unique( @@ -3283,7 +3283,7 @@ void SolveCpModelParallel(const CpModelProto& model_proto, // - randomly select random or no randoms? for (int i = 0; i < num_feasibility_jump; ++i) { SatParameters local_params = params; - local_params.set_random_seed(params.random_seed() + i); + local_params.set_random_seed(ValidSumSeed(params.random_seed(), i)); std::string name; switch (i) { case 0: { // Use defaults. diff --git a/ortools/sat/cp_model_solver.h b/ortools/sat/cp_model_solver.h index de66ca1eaf..d6bf4344d6 100644 --- a/ortools/sat/cp_model_solver.h +++ b/ortools/sat/cp_model_solver.h @@ -51,7 +51,7 @@ std::string CpSolverResponseStats(const CpSolverResponse& response, /** * Solves the given CpModelProto. * - * This advanced API accept a Model* which allows to access more adavanced + * This advanced API accept a Model* which allows to access more advanced * features by configuring some classes in the Model before solve. * * For instance: diff --git a/ortools/sat/csharp/CpModel.cs b/ortools/sat/csharp/CpModel.cs index 6a400111b2..bfea60efda 100644 --- a/ortools/sat/csharp/CpModel.cs +++ b/ortools/sat/csharp/CpModel.cs @@ -92,7 +92,7 @@ public class CpModel /** * - * Creates an Boolean variable with given domain. + * Creates a Boolean variable with given domain. * */ public BoolVar NewBoolVar(string name) diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index faf39162c2..3579d65ff1 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -2053,6 +2053,7 @@ bool FlowCoverCutHelper::GenerateCut(const SingleNodeFlow& data) { num_out_flow_ = 0; num_out_bin_ = 0; + // Note that we need to generate a <= version, so we negate everything. cut_builder_.Clear(); for (int i = 0; i < data.in_flow.size(); ++i) { const FlowInfo& info = data.in_flow[i]; @@ -2061,30 +2062,30 @@ bool FlowCoverCutHelper::GenerateCut(const SingleNodeFlow& data) { continue; } num_in_flow_++; - cut_builder_.AddTerm(info.flow_expr, -1); + cut_builder_.AddTerm(info.flow_expr, 1); if (info.capacity > slack) { num_in_bin_++; const IntegerValue coeff = info.capacity - slack; - cut_builder_.AddConstant(-coeff); - cut_builder_.AddTerm(info.bool_expr, coeff); + cut_builder_.AddConstant(coeff); + cut_builder_.AddTerm(info.bool_expr, -coeff); } } for (int i = 0; i < data.out_flow.size(); ++i) { const FlowInfo& info = data.out_flow[i]; if (out_cover[i]) { num_out_capa_++; - cut_builder_.AddConstant(info.capacity); + cut_builder_.AddConstant(-info.capacity); } else if (info.capacity > slack) { num_out_bin_++; - cut_builder_.AddTerm(info.bool_expr, slack); + cut_builder_.AddTerm(info.bool_expr, -slack); } else { num_out_flow_++; - cut_builder_.AddTerm(info.flow_expr, 1); + cut_builder_.AddTerm(info.flow_expr, -1); } } // TODO(user): Lift the cut. - cut_ = cut_builder_.BuildConstraint(-demand, kMaxIntegerValue); + cut_ = cut_builder_.BuildConstraint(kMinIntegerValue, demand); return true; } diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc index a3cb86922f..2f881b0671 100644 --- a/ortools/sat/diffn_util.cc +++ b/ortools/sat/diffn_util.cc @@ -16,6 +16,7 @@ #include #include +#include #include #include diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 239cd2d852..cce868c588 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include diff --git a/ortools/sat/lp_utils.h b/ortools/sat/lp_utils.h index a1d8e3119d..c6d0831155 100644 --- a/ortools/sat/lp_utils.h +++ b/ortools/sat/lp_utils.h @@ -103,7 +103,7 @@ bool MPModelProtoValidationBeforeConversion(const SatParameters& params, // To satisfy our scaling requirements, any terms that is almost zero can just // be set to zero. We need to do that before operations like -// DetectImpliedIntegers(), becauses really low coefficients can cause issues +// DetectImpliedIntegers(), because really low coefficients can cause issues // and might lead to less detection. void RemoveNearZeroTerms(const SatParameters& params, MPModelProto* mp_model, SolverLogger* logger); diff --git a/ortools/sat/precedences.cc b/ortools/sat/precedences.cc index bbf8405269..d38313998b 100644 --- a/ortools/sat/precedences.cc +++ b/ortools/sat/precedences.cc @@ -559,7 +559,7 @@ void PrecedencesPropagator::AddArc( } // TODO(user): On jobshop problems with a lot of tasks per machine (500), this -// takes up a big chunck of the running time even before we find a solution. +// takes up a big chunk of the running time even before we find a solution. // This is because, for each lower bound changed, we inspect 500 arcs even // though they will never be propagated because the other bound is still at the // horizon. Find an even sparser algorithm? diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 4881c99de3..a1e6ff5204 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -206,6 +206,18 @@ Domain PresolveContext::DomainSuperSetOf( return result; } +int64_t PresolveContext::ExpressionDivisor( + const LinearExpressionProto& expr) const { + DCHECK_LE(expr.vars_size(), 1); + if (IsFixed(expr)) return FixedValue(expr); + + const int64_t coeff = expr.coeffs(0); + const int64_t offset = expr.offset(); + return static_cast( + MathUtil::GCD64(static_cast(std::abs(coeff)), + static_cast(std::abs(offset)))); +} + bool PresolveContext::ExpressionIsAffineBoolean( const LinearExpressionProto& expr) const { if (expr.vars().size() != 1) return false; diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 0a53e2316d..9899811f27 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -138,6 +138,11 @@ class PresolveContext { int64_t MaxOf(const LinearExpressionProto& expr) const; bool IsFixed(const LinearExpressionProto& expr) const; int64_t FixedValue(const LinearExpressionProto& expr) const; + // Returns a constant factor that divides the expression. + // If the expression is fixed, then the result is equal to the fixed value of + // the expression. + // Currently, it only works with expression with 0 or 1 term. + int64_t ExpressionDivisor(const LinearExpressionProto& expr) const; // Accepts any proto with two parallel vector .vars() and .coeffs(), like // LinearConstraintProto or ObjectiveProto or LinearExpressionProto but beware diff --git a/ortools/sat/probing.cc b/ortools/sat/probing.cc index ea296254eb..72532c3bdc 100644 --- a/ortools/sat/probing.cc +++ b/ortools/sat/probing.cc @@ -597,7 +597,7 @@ bool FailedLiteralProbingRound(ProbingOptions options, Model* model) { // // Even if we fixed something at evel zero, next_decision might not be // fixed! But we can fix it. It can happen because when we propagate - // with clauses, we might have a => b but not not(b) => not(a). Like a + // with clauses, we might have a => b but not(b) => not(a). Like a // => b and clause (not(a), not(b), c), propagating a will set c, but // propagating not(c) will not do anything. // diff --git a/ortools/sat/samples/earliness_tardiness_cost_sample_sat.py b/ortools/sat/samples/earliness_tardiness_cost_sample_sat.py index 7cc619be2a..33d7fc1261 100755 --- a/ortools/sat/samples/earliness_tardiness_cost_sample_sat.py +++ b/ortools/sat/samples/earliness_tardiness_cost_sample_sat.py @@ -11,7 +11,7 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. -"""Encodes an convex piecewise linear function.""" +"""Encodes a convex piecewise linear function.""" from ortools.sat.python import cp_model diff --git a/ortools/sat/util.cc b/ortools/sat/util.cc index 023edefc49..ec6921aa56 100644 --- a/ortools/sat/util.cc +++ b/ortools/sat/util.cc @@ -494,7 +494,7 @@ void MaxBoundedSubsetSum::AddMultiples(int64_t coeff, int64_t max_value) { void MaxBoundedSubsetSum::AddChoicesInternal(absl::Span values) { // Mode 1: vector of all possible sums (with duplicates). - if (!sums_.empty() && sums_.size() <= kMaxComplexityPerAdd) { + if (!sums_.empty() && sums_.size() <= max_complexity_per_add_) { const int old_size = sums_.size(); for (int i = 0; i < old_size; ++i) { for (const int64_t value : values) { @@ -510,7 +510,7 @@ void MaxBoundedSubsetSum::AddChoicesInternal(absl::Span values) { } // Mode 2: bitset of all possible sums. - if (bound_ <= kMaxComplexityPerAdd) { + if (bound_ <= max_complexity_per_add_) { if (!sums_.empty()) { expanded_sums_.assign(bound_ + 1, false); for (const int64_t s : sums_) { diff --git a/ortools/sat/util.h b/ortools/sat/util.h index 7ff5195ae1..84bddfbe52 100644 --- a/ortools/sat/util.h +++ b/ortools/sat/util.h @@ -202,8 +202,11 @@ int MoveOneUnprocessedLiteralLast( // Precondition: Both bound and all added values must be >= 0. class MaxBoundedSubsetSum { public: - MaxBoundedSubsetSum() { Reset(0); } - explicit MaxBoundedSubsetSum(int64_t bound) { Reset(bound); } + MaxBoundedSubsetSum() : max_complexity_per_add_(/*default=*/50) { Reset(0); } + explicit MaxBoundedSubsetSum(int64_t bound, int max_complexity_per_add = 50) + : max_complexity_per_add_(max_complexity_per_add) { + Reset(bound); + } // Resets to an empty set of values. // We look for the maximum sum <= bound. @@ -230,8 +233,8 @@ class MaxBoundedSubsetSum { // This assumes filtered values. void AddChoicesInternal(absl::Span values); - static constexpr int kMaxComplexityPerAdd = 50; - + // Max_complexity we are willing to pay on each Add() call. + const int max_complexity_per_add_; int64_t gcd_; int64_t bound_; int64_t current_max_; diff --git a/ortools/sat/var_domination.cc b/ortools/sat/var_domination.cc index 692908367d..806b37d5a3 100644 --- a/ortools/sat/var_domination.cc +++ b/ortools/sat/var_domination.cc @@ -844,7 +844,7 @@ bool DualBoundStrengthening::Strengthen(PresolveContext* context) { } // Note(user): If we have enforced => var fixed, we could actually - // just have removed var from the constraint it it was implied by + // just have removed var from the constraint it was implied by // another constraint. If not, because of the new affine relation we // could remove it right away. processed[PositiveRef(enf)] = true; @@ -1088,7 +1088,7 @@ void DetectDominanceRelations( dual_bound_strengthening->Reset(num_vars); for (int var = 0; var < num_vars; ++var) { - // Ignore variables that have been substitued already or are unused. + // Ignore variables that have been substituted already or are unused. if (context.IsFixed(var) || context.VariableWasRemoved(var) || context.VariableIsNotUsedAnymore(var)) { dual_bound_strengthening->CannotMove({var});