enable symmetry detection by default

ortools/sat/cp_model_postsolve.cc

@@ -70,13 +70,12 @@ void PostsolveLinear(const ConstraintProto& ct,
   }
   if (free_vars.empty()) return;
 
-  Domain rhs =
-      ReadDomainFromProto(ct.linear()).AdditionWith(Domain(-fixed_activity));
-
   // Fast track for the most common case.
+  const Domain initial_rhs = ReadDomainFromProto(ct.linear());
   if (free_vars.size() == 1) {
     const int var = free_vars[0];
-    const Domain domain = rhs.InverseMultiplicationBy(free_coeffs[0])
+    const Domain domain = initial_rhs.AdditionWith(Domain(-fixed_activity))
+                              .InverseMultiplicationBy(free_coeffs[0])
                               .IntersectionWith((*domains)[var]);
     const int64 value = prefer_lower_value[var] ? domain.Min() : domain.Max();
     (*domains)[var] = Domain(value);
@@ -86,23 +85,29 @@ void PostsolveLinear(const ConstraintProto& ct,
   // The postsolve code is a bit involved if there is more than one free
   // variable, we have to postsolve them one by one.
   //
-  // Note that if there are some not free variable that are not assigned, the
-  // presolve should have made it sure that whatever the value of those
-  // variables the first variable can always be assigned to satisfy the
-  // constraint. In that case, the MultiplicationBy() below might be inexact,
-  // but that should be okay.
-  std::vector<Domain> to_add;
-  to_add.push_back(Domain(0));
+  // Here we recompute the same domains as during the presolve. Everything is
+  // like if we where substiting the variable one by one:
+  //    terms[i] + fixed_activity \in rhs_domains[i]
+  // In the reverse order.
+  std::vector<Domain> rhs_domains;
+  rhs_domains.push_back(initial_rhs);
   for (int i = 0; i + 1 < free_vars.size(); ++i) {
+    // Note that these should be exactly the same computation as the one done
+    // during presolve and should be exact. However, we have some tests that do
+    // not comply, so we don't check exactness here. Also, as long as we don't
+    // get empty domain below, and the complexity of the domain do not explode
+    // here, we should be fine.
     Domain term = (*domains)[free_vars[i]].MultiplicationBy(-free_coeffs[i]);
-    to_add.push_back(term.AdditionWith(to_add.back()));
+    rhs_domains.push_back(term.AdditionWith(rhs_domains.back()));
   }
   for (int i = free_vars.size() - 1; i >= 0; --i) {
-    // Choose a value for free_vars[i] that fall into rhs + to_add[i].
-    // This will crash if the intersection is empty, but it shouldn't be.
+    // Choose a value for free_vars[i] that fall into rhs_domains[i] -
+    // fixed_activity. This will crash if the intersection is empty, but it
+    // shouldn't be.
     const int var = free_vars[i];
     const int64 coeff = free_coeffs[i];
-    const Domain domain = rhs.AdditionWith(to_add[i])
+    const Domain domain = rhs_domains[i]
+                              .AdditionWith(Domain(-fixed_activity))
                               .InverseMultiplicationBy(coeff)
                               .IntersectionWith((*domains)[var]);
 
@@ -112,12 +117,10 @@ void PostsolveLinear(const ConstraintProto& ct,
     CHECK(!domain.IsEmpty()) << ct.ShortDebugString();
     const int64 value = prefer_lower_value[var] ? domain.Min() : domain.Max();
     (*domains)[var] = Domain(value);
-    rhs = rhs.AdditionWith(Domain(-coeff * value));
 
-    // Only needed in debug.
     fixed_activity += coeff * value;
   }
-  DCHECK(ReadDomainFromProto(ct.linear()).Contains(fixed_activity));
+  DCHECK(initial_rhs.Contains(fixed_activity));
 }
 
 // We assign any non fixed lhs variables to their minimum value. Then we assign

ortools/sat/cp_model_presolve.cc

@@ -1237,7 +1237,9 @@ bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) {
   const int num_vars = ct->linear().vars().size();
   Domain rhs = ReadDomainFromProto(ct->linear());
 
-  // First pass. Process singleton column that are not in the objective.
+  // First pass. Process singleton column that are not in the objective. Note
+  // that for postsolve, it is important that we process them in the same order
+  // in which they will be removed.
   for (int i = 0; i < num_vars; ++i) {
     const int var = ct->linear().vars(i);
     const int64 coeff = ct->linear().coeffs(i);

ortools/sat/cp_model_solver.cc

@@ -1202,7 +1202,8 @@ void LoadBaseModel(const CpModelProto& model_proto,
   // TODO(user): The core algo and symmetries seems to be problematic in some
   // cases. See for instance: neos-691058.mps.gz. This is probably because as
   // we modify the model, our symmetry might be wrong? investigate.
-  if (!parameters.optimize_with_core() && parameters.symmetry_level() > 1) {
+  if (!parameters.optimize_with_core() && parameters.symmetry_level() > 1 &&
+      !parameters.enumerate_all_solutions()) {
     mapping->LoadBooleanSymmetries(model_proto, model);
   }
 

ortools/sat/cp_model_symmetries.cc

@@ -490,6 +490,7 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
   const int initial_ct_index = proto.constraints().size();
   for (int var = 0; var < num_vars; ++var) {
     if (context->IsFixed(var)) continue;
+    if (context->VariableWasRemoved(var)) continue;
     if (context->VariableIsNotUsedAnymore(var)) continue;
 
     const AffineRelation::Relation r = context->GetAffineRelation(var);
@@ -541,10 +542,12 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
   {
     const std::vector<int> orbits = GetOrbits(num_vars, generators);
     std::vector<int> orbit_sizes;
+    int max_orbit_size = 0;
     for (const int rep : orbits) {
       if (rep == -1) continue;
       if (rep >= orbit_sizes.size()) orbit_sizes.resize(rep + 1, 0);
       orbit_sizes[rep]++;
+      max_orbit_size = std::max(max_orbit_size, orbit_sizes[rep]);
     }
 
     std::vector<int> tmp_to_clear;
@@ -592,7 +595,8 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
 
     if (log_info) {
       LOG(INFO) << "Num fixable by intersecting at_most_one with orbits: "
-                << can_be_fixed_to_false.size();
+                << can_be_fixed_to_false.size()
+                << " largest_orbit: " << max_orbit_size;
     }
   }
 
@@ -856,6 +860,24 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
     }
   }
 
+  // If we are left with a set of variable than can all be permuted, lets
+  // break the symmetry by ordering them.
+  if (orbitope.size() == 1) {
+    const int num_cols = orbitope[0].size();
+    for (int i = 0; i + 1 < num_cols; ++i) {
+      // Add orbitope[0][i] >= orbitope[0][i+1].
+      ConstraintProto* ct = context->working_model->add_constraints();
+      ct->mutable_linear()->add_coeffs(1);
+      ct->mutable_linear()->add_vars(orbitope[0][i]);
+      ct->mutable_linear()->add_coeffs(-1);
+      ct->mutable_linear()->add_vars(orbitope[0][i + 1]);
+      ct->mutable_linear()->add_domain(0);
+      ct->mutable_linear()->add_domain(kint64max);
+      context->UpdateRuleStats("symmetry: added symmetry breaking inequality");
+    }
+    context->UpdateNewConstraintsVariableUsage();
+  }
+
   return true;
 }
 

ortools/sat/sat_parameters.proto

@@ -969,7 +969,7 @@ message SatParameters {
   // exploit them. Currently, at level 1 we detect them in presolve and try
   // to fix Booleans. At level 2, we also do some form of dynamic symmetry
   // breaking during search.
-  optional int32 symmetry_level = 183 [default = 0];
+  optional int32 symmetry_level = 183 [default = 2];
 
   // ==========================================================================
   // MIP -> CP-SAT (i.e. IP with integer coeff) conversion parameters that are

ortools/sat/symmetry_util.cc

@@ -104,7 +104,7 @@ std::vector<std::vector<int>> BasicOrbitopeExtraction(
       // not appear at all.
       int num_matches_a = 0;
       int num_matches_b = 0;
-      int last_match_index;
+      int last_match_index = -1;
       for (int j = 0; j < orbitope[i].size(); ++j) {
         if (orbitope[i][j] == a) {
           ++num_matches_a;
@@ -114,6 +114,7 @@ std::vector<std::vector<int>> BasicOrbitopeExtraction(
           last_match_index = j;
         }
       }
+      if (last_match_index == -1) break;
       if (matching_column_index == -1) {
         matching_column_index = last_match_index;
       }