[CP-SAT] more work on routing cuts; fix #4582

Author: lperron

ortools/sat/BUILD.bazel

@@ -147,6 +147,31 @@ cc_test(
     ],
 )
 
+cc_library(
+    name = "primary_variables",
+    srcs = ["primary_variables.cc"],
+    hdrs = ["primary_variables.h"],
+    deps = [
+        ":cp_model_cc_proto",
+        ":cp_model_utils",
+        "//ortools/base:stl_util",
+        "//ortools/util:bitset",
+        "@com_google_absl//absl/algorithm:container",
+        "@com_google_absl//absl/container:btree",
+        "@com_google_absl//absl/log:check",
+    ],
+)
+
+cc_test(
+    name = "primary_variables_test",
+    srcs = ["primary_variables_test.cc"],
+    deps = [
+        ":primary_variables",
+        "//ortools/base:gmock_main",
+        "//ortools/base:parse_test_proto",
+    ],
+)
+
 cc_proto_library(
     name = "cp_model_cc_proto",
     deps = [":cp_model_proto"],
@@ -245,6 +270,7 @@ cc_library(
         ":cp_model_cc_proto",
         ":cp_model_utils",
         ":diffn_util",
+        ":primary_variables",
         ":sat_parameters_cc_proto",
         "//ortools/base",
         "//ortools/port:proto_utils",
@@ -254,6 +280,7 @@ cc_library(
         "@com_google_absl//absl/container:btree",
         "@com_google_absl//absl/container:flat_hash_map",
         "@com_google_absl//absl/container:flat_hash_set",
+        "@com_google_absl//absl/flags:flag",
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/meta:type_traits",
         "@com_google_absl//absl/strings",
@@ -623,6 +650,7 @@ cc_library(
         ":parameters_validation",
         ":precedences",
         ":presolve_context",
+        ":primary_variables",
         ":probing",
         ":rins",
         ":routing_cuts",
@@ -680,8 +708,8 @@ cc_library(
     deps = [
         ":cp_model_cc_proto",
         ":cp_model_utils",
+        ":integer",
         ":integer_base",
-        ":intervals",
         ":linear_constraint",
         ":model",
         ":sat_base",
@@ -1875,7 +1903,6 @@ cc_library(
         "@com_google_absl//absl/container:inlined_vector",
         "@com_google_absl//absl/log",
         "@com_google_absl//absl/log:check",
-        "@com_google_absl//absl/meta:type_traits",
         "@com_google_absl//absl/types:span",
     ],
 )
@@ -1885,6 +1912,8 @@ cc_test(
     size = "small",
     srcs = ["precedences_test.cc"],
     deps = [
+        ":cp_model_mapping",
+        ":cp_model_solver_helpers",
         ":integer",
         ":integer_base",
         ":integer_search",
@@ -1893,9 +1922,9 @@ cc_test(
         ":sat_base",
         ":sat_solver",
         "//ortools/base:gmock_main",
+        "//ortools/base:parse_test_proto",
         "//ortools/util:sorted_interval_list",
         "@com_google_absl//absl/container:flat_hash_map",
-        "@com_google_absl//absl/types:span",
     ],
 )
 
@@ -2042,6 +2071,7 @@ cc_library(
         ":sat_solver",
         "//ortools/base",
         "//ortools/graph:strongly_connected_components",
+        "//ortools/util:bitset",
         "//ortools/util:sort",
         "//ortools/util:strong_integers",
         "@com_google_absl//absl/algorithm:container",
@@ -2661,6 +2691,7 @@ cc_library(
     srcs = ["routing_cuts.cc"],
     hdrs = ["routing_cuts.h"],
     deps = [
+        ":clause",
         ":cp_model_cc_proto",
         ":cp_model_utils",
         ":cuts",
@@ -2698,6 +2729,7 @@ cc_test(
     name = "routing_cuts_test",
     srcs = ["routing_cuts_test.cc"],
     deps = [
+        ":clause",
         ":cp_model",
         ":cuts",
         ":integer",

ortools/sat/all_different.cc

@@ -424,7 +424,7 @@ AllDifferentBoundsPropagator::AllDifferentBoundsPropagator(
   const int capacity = expressions.size() + 2;
   index_to_start_index_.resize(capacity);
   index_to_end_index_.resize(capacity);
-  index_is_present_.resize(capacity, false);
+  index_is_present_.Resize(capacity);
   index_to_expr_.resize(capacity, kNoIntegerVariable);
 
   for (int i = 0; i < expressions.size(); ++i) {
@@ -535,11 +535,7 @@ bool AllDifferentBoundsPropagator::PropagateLowerBoundsInternal(
   // Make sure we change our base_ so that GetIndex() fit in our buffers.
   base_ = min_lb - IntegerValue(1);
 
-  // Sparse cleaning of index_is_present_.
-  for (const int i : indices_to_clear_) {
-    index_is_present_[i] = false;
-  }
-  indices_to_clear_.clear();
+  index_is_present_.ResetAllToFalse();
 
   // Sort bounds by increasing ub.
   std::sort(bounds.begin(), bounds.end(),
@@ -600,15 +596,14 @@ bool AllDifferentBoundsPropagator::PropagateLowerBoundsInternal(
     {
       index_to_start_index_[new_index] = start_index;
       index_to_expr_[new_index] = expr;
-      index_is_present_[new_index] = true;
-      indices_to_clear_.push_back(new_index);
+      index_is_present_.Set(new_index);
     }
 
     // In most situation, we cannot have a conflict now, because it should have
     // been detected before by pushing an interval lower bound past its upper
     // bound. However, it is possible that when we push one bound, other bounds
     // change. So if the upper bound is smaller than the current interval end,
-    // we abort so that the conflit reason will be better on the next call to
+    // we abort so that the conflict reason will be better on the next call to
     // the propagator.
     const IntegerValue end = GetValue(end_index);
     if (end > integer_trail_->UpperBound(expr)) return true;

ortools/sat/all_different.h

@@ -25,6 +25,7 @@
 #include "ortools/sat/integer_base.h"
 #include "ortools/sat/model.h"
 #include "ortools/sat/sat_base.h"
+#include "ortools/util/bitset.h"
 #include "ortools/util/strong_integers.h"
 
 namespace operations_research {
@@ -218,10 +219,9 @@ class AllDifferentBoundsPropagator : public PropagatorInterface {
 
   // Non-consecutive intervals related data-structures.
   IntegerValue base_;
-  std::vector<int> indices_to_clear_;
   std::vector<int> index_to_start_index_;
   std::vector<int> index_to_end_index_;
-  std::vector<bool> index_is_present_;
+  SparseBitset<int> index_is_present_;
   std::vector<AffineExpression> index_to_expr_;
 
   // Temporary integer reason.

ortools/sat/clause.cc

@@ -346,7 +346,7 @@ void ClauseManager::AttachAllClauses() {
   if (all_clauses_are_attached_) return;
   all_clauses_are_attached_ = true;
 
-  needs_cleaning_.ClearAll();  // This doesn't resize it.
+  needs_cleaning_.ResetAllToFalse();  // This doesn't resize it.
   watchers_on_false_.resize(needs_cleaning_.size().value());
 
   DeleteRemovedClauses();

ortools/sat/cp_model_checker.cc

@@ -27,6 +27,7 @@
 #include "absl/container/btree_map.h"
 #include "absl/container/flat_hash_map.h"
 #include "absl/container/flat_hash_set.h"
+#include "absl/flags/flag.h"
 #include "absl/log/check.h"
 #include "absl/meta/type_traits.h"
 #include "absl/strings/str_cat.h"
@@ -36,10 +37,15 @@
 #include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/cp_model_utils.h"
 #include "ortools/sat/diffn_util.h"
+#include "ortools/sat/primary_variables.h"
 #include "ortools/sat/sat_parameters.pb.h"
 #include "ortools/util/saturated_arithmetic.h"
 #include "ortools/util/sorted_interval_list.h"
 
+ABSL_FLAG(bool, cp_model_check_dependent_variables, false,
+          "When true, check that solutions can be computed only from their "
+          "free variables.");
+
 namespace operations_research {
 namespace sat {
 namespace {
@@ -1957,6 +1963,18 @@ bool SolutionIsFeasible(const CpModelProto& model,
     VLOG(2) << "Checker inner objective = " << inner_objective;
     VLOG(2) << "Checker scaled objective = " << scaled_objective;
   }
+  if (absl::GetFlag(FLAGS_cp_model_check_dependent_variables)) {
+    const VariableRelationships relationships =
+        ComputeVariableRelationships(model);
+    std::vector<int64_t> all_variables(variable_values.begin(),
+                                       variable_values.end());
+    for (const int var : relationships.secondary_variables) {
+      all_variables[var] = -999999;  // Those values should be overwritten.
+    }
+    CHECK(ComputeAllVariablesFromPrimaryVariables(model, relationships,
+                                                  &all_variables));
+    CHECK(absl::MakeSpan(all_variables) == variable_values);
+  }
 
   return true;
 }

ortools/sat/cp_model_lns.cc

@@ -2426,6 +2426,10 @@ Neighborhood RectanglesPackingRelaxOneNeighborhoodGenerator::Generate(
   Neighborhood neighborhood =
       helper_.FixGivenVariables(initial_solution, variables_to_freeze);
 
+  neighborhood.is_simple = false;
+  neighborhood.is_reduced = true;
+  neighborhood.variables_that_can_be_fixed_to_local_optimum.clear();
+
   // The call above add the relaxed variables to the neighborhood using the
   // current bounds at level 0. For big problems, this might create a hard model
   // with a large complicated landscape of fixed boxes with a lot of potential

ortools/sat/cp_model_lns.h

@@ -59,8 +59,8 @@ struct Neighborhood {
   // True if neighborhood generator was able to generate a neighborhood.
   bool is_generated = false;
 
-  // True if an optimal solution to the neighborhood is also an optimal solution
-  // to the original model.
+  // False if an optimal solution to the neighborhood is also an optimal
+  // solution to the original model.
   bool is_reduced = false;
 
   // True if this neighborhood was just obtained by fixing some variables.

ortools/sat/cp_model_loader.cc

@@ -74,27 +74,6 @@ std::vector<int64_t> ValuesFromProto(const Values& values) {
   return std::vector<int64_t>(values.begin(), values.end());
 }
 
-void ComputeLinearBounds(const LinearConstraintProto& proto,
-                         CpModelMapping* mapping, IntegerTrail* integer_trail,
-                         int64_t* sum_min, int64_t* sum_max) {
-  *sum_min = 0;
-  *sum_max = 0;
-
-  for (int i = 0; i < proto.vars_size(); ++i) {
-    const int64_t coeff = proto.coeffs(i);
-    const IntegerVariable var = mapping->Integer(proto.vars(i));
-    const int64_t lb = integer_trail->LowerBound(var).value();
-    const int64_t ub = integer_trail->UpperBound(var).value();
-    if (coeff >= 0) {
-      (*sum_min) += coeff * lb;
-      (*sum_max) += coeff * ub;
-    } else {
-      (*sum_min) += coeff * ub;
-      (*sum_max) += coeff * lb;
-    }
-  }
-}
-
 // We check if the constraint is a sum(ax * xi) == value.
 bool ConstraintIsEq(const LinearConstraintProto& proto) {
   return proto.domain_size() == 2 && proto.domain(0) == proto.domain(1);
@@ -104,9 +83,8 @@ bool ConstraintIsEq(const LinearConstraintProto& proto) {
 bool ConstraintIsNEq(const LinearConstraintProto& proto,
                      CpModelMapping* mapping, IntegerTrail* integer_trail,
                      int64_t* single_value) {
-  int64_t sum_min = 0;
-  int64_t sum_max = 0;
-  ComputeLinearBounds(proto, mapping, integer_trail, &sum_min, &sum_max);
+  const auto [sum_min, sum_max] =
+      mapping->ComputeMinMaxActivity(proto, integer_trail);
 
   const Domain complement =
       Domain(sum_min, sum_max)
@@ -1280,32 +1258,6 @@ void LoadLinearConstraint(const ConstraintProto& ct, Model* m) {
     max_sum += std::max(term_a, term_b);
   }
 
-  // Load conditional precedences and always true binary relations.
-  const SatParameters& params = *m->GetOrCreate<SatParameters>();
-  if (ct.enforcement_literal().size() <= 1 && vars.size() <= 2) {
-    // To avoid overflow in the code below, we tighten the bounds.
-    int64_t rhs_min = ct.linear().domain(0);
-    int64_t rhs_max = ct.linear().domain(ct.linear().domain().size() - 1);
-    rhs_min = std::max(rhs_min, min_sum.value());
-    rhs_max = std::min(rhs_max, max_sum.value());
-
-    auto* repository = m->GetOrCreate<BinaryRelationRepository>();
-    if (ct.enforcement_literal().empty()) {
-      if (vars.size() == 2) {
-        repository->Add(Literal(kNoLiteralIndex), {vars[0], coeffs[0]},
-                        {vars[1], coeffs[1]}, rhs_min, rhs_max);
-      }
-    } else {
-      const Literal lit = mapping->Literal(ct.enforcement_literal(0));
-      if (vars.size() == 1) {
-        repository->Add(lit, {vars[0], coeffs[0]}, {}, rhs_min, rhs_max);
-      } else if (vars.size() == 2) {
-        repository->Add(lit, {vars[0], coeffs[0]}, {vars[1], coeffs[1]},
-                        rhs_min, rhs_max);
-      }
-    }
-  }
-
   // Load precedences.
   if (!HasEnforcementLiteral(ct)) {
     auto* precedences = m->GetOrCreate<PrecedenceRelations>();
@@ -1366,6 +1318,7 @@ void LoadLinearConstraint(const ConstraintProto& ct, Model* m) {
     }
   }
 
+  const SatParameters& params = *m->GetOrCreate<SatParameters>();
   const IntegerValue domain_size_limit(
       params.max_domain_size_when_encoding_eq_neq_constraints());
   if (ct.linear().vars_size() == 2 && !integer_trail->IsFixed(vars[0]) &&

ortools/sat/cp_model_mapping.h

@@ -14,6 +14,8 @@
 #ifndef OR_TOOLS_SAT_CP_MODEL_MAPPING_H_
 #define OR_TOOLS_SAT_CP_MODEL_MAPPING_H_
 
+#include <cstdint>
+#include <utility>
 #include <vector>
 
 #include "absl/container/flat_hash_set.h"
@@ -22,8 +24,8 @@
 #include "ortools/base/strong_vector.h"
 #include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/cp_model_utils.h"
+#include "ortools/sat/integer.h"
 #include "ortools/sat/integer_base.h"
-#include "ortools/sat/intervals.h"
 #include "ortools/sat/linear_constraint.h"
 #include "ortools/sat/model.h"
 #include "ortools/sat/sat_base.h"
@@ -189,6 +191,29 @@ class CpModelMapping {
     return CanonicalizeExpr(expr);
   }
 
+  // Returns the min/max activity of the linear constraint under the current
+  // integer_trail bounds.
+  std::pair<int64_t, int64_t> ComputeMinMaxActivity(
+      const LinearConstraintProto& proto, IntegerTrail* integer_trail) {
+    int64_t sum_min = 0;
+    int64_t sum_max = 0;
+
+    for (int i = 0; i < proto.vars_size(); ++i) {
+      const int64_t coeff = proto.coeffs(i);
+      const IntegerVariable var = this->Integer(proto.vars(i));
+      const int64_t lb = integer_trail->LowerBound(var).value();
+      const int64_t ub = integer_trail->UpperBound(var).value();
+      if (coeff >= 0) {
+        sum_min += coeff * lb;
+        sum_max += coeff * ub;
+      } else {
+        sum_min += coeff * ub;
+        sum_max += coeff * lb;
+      }
+    }
+    return {sum_min, sum_max};
+  }
+
   // For logging only, these are not super efficient.
   int NumIntegerVariables() const {
     int result = 0;