[CP-SAT] fix corner cases in presolve; rewrite part of the cut management; automatic fixes

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<uint64_t>(std::abs(coeff)),
-                          static_cast<uint64_t>(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<int64_t>::min() ||
-      constant_factor == std::numeric_limits<int64_t>::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<uint64_t>(std::abs(constant_factor)),
+        static_cast<uint64_t>(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) {

ortools/sat/cp_model_search.cc

@@ -417,8 +417,6 @@ std::function<BooleanOrIntegerLiteral()> 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<int>(result);
 }
 
-}  // namespace
-
 // Note: in flatzinc setting, we know we always have a fixed search defined.
 //
 // Things to try:

ortools/sat/cp_model_search.h

@@ -112,6 +112,10 @@ std::vector<SatParameters> 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
 

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<FeasibilityJumpSolver>(
@@ -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.

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:

ortools/sat/csharp/CpModel.cs

@@ -92,7 +92,7 @@ public class CpModel
 
     /**
      * <summary>
-     * Creates an Boolean variable with given domain.
+     * Creates a Boolean variable with given domain.
      * </summary>
      */
     public BoolVar NewBoolVar(string name)

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;
 }
 

ortools/sat/diffn_util.cc

@@ -16,6 +16,7 @@
 #include <stddef.h>
 
 #include <algorithm>
+#include <ostream>
 #include <utility>
 #include <vector>
 

ortools/sat/integer.cc

@@ -18,6 +18,7 @@
 #include <deque>
 #include <functional>
 #include <limits>
+#include <ostream>
 #include <string>
 #include <utility>
 #include <vector>

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);

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?

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<int64_t>(
+      MathUtil::GCD64(static_cast<uint64_t>(std::abs(coeff)),
+                      static_cast<uint64_t>(std::abs(offset))));
+}
+
 bool PresolveContext::ExpressionIsAffineBoolean(
     const LinearExpressionProto& expr) const {
   if (expr.vars().size() != 1) return false;

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

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.
       //

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
 

ortools/sat/util.cc

@@ -494,7 +494,7 @@ void MaxBoundedSubsetSum::AddMultiples(int64_t coeff, int64_t max_value) {
 
 void MaxBoundedSubsetSum::AddChoicesInternal(absl::Span<const int64_t> 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<const int64_t> 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_) {

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<const int64_t> 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_;

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});