diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 05c80498d1..84544ae2b5 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -1137,6 +1137,114 @@ void IntegerRoundingCutHelper::ComputeCut( DivideByGCD(cut); } +bool CoverCutHelper::TrySimpleKnapsack( + const LinearConstraint base_ct, const std::vector& lp_values, + const std::vector& lower_bounds, + const std::vector& upper_bounds) { + const int base_size = lp_values.size(); + + // Fill terms with a rewrite of the base constraint where all coeffs & + // variables are positive by using either (X - LB) or (UB - X) as new + // variables. + terms_.clear(); + IntegerValue rhs = base_ct.ub; + IntegerValue sum_of_diff(0); + for (int i = 0; i < base_size; ++i) { + const IntegerValue coeff = base_ct.coeffs[i]; + const IntegerValue bound_diff = upper_bounds[i] - lower_bounds[i]; + if (!AddProductTo(IntTypeAbs(coeff), bound_diff, &sum_of_diff)) { + return false; + } + const IntegerValue diff = IntTypeAbs(coeff) * bound_diff; + if (coeff > 0) { + if (!AddProductTo(-coeff, lower_bounds[i], &rhs)) return false; + terms_.push_back({i, ToDouble(upper_bounds[i]) - lp_values[i], diff}); + } else { + if (!AddProductTo(-coeff, upper_bounds[i], &rhs)) return false; + terms_.push_back({i, lp_values[i] - ToDouble(lower_bounds[i]), diff}); + } + } + + // Try simple cover heuristic. + // Look for violated CUT if the form sum (UB - X) or (X - LB) >= 1. + double activity = 0.0; + int new_size = 0; + std::sort(terms_.begin(), terms_.end(), [](const Term& a, const Term& b) { + return a.dist_to_max_value < b.dist_to_max_value; + }); + for (int i = 0; i < terms_.size(); ++i) { + const Term& term = terms_[i]; + activity += term.dist_to_max_value; + if (activity > 0.99) { + new_size = i; // before this entry. + break; + } + + rhs -= term.diff; + } + + // If the rhs is now negative, we have a cut. + if (rhs >= 0) return false; + if (new_size == 0) return false; + + // Transfrom to a minimal cover. We want to greedily remove the largest coeff + // first, so we have more chance for the "lifting" below which can increase + // the cut violation. + terms_.resize(new_size); + std::sort(terms_.begin(), terms_.end(), + [&base_ct](const Term& a, const Term& b) { + return IntTypeAbs(base_ct.coeffs[a.index]) > + IntTypeAbs(base_ct.coeffs[b.index]); + }); + in_cut_.assign(base_ct.vars.size(), false); + + // Remove terms greedily to get a minimal cover. + // Compute the cut at the same time. + cut_.ClearTerms(); + cut_.lb = kMinIntegerValue; + cut_.ub = IntegerValue(-1); + IntegerValue max_coeff(0); + for (const Term term : terms_) { + if (term.diff + rhs < 0) { + rhs += term.diff; + continue; + } + in_cut_[term.index] = true; + max_coeff = std::max(max_coeff, IntTypeAbs(base_ct.coeffs[term.index])); + cut_.vars.push_back(base_ct.vars[term.index]); + if (base_ct.coeffs[term.index] > 0) { + cut_.coeffs.push_back(IntegerValue(1)); + cut_.ub += upper_bounds[term.index]; + } else { + cut_.coeffs.push_back(IntegerValue(-1)); + cut_.ub -= lower_bounds[term.index]; + } + } + + // Basic Lifting. A move in activity of 1 in the cut correspond to a move of + // max_coeff in the original constraint. + num_lifting_ = 0; + for (int i = 0; i < base_size; ++i) { + if (in_cut_[i]) continue; + const IntegerValue magnitude = IntTypeAbs(base_ct.coeffs[i]); + if (magnitude < max_coeff) continue; + + ++num_lifting_; + const IntegerValue f = FloorRatio(magnitude, max_coeff); + if (base_ct.coeffs[i] > 0) { // Add f * (X - LB) + cut_.coeffs.push_back(IntegerValue(f)); + cut_.vars.push_back(base_ct.vars[i]); + cut_.ub += lower_bounds[i] * f; + } else { // Add f * (UB - X) + cut_.coeffs.push_back(IntegerValue(-f)); + cut_.vars.push_back(base_ct.vars[i]); + cut_.ub -= upper_bounds[i] * f; + } + } + + return true; +} + CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z, IntegerVariable x, IntegerVariable y, diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 83bc4386ba..591f0d023a 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -226,6 +226,36 @@ class IntegerRoundingCutHelper { std::vector> tmp_terms_; }; +// Helper to find knapsack or flow cover cuts (not yet implemented). +class CoverCutHelper { + public: + // Try to find a cut with a knapsack heuristic. + // If this returns true, you can get the cut via cut(). + bool TrySimpleKnapsack(const LinearConstraint base_ct, + const std::vector& lp_values, + const std::vector& lower_bounds, + const std::vector& upper_bounds); + + // If successful, info about the last generated cut. + LinearConstraint* mutable_cut() { return &cut_; } + const LinearConstraint& cut() const { return cut_; } + + // Single line of text that we append to the cut log line. + const std::string Info() { return absl::StrCat("lift=", num_lifting_); } + + private: + struct Term { + int index; + double dist_to_max_value; + IntegerValue diff; + }; + std::vector terms_; + std::vector in_cut_; + + LinearConstraint cut_; + int num_lifting_; +}; + // If a variable is away from its upper bound by more than value 1.0, then it // cannot be part of a cover that will violate the lp solution. This method // returns a reduced constraint by removing such variables from the given diff --git a/ortools/sat/linear_constraint.h b/ortools/sat/linear_constraint.h index d8940404cc..3e68b5731c 100644 --- a/ortools/sat/linear_constraint.h +++ b/ortools/sat/linear_constraint.h @@ -50,6 +50,11 @@ struct LinearConstraint { coeffs.push_back(coeff); } + void ClearTerms() { + vars.clear(); + coeffs.clear(); + } + std::string DebugString() const { std::string result; if (lb.value() > kMinIntegerValue) { diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 577586c11b..fd14340952 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -701,27 +701,53 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( } } - // Get the cut using some integer rounding heuristic. - RoundingOptions options; - options.max_scaling = sat_parameters_.max_integer_rounding_scaling(); - integer_rounding_cut_helper_.ComputeCut(options, tmp_lp_values_, tmp_var_lbs_, - tmp_var_ubs_, - &implied_bounds_processor_, &cut_); + bool at_least_one_added = false; + // Try cover appraoch to find cut. + { + if (cover_cut_helper_.TrySimpleKnapsack(cut_, tmp_lp_values_, tmp_var_lbs_, + tmp_var_ubs_)) { + at_least_one_added |= PostprocessAndAddCut( + absl::StrCat(name, "_K"), cover_cut_helper_.Info(), first_new_var, + first_slack, ib_slack_infos, cover_cut_helper_.mutable_cut()); + } + } + + // Try integer rounding heuristic to find cut. + { + RoundingOptions options; + options.max_scaling = sat_parameters_.max_integer_rounding_scaling(); + integer_rounding_cut_helper_.ComputeCut(options, tmp_lp_values_, + tmp_var_lbs_, tmp_var_ubs_, + &implied_bounds_processor_, &cut_); + at_least_one_added |= PostprocessAndAddCut( + name, + absl::StrCat("num_lifted_booleans=", + integer_rounding_cut_helper_.NumLiftedBooleans()), + first_new_var, first_slack, ib_slack_infos, &cut_); + } + return at_least_one_added; +} + +bool LinearProgrammingConstraint::PostprocessAndAddCut( + const std::string& name, const std::string& info, + IntegerVariable first_new_var, IntegerVariable first_slack, + const std::vector& ib_slack_infos, + LinearConstraint* cut) { // Compute the activity. Warning: the cut no longer have the same size so we // cannot use tmp_lp_values_. Note that the substitution below shouldn't // change the activity by definition. double activity = 0.0; - for (int i = 0; i < cut_.vars.size(); ++i) { - if (cut_.vars[i] < first_new_var) { + for (int i = 0; i < cut->vars.size(); ++i) { + if (cut->vars[i] < first_new_var) { activity += - ToDouble(cut_.coeffs[i]) * expanded_lp_solution_[cut_.vars[i]]; + ToDouble(cut->coeffs[i]) * expanded_lp_solution_[cut->vars[i]]; } } const double kMinViolation = 1e-4; - const double violation = activity - ToDouble(cut_.ub); + const double violation = activity - ToDouble(cut->ub); if (violation < kMinViolation) { - VLOG(3) << "Bad cut " << activity << " <= " << ToDouble(cut_.ub); + VLOG(3) << "Bad cut " << activity << " <= " << ToDouble(cut->ub); return false; } @@ -729,26 +755,26 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( { int num_slack = 0; tmp_dense_vector_.assign(integer_variables_.size(), IntegerValue(0)); - IntegerValue cut_ub = cut_.ub; + IntegerValue cut_ub = cut->ub; bool overflow = false; - for (int i = 0; i < cut_.vars.size(); ++i) { - const IntegerVariable var = cut_.vars[i]; + for (int i = 0; i < cut->vars.size(); ++i) { + const IntegerVariable var = cut->vars[i]; // Simple copy for non-slack variables. if (var < first_new_var) { const glop::ColIndex col = gtl::FindOrDie(mirror_lp_variable_, PositiveVariable(var)); if (VariableIsPositive(var)) { - tmp_dense_vector_[col] += cut_.coeffs[i]; + tmp_dense_vector_[col] += cut->coeffs[i]; } else { - tmp_dense_vector_[col] -= cut_.coeffs[i]; + tmp_dense_vector_[col] -= cut->coeffs[i]; } continue; } // Replace slack from bound substitution. if (var < first_slack) { - const IntegerValue multiplier = cut_.coeffs[i]; + const IntegerValue multiplier = cut->coeffs[i]; const int index = (var.value() - first_new_var.value()) / 2; CHECK_LT(index, ib_slack_infos.size()); @@ -776,7 +802,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( ++num_slack; const int slack_index = (var.value() - first_slack.value()) / 2; const glop::RowIndex row = tmp_slack_rows_[slack_index]; - const IntegerValue multiplier = -cut_.coeffs[i]; + const IntegerValue multiplier = -cut->coeffs[i]; if (!AddLinearExpressionMultiple(multiplier, integer_lp_[row].terms, &tmp_dense_vector_)) { overflow = true; @@ -796,23 +822,22 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( } VLOG(3) << " num_slack: " << num_slack; - ConvertToLinearConstraint(tmp_dense_vector_, cut_ub, &cut_); + ConvertToLinearConstraint(tmp_dense_vector_, cut_ub, cut); } // Display some stats used for investigation of cut generation. - const std::string extra_info = absl::StrCat( - "num_ib_substitutions=", ib_slack_infos.size(), " num_lifted_booleans=", - integer_rounding_cut_helper_.NumLiftedBooleans()); + const std::string extra_info = + absl::StrCat(info, " num_ib_substitutions=", ib_slack_infos.size()); const double new_violation = - ComputeActivity(cut_, expanded_lp_solution_) - ToDouble(cut_.ub); + ComputeActivity(*cut, expanded_lp_solution_) - ToDouble(cut_.ub); if (std::abs(violation - new_violation) >= 1e-4) { VLOG(1) << "Violation discrepancy after slack removal. " << " before = " << violation << " after = " << new_violation; } - DivideByGCD(&cut_); - return constraint_manager_.AddCut(cut_, name, expanded_lp_solution_, + DivideByGCD(cut); + return constraint_manager_.AddCut(*cut, name, expanded_lp_solution_, extra_info); } diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index afa2df5a40..420a3854e7 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -197,6 +197,13 @@ class LinearProgrammingConstraint : public PropagatorInterface, const std::vector>& integer_multipliers); + // Second half of AddCutFromConstraints(). + bool PostprocessAndAddCut( + const std::string& name, const std::string& info, + IntegerVariable first_new_var, IntegerVariable first_slack, + const std::vector& ib_slack_infos, + LinearConstraint* cut); + // Computes and adds the corresponding type of cuts. // This can currently only be called at the root node. void AddCGCuts(); @@ -341,6 +348,7 @@ class LinearProgrammingConstraint : public PropagatorInterface, // Temporary data for cuts. ZeroHalfCutHelper zero_half_cut_helper_; + CoverCutHelper cover_cut_helper_; IntegerRoundingCutHelper integer_rounding_cut_helper_; LinearConstraint cut_; gtl::ITIVector tmp_dense_vector_; diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index fb818552b5..5390af10b8 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -215,37 +215,6 @@ class FuMalikSymmetryBreaker { } // namespace -void MinimizeCore(SatSolver* solver, std::vector* core) { - std::vector temp = *core; - std::reverse(temp.begin(), temp.end()); - solver->Backtrack(0); - solver->SetAssumptionLevel(0); - - // Note that this Solve() is really fast, since the solver should detect that - // the assumptions are unsat with unit propagation only. This is just a - // convenient way to remove assumptions that are propagated by the one before - // them. - const SatSolver::Status status = - solver->ResetAndSolveWithGivenAssumptions(temp); - if (status != SatSolver::ASSUMPTIONS_UNSAT) { - if (status != SatSolver::LIMIT_REACHED) { - CHECK_NE(status, SatSolver::FEASIBLE); - // This should almost never happen, but it is not impossible. The reason - // is that the solver may delete some learned clauses required by the unit - // propagation to show that the core is unsat. - LOG(WARNING) << "This should only happen rarely! otherwise, investigate. " - << "Returned status is " << SatStatusString(status); - } - return; - } - temp = solver->GetLastIncompatibleDecisions(); - if (temp.size() < core->size()) { - VLOG(1) << "minimization " << core->size() << " -> " << temp.size(); - std::reverse(temp.begin(), temp.end()); - *core = temp; - } -} - void MinimizeCoreWithPropagation(SatSolver* solver, std::vector* core) { if (solver->IsModelUnsat()) return; diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h index bdaade5da5..288cd1de11 100644 --- a/ortools/sat/optimization.h +++ b/ortools/sat/optimization.h @@ -28,15 +28,6 @@ namespace operations_research { namespace sat { -// Tries to minimize the given UNSAT core with a really simple heuristic. -// The idea is to remove literals that are consequences of others in the core. -// We already know that in the initial order, no literal is propagated by the -// one before it, so we just look for propagation in the reverse order. -// -// Important: The given SatSolver must be the one that just produced the given -// core. -void MinimizeCore(SatSolver* solver, std::vector* core); - // Like MinimizeCore() with a slower but strictly better heuristic. This // algorithm should produce a minimal core with respect to propagation. We put // each literal of the initial core "last" at least once, so if such literal can diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc index 9fd5b92452..4dda0e820f 100644 --- a/ortools/sat/sat_solver.cc +++ b/ortools/sat/sat_solver.cc @@ -2539,5 +2539,36 @@ std::string SatStatusString(SatSolver::Status status) { return "UNKNOWN"; } +void MinimizeCore(SatSolver* solver, std::vector* core) { + std::vector temp = *core; + std::reverse(temp.begin(), temp.end()); + solver->Backtrack(0); + solver->SetAssumptionLevel(0); + + // Note that this Solve() is really fast, since the solver should detect that + // the assumptions are unsat with unit propagation only. This is just a + // convenient way to remove assumptions that are propagated by the one before + // them. + const SatSolver::Status status = + solver->ResetAndSolveWithGivenAssumptions(temp); + if (status != SatSolver::ASSUMPTIONS_UNSAT) { + if (status != SatSolver::LIMIT_REACHED) { + CHECK_NE(status, SatSolver::FEASIBLE); + // This should almost never happen, but it is not impossible. The reason + // is that the solver may delete some learned clauses required by the unit + // propagation to show that the core is unsat. + LOG(WARNING) << "This should only happen rarely! otherwise, investigate. " + << "Returned status is " << SatStatusString(status); + } + return; + } + temp = solver->GetLastIncompatibleDecisions(); + if (temp.size() < core->size()) { + VLOG(1) << "minimization " << core->size() << " -> " << temp.size(); + std::reverse(temp.begin(), temp.end()); + *core = temp; + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/sat_solver.h b/ortools/sat/sat_solver.h index a3114a8a5a..5d52525967 100644 --- a/ortools/sat/sat_solver.h +++ b/ortools/sat/sat_solver.h @@ -823,6 +823,15 @@ class SatSolver { DISALLOW_COPY_AND_ASSIGN(SatSolver); }; +// Tries to minimize the given UNSAT core with a really simple heuristic. +// The idea is to remove literals that are consequences of others in the core. +// We already know that in the initial order, no literal is propagated by the +// one before it, so we just look for propagation in the reverse order. +// +// Important: The given SatSolver must be the one that just produced the given +// core. +void MinimizeCore(SatSolver* solver, std::vector* core); + // ============================================================================ // Model based functions. //