diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc index 9eb6a53f2c..0d7b7bc0a5 100644 --- a/ortools/sat/cp_model_postsolve.cc +++ b/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 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 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 diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 587397843b..002b939ffa 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/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); diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index d79e435537..f9a7a5bbdb 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/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); } diff --git a/ortools/sat/cp_model_symmetries.cc b/ortools/sat/cp_model_symmetries.cc index 522c1a04d0..41b9913644 100644 --- a/ortools/sat/cp_model_symmetries.cc +++ b/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 orbits = GetOrbits(num_vars, generators); std::vector 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 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; } diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 82ae521ebc..f2787e15e5 100644 --- a/ortools/sat/sat_parameters.proto +++ b/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 diff --git a/ortools/sat/symmetry_util.cc b/ortools/sat/symmetry_util.cc index 9c77f1f079..34b9f8e38c 100644 --- a/ortools/sat/symmetry_util.cc +++ b/ortools/sat/symmetry_util.cc @@ -104,7 +104,7 @@ std::vector> 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> BasicOrbitopeExtraction( last_match_index = j; } } + if (last_match_index == -1) break; if (matching_column_index == -1) { matching_column_index = last_match_index; }