speedup sat internals; fix bug with model with one trivially false XOr constraint

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 \

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<SatSolver>();
-  sat_solver->Propagate();
+  if (!sat_solver->Propagate()) {
+    context->is_unsat = true;
+    return;
+  }
 
   // Probe.
   //

ortools/sat/cp_model_solver.cc

@@ -1496,7 +1496,18 @@ CpSolverResponse SolveCpModelInternal(
       CHECK(SolutionIsFeasible(model_proto,
                                std::vector<int64>(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<SatParameters>()->set_max_number_of_conflicts(
         old_conflict_limit);
   }

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<double>& 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<IntegerValue(IntegerValue)> 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<double> lp_values,
+                        std::vector<IntegerValue> lower_bounds,
+                        std::vector<IntegerValue> 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<bool> 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<IntegerValue> shifts(size, IntegerValue(0));
+  std::vector<bool> 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

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<IntegerValue(IntegerValue)> 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<double> lp_values,
+                        std::vector<IntegerValue> lower_bounds,
+                        std::vector<IntegerValue> 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

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<const IntegerValue> coeffs,
     std::vector<IntegerLiteral>* 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<const IntegerValue> coeffs,
                                      std::vector<int>* 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(

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<IntegerVariable, int> tmp_var_to_trail_index_in_queue_;
   mutable SparseBitset<BooleanVariable> 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<RelaxHeapEntry> relax_heap_;
+
   // Temporary data used by AppendNewBounds().
   mutable SparseBitset<IntegerVariable> tmp_marked_;