From 5827eadbcc4484a8e49e36d46c196d5b40835b6c Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Wed, 9 Jan 2019 11:25:37 +0100 Subject: [PATCH] speedup sat internals; fix bug with model with one trivially false XOr constraint --- makefiles/Makefile.gen.mk | 4 +- ortools/sat/cp_model_presolve.cc | 5 +- ortools/sat/cp_model_solver.cc | 11 ++ ortools/sat/cuts.cc | 206 +++++++++++++++++++++++++++++++ ortools/sat/cuts.h | 59 +++++++++ ortools/sat/integer.cc | 136 +++++++++++--------- ortools/sat/integer.h | 19 +++ 7 files changed, 381 insertions(+), 59 deletions(-) diff --git a/makefiles/Makefile.gen.mk b/makefiles/Makefile.gen.mk index f2c0e3d81b..b86ed23410 100644 --- a/makefiles/Makefile.gen.mk +++ b/makefiles/Makefile.gen.mk @@ -1797,7 +1797,9 @@ objs/sat/lp_utils.$O: ortools/sat/lp_utils.cc ortools/sat/lp_utils.h \ ortools/lp_data/lp_print_utils.h ortools/lp_data/sparse_row.h \ ortools/lp_data/matrix_scaler.h ortools/sat/boolean_problem.h \ ortools/algorithms/sparse_permutation.h ortools/sat/simplification.h \ - ortools/base/adjustable_priority_queue.h | $(OBJ_DIR)/sat + ortools/base/adjustable_priority_queue.h ortools/sat/integer.h \ + ortools/util/rev.h ortools/util/saturated_arithmetic.h \ + ortools/util/sorted_interval_list.h | $(OBJ_DIR)/sat $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Slp_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Slp_utils.$O objs/sat/optimization.$O: ortools/sat/optimization.cc \ diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 6ea98502e6..1452cf41d1 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2031,7 +2031,10 @@ void Probe(TimeLimit* global_time_limit, PresolveContext* context) { } encoder->AddAllImplicationsBetweenAssociatedLiterals(); auto* sat_solver = model.GetOrCreate(); - sat_solver->Propagate(); + if (!sat_solver->Propagate()) { + context->is_unsat = true; + return; + } // Probe. // diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 6eddab00ec..c8b3bb1a1f 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -1496,7 +1496,18 @@ CpSolverResponse SolveCpModelInternal( CHECK(SolutionIsFeasible(model_proto, std::vector(response.solution().begin(), response.solution().end()))); + if (!model_proto.has_objective()) { + if (parameters.enumerate_all_solutions()) { + model->Add( + ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack()); + } else { + return response; + } + } + // TODO(user): Collect objective value in optimization model and + // constrain it in the following search. } + // TODO(user): Remove conflicts used during hint search. model->GetOrCreate()->set_max_number_of_conflicts( old_conflict_limit); } diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 987b2e9c9c..2e7b503c2c 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -552,5 +552,211 @@ CutGenerator CreateKnapsackCoverCutGenerator( return result; } +namespace { + +// Returns 0 if there is none. +// Note that we normalize the fractionality by its coefficient. +IntegerValue MagnitudeOfMostFractionalVariable( + const std::vector& lp_values, const LinearConstraint& cut) { + double best_score = 0; + IntegerValue best_magnitude(0); + int num_fractional_vars = 0; + const int size = lp_values.size(); + for (int i = 0; i < size; ++i) { + const IntegerValue coeff = cut.coeffs[i]; + const double value = lp_values[i]; + const double score = + std::abs(value - std::round(value)) * ToDouble(IntTypeAbs(coeff)); + if (score > best_score) { + best_score = score; + best_magnitude = IntTypeAbs(coeff); + } + if (std::abs(value - std::round(value)) > 0.01) { + ++num_fractional_vars; + VLOG(3) << "value: " << value << " coeff: " << coeff + << " score:" << score; + } + } + VLOG(2) << "num_fractional_vars: " << num_fractional_vars << "/" << size; + return best_magnitude; +} + +} // namespace + +std::function GetSuperAdditiveRoundingFunction( + IntegerValue remainder, IntegerValue divisor, IntegerValue max_scaling) { + const IntegerValue target = CeilRatio(divisor, IntegerValue(2)) - 1; + const IntegerValue t = std::max( + IntegerValue(1), + remainder == 0 ? IntegerValue(1) + : std::min(max_scaling, CeilRatio(target, remainder))); + const IntegerValue threshold = std::max(target, t * remainder); + return [t, threshold, divisor](IntegerValue coeff) { + const IntegerValue ratio = FloorRatio(t * coeff, divisor); + const IntegerValue remainder = t * coeff - ratio * divisor; + return 2 * ratio + (remainder > threshold ? 1 : 0); + }; +} + +void IntegerRoundingCut(std::vector lp_values, + std::vector lower_bounds, + std::vector upper_bounds, + LinearConstraint* cut) { + const int size = lp_values.size(); + CHECK_EQ(lower_bounds.size(), size); + CHECK_EQ(upper_bounds.size(), size); + CHECK_EQ(cut->vars.size(), size); + CHECK_EQ(cut->coeffs.size(), size); + CHECK_EQ(cut->lb, kMinIntegerValue); + + // Test the tighteness precondition. Note that we use a big tolerance. + // This is not really needed, but if the constraint is not tight, there is + // little chance to generate a cut that violate the LP. And given how this + // is currently used, it should always be tight. + { + double activity = 0.0; + for (int i = 0; i < size; ++i) { + activity += lp_values[i] * ToDouble(cut->coeffs[i]); + } + if (std::abs(activity - ToDouble(cut->ub)) > 1.0) { + VLOG(1) << "Issue, base constraint not tight " << activity + << " <= " << ToDouble(cut->ub); + *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); + return; + } + } + + // Find the magnitude of the most fractional variable, note that we normalize + // the fractionality by its coefficient. + const IntegerValue divisor = + MagnitudeOfMostFractionalVariable(lp_values, *cut); + if (divisor == 0) { + VLOG(1) << "Issue, no fractional variables."; + *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); + return; + } + + // To simplify the code below, we make all coefficients positive. + std::vector change_sign_at_postprocessing(size, false); + for (int i = 0; i < size; ++i) { + if (cut->coeffs[i] > 0) continue; + + change_sign_at_postprocessing[i] = true; + + cut->coeffs[i] = -cut->coeffs[i]; + lp_values[i] = -lp_values[i]; + + std::swap(lower_bounds[i], upper_bounds[i]); + lower_bounds[i] = -lower_bounds[i]; + upper_bounds[i] = -upper_bounds[i]; + } + + // Shift each variable using its lower/upper bound so that no variable can + // change sign. + bool overflow = false; + std::vector shifts(size, IntegerValue(0)); + std::vector var_is_positive_or_zero(size, true); + for (int i = 0; i < size && !overflow; ++i) { + const IntegerValue coeff = cut->coeffs[i]; + if (coeff == 0) continue; + + // Note that since we use ToDouble() this code works fine with lb/ub at + // min/max integer value. + const double value = lp_values[i]; + const IntegerValue lb = lower_bounds[i]; + const IntegerValue ub = upper_bounds[i]; + if (std::abs(value - ToDouble(lb)) < std::abs(value - ToDouble(ub))) { + // We want coeff * (X - lb) so the new var is >= 0. + var_is_positive_or_zero[i] = true; + shifts[i] = lb; + } else { + // We want coeff * (X - ub) so the new var is <= 0. + var_is_positive_or_zero[i] = false; + shifts[i] = ub; + } + + // coeff * X = coeff * (X - shift) + coeff * shift. + if (!AddProductTo(-coeff, shifts[i], &cut->ub)) { + overflow = true; + break; + } + + // Deal with fixed variable, no need to shift back in this case, we can + // just remove the term. + if (lb == ub) { + shifts[i] = IntegerValue(0); + cut->coeffs[i] = IntegerValue(0); + } + } + if (overflow) { + VLOG(1) << "Issue, overflow."; + *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); + return; + } + + // We will adjust coefficient that are close to an exact multiple of divisor + // to an exact multiple. This is meant to get rid of small errors that appears + // due to rounding error in our exact computation of the initial constraint + // given to this class. + // + // TODO(user): Tune the threshold or use a parameter. Maybe it should depend + // on the number of term in the constraint? But the basic idea is that we do + // not want to change the rhs_remainder (see below) by too much. So here we + // change it at most by: num_terms * 0.0002. Note that in practice we don't + // except a lot of terms to be close to a multiple of divisor. + const IntegerValue adjust_threshold = divisor / IntegerValue(5000); + for (int i = 0; i < size; ++i) { + const IntegerValue coeff = cut->coeffs[i]; + const IntegerValue diff( + CapSub(upper_bounds[i].value(), lower_bounds[i].value())); + if (var_is_positive_or_zero[i]) { + // Adjust coeff of the form k * divisor - epsilon. + const IntegerValue remainder = + CeilRatio(coeff, divisor) * divisor - coeff; + if (CapProd(diff.value(), remainder.value()) > adjust_threshold) continue; + cut->ub += remainder * diff; + cut->coeffs[i] += remainder; + } else { + // Adjust coeff of the form k * divisor + epsilon. + const IntegerValue remainder = + coeff - FloorRatio(coeff, divisor) * divisor; + if (CapProd(diff.value(), remainder.value()) > adjust_threshold) continue; + cut->ub += remainder * diff; + cut->coeffs[i] -= remainder; + } + } + + // Create the super-additive function f(). + const IntegerValue rhs_remainder = + cut->ub - FloorRatio(cut->ub, divisor) * divisor; + const IntegerValue kMaxScaling(4); + const auto f = + GetSuperAdditiveRoundingFunction(rhs_remainder, divisor, kMaxScaling); + + // Apply f() to the cut. + cut->ub = f(cut->ub); + for (int i = 0; i < cut->coeffs.size(); ++i) { + const IntegerValue coeff = cut->coeffs[i]; + if (coeff == 0) continue; + if (var_is_positive_or_zero[i]) { + cut->coeffs[i] = f(coeff); + } else { + cut->coeffs[i] = -f(-coeff); + } + } + + // Remove the bound shifts so the constraint is expressed in the original + // variables and do some basic post-processing. + for (int i = 0; i < size; ++i) { + cut->ub = IntegerValue( + CapAdd((cut->coeffs[i] * shifts[i]).value(), cut->ub.value())); + if (change_sign_at_postprocessing[i]) { + cut->coeffs[i] = -cut->coeffs[i]; + } + } + RemoveZeroTerms(cut); + DivideByGCD(cut); +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index b2aa0fe4f6..80fd473ecb 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -41,6 +41,65 @@ struct CutGenerator { generate_cuts; }; +// Visible for testing. Returns a function f on integers such that: +// - f is non-decreasing. +// - f is super-additive: f(a) + f(b) <= f(a + b) +// - 1 <= f(divisor) <= 2 * max_scaling +// - For all x, f(x * divisor) = x * f(divisor) +// - For all x, f(x * divisor + remainder) = x * f(divisor) +// +// Preconditions: +// - 0 <= remainder < divisor. +// - 1 <= max_scaling. +// +// This is used in IntegerRoundingCut() and is responsible for "strengthening" +// the cut. Just taking f(x) = x / divisor result in the non-strengthened cut +// and using any function that stricly dominate this one is better. +// +// Note(user): I could have used the MIR (Mixed integer rounding function), but +// I prefered to use the function described in "Strenghtening Chvatal-Gomory +// cuts and Gomory fractional cuts", Adam N. Letchfrod, Andrea Lodi. This is +// because it gives better result for small value of max_scaling, and we do not +// want to have large integer coefficient. +// +// TODO(user): Still not clear to me if this can be improved or not. +std::function GetSuperAdditiveRoundingFunction( + IntegerValue remainder, IntegerValue divisor, IntegerValue max_scaling); + +// Given an upper bounded linear constraint, this function tries to transform it +// to a valid cut that violate the given LP solution using integer rounding. +// Note that the returned cut might not always violate the LP solution, in which +// case it can be discarded. +// +// What this does is basically take the integer division of the constraint by an +// integer. If the coefficients where doubles, this would be the same as scaling +// the constraint and then rounding. We choose the coefficient of the most +// fractional variable (rescaled by its coefficient) as the divisor, but there +// are other possible alternatives. +// +// Note that if the constraint is tight under the given lp solution, and if +// there is a unique variable not at one of its bounds and fractional, then we +// are guaranteed to generate a cut that violate the current LP solution. This +// should be the case for Chvatal-Gomory base constraints modulo our loss of +// precision while doing exact integer computations. +// +// Precondition: +// - We assumes that the given initial constraint is tight using the given lp +// values. This could be relaxed, but for now it should always be the case, so +// we log a message and abort if not, to ease debugging. +// - The IntegerVariable of the cuts are not used here. We assumes that the +// first three vectors are in one to one correspondance with the initial order +// of the variable in the cut. +// +// TODO(user): There is a bunch of heuristic involved here, and we could spend +// more effort tunning them. In particular, one can try many heuristics and keep +// the best looking cut (or more than one). This is not on the critical code +// path, so we can spend more effort in finding good cuts. +void IntegerRoundingCut(std::vector lp_values, + std::vector lower_bounds, + std::vector upper_bounds, + LinearConstraint* cut); + // 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/integer.cc b/ortools/sat/integer.cc index b32970e8ca..68377073e7 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -709,10 +709,7 @@ int IntegerTrail::FindLowestTrailIndexThatExplainBound( } } -// We try to relax the reason in a smart way here by minimizing the maximum -// trail indices of the literals appearing in reason. -// -// TODO(user): use priority queue instead of O(n^2) algo. +// TODO(user): Get rid of this function and only keep the trail index one? void IntegerTrail::RelaxLinearReason( IntegerValue slack, absl::Span coeffs, std::vector* reason) const { @@ -726,73 +723,98 @@ void IntegerTrail::RelaxLinearReason( indices[i] = vars_[(*reason)[i].var].current_trail_index; } - const int num_vars = vars_.size(); - while (slack > 0) { - int best_i = -1; - for (int i = 0; i < size; ++i) { - if (indices[i] < num_vars) continue; // level zero. - if (best_i != -1 && indices[i] < indices[best_i]) continue; - const TrailEntry& entry = integer_trail_[indices[i]]; - const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index]; + RelaxLinearReason(slack, coeffs, &indices); - // Note that both terms of the product are positive. - if (CapProd(coeffs[i].value(), - (entry.bound - previous_entry.bound).value()) > slack) { - continue; - } - best_i = i; - } - if (best_i == -1) return; - - const TrailEntry& entry = integer_trail_[indices[best_i]]; - const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index]; - indices[best_i] = entry.prev_trail_index; - (*reason)[best_i].bound = previous_entry.bound; - slack -= coeffs[best_i] * (entry.bound - previous_entry.bound); + reason->clear(); + for (const int i : indices) { + reason->push_back(IntegerLiteral::GreaterOrEqual(integer_trail_[i].var, + integer_trail_[i].bound)); } } +// TODO(user): When this is called during a reason computation, we can use +// the term already part of the reason we are constructed to optimize this +// further. void IntegerTrail::RelaxLinearReason(IntegerValue slack, absl::Span coeffs, std::vector* trail_indices) const { DCHECK_GT(slack, 0); - const int size = trail_indices->size(); - const int num_vars = vars_.size(); - while (slack > 0) { - int best_i = -1; - for (int i = 0; i < size; ++i) { - if ((*trail_indices)[i] < num_vars) continue; // level zero. - if (coeffs[i] > slack) continue; - if (best_i != -1 && (*trail_indices)[i] < (*trail_indices)[best_i]) { - continue; - } + DCHECK(relax_heap_.empty()); - const TrailEntry& entry = integer_trail_[(*trail_indices)[i]]; - const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index]; - - // Note that both terms of the product are positive. - if (CapProd(coeffs[i].value(), - (entry.bound - previous_entry.bound).value()) > slack) { - continue; - } - best_i = i; - } - if (best_i == -1) break; - - const TrailEntry& entry = integer_trail_[(*trail_indices)[best_i]]; - const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index]; - (*trail_indices)[best_i] = entry.prev_trail_index; - slack -= coeffs[best_i] * (entry.bound - previous_entry.bound); - } - - // Remove level zero indices. + // We start by filtering *trail_indices: + // - remove all level zero entries. + // - keep the one that cannot be relaxed. + // - move the other one the the relax_heap_ (and creating the heap). int new_size = 0; - for (const int index : *trail_indices) { - if (index >= num_vars) { + const int size = coeffs.size(); + const int num_vars = vars_.size(); + for (int i = 0; i < size; ++i) { + const int index = (*trail_indices)[i]; + + // We ignore level zero entries. + if (index < num_vars) continue; + + // If the coeff is too large, we cannot relax this entry. + const IntegerValue coeff = coeffs[i]; + if (coeff > slack) { (*trail_indices)[new_size++] = index; + continue; } + + // Note that both terms of the product are positive. + const TrailEntry& entry = integer_trail_[index]; + const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index]; + const int64 diff = + CapProd(coeff.value(), (entry.bound - previous_entry.bound).value()); + if (diff > slack) { + (*trail_indices)[new_size++] = index; + continue; + } + + relax_heap_.push_back({index, coeff, diff}); } trail_indices->resize(new_size); + std::make_heap(relax_heap_.begin(), relax_heap_.end()); + + while (slack > 0 && !relax_heap_.empty()) { + const RelaxHeapEntry heap_entry = relax_heap_.front(); + std::pop_heap(relax_heap_.begin(), relax_heap_.end()); + relax_heap_.pop_back(); + + // The slack might have changed since the entry was added. + if (heap_entry.diff > slack) { + trail_indices->push_back(heap_entry.index); + continue; + } + + // Relax, and decide what to do with the new value of index. + slack -= heap_entry.diff; + const int index = integer_trail_[heap_entry.index].prev_trail_index; + + // Same code as in the first block. + if (index < num_vars) continue; + if (heap_entry.coeff > slack) { + trail_indices->push_back(index); + continue; + } + const TrailEntry& entry = integer_trail_[index]; + const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index]; + const int64 diff = CapProd(heap_entry.coeff.value(), + (entry.bound - previous_entry.bound).value()); + if (diff > slack) { + trail_indices->push_back(index); + continue; + } + relax_heap_.push_back({index, heap_entry.coeff, diff}); + std::push_heap(relax_heap_.begin(), relax_heap_.end()); + } + + // If we aborted early because of the slack, we need to push all remaining + // indices back into the reason. + for (const RelaxHeapEntry& entry : relax_heap_) { + trail_indices->push_back(entry.index); + } + relax_heap_.clear(); } void IntegerTrail::RemoveLevelZeroBounds( diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index 5d98f53f78..d37131948c 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -96,6 +96,16 @@ inline IntegerValue FloorRatio(IntegerValue dividend, return result - adjust; } +// Computes result += a * b, and return false iff there is an overflow. +inline bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue* result) { + const int64 prod = CapProd(a.value(), b.value()); + if (prod == kint64min || prod == kint64max) return false; + const int64 add = CapAdd(prod, result->value()); + if (add == kint64min || add == kint64max) return false; + *result = IntegerValue(add); + return true; +} + // Index of an IntegerVariable. // // Each time we create an IntegerVariable we also create its negation. This is @@ -828,6 +838,15 @@ class IntegerTrail : public SatPropagator { mutable gtl::ITIVector tmp_var_to_trail_index_in_queue_; mutable SparseBitset added_variables_; + // Temporary heap used by RelaxLinearReason(); + struct RelaxHeapEntry { + int index; + IntegerValue coeff; + int64 diff; + bool operator<(const RelaxHeapEntry& o) const { return index < o.index; } + }; + mutable std::vector relax_heap_; + // Temporary data used by AppendNewBounds(). mutable SparseBitset tmp_marked_;