diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 2b7d160e28..5456ddd447 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -1117,15 +1117,15 @@ int CoverCutHelper::GetCoverSizeForBooleans(int relevant_size) { // Sorting can be slow, so we start by splitting the vector in 3 parts // [can always be in cover, candidates, can never be in cover]. int part1 = 0; - const double threshold = 1.0 / static_cast(relevant_size); + const double threshold = 1.0 - 1.0 / static_cast(relevant_size); for (int i = 0; i < relevant_size;) { - const double dist = base_ct_.terms[i].LpDistToMaxValue(); - if (dist < threshold) { + const double lp_value = base_ct_.terms[i].lp_value; + if (lp_value >= threshold) { // Move to part 1. std::swap(base_ct_.terms[i], base_ct_.terms[part1]); ++i; ++part1; - } else if (dist < 0.9999) { + } else if (lp_value >= 0.001) { // Keep in part 2. ++i; } else { @@ -1134,16 +1134,16 @@ int CoverCutHelper::GetCoverSizeForBooleans(int relevant_size) { std::swap(base_ct_.terms[i], base_ct_.terms[relevant_size]); } } + + // Sort by decreasing Lp value. std::sort(base_ct_.terms.begin() + part1, base_ct_.terms.begin() + relevant_size, [](const CutTerm& a, const CutTerm& b) { - const double dist_a = a.LpDistToMaxValue(); - const double dist_b = b.LpDistToMaxValue(); - if (dist_a == dist_b) { + if (a.lp_value == b.lp_value) { // Prefer low coefficients if the distance is the same. return a.coeff < b.coeff; } - return dist_a < dist_b; + return a.lp_value > b.lp_value; }); double activity = 0.0; diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 46eaabf679..3ff17185b2 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -276,7 +276,7 @@ bool LinearConstraintManager::AddCut( const double l2_norm = ComputeL2Norm(ct); // Only add cut with sufficient efficacy. - if (violation / l2_norm < 1e-5) { + if (violation / l2_norm < 1e-4) { VLOG(3) << "BAD Cut '" << type_name << "'" << " size=" << ct.vars.size() << " max_magnitude=" << ComputeInfinityNorm(ct) diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 4e8c7569ab..e600c51b6c 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -165,6 +165,35 @@ void ScatteredIntegerVector::ConvertToLinearConstraint( result->ub = upper_bound; } +void ScatteredIntegerVector::ConvertToCutData( + absl::int128 rhs, const std::vector& integer_variables, + const std::vector& lp_solution, IntegerTrail* integer_trail, + CutData* result) { + result->terms.clear(); + result->rhs = rhs; + if (is_sparse_) { + std::sort(non_zeros_.begin(), non_zeros_.end()); + for (const glop::ColIndex col : non_zeros_) { + const IntegerValue coeff = dense_vector_[col]; + if (coeff == 0) continue; + const IntegerVariable var = integer_variables[col.value()]; + CHECK(result->AppendOneTerm(var, coeff, lp_solution[col.value()], + integer_trail->LevelZeroLowerBound(var), + integer_trail->LevelZeroUpperBound(var))); + } + } else { + const int size = dense_vector_.size(); + for (glop::ColIndex col(0); col < size; ++col) { + const IntegerValue coeff = dense_vector_[col]; + if (coeff == 0) continue; + const IntegerVariable var = integer_variables[col.value()]; + CHECK(result->AppendOneTerm(var, coeff, lp_solution[col.value()], + integer_trail->LevelZeroLowerBound(var), + integer_trail->LevelZeroUpperBound(var))); + } + } +} + std::vector> ScatteredIntegerVector::GetTerms() { std::vector> result; @@ -356,12 +385,14 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() { for (const LinearConstraintInternal& ct : integer_lp_) { const ConstraintIndex row = lp_data_.CreateNewConstraint(); - // Note that we use the bounds even if they are trivial. - // - // TODO(user): This seems good for things like sum bool <= 1 since setting - // the slack in [0, 1] can lead to bound flip in the simplex. However if the - // bound is large, maybe it make more sense to use +/- infinity. - lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub)); + // TODO(user): Using trivial bound might be good for things like + // sum bool <= 1 since setting the slack in [0, 1] can lead to bound flip in + // the simplex. However if the bound is large, maybe it make more sense to + // use +/- infinity. + const double infinity = std::numeric_limits::infinity(); + lp_data_.SetConstraintBounds( + row, ct.lb_is_trivial ? -infinity : ToDouble(ct.lb), + ct.ub_is_trivial ? +infinity : ToDouble(ct.ub)); for (const auto& term : ct.terms) { lp_data_.SetCoefficient(row, term.first, ToDouble(term.second)); } @@ -831,7 +862,12 @@ bool LinearProgrammingConstraint::AnalyzeLp() { return true; } -// TODO(user): We could also strengthen the coefficients. +// Note that since we call this on the constraint with slack, we actually have +// linear expression == rhs, we can use this to propagate more! +// +// TODO(user): Also propagate on -cut ? in practice we already do that in many +// places were we try to generate the cut on -cut... But we coould do it sooner +// and more cleanly here. bool LinearProgrammingConstraint::PreprocessCut(IntegerVariable first_slack, CutData* cut) { // Because of complement, all coeffs and all terms are positive after this. @@ -841,15 +877,47 @@ bool LinearProgrammingConstraint::PreprocessCut(IntegerVariable first_slack, return false; } + // Limited DP to compute first few reachable values. + // Note that all coeff are positive. + reachable_.Reset(); + for (const CutTerm& term : cut->terms) { + reachable_.Add(term.coeff.value()); + } + + // Extra propag since we know it is actually an equality constraint. + if (cut->rhs < absl::int128(reachable_.LastValue()) && + !reachable_.MightBeReachable(static_cast(cut->rhs))) { + problem_proven_infeasible_by_cuts_ = true; + return false; + } + bool some_fixed_terms = false; bool some_relevant_positions = false; for (CutTerm& term : cut->terms) { const absl::int128 magnitude128 = term.coeff.value(); const absl::int128 range = absl::int128(term.bound_diff.value()) * magnitude128; + + IntegerValue new_diff = term.bound_diff; if (range > cut->rhs) { - // Update the term bound_diff. - term.bound_diff = static_cast(cut->rhs / magnitude128); + new_diff = static_cast(cut->rhs / magnitude128); + } + { + // Extra propag since we know it is actually an equality constraint. + absl::int128 rest128 = + cut->rhs - absl::int128(new_diff.value()) * magnitude128; + while (rest128 < absl::int128(reachable_.LastValue()) && + !reachable_.MightBeReachable(static_cast(rest128))) { + ++total_num_eq_propagations_; + CHECK_GT(new_diff, 0); + --new_diff; + rest128 += magnitude128; + } + } + + if (new_diff < term.bound_diff) { + term.bound_diff = new_diff; + const IntegerVariable var = term.expr_vars[0]; if (var < first_slack) { // Normal variable. @@ -940,19 +1008,16 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( return false; } - // Important: because we use integer_multipliers below, we cannot just - // divide by GCD or call PreventOverflow() here. + // Important: because we use integer_multipliers below to create the slack, + // we cannot try to adjust this linear combination easily. // // TODO(user): the conversion col_index -> IntegerVariable is slow and could // in principle be removed. Easy for cuts, but not so much for // implied_bounds_processor_. Note that in theory this could allow us to // use Literal directly without the need to have an IntegerVariable for them. - tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, cut_ub, - &cut_); - - // TODO(user): Remove one conversion step. - CHECK(base_ct_.FillFromLinearConstraint(cut_, expanded_lp_solution_, - integer_trail_)); + tmp_scattered_vector_.ConvertToCutData(cut_ub.value(), integer_variables_, + lp_solution_, integer_trail_, + &base_ct_); // If there are no integer (all Booleans), no need to try implied bounds // heurititics. By setting this to nullptr, we are a bit faster. @@ -1006,6 +1071,9 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( // This also make all coefficients positive. if (!PreprocessCut(first_slack, &base_ct_)) return false; + // We cannot cut sum Bool <= 1. + if (base_ct_.rhs == 1) return false; + // Constraint with slack should be tight. if (DEBUG_MODE) { double activity = 0.0; @@ -1033,6 +1101,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints( } // Try cover approach to find cut. + // TODO(user): Share common computation between two kinds. { DCHECK(base_ct_.AllCoefficientsArePositive()); @@ -1206,7 +1275,7 @@ void LinearProgrammingConstraint::AddCGCuts() { tmp_integer_multipliers_ = ScaleLpMultiplier( /*take_objective_into_account=*/false, /*ignore_trivial_constraints=*/old_gomory, tmp_lp_multipliers_, - &scaling, int64_t{1} << 52); + &scaling); if (scaling != 0) { AddCutFromConstraints("CG", tmp_integer_multipliers_); } @@ -2633,6 +2702,8 @@ std::string LinearProgrammingConstraint::Statistics() const { FormatCounter(total_num_simplex_iterations_), "\n"); absl::StrAppend(&result, " total num cut propagation: ", FormatCounter(total_num_cut_propagations_), "\n"); + absl::StrAppend(&result, " total num eq cut propagation: ", + FormatCounter(total_num_eq_propagations_), "\n"); absl::StrAppend(&result, " total num cut overflow: ", FormatCounter(num_cut_overflows_), "\n"); absl::StrAppend(&result, " total num adjust: ", FormatCounter(num_adjusts_), diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 86a919d8fe..54a4af875e 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -97,6 +97,11 @@ class ScatteredIntegerVector { const std::vector& integer_variables, IntegerValue upper_bound, LinearConstraint* result); + void ConvertToCutData(absl::int128 rhs, + const std::vector& integer_variables, + const std::vector& lp_solution, + IntegerTrail* integer_trail, CutData* result); + // Similar to ConvertToLinearConstraint(). std::vector> GetTerms(); @@ -575,7 +580,9 @@ class LinearProgrammingConstraint : public PropagatorInterface, // As we form candidate form cuts, sometimes we can propagate level zero // bounds with them. + FirstFewValues<10> reachable_; int64_t total_num_cut_propagations_ = 0; + int64_t total_num_eq_propagations_ = 0; // Some stats on the LP statuses encountered. int64_t num_solves_ = 0; diff --git a/ortools/sat/util.h b/ortools/sat/util.h index 5f78d02708..bdf0abe192 100644 --- a/ortools/sat/util.h +++ b/ortools/sat/util.h @@ -34,6 +34,7 @@ #include "ortools/sat/sat_base.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/util/random_engine.h" +#include "ortools/util/saturated_arithmetic.h" #include "ortools/util/sorted_interval_list.h" #include "ortools/util/time_limit.h" @@ -236,6 +237,65 @@ class MaxBoundedSubsetSum { std::vector filtered_values_; }; +// Simple DP to keep the set of the first n reachable value (n > 1). +// +// TODO(user): Maybe modulo some prime number we can keep more info. +// TODO(user): Another common case is a bunch of really small values and larger +// ones, so we could bound the sum of the small values and keep the first few +// reachable by the big ones. This is similar to some presolve transformations. +template +class FirstFewValues { + public: + FirstFewValues() { Reset(); } + + void Reset() { + reachable_.fill(std::numeric_limits::max()); + reachable_[0] = 0; + new_reachable_[0] = 0; + } + + // We assume the given positive value can be added as many time as wanted. + // + // TODO(user): Implement Add() with an upper bound on the multiplicity. + void Add(const int64_t positive_value) { + DCHECK_GT(positive_value, 0); + if (positive_value >= reachable_.back()) return; + + // We copy from reachable_[i] to new_reachable_[j]. + // The position zero is already copied. + int i = 1; + int j = 1; + for (int base = 0; j < n && base < n; ++base) { + const int64_t candidate = CapAdd(new_reachable_[base], positive_value); + while (j < n && i < n && reachable_[i] < candidate) { + new_reachable_[j++] = reachable_[i++]; + } + if (j < n) { + // Eliminate duplicates. + while (i < n && reachable_[i] == candidate) i++; + + // insert candidate in its final place. + new_reachable_[j++] = candidate; + } + } + std::swap(reachable_, new_reachable_); + } + + // Returns true iff sum might be expressible as a weighted sum of the added + // value. Any sum >= LastValue() is always considered potentially reachable. + bool MightBeReachable(int64_t sum) const { + if (sum >= reachable_.back()) return true; + return std::binary_search(reachable_.begin(), reachable_.end(), sum); + } + + const std::array& reachable() const { return reachable_; } + int64_t LastValue() const { return reachable_.back(); } + + private: + std::array reachable_; + std::array new_reachable_; +}; + // Use Dynamic programming to solve a single knapsack. This is used by the // presolver to simplify variables appearing in a single linear constraint. //