18#include "absl/container/flat_hash_map.h"
19#include "absl/memory/memory.h"
20#include "google/protobuf/repeated_field.h"
32 std::size_t operator()(
const std::vector<int64>& values)
const {
49 int GetId(
const std::vector<int64>& color) {
53 int NextFreeId()
const {
return id_map_.size(); }
56 absl::flat_hash_map<std::vector<int64>, int, VectorHash> id_map_;
62template <
typename FieldInt64Type>
64 const google::protobuf::RepeatedField<FieldInt64Type>& repeated_field,
65 std::vector<int64>* vector) {
66 CHECK(vector !=
nullptr);
67 for (
const FieldInt64Type
value : repeated_field) {
68 vector->push_back(
value);
86template <
typename Graph>
88 bool log_info,
const CpModelProto& problem,
89 std::vector<int>* initial_equivalence_classes) {
90 CHECK(initial_equivalence_classes !=
nullptr);
92 const int num_variables = problem.variables_size();
93 auto graph = absl::make_unique<Graph>();
103 VAR_COEFFICIENT_NODE,
106 IdGenerator color_id_generator;
107 initial_equivalence_classes->clear();
108 auto new_node = [&initial_equivalence_classes, &graph,
109 &color_id_generator](
const std::vector<int64>& color) {
112 const int node = initial_equivalence_classes->size();
113 initial_equivalence_classes->push_back(color_id_generator.GetId(color));
117 graph->AddNode(node);
127 std::vector<int64> objective_by_var(num_variables, 0);
128 for (
int i = 0; i < problem.objective().vars_size(); ++i) {
129 const int ref = problem.objective().vars(i);
131 const int64 coeff = problem.objective().coeffs(i);
137 std::vector<int64> tmp_color;
138 for (
int v = 0; v < num_variables; ++v) {
139 tmp_color = {VARIABLE_NODE, objective_by_var[v]};
140 Append(problem.variables(v).domain(), &tmp_color);
146 absl::flat_hash_map<std::pair<int64, int64>,
int> coefficient_nodes;
147 auto get_coefficient_node = [&new_node, &graph, &coefficient_nodes,
149 const int var_node =
var;
155 if (coeff == 1)
return var_node;
158 coefficient_nodes.insert({std::make_pair(
var, coeff), 0});
159 if (!insert.second)
return insert.first->second;
161 tmp_color = {VAR_COEFFICIENT_NODE, coeff};
162 const int secondary_node = new_node(tmp_color);
163 graph->AddArc(var_node, secondary_node);
164 insert.first->second = secondary_node;
165 return secondary_node;
171 auto get_literal_node = [&get_coefficient_node](
int ref) {
187 absl::flat_hash_set<std::pair<int, int>> implications;
188 auto get_implication_node = [&new_node, &graph, &coefficient_nodes,
189 &tmp_color](
int ref) {
193 coefficient_nodes.insert({std::make_pair(
var, coeff), 0});
194 if (!insert.second)
return insert.first->second;
195 tmp_color = {VAR_COEFFICIENT_NODE, coeff};
196 const int secondary_node = new_node(tmp_color);
197 graph->AddArc(
var, secondary_node);
198 insert.first->second = secondary_node;
199 return secondary_node;
201 auto add_implication = [&get_implication_node, &graph, &implications](
202 int ref_a,
int ref_b) {
203 const auto insert = implications.insert({ref_a, ref_b});
204 if (!insert.second)
return;
205 graph->AddArc(get_implication_node(ref_a), get_implication_node(ref_b));
209 graph->AddArc(get_implication_node(
NegatedRef(ref_b)),
214 for (
const ConstraintProto& constraint : problem.constraints()) {
215 const int constraint_node = initial_equivalence_classes->size();
216 std::vector<int64> color = {CONSTRAINT_NODE, constraint.constraint_case()};
218 switch (constraint.constraint_case()) {
219 case ConstraintProto::CONSTRAINT_NOT_SET:
221 case ConstraintProto::kLinear: {
226 Append(constraint.linear().domain(), &color);
227 CHECK_EQ(constraint_node, new_node(color));
228 for (
int i = 0; i < constraint.linear().vars_size(); ++i) {
229 const int ref = constraint.linear().vars(i);
232 ? constraint.linear().coeffs(i)
233 : -constraint.linear().coeffs(i);
234 graph->AddArc(get_coefficient_node(variable_node, coeff),
239 case ConstraintProto::kBoolOr: {
240 CHECK_EQ(constraint_node, new_node(color));
241 for (
const int ref : constraint.bool_or().literals()) {
242 graph->AddArc(get_literal_node(ref), constraint_node);
246 case ConstraintProto::kAtMostOne: {
247 if (constraint.at_most_one().literals().size() == 2) {
249 add_implication(constraint.at_most_one().literals(0),
250 NegatedRef(constraint.at_most_one().literals(1)));
254 CHECK_EQ(constraint_node, new_node(color));
255 for (
const int ref : constraint.at_most_one().literals()) {
256 graph->AddArc(get_literal_node(ref), constraint_node);
260 case ConstraintProto::kExactlyOne: {
261 CHECK_EQ(constraint_node, new_node(color));
262 for (
const int ref : constraint.exactly_one().literals()) {
263 graph->AddArc(get_literal_node(ref), constraint_node);
267 case ConstraintProto::kBoolXor: {
268 CHECK_EQ(constraint_node, new_node(color));
269 for (
const int ref : constraint.bool_xor().literals()) {
270 graph->AddArc(get_literal_node(ref), constraint_node);
274 case ConstraintProto::kBoolAnd: {
277 if (constraint.enforcement_literal_size() != 1) {
279 LOG(
INFO) <<
"BoolAnd with multiple enforcement literal are not "
280 "supported in symmetry code:"
281 << constraint.ShortDebugString();
286 CHECK_EQ(constraint.enforcement_literal_size(), 1);
287 const int ref_a = constraint.enforcement_literal(0);
288 for (
const int ref_b : constraint.bool_and().literals()) {
289 add_implication(ref_a, ref_b);
301 LOG(
INFO) <<
"Unsupported constraint type "
312 if (constraint.constraint_case() != ConstraintProto::kBoolAnd) {
313 for (
const int ref : constraint.enforcement_literal()) {
314 graph->AddArc(constraint_node, get_literal_node(ref));
320 DCHECK_EQ(graph->num_nodes(), initial_equivalence_classes->size());
327 const int num_nodes = graph->num_nodes();
328 std::vector<int> in_degree(num_nodes, 0);
329 std::vector<int> out_degree(num_nodes, 0);
330 for (
int i = 0; i < num_nodes; ++i) {
331 out_degree[i] = graph->OutDegree(i);
332 for (
const int head : (*graph)[i]) {
336 for (
int i = 0; i < num_nodes; ++i) {
337 if (in_degree[i] >= num_nodes || out_degree[i] >= num_nodes) {
338 if (log_info)
LOG(
INFO) <<
"Too many multi-arcs in symmetry code.";
351 int next_id = color_id_generator.NextFreeId();
352 for (
int i = 0; i < num_variables; ++i) {
353 if ((*graph)[i].empty()) {
354 (*initial_equivalence_classes)[i] = next_id++;
360 std::vector<int> mapping(next_id, -1);
361 for (
int& ref : *initial_equivalence_classes) {
362 if (mapping[ref] == -1) {
363 ref = mapping[ref] =
id++;
374 const SatParameters& params,
const CpModelProto& problem,
375 std::vector<std::unique_ptr<SparsePermutation>>* generators,
376 double deterministic_limit) {
377 const bool log_info = params.log_search_progress() ||
VLOG_IS_ON(1);
378 CHECK(generators !=
nullptr);
383 std::vector<int> equivalence_classes;
384 std::unique_ptr<Graph> graph(GenerateGraphForSymmetryDetection<Graph>(
385 log_info, problem, &equivalence_classes));
386 if (graph ==
nullptr)
return;
389 LOG(
INFO) <<
"Graph for symmetry has " << graph->num_nodes()
390 <<
" nodes and " << graph->num_arcs() <<
" arcs.";
392 if (graph->num_nodes() == 0)
return;
395 std::vector<int> factorized_automorphism_group_size;
399 &equivalence_classes, generators, &factorized_automorphism_group_size,
404 if (log_info && !status.ok()) {
405 LOG(
INFO) <<
"GraphSymmetryFinder error: " << status.message();
411 double average_support_size = 0.0;
412 int num_generators = 0;
413 int num_duplicate_constraints = 0;
414 for (
int i = 0; i < generators->size(); ++i) {
416 std::vector<int> to_delete;
417 for (
int j = 0; j < permutation->
NumCycles(); ++j) {
421 if (*(permutation->
Cycle(j).
begin()) >= problem.variables_size()) {
422 to_delete.push_back(j);
425 for (
const int node : permutation->
Cycle(j)) {
426 DCHECK_GE(node, problem.variables_size());
433 if (!permutation->
Support().empty()) {
434 average_support_size += permutation->
Support().size();
435 swap((*generators)[num_generators], (*generators)[i]);
438 ++num_duplicate_constraints;
441 generators->resize(num_generators);
442 average_support_size /= num_generators;
444 LOG(
INFO) <<
"Symmetry computation done. time: "
446 <<
" dtime: " <<
time_limit->GetElapsedDeterministicTime();
447 if (num_generators > 0) {
448 LOG(
INFO) <<
"# of generators: " << num_generators;
449 LOG(
INFO) <<
"Average support size: " << average_support_size;
450 if (num_duplicate_constraints > 0) {
451 LOG(
INFO) <<
"The model contains " << num_duplicate_constraints
452 <<
" duplicate constraints !";
459 CpModelProto*
proto) {
460 const bool log_info = params.log_search_progress() ||
VLOG_IS_ON(1);
461 SymmetryProto* symmetry =
proto->mutable_symmetry();
464 std::vector<std::unique_ptr<SparsePermutation>> generators;
467 if (generators.empty())
return;
469 for (
const std::unique_ptr<SparsePermutation>& perm : generators) {
470 SparsePermutationProto* perm_proto = symmetry->add_permutations();
471 const int num_cycle = perm->NumCycles();
472 for (
int i = 0; i < num_cycle; ++i) {
473 const int old_size = perm_proto->support().size();
474 for (
const int var : perm->Cycle(i)) {
475 perm_proto->add_support(
var);
477 perm_proto->add_cycle_sizes(perm_proto->support().size() - old_size);
482 if (orbitope.empty())
return;
484 LOG(
INFO) <<
"Found orbitope of size " << orbitope.size() <<
" x "
485 << orbitope[0].size();
487 DenseMatrixProto* matrix = symmetry->add_orbitopes();
488 matrix->set_num_rows(orbitope.size());
489 matrix->set_num_cols(orbitope[0].size());
490 for (
const std::vector<int>&
row : orbitope) {
491 for (
const int entry :
row) {
492 matrix->add_entries(entry);
498 const SatParameters& params =
context->params();
499 const bool log_info = params.log_search_progress() ||
VLOG_IS_ON(1);
503 if (
context->working_model->has_objective()) {
504 context->WriteObjectiveToProto();
506 const int num_vars =
proto.variables_size();
507 for (
int i = 0; i < num_vars; ++i) {
509 context->working_model->mutable_variables(i));
515 const int initial_ct_index =
proto.constraints().size();
516 for (
int var = 0;
var < num_vars; ++
var) {
518 if (
context->VariableWasRemoved(
var))
continue;
519 if (
context->VariableIsNotUsedAnymore(
var))
continue;
525 ConstraintProto*
ct =
context->working_model->add_constraints();
526 auto* arg =
ct->mutable_linear();
530 arg->add_coeffs(-r.
coeff);
531 arg->add_domain(r.
offset);
532 arg->add_domain(r.
offset);
535 std::vector<std::unique_ptr<SparsePermutation>> generators;
540 context->working_model->mutable_constraints()->DeleteSubrange(
541 initial_ct_index, num_added);
543 if (generators.empty())
return true;
550 std::vector<const google::protobuf::RepeatedField<int32>*> at_most_ones;
551 for (
int i = 0; i <
proto.constraints_size(); ++i) {
552 if (
proto.constraints(i).constraint_case() == ConstraintProto::kAtMostOne) {
553 at_most_ones.push_back(&
proto.constraints(i).at_most_one().literals());
555 if (
proto.constraints(i).constraint_case() ==
556 ConstraintProto::kExactlyOne) {
557 at_most_ones.push_back(&
proto.constraints(i).exactly_one().literals());
566 std::vector<int> can_be_fixed_to_false;
568 const std::vector<int> orbits =
GetOrbits(num_vars, generators);
569 std::vector<int> orbit_sizes;
570 int max_orbit_size = 0;
571 for (
const int rep : orbits) {
572 if (rep == -1)
continue;
573 if (rep >= orbit_sizes.size()) orbit_sizes.resize(rep + 1, 0);
575 max_orbit_size =
std::max(max_orbit_size, orbit_sizes[rep]);
578 std::vector<int> tmp_to_clear;
579 std::vector<int> tmp_sizes(num_vars, 0);
580 for (
const google::protobuf::RepeatedField<int32>* literals :
582 tmp_to_clear.clear();
586 for (
const int literal : *literals) {
591 const int rep = orbits[
var];
592 if (rep == -1)
continue;
595 if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
596 if (tmp_sizes[rep] > 0) ++num_fixable;
601 if (num_fixable > can_be_fixed_to_false.size()) {
602 can_be_fixed_to_false.clear();
603 for (
const int literal : *literals) {
608 const int rep = orbits[
var];
609 if (rep == -1)
continue;
612 if (tmp_sizes[rep] == 0) can_be_fixed_to_false.push_back(
var);
617 for (
const int rep : tmp_to_clear) tmp_sizes[rep] = 0;
622 LOG(
INFO) <<
"Num fixable by intersecting at_most_one with orbits: "
623 << can_be_fixed_to_false.size()
624 <<
" largest_orbit: " << max_orbit_size;
644 if (!orbitope.empty() && log_info) {
645 LOG(
INFO) <<
"Found orbitope of size " << orbitope.size() <<
" x "
646 << orbitope[0].size();
656 int max_num_fixed_in_orbitope = 0;
657 if (!orbitope.empty()) {
658 const int num_rows = orbitope[0].size();
659 int size_left = num_rows;
660 for (
int col = 0; size_left > 1 &&
col < orbitope.size(); ++
col) {
661 max_num_fixed_in_orbitope += size_left - 1;
665 if (max_num_fixed_in_orbitope < can_be_fixed_to_false.size()) {
666 for (
int i = 0; i < can_be_fixed_to_false.size(); ++i) {
667 const int var = can_be_fixed_to_false[i];
668 context->UpdateRuleStats(
"symmetry: fixed to false in general orbit");
669 if (!
context->SetLiteralToFalse(
var))
return false;
673 if (orbitope.empty())
return true;
676 std::vector<int> tmp_to_clear;
677 std::vector<int> tmp_sizes(num_vars, 0);
678 std::vector<int> tmp_num_positive(num_vars, 0);
683 for (
const google::protobuf::RepeatedField<int32>* literals : at_most_ones) {
684 for (
const int lit : *literals) {
689 for (
const int lit : *literals) {
694 while (!orbitope.empty() && orbitope[0].size() > 1) {
695 const int num_cols = orbitope[0].size();
726 std::vector<bool> all_equivalent_rows(orbitope.size(),
false);
732 bool at_most_one_in_best_rows;
733 int64 best_score = 0;
734 std::vector<int> best_rows;
736 std::vector<int> rows_in_at_most_one;
737 for (
const google::protobuf::RepeatedField<int32>* literals :
739 tmp_to_clear.clear();
740 for (
const int literal : *literals) {
743 const int rep = orbits[
var];
744 if (rep == -1)
continue;
746 if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
751 int num_positive_direction = 0;
752 int num_negative_direction = 0;
757 bool possible_extension =
false;
759 rows_in_at_most_one.clear();
760 for (
const int row : tmp_to_clear) {
761 const int size = tmp_sizes[
row];
762 const int num_positive = tmp_num_positive[
row];
763 const int num_negative = tmp_sizes[
row] - tmp_num_positive[
row];
765 tmp_num_positive[
row] = 0;
767 if (num_positive > 1 && num_negative == 0) {
768 if (size < num_cols) possible_extension =
true;
769 rows_in_at_most_one.push_back(
row);
770 ++num_positive_direction;
771 }
else if (num_positive == 0 && num_negative > 1) {
772 if (size < num_cols) possible_extension =
true;
773 rows_in_at_most_one.push_back(
row);
774 ++num_negative_direction;
775 }
else if (num_positive > 0 && num_negative > 0) {
776 all_equivalent_rows[
row] =
true;
780 if (possible_extension) {
782 "TODO symmetry: possible at most one extension.");
785 if (num_positive_direction > 0 && num_negative_direction > 0) {
786 return context->NotifyThatModelIsUnsat(
"Symmetry and at most ones");
788 const bool direction = num_positive_direction > 0;
800 for (
const int row : rows_in_at_most_one) {
804 if (score > best_score) {
805 at_most_one_in_best_rows = direction;
807 best_rows = rows_in_at_most_one;
816 for (
int i = 0; i < all_equivalent_rows.size(); ++i) {
817 if (all_equivalent_rows[i]) {
818 for (
int j = 1; j < num_cols; ++j) {
819 context->StoreBooleanEqualityRelation(orbitope[i][0], orbitope[i][j]);
820 context->UpdateRuleStats(
"symmetry: all equivalent in orbit");
821 if (
context->ModelIsUnsat())
return false;
834 if (best_score == 0) {
836 "TODO symmetry: add symmetry breaking inequalities?");
843 for (
const int i : best_rows) {
844 for (
int j = 0; j < num_cols; ++j) {
845 const int var = orbitope[i][j];
846 if ((at_most_one_in_best_rows &&
context->LiteralIsTrue(
var)) ||
847 (!at_most_one_in_best_rows &&
context->LiteralIsFalse(
var))) {
848 return context->NotifyThatModelIsUnsat(
"Symmetry and at most one");
857 const int best_col = 0;
858 for (
const int i : best_rows) {
859 for (
int j = 0; j < num_cols; ++j) {
860 if (j == best_col)
continue;
861 const int var = orbitope[i][j];
862 if (at_most_one_in_best_rows) {
863 context->UpdateRuleStats(
"symmetry: fixed to false");
864 if (!
context->SetLiteralToFalse(
var))
return false;
866 context->UpdateRuleStats(
"symmetry: fixed to true");
867 if (!
context->SetLiteralToTrue(
var))
return false;
873 for (
const int i : best_rows) orbitope[i].clear();
875 for (
int i = 0; i < orbitope.size(); ++i) {
876 if (!orbitope[i].empty()) orbitope[new_size++] = orbitope[i];
878 CHECK_LT(new_size, orbitope.size());
879 orbitope.resize(new_size);
882 for (
int i = 0; i < orbitope.size(); ++i) {
883 std::swap(orbitope[i][best_col], orbitope[i].back());
884 orbitope[i].pop_back();
890 if (orbitope.size() == 1) {
891 const int num_cols = orbitope[0].size();
892 for (
int i = 0; i + 1 < num_cols; ++i) {
894 ConstraintProto*
ct =
context->working_model->add_constraints();
895 ct->mutable_linear()->add_coeffs(1);
896 ct->mutable_linear()->add_vars(orbitope[0][i]);
897 ct->mutable_linear()->add_coeffs(-1);
898 ct->mutable_linear()->add_vars(orbitope[0][i + 1]);
899 ct->mutable_linear()->add_domain(0);
901 context->UpdateRuleStats(
"symmetry: added symmetry breaking inequality");
903 context->UpdateNewConstraintsVariableUsage();
#define CHECK_LT(val1, val2)
#define CHECK_EQ(val1, val2)
#define DCHECK_GE(val1, val2)
#define CHECK_NE(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
absl::Status FindSymmetries(std::vector< int > *node_equivalence_classes_io, std::vector< std::unique_ptr< SparsePermutation > > *generators, std::vector< int > *factorized_automorphism_group_size, TimeLimit *time_limit=nullptr)
const std::vector< int > & Support() const
void RemoveCycles(const std::vector< int > &cycle_indices)
Iterator Cycle(int i) const
static std::unique_ptr< TimeLimit > FromDeterministicTime(double deterministic_limit)
Creates a time limit object that puts limit only on the deterministic time.
SharedTimeLimit * time_limit
GurobiMPCallbackContext * context
static const int64 kint64max
Collection::value_type::second_type & LookupOrInsert(Collection *const collection, const typename Collection::value_type::first_type &key, const typename Collection::value_type::second_type &value)
bool DetectAndExploitSymmetriesInPresolve(PresolveContext *context)
std::vector< int > GetOrbitopeOrbits(int n, const std::vector< std::vector< int > > &orbitope)
void FindCpModelSymmetries(const SatParameters ¶ms, const CpModelProto &problem, std::vector< std::unique_ptr< SparsePermutation > > *generators, double deterministic_limit)
Graph * GenerateGraphForSymmetryDetection(const LinearBooleanProblem &problem, std::vector< int > *initial_equivalence_classes)
void DetectAndAddSymmetryToProto(const SatParameters ¶ms, CpModelProto *proto)
std::string ConstraintCaseName(ConstraintProto::ConstraintCase constraint_case)
std::vector< int > GetOrbits(int n, const std::vector< std::unique_ptr< SparsePermutation > > &generators)
std::vector< std::vector< int > > BasicOrbitopeExtraction(const std::vector< std::unique_ptr< SparsePermutation > > &generators)
bool RefIsPositive(int ref)
void FillDomainInProto(const Domain &domain, ProtoWithDomain *proto)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
uint64 Hash(uint64 num, uint64 c)
std::vector< int >::const_iterator begin() const
#define VLOG_IS_ON(verboselevel)