always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

more work on symmetry detection and usage in CP-SAT

Browse Source
This commit is contained in:
Laurent Perron
2021-01-28 14:41:20 +01:00
parent 7b9c099453
commit 8abff55ab6
9 changed files with 755 additions and 116 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
makefiles/Makefile.gen.mk
View File

@@ -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

14
ortools/sat/BUILD
View File

@@ -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",

206
ortools/sat/cp_model_presolve.cc
View File

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

3
ortools/sat/cp_model_presolve.h
View File

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

336
ortools/sat/cp_model_symmetries.cc
View File

@@ -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<std::vector<int64>, int, VectorHash> id_map_;
};
@@ -116,9 +119,11 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
std::vector<int64> 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<Graph> GenerateGraphForSymmetryDetection(
std::vector<int64> 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<Graph> 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<int> 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<std::unique_ptr<SparsePermutation>> 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<std::vector<int>> 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<const google::protobuf::RepeatedField<int32>*> 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<int> tmp_to_clear;
std::vector<int> tmp_sizes(num_vars, 0);
std::vector<int> 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<int32>* 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<int> 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<bool> all_equivalent_rows(orbitope.size(), false);
std::vector<bool> at_most_one_one(orbitope.size(), false);
std::vector<bool> at_most_one_zero(orbitope.size(), false);
for (const google::protobuf::RepeatedField<int32>* 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

12
ortools/sat/cp_model_symmetries.h
View File

@@ -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<std::unique_ptr<SparsePermutation>>* generators,
double time_limit_seconds = std::numeric_limits<double>::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

16
ortools/sat/intervals.h
View File

@@ -710,7 +710,7 @@ NewOptionalIntervalWithVariableSize(int64 min_start, int64 max_end,
// This requires that all the alternatives are optional tasks.
inline std::function<void(Model*)> IntervalWithAlternatives(
IntervalVariable main, const std::vector<IntervalVariable>& members) {
IntervalVariable parent, const std::vector<IntervalVariable>& members) {
return [=](Model* model) {
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
auto* intervals = model->GetOrCreate<IntervalsRepository>();
@@ -725,9 +725,9 @@ inline std::function<void(Model*)> 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<void(Model*)> 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<IntegerVariable> starts;
starts.reserve(members.size());
@@ -748,7 +748,7 @@ inline std::function<void(Model*)> IntervalWithAlternatives(
starts.push_back(intervals->StartVar(member));
}
model->Add(
PartialIsOneOfVar(intervals->StartVar(main), starts, presences));
PartialIsOneOfVar(intervals->StartVar(parent), starts, presences));
}
{
std::vector<IntegerVariable> ends;
@@ -756,7 +756,7 @@ inline std::function<void(Model*)> 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));
}
};
}

188
ortools/sat/symmetry_util.cc Normal file
View File

@@ -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<std::vector<int>> BasicOrbitopeExtraction(
const std::vector<std::unique_ptr<SparsePermutation>>& generators) {
// Count the number of permutations that are compositions of 2-cycle and
// regroup them according to the number of cycles.
std::vector<std::vector<int>> num_cycles_to_2cyclers;
for (int g = 0; g < generators.size(); ++g) {
const std::unique_ptr<SparsePermutation>& 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<std::vector<int>> orbitope;
if (best == -1) return orbitope;
// We will track the element already added so we never have duplicates.
std::vector<bool> 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<SparsePermutation>& 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<int> grow;
int matching_column_index = -1;
const std::unique_ptr<SparsePermutation>& 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<int> 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<int> GetOrbits(
int n, const std::vector<std::unique_ptr<SparsePermutation>>& generators) {
MergingPartition union_find;
union_find.Reset(n);
for (const std::unique_ptr<SparsePermutation>& 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<int> 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<int> GetOrbitopeOrbits(
int n, const std::vector<std::vector<int>>& orbitope) {
std::vector<int> 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

82
ortools/sat/symmetry_util.h Normal file
View File

@@ -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<std::vector<int>> BasicOrbitopeExtraction(
const std::vector<std::unique_ptr<SparsePermutation>>& 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<int> GetOrbits(
int n, const std::vector<std::unique_ptr<SparsePermutation>>& 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<int> GetOrbitopeOrbits(
int n, const std::vector<std::vector<int>>& 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<std::unique_ptr<SparsePermutation>>* generators) {}
} // namespace sat
} // namespace operations_research
#endif // OR_TOOLS_SAT_SYMMETRY_UTIL_H_