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

[CP-SAT] compute symmetries during presolve; save them in the cp_model.proto

Browse Source
This commit is contained in:
Laurent Perron
2021-02-07 00:46:28 +01:00
parent 3ea59976a7
commit dfe40729e2
8 changed files with 284 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
ortools/sat/cp_model.proto
View File

@@ -230,6 +230,12 @@ message RoutesConstraintProto {
// Experimental. The demands for each node, and the maximum capacity for each
// route. Note that this is currently only used for the LP relaxation and one
// need to add the corresponding constraint to enforce this outside of the LP.
//
// TODO(user): Ideally, we should be able to extract any dimension like these
// (i.e. capacity, route_length, etc..) automatically from the encoding. The
// classical way to encode that is to have "current_capacity" variables along
// the route and linear equations of the form:
// arc_literal => (current_capacity_tail + demand <= current_capacity_head)
repeated int32 demands = 4;
int64 capacity = 5;
}
@@ -515,6 +521,53 @@ message PartialVariableAssignment {
repeated int64 values = 2;
}
// A permutation of integers encoded as a list of cycles, hence the "sparse"
// format. The image of an element cycle[i] is cycle[(i + 1) % cycle_length].
message SparsePermutationProto {
// Each cycle is listed one after the other in the support field.
// The size of each cycle is given (in order) in the cycle_sizes field.
repeated int32 support = 1;
repeated int32 cycle_sizes = 2;
}
// A dense matrix of numbers encoded in a flat way, row by row.
// That is matrix[i][j] = entries[i * num_cols + j];
message DenseMatrixProto {
int32 num_rows = 1;
int32 num_cols = 2;
repeated int32 entries = 3;
}
// Experimental. For now, this is meant to be used by the solver and not filled
// by clients.
//
// Hold symmetry information about the set of feasible solutions. If we permute
// the variable values of any feasible solution using one of the permutation
// described here, we should always get another feasible solution.
//
// We usually also enforce that the objective of the new solution is the same.
//
// The group of permutations encoded here is usually computed from the encoding
// of the model, so it is not meant to be a complete representation of the
// feasible solution symmetries, just a valid subgroup.
message SymmetryProto {
// A list of variable indices permutations that leave the feasible space of
// solution invariant. Usually, we only encode a set of generators of the
// group.
repeated SparsePermutationProto permutations = 1;
// An orbitope is a special symmetry structure of the solution space. If the
// variable indices are arranged in a matrix (with no duplicates), then any
// permutation of the columns will be a valid permutation of the feasible
// space.
//
// This arise quite often. The typical example is a graph coloring problem
// where for each node i, you have j booleans to indicate its color. If the
// variables color_of_i_is_j are arranged in a matrix[i][j], then any columns
// permutations leave the problem invariant.
repeated DenseMatrixProto orbitopes = 2;
}
// A constraint programming problem.
message CpModelProto {
// For debug/logging only. Can be empty.
@@ -566,6 +619,13 @@ message CpModelProto {
// Such infeasibility explanation will be available in the
// sufficient_assumptions_for_infeasibility response field.
repeated int32 assumptions = 7;
// For now, this is not meant to be filled by a client writing a model, but
// by our preprocessing step.
//
// Information about the symmetries of the feasible solution space.
// These usually leaves the objective invariant.
SymmetryProto symmetry = 8;
}
// The status returned by a solver trying to solve a CpModelProto.

30
ortools/sat/cp_model_loader.cc
View File

@@ -32,7 +32,6 @@
#include "ortools/sat/circuit.h"
#include "ortools/sat/cp_constraints.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_symmetries.h"
#include "ortools/sat/cp_model_utils.h"
#include "ortools/sat/cumulative.h"
#include "ortools/sat/diffn.h"
@@ -271,25 +270,20 @@ void CpModelMapping::CreateVariables(const CpModelProto& model_proto,
}
}
void CpModelMapping::DetectAndLoadBooleanSymmetries(
const CpModelProto& model_proto, Model* m) {
void CpModelMapping::LoadBooleanSymmetries(const CpModelProto& model_proto,
Model* m) {
const SatParameters& params = *m->GetOrCreate<SatParameters>();
if (!params.detect_symmetries()) return;
// TODO(user): Use a deterministic time limit.
std::vector<std::unique_ptr<SparsePermutation>> generators;
FindCpModelSymmetries(params, model_proto, &generators,
/*time_limit_seconds=*/1.0);
if (generators.empty()) return;
const SymmetryProto symmetry = model_proto.symmetry();
if (symmetry.permutations().empty()) return;
auto* sat_solver = m->GetOrCreate<SatSolver>();
auto* symmetry_handler = m->GetOrCreate<SymmetryPropagator>();
sat_solver->AddPropagator(symmetry_handler);
const int num_literals = 2 * sat_solver->NumVariables();
for (const auto& perm : generators) {
for (const SparsePermutationProto& perm : symmetry.permutations()) {
bool all_bool = true;
for (const int var : perm->Support()) {
for (const int var : perm.support()) {
if (!IsBoolean(var)) {
all_bool = false;
break;
@@ -300,16 +294,22 @@ void CpModelMapping::DetectAndLoadBooleanSymmetries(
// Convert the variable symmetry to a "literal" one.
auto literal_permutation =
absl::make_unique<SparsePermutation>(num_literals);
const int num_cycle = perm->NumCycles();
int support_index = 0;
const int num_cycle = perm.cycle_sizes().size();
for (int i = 0; i < num_cycle; ++i) {
for (const int var : perm->Cycle(i)) {
const int size = perm.cycle_sizes(i);
const int saved_support_index = support_index;
for (int j = 0; j < size; ++j) {
const int var = perm.support(support_index++);
literal_permutation->AddToCurrentCycle(Literal(var).Index().value());
}
literal_permutation->CloseCurrentCycle();
// Note that we also need to add the corresponding cycle for the negated
// literals.
for (const int var : perm->Cycle(i)) {
support_index = saved_support_index;
for (int j = 0; j < size; ++j) {
const int var = perm.support(support_index++);
literal_permutation->AddToCurrentCycle(
Literal(var).NegatedIndex().value());
}

12
ortools/sat/cp_model_loader.h
View File

@@ -80,12 +80,12 @@ class CpModelMapping {
// Automatically detect optional variables.
void DetectOptionalVariables(const CpModelProto& model_proto, Model* m);
// Experimental. Detects symmetries of the model, and if they only involve
// Booleans, add the corresponding symmetries to the solver. We only have the
// code for Booleans, it is why we currently ignore symmetries involving
// integer variables.
void DetectAndLoadBooleanSymmetries(const CpModelProto& model_proto,
Model* m);
// Experimental. Loads the symmetry form the proto symmetry field, as long as
// they only involve Booleans.
//
// TODO(user): We currently only have the code for Booleans, it is why we
// currently ignore symmetries involving integer variables.
void LoadBooleanSymmetries(const CpModelProto& model_proto, Model* m);
// Extract the encodings (IntegerVariable <-> Booleans) present in the model.
// This effectively load some linear constraints of size 1 that will be marked

2
ortools/sat/cp_model_presolve.cc
View File

@@ -5255,7 +5255,7 @@ bool CpModelPresolver::Presolve() {
// 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() &&
if (context_->params().symmetry_level() > 0 && !context_->ModelIsUnsat() &&
!context_->time_limit()->LimitReached() &&
!context_->keep_all_feasible_solutions) {
DetectAndExploitSymmetriesInPresolve(context_);

10
ortools/sat/cp_model_solver.cc
View File

@@ -61,6 +61,7 @@
#include "ortools/sat/cp_model_postsolve.h"
#include "ortools/sat/cp_model_presolve.h"
#include "ortools/sat/cp_model_search.h"
#include "ortools/sat/cp_model_symmetries.h"
#include "ortools/sat/cp_model_utils.h"
#include "ortools/sat/cuts.h"
#include "ortools/sat/drat_checker.h"
@@ -1201,8 +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()) {
mapping->DetectAndLoadBooleanSymmetries(model_proto, model);
if (!parameters.optimize_with_core() && parameters.symmetry_level() > 1) {
mapping->LoadBooleanSymmetries(model_proto, model);
}
mapping->ExtractEncoding(model_proto, model);
@@ -2469,6 +2470,7 @@ class LnsSolver : public SubSolver {
local_params.set_stop_after_first_solution(false);
local_params.set_log_search_progress(false);
local_params.set_cp_model_probing_level(0);
local_params.set_symmetry_level(0);
if (absl::GetFlag(FLAGS_cp_model_dump_lns)) {
const std::string name =
@@ -3108,6 +3110,10 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) {
// Delete the context.
context.reset(nullptr);
if (params.symmetry_level() > 1) {
DetectAndAddSymmetryToProto(params, &new_cp_model_proto);
}
SharedResponseManager shared_response_manager(
log_search, params.enumerate_all_solutions(), &new_cp_model_proto,
&wall_timer, &shared_time_limit);

219
ortools/sat/cp_model_symmetries.cc
View File

@@ -105,12 +105,16 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
};
IdGenerator color_id_generator;
initial_equivalence_classes->clear();
auto new_node = [&initial_equivalence_classes,
auto new_node = [&initial_equivalence_classes, &graph,
&color_id_generator](const std::vector<int64>& color) {
// Since we add nodes one by one, initial_equivalence_classes->size() gives
// the number of nodes at any point, which we use as the next node index.
const int node = initial_equivalence_classes->size();
initial_equivalence_classes->push_back(color_id_generator.GetId(color));
// In some corner cases, we create a node but never uses it. We still
// want it to be there.
graph->AddNode(node);
return node;
};
@@ -122,8 +126,10 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
// with the smallest cost?
std::vector<int64> objective_by_var(num_variables, 0);
for (int i = 0; i < problem.objective().vars_size(); ++i) {
objective_by_var[problem.objective().vars(i)] =
problem.objective().coeffs(i);
const int ref = problem.objective().vars(i);
const int var = PositiveRef(ref);
const int64 coeff = problem.objective().coeffs(i);
objective_by_var[var] = RefIsPositive(ref) ? coeff : -coeff;
}
// Create one node for each variable. Note that the code rely on the fact that
@@ -133,10 +139,6 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
tmp_color = {VARIABLE_NODE, objective_by_var[v]};
Append(problem.variables(v).domain(), &tmp_color);
CHECK_EQ(v, new_node(tmp_color));
// Make sure the graph contains all the variable nodes, even if no edges are
// attached to them through constraints.
graph->AddNode(v);
}
// We will lazily create "coefficient nodes" that correspond to a variable
@@ -170,21 +172,42 @@ std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
return get_coefficient_node(PositiveRef(ref), RefIsPositive(ref) ? 1 : -1);
};
// Because the implication can be numerous, we encode them without constraints
// node by using an arc from the lhs to the rhs. Note that we also always add
// the other direction. We use a set to remove duplicates both for efficiency
// and to not artificially break symmetries by using multi-arcs.
// Because the implications can be numerous, we encode them without
// constraints node by using an arc from the lhs to the rhs. Note that we also
// always add the other direction. We use a set to remove duplicates both for
// efficiency and to not artificially break symmetries by using multi-arcs.
//
// Tricky: We cannot use the base variable node here to avoid situation like
// both a variable a and b having the same children (not(a), not(b)) in the
// graph. Because if that happen, we can permute a and b without permuting
// their associated not(a) and not(b) node! To be sure this cannot happen, a
// variable node can not have as children a VAR_COEFFICIENT_NODE from another
// node. This makes sure that any permutation that touch a variable, must
// permute its coefficient nodes accordingly.
absl::flat_hash_set<std::pair<int, int>> implications;
auto add_implication = [&get_literal_node, &graph, &implications](int ref_a,
int ref_b) {
auto get_implication_node = [&new_node, &graph, &coefficient_nodes,
&tmp_color](int ref) {
const int var = PositiveRef(ref);
const int64 coeff = RefIsPositive(ref) ? 1 : -1;
const auto insert =
coefficient_nodes.insert({std::make_pair(var, coeff), 0});
if (!insert.second) return insert.first->second;
tmp_color = {VAR_COEFFICIENT_NODE, coeff};
const int secondary_node = new_node(tmp_color);
graph->AddArc(var, secondary_node);
insert.first->second = secondary_node;
return secondary_node;
};
auto add_implication = [&get_implication_node, &graph, &implications](
int ref_a, int ref_b) {
const auto insert = implications.insert({ref_a, ref_b});
if (!insert.second) return;
graph->AddArc(get_literal_node(ref_a), get_literal_node(ref_b));
graph->AddArc(get_implication_node(ref_a), get_implication_node(ref_b));
// Always add the other side.
implications.insert({NegatedRef(ref_b), NegatedRef(ref_a)});
graph->AddArc(get_literal_node(NegatedRef(ref_b)),
get_literal_node(NegatedRef(ref_a)));
graph->AddArc(get_implication_node(NegatedRef(ref_b)),
get_implication_node(NegatedRef(ref_a)));
};
// Add constraints to the graph.
@@ -407,13 +430,54 @@ void FindCpModelSymmetries(
}
}
void DetectAndAddSymmetryToProto(const SatParameters& params,
CpModelProto* proto) {
const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
SymmetryProto* symmetry = proto->mutable_symmetry();
symmetry->Clear();
std::vector<std::unique_ptr<SparsePermutation>> generators;
FindCpModelSymmetries(params, *proto, &generators,
/*deterministic_limit=*/1.0);
if (generators.empty()) return;
for (const std::unique_ptr<SparsePermutation>& perm : generators) {
SparsePermutationProto* perm_proto = symmetry->add_permutations();
const int num_cycle = perm->NumCycles();
for (int i = 0; i < num_cycle; ++i) {
const int old_size = perm_proto->support().size();
for (const int var : perm->Cycle(i)) {
perm_proto->add_support(var);
}
perm_proto->add_cycle_sizes(perm_proto->support().size() - old_size);
}
}
std::vector<std::vector<int>> orbitope = BasicOrbitopeExtraction(generators);
if (orbitope.empty()) return;
if (log_info) {
LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x "
<< orbitope[0].size();
}
DenseMatrixProto* matrix = symmetry->add_orbitopes();
matrix->set_num_rows(orbitope.size());
matrix->set_num_cols(orbitope[0].size());
for (const std::vector<int>& row : orbitope) {
for (const int entry : row) {
matrix->add_entries(entry);
}
}
}
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();
if (context->working_model->has_objective()) {
context->WriteObjectiveToProto();
}
const int num_vars = proto.variables_size();
for (int i = 0; i < num_vars; ++i) {
FillDomainInProto(context->DomainOf(i),
@@ -443,7 +507,8 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
}
std::vector<std::unique_ptr<SparsePermutation>> generators;
FindCpModelSymmetries(params, proto, &generators, /*time_limit_seconds=*/1.0);
FindCpModelSymmetries(params, proto, &generators,
/*deterministic_limit=*/1.0);
// Remove temporary affine relation.
context->working_model->mutable_constraints()->DeleteSubrange(
@@ -451,6 +516,86 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
if (generators.empty()) return true;
// Collect the at most ones.
//
// Note(user): This relies on the fact that the pointers 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());
}
}
// Experimental. Generic approach. First step.
//
// If an at most one intersect with one or more orbit, in each intersection,
// we can fix all but one variable to zero. For now we only test positive
// literal, and maximize the number of fixing.
std::vector<int> can_be_fixed_to_false;
{
const std::vector<int> orbits = GetOrbits(num_vars, generators);
std::vector<int> orbit_sizes;
for (const int rep : orbits) {
if (rep == -1) continue;
if (rep >= orbit_sizes.size()) orbit_sizes.resize(rep + 1, 0);
orbit_sizes[rep]++;
}
std::vector<int> tmp_to_clear;
std::vector<int> tmp_sizes(num_vars, 0);
for (const google::protobuf::RepeatedField<int32>* literals :
at_most_ones) {
tmp_to_clear.clear();
// Compute how many variables we can fix with this at most one.
int num_fixable = 0;
for (const int literal : *literals) {
if (!RefIsPositive(literal)) continue;
if (context->IsFixed(literal)) continue;
const int var = PositiveRef(literal);
const int rep = orbits[var];
if (rep == -1) continue;
// We count all but the first one in each orbit.
if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
if (tmp_sizes[rep] > 0) ++num_fixable;
tmp_sizes[rep]++;
}
// Redo a pass to copy the intersection.
if (num_fixable > can_be_fixed_to_false.size()) {
can_be_fixed_to_false.clear();
for (const int literal : *literals) {
if (!RefIsPositive(literal)) continue;
if (context->IsFixed(literal)) continue;
const int var = PositiveRef(literal);
const int rep = orbits[var];
if (rep == -1) continue;
// We push all but the first one in each orbit.
if (tmp_sizes[rep] == 0) can_be_fixed_to_false.push_back(var);
tmp_sizes[rep] = 0;
}
} else {
// Sparse clean up.
for (const int rep : tmp_to_clear) tmp_sizes[rep] = 0;
}
}
if (log_info) {
LOG(INFO) << "Num fixable by intersecting at_most_one with orbits: "
<< can_be_fixed_to_false.size();
}
}
// Orbitope approach.
//
// This is basically the same as the generic approach, but because of the
@@ -467,28 +612,36 @@ bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
//
// 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) {
if (!orbitope.empty() && log_info) {
LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x "
<< orbitope[0].size();
}
// Collect the at most ones.
// Supper simple heuristic to use the orbitope or not.
//
// Note(user): This relies on the fact that the pointers 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());
// In an orbitope with an at most one on each row, we can fix the upper right
// triangle. We could use a formula, but the loop is fast enough.
//
// TODO(user): Compute the stabilizer under the only non-fixed element and
// iterate!
int max_num_fixed_in_orbitope = 0;
if (!orbitope.empty()) {
const int num_rows = orbitope[0].size();
int size_left = num_rows;
for (int col = 0; size_left > 1 && col < orbitope.size(); ++col) {
max_num_fixed_in_orbitope += size_left - 1;
--size_left;
}
}
if (max_num_fixed_in_orbitope < can_be_fixed_to_false.size()) {
for (int i = 0; i < can_be_fixed_to_false.size(); ++i) {
const int var = can_be_fixed_to_false[i];
context->UpdateRuleStats("symmetry: fixed to false in general orbit");
if (!context->SetLiteralToFalse(var)) return false;
}
return true;
}
if (orbitope.empty()) return true;
// This will always be kept all zero after usage.
std::vector<int> tmp_to_clear;

4
ortools/sat/cp_model_symmetries.h
View File

@@ -46,6 +46,10 @@ void FindCpModelSymmetries(
std::vector<std::unique_ptr<SparsePermutation>>* generators,
double deterministic_limit = std::numeric_limits<double>::infinity());
// Detects symmetries and fill the symmetry field.
void DetectAndAddSymmetryToProto(const SatParameters& params,
CpModelProto* proto);
// Basic implementation of some symmetry breaking during presolve.
//
// Currently this just try to fix variables by detecting symmetries between

6
ortools/sat/sat_parameters.proto
View File

@@ -966,8 +966,10 @@ message SatParameters {
optional bool convert_intervals = 177 [default = false];
// Whether we try to automatically detect the symmetries in a model and
// exploit them.
optional bool detect_symmetries = 183 [default = false];
// 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];
// ==========================================================================
// MIP -> CP-SAT (i.e. IP with integer coeff) conversion parameters that are