diff --git a/makefiles/Makefile.gen.mk b/makefiles/Makefile.gen.mk index fa0701613e..a0d2489b83 100644 --- a/makefiles/Makefile.gen.mk +++ b/makefiles/Makefile.gen.mk @@ -1320,6 +1320,7 @@ SAT_DEPS = \ $(SRC_DIR)/ortools/sat/subsolver.h \ $(SRC_DIR)/ortools/sat/swig_helper.h \ $(SRC_DIR)/ortools/sat/symmetry.h \ + $(SRC_DIR)/ortools/sat/symmetry_util.h \ $(SRC_DIR)/ortools/sat/synchronization.h \ $(SRC_DIR)/ortools/sat/table.h \ $(SRC_DIR)/ortools/sat/theta_tree.h \ @@ -1386,6 +1387,7 @@ SAT_LIB_OBJS = \ $(OBJ_DIR)/sat/simplification.$O \ $(OBJ_DIR)/sat/subsolver.$O \ $(OBJ_DIR)/sat/symmetry.$O \ + $(OBJ_DIR)/sat/symmetry_util.$O \ $(OBJ_DIR)/sat/synchronization.$O \ $(OBJ_DIR)/sat/table.$O \ $(OBJ_DIR)/sat/theta_tree.$O \ @@ -2634,6 +2636,17 @@ objs/sat/symmetry.$O: ortools/sat/symmetry.cc ortools/sat/symmetry.h \ ortools/base/timer.h ortools/base/basictypes.h | $(OBJ_DIR)/sat $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Ssymmetry.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Ssymmetry.$O +objs/sat/symmetry_util.$O: ortools/sat/symmetry_util.cc ortools/sat/symmetry_util.h \ + ortools/algorithms/sparse_permutation.h ortools/base/logging.h \ + ortools/base/commandlineflags.h ortools/base/integral_types.h \ + ortools/base/log_severity.h ortools/base/logging_export.h \ + ortools/base/macros.h ortools/base/vlog_is_on.h \ + ortools/base/strong_vector.h ortools/base/int_type.h \ + ortools/sat/sat_base.h ortools/sat/model.h ortools/base/map_util.h \ + ortools/base/typeid.h ortools/util/bitset.h ortools/util/stats.h \ + ortools/base/timer.h ortools/base/basictypes.h | $(OBJ_DIR)/sat + $(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Ssymmetry_util.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Ssymmetry_util.$O + objs/sat/synchronization.$O: ortools/sat/synchronization.cc \ ortools/sat/synchronization.h ortools/base/integral_types.h \ ortools/base/logging.h ortools/base/commandlineflags.h \ @@ -4950,4 +4963,3 @@ $(GEN_DIR)/ortools/gscip/gscip.pb.h: \ $(OBJ_DIR)/gscip/gscip.pb.$O: \ $(GEN_DIR)/ortools/gscip/gscip.pb.cc | $(OBJ_DIR)/gscip $(CCC) $(CFLAGS) -c $(GEN_PATH)$Sortools$Sgscip$Sgscip.pb.cc $(OBJ_OUT)$(OBJ_DIR)$Sgscip$Sgscip.pb.$O - diff --git a/ortools/sat/BUILD b/ortools/sat/BUILD index d18c0c5b7f..a14406b9f4 100644 --- a/ortools/sat/BUILD +++ b/ortools/sat/BUILD @@ -305,6 +305,7 @@ cc_library( ":cp_model_expand", ":cp_model_loader", ":cp_model_objective", + ":cp_model_symmetries", ":cp_model_utils", ":presolve_context", ":presolve_util", @@ -565,6 +566,17 @@ cc_library( ], ) +cc_library( + name = "symmetry_util", + srcs = ["symmetry_util.cc"], + hdrs = ["symmetry_util.h"], + deps = [ + "//ortools/base", + "//ortools/algorithms:dynamic_partition", + "//ortools/algorithms:sparse_permutation", + ], +) + cc_library( name = "var_domination", srcs = ["var_domination.cc"], @@ -1307,7 +1319,9 @@ cc_library( deps = [ ":cp_model_cc_proto", ":cp_model_utils", + ":presolve_context", ":sat_parameters_cc_proto", + ":symmetry_util", "//ortools/algorithms:find_graph_symmetries", "//ortools/algorithms:sparse_permutation", "//ortools/base:hash", diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index b9476d0815..194cd3aa8b 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -43,6 +43,7 @@ #include "ortools/sat/cp_model_expand.h" #include "ortools/sat/cp_model_loader.h" #include "ortools/sat/cp_model_objective.h" +#include "ortools/sat/cp_model_symmetries.h" #include "ortools/sat/cp_model_utils.h" #include "ortools/sat/presolve_util.h" #include "ortools/sat/probing.h" @@ -341,34 +342,57 @@ bool CpModelPresolver::PresolveBoolAnd(ConstraintProto* ct) { return changed; } -bool CpModelPresolver::PresolveAtMostOne(ConstraintProto* ct) { - if (context_->ModelIsUnsat()) return false; - CHECK(!HasEnforcementLiteral(*ct)); +bool CpModelPresolver::PresolveAtMostOrExactlyOne(ConstraintProto* ct) { + const bool is_at_most_one = + ct->constraint_case() == ConstraintProto::kAtMostOne; + const std::string name = is_at_most_one ? "at_most_one: " : "exactly_one: "; + auto* literals = is_at_most_one + ? ct->mutable_at_most_one()->mutable_literals() + : ct->mutable_exactly_one()->mutable_literals(); - // An at most one with just one literal is always satisfied. - if (ct->at_most_one().literals_size() == 1) { - context_->UpdateRuleStats("at_most_one: size one"); - return RemoveConstraint(ct); - } - - // Fix to false any duplicate literals. - std::sort(ct->mutable_at_most_one()->mutable_literals()->begin(), - ct->mutable_at_most_one()->mutable_literals()->end()); - int previous = kint32max; - for (const int literal : ct->at_most_one().literals()) { - if (literal == previous) { + // Deal with duplicate variable reference. + context_->tmp_literal_set.clear(); + for (const int literal : *literals) { + if (context_->tmp_literal_set.contains(literal)) { if (!context_->SetLiteralToFalse(literal)) return true; - context_->UpdateRuleStats("at_most_one: duplicate literals"); + context_->UpdateRuleStats(absl::StrCat(name, "duplicate literals")); } - previous = literal; + if (context_->tmp_literal_set.contains(NegatedRef(literal))) { + int num_positive = 0; + int num_negative = 0; + for (const int other : *literals) { + if (PositiveRef(other) != PositiveRef(literal)) { + if (!context_->SetLiteralToFalse(other)) return true; + context_->UpdateRuleStats(absl::StrCat(name, "x and not(x)")); + } else { + if (other == literal) { + ++num_positive; + } else { + ++num_negative; + } + } + } + + // This is tricky for the case where the at most one reduce to (lit, + // not(lit), not(lit)) for instance. + if (num_positive > 1 && !context_->SetLiteralToFalse(literal)) { + return true; + } + if (num_negative > 1 && !context_->SetLiteralToTrue(literal)) { + return true; + } + return RemoveConstraint(ct); + } + context_->tmp_literal_set.insert(literal); } + // Remove fixed variables. bool changed = false; context_->tmp_literals.clear(); - for (const int literal : ct->at_most_one().literals()) { + for (const int literal : *literals) { if (context_->LiteralIsTrue(literal)) { - context_->UpdateRuleStats("at_most_one: satisfied"); - for (const int other : ct->at_most_one().literals()) { + context_->UpdateRuleStats(absl::StrCat(name, "satisfied")); + for (const int other : *literals) { if (other != literal) { if (!context_->SetLiteralToFalse(other)) return true; } @@ -385,102 +409,71 @@ bool CpModelPresolver::PresolveAtMostOne(ConstraintProto* ct) { // code in var_domination.cc, so maybe we don't need it here. if (context_->VariableIsUniqueAndRemovable(literal) || context_->VariableWithCostIsUniqueAndRemovable(literal)) { - context_->UpdateRuleStats("TODO at_most_one: singleton"); + context_->UpdateRuleStats("TODO [exactly/at_most]_one: singleton"); } context_->tmp_literals.push_back(literal); } - if (context_->tmp_literals.empty()) { - context_->UpdateRuleStats("at_most_one: all false"); - return RemoveConstraint(ct); - } if (changed) { - ct->mutable_at_most_one()->mutable_literals()->Clear(); + literals->Clear(); for (const int lit : context_->tmp_literals) { - ct->mutable_at_most_one()->add_literals(lit); + literals->Add(lit); } - context_->UpdateRuleStats("at_most_one: removed literals"); + context_->UpdateRuleStats(absl::StrCat(name, "removed literals")); } return changed; } -// TODO(user): This is really close to the at most one code. Try to remove -// duplication. +bool CpModelPresolver::PresolveAtMostOne(ConstraintProto* ct) { + if (context_->ModelIsUnsat()) return false; + CHECK(!HasEnforcementLiteral(*ct)); + const bool changed = PresolveAtMostOrExactlyOne(ct); + if (ct->constraint_case() != ConstraintProto::kAtMostOne) return changed; + + // Size zero: ok. + const auto& literals = ct->at_most_one().literals(); + if (literals.empty()) { + context_->UpdateRuleStats("at_most_one: empty or all false"); + return RemoveConstraint(ct); + } + + // Size one: always satisfied. + if (literals.size() == 1) { + context_->UpdateRuleStats("at_most_one: size one"); + return RemoveConstraint(ct); + } + + return changed; +} + bool CpModelPresolver::PresolveExactlyOne(ConstraintProto* ct) { if (context_->ModelIsUnsat()) return false; CHECK(!HasEnforcementLiteral(*ct)); + const bool changed = PresolveAtMostOrExactlyOne(ct); + if (ct->constraint_case() != ConstraintProto::kExactlyOne) return changed; - { - // Size one: fix variable. - const auto& literals = ct->exactly_one().literals(); - if (literals.size() == 1) { - context_->UpdateRuleStats("exactly_one: size one"); - if (!context_->SetLiteralToTrue(literals[0])) return false; - return RemoveConstraint(ct); - } - - // Size two: Equivalence. - if (literals.size() == 2) { - context_->UpdateRuleStats("exactly_one: size two"); - context_->StoreBooleanEqualityRelation(literals[0], - NegatedRef(literals[1])); - return RemoveConstraint(ct); - } - } - - // Fix to false any duplicate literals. - std::sort(ct->mutable_exactly_one()->mutable_literals()->begin(), - ct->mutable_exactly_one()->mutable_literals()->end()); - int previous = kint32max; - for (const int literal : ct->exactly_one().literals()) { - if (literal == previous) { - if (!context_->SetLiteralToFalse(literal)) return true; - context_->UpdateRuleStats("exactly_one: duplicate literals"); - } - previous = literal; - } - - // Remove fixed variables. - bool changed = false; - context_->tmp_literals.clear(); - for (const int literal : ct->exactly_one().literals()) { - if (context_->LiteralIsTrue(literal)) { - context_->UpdateRuleStats("exactly_one: satisfied"); - for (const int other : ct->exactly_one().literals()) { - if (other != literal) { - if (!context_->SetLiteralToFalse(other)) return true; - } - } - return RemoveConstraint(ct); - } - - if (context_->LiteralIsFalse(literal)) { - changed = true; - continue; - } - - // TODO(user): We could remove the variable, and leave an at most one on - // the other. Even if there is an objective on this variable, we could - // transfer the cost to the other. - if (context_->VariableIsUniqueAndRemovable(literal) || - context_->VariableWithCostIsUniqueAndRemovable(literal)) { - context_->UpdateRuleStats("TODO exactly_one: singleton"); - } - - context_->tmp_literals.push_back(literal); - } - if (context_->tmp_literals.empty()) { + // Size zero: UNSAT. + const auto& literals = ct->exactly_one().literals(); + if (literals.empty()) { return context_->NotifyThatModelIsUnsat("exactly_one: empty or all false"); } - if (changed) { - ct->mutable_exactly_one()->mutable_literals()->Clear(); - for (const int lit : context_->tmp_literals) { - ct->mutable_exactly_one()->add_literals(lit); - } - context_->UpdateRuleStats("exactly_one: removed literals"); + // Size one: fix variable. + if (literals.size() == 1) { + context_->UpdateRuleStats("exactly_one: size one"); + if (!context_->SetLiteralToTrue(literals[0])) return false; + return RemoveConstraint(ct); } + + // Size two: Equivalence. + if (literals.size() == 2) { + context_->UpdateRuleStats("exactly_one: size two"); + context_->StoreBooleanEqualityRelation(literals[0], + NegatedRef(literals[1])); + return RemoveConstraint(ct); + } + return changed; } @@ -4348,6 +4341,7 @@ void CpModelPresolver::TransformIntoMaxCliques() { if (ct->enforcement_literal().size() != 1) continue; const Literal enforcement = convert(ct->enforcement_literal(0)); for (const int ref : ct->bool_and().literals()) { + if (ref == ct->enforcement_literal(0)) continue; cliques.push_back({enforcement, convert(ref).Negated()}); } if (RemoveConstraint(ct)) { @@ -4401,6 +4395,9 @@ void CpModelPresolver::TransformIntoMaxCliques() { NegatedRef(literal.Variable().value())); } } + + // Make sure we do not have duplicate variable reference. + PresolveAtMostOne(ct); } context_->UpdateNewConstraintsVariableUsage(); if (num_new_cliques != num_old_cliques) { @@ -5219,8 +5216,7 @@ bool CpModelPresolver::Presolve() { // TODO(user): do that and the pure-SAT part below more than once. if (context_->params().cp_model_probing_level() > 0) { - if (context_->time_limit() == nullptr || - !context_->time_limit()->LimitReached()) { + if (!context_->time_limit()->LimitReached()) { Probe(); PresolveToFixPoint(); } @@ -5230,8 +5226,7 @@ bool CpModelPresolver::Presolve() { // Note that because this can only remove/fix variable not used in the other // part of the problem, there is no need to redo more presolve afterwards. if (context_->params().cp_model_use_sat_presolve()) { - if (context_->time_limit() == nullptr || - !context_->time_limit()->LimitReached()) { + if (!context_->time_limit()->LimitReached()) { PresolvePureSatPart(); } } @@ -5258,9 +5253,16 @@ bool CpModelPresolver::Presolve() { if (iter == 0) TransformIntoMaxCliques(); + // TODO(user): Decide where is the best place for this. Fow now we do it + // after max clique to get all the bool_and converted to at most ones. + if (context_->params().detect_symmetries() && !context_->ModelIsUnsat() && + !context_->time_limit()->LimitReached() && + !context_->keep_all_feasible_solutions) { + DetectAndExploitSymmetriesInPresolve(context_); + } + // Process set packing, partitioning and covering constraint. - if (context_->time_limit() == nullptr || - !context_->time_limit()->LimitReached()) { + if (!context_->time_limit()->LimitReached()) { ProcessSetPPC(); ExtractBoolAnd(); diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h index 8234afe164..22f2d6298d 100644 --- a/ortools/sat/cp_model_presolve.h +++ b/ortools/sat/cp_model_presolve.h @@ -112,8 +112,11 @@ class CpModelPresolver { bool PresolveLinMax(ConstraintProto* ct); bool PresolveIntAbs(ConstraintProto* ct); bool PresolveBoolXor(ConstraintProto* ct); + + bool PresolveAtMostOrExactlyOne(ConstraintProto* ct); bool PresolveAtMostOne(ConstraintProto* ct); bool PresolveExactlyOne(ConstraintProto* ct); + bool PresolveBoolAnd(ConstraintProto* ct); bool PresolveBoolOr(ConstraintProto* ct); bool PresolveEnforcementLiteral(ConstraintProto* ct); diff --git a/ortools/sat/cp_model_symmetries.cc b/ortools/sat/cp_model_symmetries.cc index 2deab4028f..d9af14ae38 100644 --- a/ortools/sat/cp_model_symmetries.cc +++ b/ortools/sat/cp_model_symmetries.cc @@ -22,6 +22,7 @@ #include "ortools/base/hash.h" #include "ortools/base/map_util.h" #include "ortools/sat/cp_model_utils.h" +#include "ortools/sat/symmetry_util.h" namespace operations_research { namespace sat { @@ -49,6 +50,8 @@ class IdGenerator { return gtl::LookupOrInsert(&id_map_, key, id_map_.size()); } + int NextFreeId() const { return id_map_.size(); } + private: absl::flat_hash_map, int, VectorHash> id_map_; }; @@ -116,9 +119,11 @@ std::unique_ptr GenerateGraphForSymmetryDetection( std::vector tmp_key; for (int v = 0; v < num_variables; ++v) { - const IntegerVariableProto& variable = problem.variables(v); tmp_key = {VARIABLE_NODE, objective_by_var[v]}; - Append(variable.domain(), &tmp_key); + Append(problem.variables(v).domain(), &tmp_key); + + // Note that the code rely on the fact that the index of a VARIABLE_NODE + // type is the same as the variable index. CHECK_EQ(v, new_node(tmp_key)); // Make sure the graph contains all the variable nodes, even if no edges are @@ -171,6 +176,8 @@ std::unique_ptr GenerateGraphForSymmetryDetection( std::vector key = {CONSTRAINT_NODE, constraint.constraint_case()}; switch (constraint.constraint_case()) { + case ConstraintProto::CONSTRAINT_NOT_SET: + break; case ConstraintProto::kLinear: { Append(constraint.linear().domain(), &key); CHECK_EQ(constraint_node, new_node(key)); @@ -278,6 +285,33 @@ std::unique_ptr GenerateGraphForSymmetryDetection( graph->Build(); DCHECK_EQ(graph->num_nodes(), initial_equivalence_classes->size()); + + // Because this code is running during presolve, a lot a variable might have + // no edges. We do not want to detect symmetries between these. + // + // Note that this code forces us to "densify" the ids aftewards because the + // symmetry detection code relies on that. + // + // TODO(user): It will probably be more efficient to not even create these + // nodes, but we will need a mapping to know the variable <-> node index. + int next_id = id_generator.NextFreeId(); + for (int i = 0; i < num_variables; ++i) { + if ((*graph)[i].empty()) { + (*initial_equivalence_classes)[i] = next_id++; + } + } + + // Densify ids. + int id = 0; + std::vector mapping(next_id, -1); + for (int& ref : *initial_equivalence_classes) { + if (mapping[ref] == -1) { + ref = mapping[ref] = id++; + } else { + ref = mapping[ref]; + } + } + return graph; } } // namespace @@ -298,8 +332,8 @@ void FindCpModelSymmetries( if (graph == nullptr) return; if (log_info) { - LOG(INFO) << "Graph has " << graph->num_nodes() << " nodes and " - << graph->num_arcs() / 2 << " edges."; + LOG(INFO) << "Graph for symmetry has " << graph->num_nodes() + << " nodes and " << graph->num_arcs() / 2 << " edges."; } if (graph->num_nodes() == 0) return; @@ -350,7 +384,7 @@ void FindCpModelSymmetries( } generators->resize(num_generators); average_support_size /= num_generators; - if (log_info) { + if (log_info && num_generators > 0) { LOG(INFO) << "# of generators: " << num_generators; LOG(INFO) << "Average support size: " << average_support_size; if (num_duplicate_constraints > 0) { @@ -360,5 +394,297 @@ void FindCpModelSymmetries( } } +bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) { + const SatParameters& params = context->params(); + const bool log_info = params.log_search_progress() || VLOG_IS_ON(1); + const CpModelProto& proto = *context->working_model; + + // We need to make sure the proto is up to date before computing symmetries! + context->WriteObjectiveToProto(); + const int num_vars = proto.variables_size(); + for (int i = 0; i < num_vars; ++i) { + FillDomainInProto(context->DomainOf(i), + context->working_model->mutable_variables(i)); + } + + // Tricky: the equivalence relation are not part of the proto. + // We thus add them temporarily to compute the symmetry. + int64 num_added = 0; + const int initial_ct_index = proto.constraints().size(); + for (int var = 0; var < num_vars; ++var) { + if (context->IsFixed(var)) continue; + if (context->VariableIsNotUsedAnymore(var)) continue; + + const AffineRelation::Relation r = context->GetAffineRelation(var); + if (r.representative == var) continue; + + ++num_added; + ConstraintProto* ct = context->working_model->add_constraints(); + auto* arg = ct->mutable_linear(); + arg->add_vars(var); + arg->add_coeffs(1); + arg->add_vars(r.representative); + arg->add_coeffs(-r.coeff); + arg->add_domain(r.offset); + arg->add_domain(r.offset); + } + + std::vector> generators; + FindCpModelSymmetries(params, proto, &generators, /*time_limit_seconds=*/1.0); + + // Remove temporary affine relation. + context->working_model->mutable_constraints()->DeleteSubrange( + initial_ct_index, num_added); + + if (generators.empty()) return true; + + // Orbitope approach. + // + // This is basically the same as the generic approach, but because of the + // extra structure, computing the orbit of any stabilizer subgroup is easy. + // We look for orbits intersecting at most one constraints, so we can break + // symmetry by fixing variables. + // + // TODO(user): The same effect could be achieved by adding symmetry breaking + // constraints of the form "a >= b " between Booleans and let the presolve do + // the reduction. This might be less code, but it is also less efficient. + // Similarly, when we cannot just fix variables to break symmetries, we could + // add these constraints, but it is unclear if we should do it all the time or + // not. + // + // TODO(user): code the generic approach with orbits and stabilizer. + std::vector> orbitope = BasicOrbitopeExtraction(generators); + if (orbitope.empty()) return true; + + if (log_info) { + LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x " + << orbitope[0].size(); + } + + // Collect the at most ones. + // + // Note(user): This relies on the fact that the pointer remain stable when + // we adds new constraints. It should be the case, but it is a bit unsafe. + // On the other hand it is annoying to deal with both cases below. + std::vector*> at_most_ones; + for (int i = 0; i < proto.constraints_size(); ++i) { + if (proto.constraints(i).constraint_case() == ConstraintProto::kAtMostOne) { + at_most_ones.push_back(&proto.constraints(i).at_most_one().literals()); + } + if (proto.constraints(i).constraint_case() == + ConstraintProto::kExactlyOne) { + at_most_ones.push_back(&proto.constraints(i).exactly_one().literals()); + } + } + + // This will always be kept all zero after usage. + std::vector tmp_to_clear; + std::vector tmp_sizes(num_vars, 0); + std::vector tmp_num_positive(num_vars, 0); + + // TODO(user): The code below requires that no variable appears twice in the + // same at most one. In particular lit and not(lit) cannot appear in the same + // at most one. + for (const google::protobuf::RepeatedField* literals : at_most_ones) { + for (const int lit : *literals) { + const int var = PositiveRef(lit); + CHECK_NE(tmp_sizes[var], 1); + tmp_sizes[var] = 1; + } + for (const int lit : *literals) { + tmp_sizes[PositiveRef(lit)] = 0; + } + } + + while (!orbitope.empty() && !orbitope[0].empty()) { + const std::vector orbits = GetOrbitopeOrbits(num_vars, orbitope); + + // Because in the orbitope case, we have a full symmetry group of the + // columns, we can infer more than just using the orbits under a general + // permutation group. If an at most one contains two variables from the + // orbit, we can infer: + // 1/ If the two variables appear positively, then there is an at most one + // on the full orbit, and we can set n - 1 variables to zero to break the + // symmetry. + // 2/ If the two variables appear negatively, then the opposite situation + // arise and there is at most one zero on the orbit, we can set n - 1 + // variables to one. + // 3/ If two literals of opposite sign appear, then the only possibility + // for the orbit are all at one or all at zero, thus we can mark all + // variables as equivalent. + // + // These property comes from the fact that when we permute a line of the + // orbitope in any way, then the position than ends up in the at most one + // must never be both at one. + // + // Note that 1/ can be done without breaking any symmetry, but for 2/ and 3/ + // by choosing which variable is not fixed, we will break some symmetry, and + // we will need to update the orbitope to stabilize this choice before + // continuing. + // + // TODO(user): for 2/ and 3/ we could add an at most one constraint on the + // full orbit if it is not already there! + // + // Note(user): On the miplib, only 1/ happens currently. Not sure with LNS + // though. + std::vector all_equivalent_rows(orbitope.size(), false); + std::vector at_most_one_one(orbitope.size(), false); + std::vector at_most_one_zero(orbitope.size(), false); + + for (const google::protobuf::RepeatedField* literals : + at_most_ones) { + tmp_to_clear.clear(); + int num_in_intersections = 0; + for (const int literal : *literals) { + if (context->IsFixed(literal)) continue; + const int var = PositiveRef(literal); + const int rep = orbits[var]; + if (rep == -1) continue; + + ++num_in_intersections; + if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep); + tmp_sizes[rep]++; + if (RefIsPositive(literal)) tmp_num_positive[rep]++; + } + + for (const int row : tmp_to_clear) { + const int size = tmp_sizes[row]; + const int num_positive = tmp_num_positive[row]; + const int num_negative = tmp_sizes[row] - tmp_num_positive[row]; + tmp_sizes[row] = 0; + tmp_num_positive[row] = 0; + + if (num_positive > 1 && num_negative == 0) { + at_most_one_one[row] = true; + } else if (num_positive == 0 && num_negative > 1) { + at_most_one_zero[row] = true; + } else if (num_positive > 0 && num_negative > 0) { + all_equivalent_rows[row] = true; + } + + // We might be able to presolve more in these cases. + if (at_most_one_zero[row] || at_most_one_one[row] || + all_equivalent_rows[row]) { + if (tmp_to_clear.size() > 1) { + context->UpdateRuleStats( + "TODO symmetry: at most one across orbits."); + } else if (size < orbitope[0].size()) { + context->UpdateRuleStats( + "TODO symmetry: at most one can be extended"); + } + } + } + } + + // Heuristically choose a "best" row/col to "distinguish" and break the + // symmetry on. + int best_row = 0; + int best_col = 0; + int best_score = 0; + bool fix_others_to_zero = true; + for (int i = 0; i < all_equivalent_rows.size(); ++i) { + const int num_cols = orbitope[i].size(); + + // Note that this operation do not change the symmetry group. + if (all_equivalent_rows[i]) { + for (int j = 1; j < num_cols; ++j) { + context->StoreBooleanEqualityRelation(orbitope[i][0], orbitope[i][j]); + context->UpdateRuleStats("symmetry: all equivalent in orbit"); + if (context->ModelIsUnsat()) return false; + } + } + + // Because of symmetry, the choice of the column shouldn't matter (they + // will all appear in the same number of constraints of the same types), + // however we prefer to fix a variable that seems to touch more + // constraints. + // + // TODO(user): maybe we should simplify the constraint using the variable + // we fix before choosing the next row to break symmetry on. + const int row_score = + context->VarToConstraints(PositiveRef(orbitope[i][0])).size(); + + // TODO(user): If one variable make the line already fixed, we should just + // ignore this row. Not too important as actually this shouldn't happen + // because we never compute symmetries involving fixed variables. But in + // the future, fixing some literal might have some side effects and fix + // others. + if (at_most_one_one[i] && row_score > best_score) { + best_col = 0; + for (int j = 0; j < num_cols; ++j) { + if (context->LiteralIsTrue(orbitope[i][j])) { + best_col = j; + break; + } + } + best_row = i; + best_score = row_score; + fix_others_to_zero = true; + } + if (at_most_one_zero[i] && row_score > best_score) { + best_col = 0; + for (int j = 0; j < num_cols; ++j) { + if (context->LiteralIsFalse(orbitope[i][j])) { + best_col = j; + break; + } + } + best_row = i; + best_score = row_score; + fix_others_to_zero = false; + } + } + + if (best_score == 0) break; + for (int j = 1; j < orbitope[best_row].size(); ++j) { + if (fix_others_to_zero) { + context->UpdateRuleStats("symmetry: fixed to false"); + if (!context->SetLiteralToFalse(orbitope[best_row][j])) return false; + } else { + context->UpdateRuleStats("symmetry: fixed to true"); + if (!context->SetLiteralToTrue(orbitope[best_row][j])) return false; + } + } + + // We add the symmetry breaking inequalities: best_var >= all other var + // in orbit. That is not(best_var) => not(other) for Booleans. We only add + // them if we didn't fix any variable just above. + // + // TODO(user): Add the inequality for non-Boolean too? Also note that this + // code only run if the code above is disabled. It is here for testing + // alternatives. In particular, if there is no at most one, we cannot fix + // n-1 variables, but we can still add inequalities. + const int best_var = orbitope[best_row][best_col]; + const bool maximize_best_var = fix_others_to_zero; + if (context->CanBeUsedAsLiteral(best_var) && !context->IsFixed(best_var)) { + ConstraintProto* ct = context->working_model->add_constraints(); + ct->add_enforcement_literal(maximize_best_var ? NegatedRef(best_var) + : best_var); + for (const int other : orbitope[best_row]) { + if (other == best_var) continue; + if (context->IsFixed(other)) continue; + ct->mutable_bool_and()->add_literals( + maximize_best_var ? NegatedRef(other) : other); + context->UpdateRuleStats("symmetry: added implication"); + } + context->UpdateNewConstraintsVariableUsage(); + } + + // Remove the column of best_var. + CHECK_NE(best_col, -1); + for (int i = 0; i < orbitope.size(); ++i) { + std::swap(orbitope[i][best_col], orbitope[i].back()); + orbitope[i].pop_back(); + } + + // We also remove the line of best_var since heuristicially, it is better to + // not add symmetries involving any of the variable on this line. + std::swap(orbitope[best_row], orbitope.back()); + orbitope.pop_back(); + } + + return true; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_symmetries.h b/ortools/sat/cp_model_symmetries.h index 61ef00f2ec..40705f2707 100644 --- a/ortools/sat/cp_model_symmetries.h +++ b/ortools/sat/cp_model_symmetries.h @@ -20,6 +20,7 @@ #include "ortools/algorithms/sparse_permutation.h" #include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/presolve_context.h" #include "ortools/sat/sat_parameters.pb.h" namespace operations_research { @@ -30,6 +31,9 @@ namespace sat { // of variables of the problem. They are permutations of the (index // representation of the) problem variables. // +// Note that we ignore the variables that appear in no constraint, instead of +// outputing the full symmetry group involving them. +// // TODO(user): On SAT problems it is more powerful to detect permutations also // involving the negation of the problem variables. So that we could find a // symmetry x <-> not(y) for instance. @@ -37,11 +41,19 @@ namespace sat { // TODO(user): As long as we only exploit symmetry involving only Boolean // variables we can make this code more efficient by not detecting symmetries // involing integer variable. +// +// TODO(user): Make the limit deterministic. void FindCpModelSymmetries( const SatParameters& params, const CpModelProto& problem, std::vector>* generators, double time_limit_seconds = std::numeric_limits::infinity()); +// Basic implementation of some symmetry breaking during presolve. +// +// Currently this just try to fix variables by detecting symmetries between +// Booleans in bool_and, at_most_one or exactly_one constraints. +bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h index e99ac7427d..95b3abe7e8 100644 --- a/ortools/sat/intervals.h +++ b/ortools/sat/intervals.h @@ -710,7 +710,7 @@ NewOptionalIntervalWithVariableSize(int64 min_start, int64 max_end, // This requires that all the alternatives are optional tasks. inline std::function IntervalWithAlternatives( - IntervalVariable main, const std::vector& members) { + IntervalVariable parent, const std::vector& members) { return [=](Model* model) { auto* integer_trail = model->GetOrCreate(); auto* intervals = model->GetOrCreate(); @@ -725,9 +725,9 @@ inline std::function IntervalWithAlternatives( const Literal is_present = intervals->PresenceLiteral(member); sat_ct.push_back({is_present, Coefficient(1)}); model->Add( - Equality(model->Get(StartVar(main)), model->Get(StartVar(member)))); + Equality(model->Get(StartVar(parent)), model->Get(StartVar(member)))); model->Add( - Equality(model->Get(EndVar(main)), model->Get(EndVar(member)))); + Equality(model->Get(EndVar(parent)), model->Get(EndVar(member)))); // TODO(user): IsOneOf() only work for members with fixed size. // Generalize to an "int_var_element" constraint. @@ -735,12 +735,12 @@ inline std::function IntervalWithAlternatives( presences.push_back(is_present); sizes.push_back(intervals->MinSize(member)); } - if (intervals->SizeVar(main) != kNoIntegerVariable) { - model->Add(IsOneOf(intervals->SizeVar(main), presences, sizes)); + if (intervals->SizeVar(parent) != kNoIntegerVariable) { + model->Add(IsOneOf(intervals->SizeVar(parent), presences, sizes)); } model->Add(BooleanLinearConstraint(1, 1, &sat_ct)); - // Propagate from the candidate bounds to the main interval ones. + // Propagate from the candidate bounds to the parent interval ones. { std::vector starts; starts.reserve(members.size()); @@ -748,7 +748,7 @@ inline std::function IntervalWithAlternatives( starts.push_back(intervals->StartVar(member)); } model->Add( - PartialIsOneOfVar(intervals->StartVar(main), starts, presences)); + PartialIsOneOfVar(intervals->StartVar(parent), starts, presences)); } { std::vector ends; @@ -756,7 +756,7 @@ inline std::function IntervalWithAlternatives( for (const IntervalVariable member : members) { ends.push_back(intervals->EndVar(member)); } - model->Add(PartialIsOneOfVar(intervals->EndVar(main), ends, presences)); + model->Add(PartialIsOneOfVar(intervals->EndVar(parent), ends, presences)); } }; } diff --git a/ortools/sat/symmetry_util.cc b/ortools/sat/symmetry_util.cc new file mode 100644 index 0000000000..9c77f1f079 --- /dev/null +++ b/ortools/sat/symmetry_util.cc @@ -0,0 +1,188 @@ +// Copyright 2010-2018 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// 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. + +#include "ortools/sat/symmetry_util.h" + +#include "ortools/algorithms/dynamic_partition.h" + +namespace operations_research { +namespace sat { + +std::vector> BasicOrbitopeExtraction( + const std::vector>& generators) { + // Count the number of permutations that are compositions of 2-cycle and + // regroup them according to the number of cycles. + std::vector> num_cycles_to_2cyclers; + for (int g = 0; g < generators.size(); ++g) { + const std::unique_ptr& perm = generators[g]; + bool contain_only_2cycles = true; + const int num_cycles = perm->NumCycles(); + for (int i = 0; i < num_cycles; ++i) { + if (perm->Cycle(i).size() != 2) { + contain_only_2cycles = false; + break; + } + } + if (!contain_only_2cycles) continue; + if (num_cycles >= num_cycles_to_2cyclers.size()) { + num_cycles_to_2cyclers.resize(num_cycles + 1); + } + num_cycles_to_2cyclers[num_cycles].push_back(g); + } + + // Heuristic: we try to grow the orbitope that has the most potential for + // fixing variables. + // + // TODO(user): We could grow each and keep the real maximum. + int best = -1; + int best_score = 0; + for (int i = 0; i < num_cycles_to_2cyclers.size(); ++i) { + if (num_cycles_to_2cyclers[i].size() > 1) { + const int num_perms = num_cycles_to_2cyclers[i].size() + 1; + VLOG(1) << "Potential orbitope: " << i << " x " << num_perms; + const int64 score = std::min(i, num_perms); + if (score > best_score) { + best = i; + best_score = score; + } + } + } + + std::vector> orbitope; + if (best == -1) return orbitope; + + // We will track the element already added so we never have duplicates. + std::vector in_matrix; + + // Greedily grow the orbitope. + orbitope.resize(best); + for (const int g : num_cycles_to_2cyclers[best]) { + // Start using the first permutation. + if (orbitope[0].empty()) { + const std::unique_ptr& perm = generators[g]; + const int num_cycles = perm->NumCycles(); + for (int i = 0; i < num_cycles; ++i) { + for (const int x : perm->Cycle(i)) { + orbitope[i].push_back(x); + if (x >= in_matrix.size()) in_matrix.resize(x + 1, false); + in_matrix[x] = true; + } + } + continue; + } + + // We want to find a column such that g sends it to variables not already + // in the orbitope matrix. + // + // Note(user): This relies on the cycle in each permutation to be ordered by + // smaller element first. This way we don't have to account any row + // permutation of the orbitope matrix. The code that detect the symmetries + // of the problem should already return permutation in this canonical + // format. + std::vector grow; + int matching_column_index = -1; + const std::unique_ptr& perm = generators[g]; + const int num_cycles = perm->NumCycles(); + for (int i = 0; i < num_cycles; ++i) { + // Extract the two elements of this transposition. + std::vector tmp; + for (const int x : perm->Cycle(i)) tmp.push_back(x); + const int a = tmp[0]; + const int b = tmp[1]; + + // We want one element to appear in matching_column_index and the other to + // not appear at all. + int num_matches_a = 0; + int num_matches_b = 0; + int last_match_index; + for (int j = 0; j < orbitope[i].size(); ++j) { + if (orbitope[i][j] == a) { + ++num_matches_a; + last_match_index = j; + } else if (orbitope[i][j] == b) { + ++num_matches_b; + last_match_index = j; + } + } + if (matching_column_index == -1) { + matching_column_index = last_match_index; + } + if (matching_column_index != last_match_index) break; + if (num_matches_a == 0 && num_matches_b == 1) { + if (a >= in_matrix.size() || !in_matrix[a]) grow.push_back(a); + } else if (num_matches_a == 1 && num_matches_b == 0) { + if (b >= in_matrix.size() || !in_matrix[b]) grow.push_back(b); + } else { + break; + } + } + + // If grow is of full size, we can extend the orbitope. + if (grow.size() == num_cycles) { + for (int i = 0; i < orbitope.size(); ++i) { + orbitope[i].push_back(grow[i]); + if (grow[i] >= in_matrix.size()) in_matrix.resize(grow[i] + 1, false); + in_matrix[grow[i]] = true; + } + } + } + + return orbitope; +} + +std::vector GetOrbits( + int n, const std::vector>& generators) { + MergingPartition union_find; + union_find.Reset(n); + for (const std::unique_ptr& perm : generators) { + const int num_cycles = perm->NumCycles(); + for (int i = 0; i < num_cycles; ++i) { + // Note that there is currently no random access api like cycle[j]. + int first; + bool is_first = true; + for (const int x : perm->Cycle(i)) { + if (is_first) { + first = x; + is_first = false; + } else { + union_find.MergePartsOf(first, x); + } + } + } + } + + int num_parts = 0; + std::vector orbits(n, -1); + for (int i = 0; i < n; ++i) { + if (union_find.NumNodesInSamePartAs(i) == 1) continue; + const int root = union_find.GetRootAndCompressPath(i); + if (orbits[root] == -1) orbits[root] = num_parts++; + orbits[i] = orbits[root]; + } + return orbits; +} + +std::vector GetOrbitopeOrbits( + int n, const std::vector>& orbitope) { + std::vector orbits(n, -1); + for (int i = 0; i < orbitope.size(); ++i) { + for (int j = 0; j < orbitope[i].size(); ++j) { + CHECK_EQ(orbits[orbitope[i][j]], -1); + orbits[orbitope[i][j]] = i; + } + } + return orbits; +} + +} // namespace sat +} // namespace operations_research diff --git a/ortools/sat/symmetry_util.h b/ortools/sat/symmetry_util.h new file mode 100644 index 0000000000..142002b445 --- /dev/null +++ b/ortools/sat/symmetry_util.h @@ -0,0 +1,82 @@ +// Copyright 2010-2018 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// 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. + +#ifndef OR_TOOLS_SAT_SYMMETRY_UTIL_H_ +#define OR_TOOLS_SAT_SYMMETRY_UTIL_H_ + +#include "ortools/algorithms/sparse_permutation.h" + +namespace operations_research { +namespace sat { + +// Given the generator for a permutation group of [0, n-1], tries to identify +// a grouping of the variables in an p x q matrix such that any permutations +// of the columns of this matrix is in the given group. +// +// The name comes from: "Packing and Partitioning Orbitopes", Volker Kaibel, +// Marc E. Pfetsch, https://arxiv.org/abs/math/0603678 . Here we just detect it, +// independently of the constraints on the variables in this matrix. We can also +// detect non-Boolean orbitope. +// +// In order to detect orbitope, this basic algorithm requires that the +// generators of the orbitope must only contain one or more 2-cyle (i.e +// transpositions). Thus they must be involutions. The list of transpositions in +// the SparsePermutation must also be listed in a canonical order. +// +// TODO(user): Detect more than one orbitope? Note that once detected, the +// structure can be exploited efficiently, but for now, a more "generic" +// algorithm based on stabilizator should achieve the same preprocessing power, +// so I don't know how hard we need to invest in orbitope detection. +// +// TODO(user): The heuristic is quite limited for now, but this works on +// graph20-20-1rand.mps.gz. I suspect the generators provided by the detection +// code follow our preconditions. +std::vector> BasicOrbitopeExtraction( + const std::vector>& generators); + +// Returns a vector of size n such that +// - orbits[i] == -1 iff i is never touched by the generators (singleton orbit). +// - orbits[i] = orbit_index, where orbits are numbered from 0 to num_orbits - 1 +// +// TODO(user): We could reuse the internal memory if needed. +std::vector GetOrbits( + int n, const std::vector>& generators); + +// Returns the orbits under the given orbitope action. +// Same results format as in GetOrbits(). Note that here, the orbit index +// is simply the row index of an element in the orbitope matrix. +std::vector GetOrbitopeOrbits( + int n, const std::vector>& orbitope); + +// Given the generators for a permutation group of [0, n-1], update it to +// a set of generators of the group stabilizing the given element. +// +// Note that one can add symmetry breaking constraints by repeatedly doing: +// 1/ Call GetOrbits() using the current set of generators. +// 2/ Choose an element x0 in a large orbit (x0, .. xi ..) , and add x0 >= xi +// for all i. +// 3/ Update the set of generators to the one stabilizing x0. +// +// This is more or less what is described in "Symmetry Breaking Inequalities +// from the Schreier-Sims Table", Domenico Salvagnin, +// https://link.springer.com/chapter/10.1007/978-3-319-93031-2_37 +// +// TODO(user): Implement! +inline void TransformToGeneratorOfStabilizer( + int to_stabilize, + std::vector>* generators) {} + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_SYMMETRY_UTIL_H_