OR-Tools  8.2
cp_model_symmetries.cc
Go to the documentation of this file.
1// Copyright 2010-2018 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <memory>
17
18#include "absl/container/flat_hash_map.h"
19#include "absl/memory/memory.h"
20#include "google/protobuf/repeated_field.h"
22#include "ortools/base/hash.h"
26
27namespace operations_research {
28namespace sat {
29
30namespace {
31struct VectorHash {
32 std::size_t operator()(const std::vector<int64>& values) const {
33 size_t hash = 0;
34 for (const int64 value : values) {
36 }
37 return hash;
38 }
39};
40
41// A simple class to generate equivalence class number for
42// GenerateGraphForSymmetryDetection().
43class IdGenerator {
44 public:
45 IdGenerator() {}
46
47 // If the color was never seen before, then generate a new id, otherwise
48 // return the previously generated id.
49 int GetId(const std::vector<int64>& color) {
50 return gtl::LookupOrInsert(&id_map_, color, id_map_.size());
51 }
52
53 int NextFreeId() const { return id_map_.size(); }
54
55 private:
56 absl::flat_hash_map<std::vector<int64>, int, VectorHash> id_map_;
57};
58
59// Appends values in `repeated_field` to `vector`.
60//
61// We use a template as proto int64 != C++ int64 in open source.
62template <typename FieldInt64Type>
63void Append(
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);
69 }
70}
71
72// Returns a graph whose automorphisms can be mapped back to the symmetries of
73// the model described in the given CpModelProto.
74//
75// Any permutation of the graph that respects the initial_equivalence_classes
76// output can be mapped to a symmetry of the given problem simply by taking its
77// restriction on the first num_variables nodes and interpreting its index as a
78// variable index. In a sense, a node with a low enough index #i is in
79// one-to-one correspondence with the variable #i (using the index
80// representation of variables).
81//
82// The format of the initial_equivalence_classes is the same as the one
83// described in GraphSymmetryFinder::FindSymmetries(). The classes must be dense
84// in [0, num_classes) and any symmetry will only map nodes with the same class
85// between each other.
86template <typename Graph>
87std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
88 bool log_info, const CpModelProto& problem,
89 std::vector<int>* initial_equivalence_classes) {
90 CHECK(initial_equivalence_classes != nullptr);
91
92 const int num_variables = problem.variables_size();
93 auto graph = absl::make_unique<Graph>();
94
95 // Each node will be created with a given color. Two nodes of different color
96 // can never be send one into another by a symmetry. The first element of
97 // the color vector will always be the NodeType.
98 //
99 // TODO(user): Using a full int64 for storing 3 values is not great. We
100 // can optimize this at the price of a bit more code.
101 enum NodeType {
102 VARIABLE_NODE,
103 VAR_COEFFICIENT_NODE,
104 CONSTRAINT_NODE,
105 };
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) {
110 // Since we add nodes one by one, initial_equivalence_classes->size() gives
111 // the number of nodes at any point, which we use as the next node index.
112 const int node = initial_equivalence_classes->size();
113 initial_equivalence_classes->push_back(color_id_generator.GetId(color));
114
115 // In some corner cases, we create a node but never uses it. We still
116 // want it to be there.
117 graph->AddNode(node);
118 return node;
119 };
120
121 // For two variables to be in the same equivalence class, they need to have
122 // the same objective coefficient, and the same possible bounds.
123 //
124 // TODO(user): We could ignore the objective coefficients, and just make sure
125 // that when we break symmetry amongst variables, we choose the possibility
126 // with the smallest cost?
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);
130 const int var = PositiveRef(ref);
131 const int64 coeff = problem.objective().coeffs(i);
132 objective_by_var[var] = RefIsPositive(ref) ? coeff : -coeff;
133 }
134
135 // Create one node for each variable. Note that the code rely on the fact that
136 // the index of a VARIABLE_NODE type is the same as the variable index.
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);
141 CHECK_EQ(v, new_node(tmp_color));
142 }
143
144 // We will lazily create "coefficient nodes" that correspond to a variable
145 // with a given coefficient.
146 absl::flat_hash_map<std::pair<int64, int64>, int> coefficient_nodes;
147 auto get_coefficient_node = [&new_node, &graph, &coefficient_nodes,
148 &tmp_color](int var, int64 coeff) {
149 const int var_node = var;
151
152 // For a coefficient of one, which are the most common, we can optimize the
153 // size of the graph by omitting the coefficient node altogether and using
154 // directly the var_node in this case.
155 if (coeff == 1) return var_node;
156
157 const auto insert =
158 coefficient_nodes.insert({std::make_pair(var, coeff), 0});
159 if (!insert.second) return insert.first->second;
160
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;
166 };
167
168 // For a literal we use the same as a coefficient 1 or -1. We can do that
169 // because literal and (var, coefficient) never appear together in the same
170 // constraint.
171 auto get_literal_node = [&get_coefficient_node](int ref) {
172 return get_coefficient_node(PositiveRef(ref), RefIsPositive(ref) ? 1 : -1);
173 };
174
175 // Because the implications can be numerous, we encode them without
176 // constraints node by using an arc from the lhs to the rhs. Note that we also
177 // always add the other direction. We use a set to remove duplicates both for
178 // efficiency and to not artificially break symmetries by using multi-arcs.
179 //
180 // Tricky: We cannot use the base variable node here to avoid situation like
181 // both a variable a and b having the same children (not(a), not(b)) in the
182 // graph. Because if that happen, we can permute a and b without permuting
183 // their associated not(a) and not(b) node! To be sure this cannot happen, a
184 // variable node can not have as children a VAR_COEFFICIENT_NODE from another
185 // node. This makes sure that any permutation that touch a variable, must
186 // permute its coefficient nodes accordingly.
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) {
190 const int var = PositiveRef(ref);
191 const int64 coeff = RefIsPositive(ref) ? 1 : -1;
192 const auto insert =
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;
200 };
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));
206
207 // Always add the other side.
208 implications.insert({NegatedRef(ref_b), NegatedRef(ref_a)});
209 graph->AddArc(get_implication_node(NegatedRef(ref_b)),
210 get_implication_node(NegatedRef(ref_a)));
211 };
212
213 // Add constraints to the graph.
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()};
217
218 switch (constraint.constraint_case()) {
219 case ConstraintProto::CONSTRAINT_NOT_SET:
220 break;
221 case ConstraintProto::kLinear: {
222 // TODO(user): We can use the same trick as for the implications to
223 // encode relations of the form coeff * var_a <= coeff * var_b without
224 // creating a constraint node by directly adding an arc between the two
225 // var coefficient nodes.
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);
230 const int variable_node = PositiveRef(ref);
231 const int64 coeff = RefIsPositive(ref)
232 ? constraint.linear().coeffs(i)
233 : -constraint.linear().coeffs(i);
234 graph->AddArc(get_coefficient_node(variable_node, coeff),
235 constraint_node);
236 }
237 break;
238 }
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);
243 }
244 break;
245 }
246 case ConstraintProto::kAtMostOne: {
247 if (constraint.at_most_one().literals().size() == 2) {
248 // Treat it as an implication to avoid creating a node.
249 add_implication(constraint.at_most_one().literals(0),
250 NegatedRef(constraint.at_most_one().literals(1)));
251 break;
252 }
253
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);
257 }
258 break;
259 }
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);
264 }
265 break;
266 }
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);
271 }
272 break;
273 }
274 case ConstraintProto::kBoolAnd: {
275 // The other cases should be presolved before this is called.
276 // TODO(user): not 100% true, this happen on rmatr200-p5, Fix.
277 if (constraint.enforcement_literal_size() != 1) {
278 if (log_info) {
279 LOG(INFO) << "BoolAnd with multiple enforcement literal are not "
280 "supported in symmetry code:"
281 << constraint.ShortDebugString();
282 }
283 return nullptr;
284 }
285
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);
290 }
291 break;
292 }
293 default: {
294 // If the model contains any non-supported constraints, return an empty
295 // graph.
296 //
297 // TODO(user): support other types of constraints. Or at least, we
298 // could associate to them an unique node so that their variables can
299 // appear in no symmetry.
300 if (log_info) {
301 LOG(INFO) << "Unsupported constraint type "
302 << ConstraintCaseName(constraint.constraint_case());
303 }
304 return nullptr;
305 }
306 }
307
308 // For enforcement, we use a similar trick than for the implications.
309 // Because all our constraint arcs are in the direction var_node to
310 // constraint_node, we just use the reverse direction for the enforcement
311 // part. This way we can reuse the same get_literal_node() function.
312 if (constraint.constraint_case() != ConstraintProto::kBoolAnd) {
313 for (const int ref : constraint.enforcement_literal()) {
314 graph->AddArc(constraint_node, get_literal_node(ref));
315 }
316 }
317 }
318
319 graph->Build();
320 DCHECK_EQ(graph->num_nodes(), initial_equivalence_classes->size());
321
322 // TODO(user): The symmetry code does not officially support multi-arcs. And
323 // we shouldn't have any as long as there is no duplicates variable in our
324 // constraints (but of course, we can't always guarantee that). That said,
325 // because the symmetry code really only look at the degree, it works as long
326 // as the maximum degree is bounded by num_nodes.
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]) {
333 in_degree[head]++;
334 }
335 }
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.";
339 return nullptr;
340 }
341 }
342
343 // Because this code is running during presolve, a lot a variable might have
344 // no edges. We do not want to detect symmetries between these.
345 //
346 // Note that this code forces us to "densify" the ids afterwards because the
347 // symmetry detection code relies on that.
348 //
349 // TODO(user): It will probably be more efficient to not even create these
350 // nodes, but we will need a mapping to know the variable <-> node index.
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++;
355 }
356 }
357
358 // Densify ids.
359 int id = 0;
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++;
364 } else {
365 ref = mapping[ref];
366 }
367 }
368
369 return graph;
370}
371} // namespace
372
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);
379 generators->clear();
380
382
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;
387
388 if (log_info) {
389 LOG(INFO) << "Graph for symmetry has " << graph->num_nodes()
390 << " nodes and " << graph->num_arcs() << " arcs.";
391 }
392 if (graph->num_nodes() == 0) return;
393
394 GraphSymmetryFinder symmetry_finder(*graph, /*is_undirected=*/false);
395 std::vector<int> factorized_automorphism_group_size;
396 std::unique_ptr<TimeLimit> time_limit =
397 TimeLimit::FromDeterministicTime(deterministic_limit);
398 const absl::Status status = symmetry_finder.FindSymmetries(
399 &equivalence_classes, generators, &factorized_automorphism_group_size,
400 time_limit.get());
401
402 // TODO(user): Change the API to not return an error when the time limit is
403 // reached.
404 if (log_info && !status.ok()) {
405 LOG(INFO) << "GraphSymmetryFinder error: " << status.message();
406 }
407
408 // Remove from the permutations the part not concerning the variables.
409 // Note that some permutations may become empty, which means that we had
410 // duplicate constraints.
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) {
415 SparsePermutation* permutation = (*generators)[i].get();
416 std::vector<int> to_delete;
417 for (int j = 0; j < permutation->NumCycles(); ++j) {
418 // Because variable nodes are in a separate equivalence class than any
419 // other node, a cycle can either contain only variable nodes or none, so
420 // we just need to check one element of the cycle.
421 if (*(permutation->Cycle(j).begin()) >= problem.variables_size()) {
422 to_delete.push_back(j);
423 if (DEBUG_MODE) {
424 // Verify that the cycle's entire support does not touch any variable.
425 for (const int node : permutation->Cycle(j)) {
426 DCHECK_GE(node, problem.variables_size());
427 }
428 }
429 }
430 }
431
432 permutation->RemoveCycles(to_delete);
433 if (!permutation->Support().empty()) {
434 average_support_size += permutation->Support().size();
435 swap((*generators)[num_generators], (*generators)[i]);
436 ++num_generators;
437 } else {
438 ++num_duplicate_constraints;
439 }
440 }
441 generators->resize(num_generators);
442 average_support_size /= num_generators;
443 if (log_info) {
444 LOG(INFO) << "Symmetry computation done. time: "
445 << time_limit->GetElapsedTime()
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 !";
453 }
454 }
455 }
456}
457
458void DetectAndAddSymmetryToProto(const SatParameters& params,
459 CpModelProto* proto) {
460 const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
461 SymmetryProto* symmetry = proto->mutable_symmetry();
462 symmetry->Clear();
463
464 std::vector<std::unique_ptr<SparsePermutation>> generators;
465 FindCpModelSymmetries(params, *proto, &generators,
466 /*deterministic_limit=*/1.0);
467 if (generators.empty()) return;
468
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);
476 }
477 perm_proto->add_cycle_sizes(perm_proto->support().size() - old_size);
478 }
479 }
480
481 std::vector<std::vector<int>> orbitope = BasicOrbitopeExtraction(generators);
482 if (orbitope.empty()) return;
483 if (log_info) {
484 LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x "
485 << orbitope[0].size();
486 }
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);
493 }
494 }
495}
496
498 const SatParameters& params = context->params();
499 const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
500 const CpModelProto& proto = *context->working_model;
501
502 // We need to make sure the proto is up to date before computing symmetries!
503 if (context->working_model->has_objective()) {
504 context->WriteObjectiveToProto();
505 }
506 const int num_vars = proto.variables_size();
507 for (int i = 0; i < num_vars; ++i) {
508 FillDomainInProto(context->DomainOf(i),
509 context->working_model->mutable_variables(i));
510 }
511
512 // Tricky: the equivalence relation are not part of the proto.
513 // We thus add them temporarily to compute the symmetry.
514 int64 num_added = 0;
515 const int initial_ct_index = proto.constraints().size();
516 for (int var = 0; var < num_vars; ++var) {
517 if (context->IsFixed(var)) continue;
518 if (context->VariableWasRemoved(var)) continue;
519 if (context->VariableIsNotUsedAnymore(var)) continue;
520
521 const AffineRelation::Relation r = context->GetAffineRelation(var);
522 if (r.representative == var) continue;
523
524 ++num_added;
525 ConstraintProto* ct = context->working_model->add_constraints();
526 auto* arg = ct->mutable_linear();
527 arg->add_vars(var);
528 arg->add_coeffs(1);
529 arg->add_vars(r.representative);
530 arg->add_coeffs(-r.coeff);
531 arg->add_domain(r.offset);
532 arg->add_domain(r.offset);
533 }
534
535 std::vector<std::unique_ptr<SparsePermutation>> generators;
536 FindCpModelSymmetries(params, proto, &generators,
537 /*deterministic_limit=*/1.0);
538
539 // Remove temporary affine relation.
540 context->working_model->mutable_constraints()->DeleteSubrange(
541 initial_ct_index, num_added);
542
543 if (generators.empty()) return true;
544
545 // Collect the at most ones.
546 //
547 // Note(user): This relies on the fact that the pointers remain stable when
548 // we adds new constraints. It should be the case, but it is a bit unsafe.
549 // On the other hand it is annoying to deal with both cases below.
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());
554 }
555 if (proto.constraints(i).constraint_case() ==
556 ConstraintProto::kExactlyOne) {
557 at_most_ones.push_back(&proto.constraints(i).exactly_one().literals());
558 }
559 }
560
561 // Experimental. Generic approach. First step.
562 //
563 // If an at most one intersect with one or more orbit, in each intersection,
564 // we can fix all but one variable to zero. For now we only test positive
565 // literal, and maximize the number of fixing.
566 std::vector<int> can_be_fixed_to_false;
567 {
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);
574 orbit_sizes[rep]++;
575 max_orbit_size = std::max(max_orbit_size, orbit_sizes[rep]);
576 }
577
578 std::vector<int> tmp_to_clear;
579 std::vector<int> tmp_sizes(num_vars, 0);
580 for (const google::protobuf::RepeatedField<int32>* literals :
581 at_most_ones) {
582 tmp_to_clear.clear();
583
584 // Compute how many variables we can fix with this at most one.
585 int num_fixable = 0;
586 for (const int literal : *literals) {
587 if (!RefIsPositive(literal)) continue;
588 if (context->IsFixed(literal)) continue;
589
590 const int var = PositiveRef(literal);
591 const int rep = orbits[var];
592 if (rep == -1) continue;
593
594 // We count all but the first one in each orbit.
595 if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
596 if (tmp_sizes[rep] > 0) ++num_fixable;
597 tmp_sizes[rep]++;
598 }
599
600 // Redo a pass to copy the intersection.
601 if (num_fixable > can_be_fixed_to_false.size()) {
602 can_be_fixed_to_false.clear();
603 for (const int literal : *literals) {
604 if (!RefIsPositive(literal)) continue;
605 if (context->IsFixed(literal)) continue;
606
607 const int var = PositiveRef(literal);
608 const int rep = orbits[var];
609 if (rep == -1) continue;
610
611 // We push all but the first one in each orbit.
612 if (tmp_sizes[rep] == 0) can_be_fixed_to_false.push_back(var);
613 tmp_sizes[rep] = 0;
614 }
615 } else {
616 // Sparse clean up.
617 for (const int rep : tmp_to_clear) tmp_sizes[rep] = 0;
618 }
619 }
620
621 if (log_info) {
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;
625 }
626 }
627
628 // Orbitope approach.
629 //
630 // This is basically the same as the generic approach, but because of the
631 // extra structure, computing the orbit of any stabilizer subgroup is easy.
632 // We look for orbits intersecting at most one constraints, so we can break
633 // symmetry by fixing variables.
634 //
635 // TODO(user): The same effect could be achieved by adding symmetry breaking
636 // constraints of the form "a >= b " between Booleans and let the presolve do
637 // the reduction. This might be less code, but it is also less efficient.
638 // Similarly, when we cannot just fix variables to break symmetries, we could
639 // add these constraints, but it is unclear if we should do it all the time or
640 // not.
641 //
642 // TODO(user): code the generic approach with orbits and stabilizer.
643 std::vector<std::vector<int>> orbitope = BasicOrbitopeExtraction(generators);
644 if (!orbitope.empty() && log_info) {
645 LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x "
646 << orbitope[0].size();
647 }
648
649 // Supper simple heuristic to use the orbitope or not.
650 //
651 // In an orbitope with an at most one on each row, we can fix the upper right
652 // triangle. We could use a formula, but the loop is fast enough.
653 //
654 // TODO(user): Compute the stabilizer under the only non-fixed element and
655 // iterate!
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;
662 --size_left;
663 }
664 }
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;
670 }
671 return true;
672 }
673 if (orbitope.empty()) return true;
674
675 // This will always be kept all zero after usage.
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);
679
680 // TODO(user): The code below requires that no variable appears twice in the
681 // same at most one. In particular lit and not(lit) cannot appear in the same
682 // at most one.
683 for (const google::protobuf::RepeatedField<int32>* literals : at_most_ones) {
684 for (const int lit : *literals) {
685 const int var = PositiveRef(lit);
686 CHECK_NE(tmp_sizes[var], 1);
687 tmp_sizes[var] = 1;
688 }
689 for (const int lit : *literals) {
690 tmp_sizes[PositiveRef(lit)] = 0;
691 }
692 }
693
694 while (!orbitope.empty() && orbitope[0].size() > 1) {
695 const int num_cols = orbitope[0].size();
696 const std::vector<int> orbits = GetOrbitopeOrbits(num_vars, orbitope);
697
698 // Because in the orbitope case, we have a full symmetry group of the
699 // columns, we can infer more than just using the orbits under a general
700 // permutation group. If an at most one contains two variables from the
701 // orbit, we can infer:
702 // 1/ If the two variables appear positively, then there is an at most one
703 // on the full orbit, and we can set n - 1 variables to zero to break the
704 // symmetry.
705 // 2/ If the two variables appear negatively, then the opposite situation
706 // arise and there is at most one zero on the orbit, we can set n - 1
707 // variables to one.
708 // 3/ If two literals of opposite sign appear, then the only possibility
709 // for the orbit are all at one or all at zero, thus we can mark all
710 // variables as equivalent.
711 //
712 // These property comes from the fact that when we permute a line of the
713 // orbitope in any way, then the position than ends up in the at most one
714 // must never be both at one.
715 //
716 // Note that 1/ can be done without breaking any symmetry, but for 2/ and 3/
717 // by choosing which variable is not fixed, we will break some symmetry, and
718 // we will need to update the orbitope to stabilize this choice before
719 // continuing.
720 //
721 // TODO(user): for 2/ and 3/ we could add an at most one constraint on the
722 // full orbit if it is not already there!
723 //
724 // Note(user): On the miplib, only 1/ happens currently. Not sure with LNS
725 // though.
726 std::vector<bool> all_equivalent_rows(orbitope.size(), false);
727
728 // The result described above can be generalized if an at most one intersect
729 // many of the orbitope rows, each in at leat two positions. We will track
730 // the set of best rows on which we have an at most one (or at most one
731 // zero) on all their entries.
732 bool at_most_one_in_best_rows; // The alternative is at most one zero.
733 int64 best_score = 0;
734 std::vector<int> best_rows;
735
736 std::vector<int> rows_in_at_most_one;
737 for (const google::protobuf::RepeatedField<int32>* literals :
738 at_most_ones) {
739 tmp_to_clear.clear();
740 for (const int literal : *literals) {
741 if (context->IsFixed(literal)) continue;
742 const int var = PositiveRef(literal);
743 const int rep = orbits[var];
744 if (rep == -1) continue;
745
746 if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
747 tmp_sizes[rep]++;
748 if (RefIsPositive(literal)) tmp_num_positive[rep]++;
749 }
750
751 int num_positive_direction = 0;
752 int num_negative_direction = 0;
753
754 // An at most one touching two positions in an orbitope row can possibly
755 // be extended, depending if it has singleton intersection swith other
756 // rows and where.
757 bool possible_extension = false;
758
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];
764 tmp_sizes[row] = 0;
765 tmp_num_positive[row] = 0;
766
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;
777 }
778 }
779
780 if (possible_extension) {
781 context->UpdateRuleStats(
782 "TODO symmetry: possible at most one extension.");
783 }
784
785 if (num_positive_direction > 0 && num_negative_direction > 0) {
786 return context->NotifyThatModelIsUnsat("Symmetry and at most ones");
787 }
788 const bool direction = num_positive_direction > 0;
789
790 // Because of symmetry, the choice of the column shouldn't matter (they
791 // will all appear in the same number of constraints of the same types),
792 // however we prefer to fix the variables that seems to touch more
793 // constraints.
794 //
795 // TODO(user): maybe we should simplify the constraint using the variable
796 // we fix before choosing the next row to break symmetry on. If there are
797 // multiple row involved, we could also take the intersection instead of
798 // probably counting the same constraints more than once.
799 int64 score = 0;
800 for (const int row : rows_in_at_most_one) {
801 score +=
802 context->VarToConstraints(PositiveRef(orbitope[row][0])).size();
803 }
804 if (score > best_score) {
805 at_most_one_in_best_rows = direction;
806 best_score = score;
807 best_rows = rows_in_at_most_one;
808 }
809 }
810
811 // Mark all the equivalence.
812 // Note that this operation do not change the symmetry group.
813 //
814 // TODO(user): We could remove these rows from the orbitope. Note that
815 // currently this never happen on the miplib (maybe in LNS though).
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;
822 }
823 }
824 }
825
826 // Break the symmetry on our set of best rows by picking one columns
827 // and setting all the other entries to zero or one. Note that the at most
828 // one applies to all entries in all rows.
829 //
830 // TODO(user): We don't have any at most one relation on this orbitope,
831 // but we could still add symmetry breaking inequality by picking any matrix
832 // entry and making it the largest/lowest value on its row. This also work
833 // for non-Booleans.
834 if (best_score == 0) {
835 context->UpdateRuleStats(
836 "TODO symmetry: add symmetry breaking inequalities?");
837 break;
838 }
839
840 // If our symmetry group is valid, they cannot be any variable already
841 // fixed to one (or zero if !at_most_one_in_best_rows). Otherwise all would
842 // be fixed to one and the problem would be unsat.
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");
849 }
850 }
851 }
852
853 // We have an at most one on a set of rows, we will pick a column, and set
854 // all other entries on these rows to zero.
855 //
856 // TODO(user): All choices should be equivalent, but double check?
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;
865 } else {
866 context->UpdateRuleStats("symmetry: fixed to true");
867 if (!context->SetLiteralToTrue(var)) return false;
868 }
869 }
870 }
871
872 // Remove all best rows.
873 for (const int i : best_rows) orbitope[i].clear();
874 int new_size = 0;
875 for (int i = 0; i < orbitope.size(); ++i) {
876 if (!orbitope[i].empty()) orbitope[new_size++] = orbitope[i];
877 }
878 CHECK_LT(new_size, orbitope.size());
879 orbitope.resize(new_size);
880
881 // Remove best_col.
882 for (int i = 0; i < orbitope.size(); ++i) {
883 std::swap(orbitope[i][best_col], orbitope[i].back());
884 orbitope[i].pop_back();
885 }
886 }
887
888 // If we are left with a set of variable than can all be permuted, lets
889 // break the symmetry by ordering them.
890 if (orbitope.size() == 1) {
891 const int num_cols = orbitope[0].size();
892 for (int i = 0; i + 1 < num_cols; ++i) {
893 // Add orbitope[0][i] >= orbitope[0][i+1].
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);
900 ct->mutable_linear()->add_domain(kint64max);
901 context->UpdateRuleStats("symmetry: added symmetry breaking inequality");
902 }
903 context->UpdateNewConstraintsVariableUsage();
904 }
905
906 return true;
907}
908
909} // namespace sat
910} // namespace operations_research
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define CHECK_LT(val1, val2)
Definition: base/logging.h:700
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define CHECK_NE(val1, val2)
Definition: base/logging.h:698
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
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)
static std::unique_ptr< TimeLimit > FromDeterministicTime(double deterministic_limit)
Creates a time limit object that puts limit only on the deterministic time.
Definition: time_limit.h:144
CpModelProto proto
SharedTimeLimit * time_limit
const Constraint * ct
int64 value
IntVar * var
Definition: expr_array.cc:1858
GurobiMPCallbackContext * context
static const int64 kint64max
int64_t int64
const int INFO
Definition: log_severity.h:31
const bool DEBUG_MODE
Definition: macros.h:24
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
int64 hash
Definition: matrix_utils.cc:60
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)
Definition: map_util.h:207
bool DetectAndExploitSymmetriesInPresolve(PresolveContext *context)
std::vector< int > GetOrbitopeOrbits(int n, const std::vector< std::vector< int > > &orbitope)
void FindCpModelSymmetries(const SatParameters &params, 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 &params, 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)
Definition: hash.h:150
ListGraph Graph
Definition: graph.h:2360
Literal literal
Definition: optimization.cc:84
int64 head
std::vector< int >::const_iterator begin() const
#define VLOG_IS_ON(verboselevel)
Definition: vlog_is_on.h:41