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"
23 #include "ortools/base/map_util.h"
26 
27 namespace operations_research {
28 namespace sat {
29 
30 namespace {
31 struct 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().
43 class 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.
62 template <typename FieldInt64Type>
63 void 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.
86 template <typename Graph>
87 std::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 
458 void 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)
void RemoveCycles(const std::vector< int > &cycle_indices)
const std::vector< int > & Support() const
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)
bool RefIsPositive(int ref)
std::vector< int > GetOrbitopeOrbits(int n, const std::vector< std::vector< int >> &orbitope)
Graph * GenerateGraphForSymmetryDetection(const LinearBooleanProblem &problem, std::vector< int > *initial_equivalence_classes)
void DetectAndAddSymmetryToProto(const SatParameters &params, CpModelProto *proto)
void FillDomainInProto(const Domain &domain, ProtoWithDomain *proto)
void FindCpModelSymmetries(const SatParameters &params, const CpModelProto &problem, std::vector< std::unique_ptr< SparsePermutation >> *generators, double deterministic_limit)
std::string ConstraintCaseName(ConstraintProto::ConstraintCase constraint_case)
std::vector< std::vector< int > > BasicOrbitopeExtraction(const std::vector< std::unique_ptr< SparsePermutation >> &generators)
std::vector< int > GetOrbits(int n, const std::vector< std::unique_ptr< SparsePermutation >> &generators)
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