text: // Copyright 2010-2018 Google LLC
  // Licensed under the Apache License, Version 2.0 (the "License");
  // you may not use this file except in compliance with the License.
  // You may obtain a copy of the License at
  //
  //     http://www.apache.org/licenses/LICENSE-2.0
  //
  // Unless required by applicable law or agreed to in writing, software
  // distributed under the License is distributed on an "AS IS" BASIS,
  // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  // See the License for the specific language governing permissions and
  // limitations under the License.
  
  #include "ortools/sat/cp_model_symmetries.h"
  
  #include <memory>
  
  #include "absl/container/flat_hash_map.h"
  #include "absl/memory/memory.h"
  #include "google/protobuf/repeated_field.h"
  #include "ortools/algorithms/find_graph_symmetries.h"
  #include "ortools/base/hash.h"
  #include "ortools/base/map_util.h"
  #include "ortools/sat/cp_model_utils.h"
  #include "ortools/sat/symmetry_util.h"
  
  namespace operations_research {
  namespace sat {
  
  namespace {
  struct VectorHash {
    std::size_t operator()(const std::vector<int64>& values) const {
      size_t hash = 0;
      for (const int64 value : values) {
        hash = util_hash::Hash(value, hash);
      }
      return hash;
    }
  };
  
  // A simple class to generate equivalence class number for
  // GenerateGraphForSymmetryDetection().
  class IdGenerator {
   public:
    IdGenerator() {}
  
    // If the color was never seen before, then generate a new id, otherwise
    // return the previously generated id.
    int GetId(const std::vector<int64>& color) {
      return gtl::LookupOrInsert(&id_map_, color, id_map_.size());
    }
  
    int NextFreeId() const { return id_map_.size(); }
  
   private:
    absl::flat_hash_map<std::vector<int64>, int, VectorHash> id_map_;
  };
  
  // Appends values in `repeated_field` to `vector`.
  //
  // We use a template as proto int64 != C++ int64 in open source.
  template <typename FieldInt64Type>
  void Append(
      const google::protobuf::RepeatedField<FieldInt64Type>& repeated_field,
      std::vector<int64>* vector) {
    CHECK(vector != nullptr);
    for (const FieldInt64Type value : repeated_field) {
      vector->push_back(value);
    }
  }
  
  // Returns a graph whose automorphisms can be mapped back to the symmetries of
  // the model described in the given CpModelProto.
  //
  // Any permutation of the graph that respects the initial_equivalence_classes
  // output can be mapped to a symmetry of the given problem simply by taking its
  // restriction on the first num_variables nodes and interpreting its index as a
  // variable index. In a sense, a node with a low enough index #i is in
  // one-to-one correspondence with the variable #i (using the index
  // representation of variables).
  //
  // The format of the initial_equivalence_classes is the same as the one
  // described in GraphSymmetryFinder::FindSymmetries(). The classes must be dense
  // in [0, num_classes) and any symmetry will only map nodes with the same class
  // between each other.
  template <typename Graph>
  std::unique_ptr<Graph> GenerateGraphForSymmetryDetection(
      bool log_info, const CpModelProto& problem,
      std::vector<int>* initial_equivalence_classes) {
    CHECK(initial_equivalence_classes != nullptr);
  
    const int num_variables = problem.variables_size();
    auto graph = absl::make_unique<Graph>();
  
    // Each node will be created with a given color. Two nodes of different color
    // can never be send one into another by a symmetry. The first element of
    // the color vector will always be the NodeType.
    //
    // TODO(user): Using a full int64 for storing 3 values is not great. We
    // can optimize this at the price of a bit more code.
    enum NodeType {
      VARIABLE_NODE,
      VAR_COEFFICIENT_NODE,
      CONSTRAINT_NODE,
    };
    IdGenerator color_id_generator;
    initial_equivalence_classes->clear();
    auto new_node = [&initial_equivalence_classes, &graph,
                     &color_id_generator](const std::vector<int64>& color) {
      // Since we add nodes one by one, initial_equivalence_classes->size() gives
      // the number of nodes at any point, which we use as the next node index.
      const int node = initial_equivalence_classes->size();
      initial_equivalence_classes->push_back(color_id_generator.GetId(color));
  
      // In some corner cases, we create a node but never uses it. We still
      // want it to be there.
      graph->AddNode(node);
      return node;
    };
  
    // For two variables to be in the same equivalence class, they need to have
    // the same objective coefficient, and the same possible bounds.
    //
    // TODO(user): We could ignore the objective coefficients, and just make sure
    // that when we break symmetry amongst variables, we choose the possibility
    // with the smallest cost?
    std::vector<int64> objective_by_var(num_variables, 0);
    for (int i = 0; i < problem.objective().vars_size(); ++i) {
      const int ref = problem.objective().vars(i);
      const int var = PositiveRef(ref);
      const int64 coeff = problem.objective().coeffs(i);
      objective_by_var[var] = RefIsPositive(ref) ? coeff : -coeff;
    }
  
    // Create one node for each variable. Note that the code rely on the fact that
    // the index of a VARIABLE_NODE type is the same as the variable index.
    std::vector<int64> tmp_color;
    for (int v = 0; v < num_variables; ++v) {
      tmp_color = {VARIABLE_NODE, objective_by_var[v]};
      Append(problem.variables(v).domain(), &tmp_color);
      CHECK_EQ(v, new_node(tmp_color));
    }
  
    // We will lazily create "coefficient nodes" that correspond to a variable
    // with a given coefficient.
    absl::flat_hash_map<std::pair<int64, int64>, int> coefficient_nodes;
    auto get_coefficient_node = [&new_node, &graph, &coefficient_nodes,
                                 &tmp_color](int var, int64 coeff) {
      const int var_node = var;
      DCHECK(RefIsPositive(var));
  
      // For a coefficient of one, which are the most common, we can optimize the
      // size of the graph by omitting the coefficient node altogether and using
      // directly the var_node in this case.
      if (coeff == 1) return var_node;
  
      const auto insert =
          coefficient_nodes.insert({std::make_pair(var, coeff), 0});
      if (!insert.second) return insert.first->second;
  
      tmp_color = {VAR_COEFFICIENT_NODE, coeff};
      const int secondary_node = new_node(tmp_color);
      graph->AddArc(var_node, secondary_node);
      insert.first->second = secondary_node;
      return secondary_node;
    };
  
    // For a literal we use the same as a coefficient 1 or -1. We can do that
    // because literal and (var, coefficient) never appear together in the same
    // constraint.
    auto get_literal_node = [&get_coefficient_node](int ref) {
      return get_coefficient_node(PositiveRef(ref), RefIsPositive(ref) ? 1 : -1);
    };
  
    // Because the implications can be numerous, we encode them without
    // constraints node by using an arc from the lhs to the rhs. Note that we also
    // always add the other direction. We use a set to remove duplicates both for
    // efficiency and to not artificially break symmetries by using multi-arcs.
    //
    // Tricky: We cannot use the base variable node here to avoid situation like
    // both a variable a and b having the same children (not(a), not(b)) in the
    // graph. Because if that happen, we can permute a and b without permuting
    // their associated not(a) and not(b) node! To be sure this cannot happen, a
    // variable node can not have as children a VAR_COEFFICIENT_NODE from another
    // node. This makes sure that any permutation that touch a variable, must
    // permute its coefficient nodes accordingly.
    absl::flat_hash_set<std::pair<int, int>> implications;
    auto get_implication_node = [&new_node, &graph, &coefficient_nodes,
                                 &tmp_color](int ref) {
      const int var = PositiveRef(ref);
      const int64 coeff = RefIsPositive(ref) ? 1 : -1;
      const auto insert =
          coefficient_nodes.insert({std::make_pair(var, coeff), 0});
      if (!insert.second) return insert.first->second;
      tmp_color = {VAR_COEFFICIENT_NODE, coeff};
      const int secondary_node = new_node(tmp_color);
      graph->AddArc(var, secondary_node);
      insert.first->second = secondary_node;
      return secondary_node;
    };
    auto add_implication = [&get_implication_node, &graph, &implications](
                               int ref_a, int ref_b) {
      const auto insert = implications.insert({ref_a, ref_b});
      if (!insert.second) return;
      graph->AddArc(get_implication_node(ref_a), get_implication_node(ref_b));
  
      // Always add the other side.
      implications.insert({NegatedRef(ref_b), NegatedRef(ref_a)});
      graph->AddArc(get_implication_node(NegatedRef(ref_b)),
                    get_implication_node(NegatedRef(ref_a)));
    };
  
    // Add constraints to the graph.
    for (const ConstraintProto& constraint : problem.constraints()) {
      const int constraint_node = initial_equivalence_classes->size();
      std::vector<int64> color = {CONSTRAINT_NODE, constraint.constraint_case()};
  
      switch (constraint.constraint_case()) {
        case ConstraintProto::CONSTRAINT_NOT_SET:
          break;
        case ConstraintProto::kLinear: {
          // TODO(user): We can use the same trick as for the implications to
          // encode relations of the form coeff * var_a <= coeff * var_b without
          // creating a constraint node by directly adding an arc between the two
          // var coefficient nodes.
          Append(constraint.linear().domain(), &color);
          CHECK_EQ(constraint_node, new_node(color));
          for (int i = 0; i < constraint.linear().vars_size(); ++i) {
            const int ref = constraint.linear().vars(i);
            const int variable_node = PositiveRef(ref);
            const int64 coeff = RefIsPositive(ref)
                                    ? constraint.linear().coeffs(i)
                                    : -constraint.linear().coeffs(i);
            graph->AddArc(get_coefficient_node(variable_node, coeff),
                          constraint_node);
          }
          break;
        }
        case ConstraintProto::kBoolOr: {
          CHECK_EQ(constraint_node, new_node(color));
          for (const int ref : constraint.bool_or().literals()) {
            graph->AddArc(get_literal_node(ref), constraint_node);
          }
          break;
        }
        case ConstraintProto::kAtMostOne: {
          if (constraint.at_most_one().literals().size() == 2) {
            // Treat it as an implication to avoid creating a node.
            add_implication(constraint.at_most_one().literals(0),
                            NegatedRef(constraint.at_most_one().literals(1)));
            break;
          }
  
          CHECK_EQ(constraint_node, new_node(color));
          for (const int ref : constraint.at_most_one().literals()) {
            graph->AddArc(get_literal_node(ref), constraint_node);
          }
          break;
        }
        case ConstraintProto::kExactlyOne: {
          CHECK_EQ(constraint_node, new_node(color));
          for (const int ref : constraint.exactly_one().literals()) {
            graph->AddArc(get_literal_node(ref), constraint_node);
          }
          break;
        }
        case ConstraintProto::kBoolXor: {
          CHECK_EQ(constraint_node, new_node(color));
          for (const int ref : constraint.bool_xor().literals()) {
            graph->AddArc(get_literal_node(ref), constraint_node);
          }
          break;
        }
        case ConstraintProto::kBoolAnd: {
          // The other cases should be presolved before this is called.
          // TODO(user): not 100% true, this happen on rmatr200-p5, Fix.
          if (constraint.enforcement_literal_size() != 1) {
            if (log_info) {
              LOG(INFO) << "BoolAnd with multiple enforcement literal are not "
                           "supported in symmetry code:"
                        << constraint.ShortDebugString();
            }
            return nullptr;
          }
  
          CHECK_EQ(constraint.enforcement_literal_size(), 1);
          const int ref_a = constraint.enforcement_literal(0);
          for (const int ref_b : constraint.bool_and().literals()) {
            add_implication(ref_a, ref_b);
          }
          break;
        }
        default: {
          // If the model contains any non-supported constraints, return an empty
          // graph.
          //
          // TODO(user): support other types of constraints. Or at least, we
          // could associate to them an unique node so that their variables can
          // appear in no symmetry.
          if (log_info) {
            LOG(INFO) << "Unsupported constraint type "
                      << ConstraintCaseName(constraint.constraint_case());
          }
          return nullptr;
        }
      }
  
      // For enforcement, we use a similar trick than for the implications.
      // Because all our constraint arcs are in the direction var_node to
      // constraint_node, we just use the reverse direction for the enforcement
      // part. This way we can reuse the same get_literal_node() function.
      if (constraint.constraint_case() != ConstraintProto::kBoolAnd) {
        for (const int ref : constraint.enforcement_literal()) {
          graph->AddArc(constraint_node, get_literal_node(ref));
        }
      }
    }
  
    graph->Build();
    DCHECK_EQ(graph->num_nodes(), initial_equivalence_classes->size());
  
    // TODO(user): The symmetry code does not officially support multi-arcs. And
    // we shouldn't have any as long as there is no duplicates variable in our
    // constraints (but of course, we can't always guarantee that). That said,
    // because the symmetry code really only look at the degree, it works as long
    // as the maximum degree is bounded by num_nodes.
    const int num_nodes = graph->num_nodes();
    std::vector<int> in_degree(num_nodes, 0);
    std::vector<int> out_degree(num_nodes, 0);
    for (int i = 0; i < num_nodes; ++i) {
      out_degree[i] = graph->OutDegree(i);
      for (const int head : (*graph)[i]) {
        in_degree[head]++;
      }
    }
    for (int i = 0; i < num_nodes; ++i) {
      if (in_degree[i] >= num_nodes || out_degree[i] >= num_nodes) {
        if (log_info) LOG(INFO) << "Too many multi-arcs in symmetry code.";
        return nullptr;
      }
    }
  
    // Because this code is running during presolve, a lot a variable might have
    // no edges. We do not want to detect symmetries between these.
    //
    // Note that this code forces us to "densify" the ids afterwards because the
    // symmetry detection code relies on that.
    //
    // TODO(user): It will probably be more efficient to not even create these
    // nodes, but we will need a mapping to know the variable <-> node index.
    int next_id = color_id_generator.NextFreeId();
    for (int i = 0; i < num_variables; ++i) {
      if ((*graph)[i].empty()) {
        (*initial_equivalence_classes)[i] = next_id++;
      }
    }
  
    // Densify ids.
    int id = 0;
    std::vector<int> mapping(next_id, -1);
    for (int& ref : *initial_equivalence_classes) {
      if (mapping[ref] == -1) {
        ref = mapping[ref] = id++;
      } else {
        ref = mapping[ref];
      }
    }
  
    return graph;
  }
  }  // namespace
  
  void FindCpModelSymmetries(
      const SatParameters& params, const CpModelProto& problem,
      std::vector<std::unique_ptr<SparsePermutation>>* generators,
      double deterministic_limit) {
    const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
    CHECK(generators != nullptr);
    generators->clear();
  
    typedef GraphSymmetryFinder::Graph Graph;
  
    std::vector<int> equivalence_classes;
    std::unique_ptr<Graph> graph(GenerateGraphForSymmetryDetection<Graph>(
        log_info, problem, &equivalence_classes));
    if (graph == nullptr) return;
  
    if (log_info) {
      LOG(INFO) << "Graph for symmetry has " << graph->num_nodes()
                << " nodes and " << graph->num_arcs() << " arcs.";
    }
    if (graph->num_nodes() == 0) return;
  
    GraphSymmetryFinder symmetry_finder(*graph, /*is_undirected=*/false);
    std::vector<int> factorized_automorphism_group_size;
    std::unique_ptr<TimeLimit> time_limit =
        TimeLimit::FromDeterministicTime(deterministic_limit);
    const absl::Status status = symmetry_finder.FindSymmetries(
        &equivalence_classes, generators, &factorized_automorphism_group_size,
        time_limit.get());
  
    // TODO(user): Change the API to not return an error when the time limit is
    // reached.
    if (log_info && !status.ok()) {
      LOG(INFO) << "GraphSymmetryFinder error: " << status.message();
    }
  
    // Remove from the permutations the part not concerning the variables.
    // Note that some permutations may become empty, which means that we had
    // duplicate constraints.
    double average_support_size = 0.0;
    int num_generators = 0;
    int num_duplicate_constraints = 0;
    for (int i = 0; i < generators->size(); ++i) {
      SparsePermutation* permutation = (*generators)[i].get();
      std::vector<int> to_delete;
      for (int j = 0; j < permutation->NumCycles(); ++j) {
        // Because variable nodes are in a separate equivalence class than any
        // other node, a cycle can either contain only variable nodes or none, so
        // we just need to check one element of the cycle.
        if (*(permutation->Cycle(j).begin()) >= problem.variables_size()) {
          to_delete.push_back(j);
          if (DEBUG_MODE) {
            // Verify that the cycle's entire support does not touch any variable.
            for (const int node : permutation->Cycle(j)) {
              DCHECK_GE(node, problem.variables_size());
            }
          }
        }
      }
  
      permutation->RemoveCycles(to_delete);
      if (!permutation->Support().empty()) {
        average_support_size += permutation->Support().size();
        swap((*generators)[num_generators], (*generators)[i]);
        ++num_generators;
      } else {
        ++num_duplicate_constraints;
      }
    }
    generators->resize(num_generators);
    average_support_size /= num_generators;
    if (log_info) {
      LOG(INFO) << "Symmetry computation done. time: "
                << time_limit->GetElapsedTime()
                << " dtime: " << time_limit->GetElapsedDeterministicTime();
      if (num_generators > 0) {
        LOG(INFO) << "# of generators: " << num_generators;
        LOG(INFO) << "Average support size: " << average_support_size;
        if (num_duplicate_constraints > 0) {
          LOG(INFO) << "The model contains " << num_duplicate_constraints
                    << " duplicate constraints !";
        }
      }
    }
  }
  
  void DetectAndAddSymmetryToProto(const SatParameters& params,
                                   CpModelProto* proto) {
    const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
    SymmetryProto* symmetry = proto->mutable_symmetry();
    symmetry->Clear();
  
    std::vector<std::unique_ptr<SparsePermutation>> generators;
    FindCpModelSymmetries(params, *proto, &generators,
                          /*deterministic_limit=*/1.0);
    if (generators.empty()) return;
  
    for (const std::unique_ptr<SparsePermutation>& perm : generators) {
      SparsePermutationProto* perm_proto = symmetry->add_permutations();
      const int num_cycle = perm->NumCycles();
      for (int i = 0; i < num_cycle; ++i) {
        const int old_size = perm_proto->support().size();
        for (const int var : perm->Cycle(i)) {
          perm_proto->add_support(var);
        }
        perm_proto->add_cycle_sizes(perm_proto->support().size() - old_size);
      }
    }
  
    std::vector<std::vector<int>> orbitope = BasicOrbitopeExtraction(generators);
    if (orbitope.empty()) return;
    if (log_info) {
      LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x "
                << orbitope[0].size();
    }
    DenseMatrixProto* matrix = symmetry->add_orbitopes();
    matrix->set_num_rows(orbitope.size());
    matrix->set_num_cols(orbitope[0].size());
    for (const std::vector<int>& row : orbitope) {
      for (const int entry : row) {
        matrix->add_entries(entry);
      }
    }
  }
  
  bool DetectAndExploitSymmetriesInPresolve(PresolveContext* context) {
    const SatParameters& params = context->params();
    const bool log_info = params.log_search_progress() || VLOG_IS_ON(1);
    const CpModelProto& proto = *context->working_model;
  
    // We need to make sure the proto is up to date before computing symmetries!
    if (context->working_model->has_objective()) {
      context->WriteObjectiveToProto();
    }
    const int num_vars = proto.variables_size();
    for (int i = 0; i < num_vars; ++i) {
      FillDomainInProto(context->DomainOf(i),
                        context->working_model->mutable_variables(i));
    }
  
    // Tricky: the equivalence relation are not part of the proto.
    // We thus add them temporarily to compute the symmetry.
    int64 num_added = 0;
    const int initial_ct_index = proto.constraints().size();
    for (int var = 0; var < num_vars; ++var) {
      if (context->IsFixed(var)) continue;
      if (context->VariableWasRemoved(var)) continue;
      if (context->VariableIsNotUsedAnymore(var)) continue;
  
      const AffineRelation::Relation r = context->GetAffineRelation(var);
      if (r.representative == var) continue;
  
      ++num_added;
      ConstraintProto* ct = context->working_model->add_constraints();
      auto* arg = ct->mutable_linear();
      arg->add_vars(var);
      arg->add_coeffs(1);
      arg->add_vars(r.representative);
      arg->add_coeffs(-r.coeff);
      arg->add_domain(r.offset);
      arg->add_domain(r.offset);
    }
  
    std::vector<std::unique_ptr<SparsePermutation>> generators;
    FindCpModelSymmetries(params, proto, &generators,
                          /*deterministic_limit=*/1.0);
  
    // Remove temporary affine relation.
    context->working_model->mutable_constraints()->DeleteSubrange(
        initial_ct_index, num_added);
  
    if (generators.empty()) return true;
  
    // Collect the at most ones.
    //
    // Note(user): This relies on the fact that the pointers remain stable when
    // we adds new constraints. It should be the case, but it is a bit unsafe.
    // On the other hand it is annoying to deal with both cases below.
    std::vector<const google::protobuf::RepeatedField<int32>*> at_most_ones;
    for (int i = 0; i < proto.constraints_size(); ++i) {
      if (proto.constraints(i).constraint_case() == ConstraintProto::kAtMostOne) {
        at_most_ones.push_back(&proto.constraints(i).at_most_one().literals());
      }
      if (proto.constraints(i).constraint_case() ==
          ConstraintProto::kExactlyOne) {
        at_most_ones.push_back(&proto.constraints(i).exactly_one().literals());
      }
    }
  
    // Experimental. Generic approach. First step.
    //
    // If an at most one intersect with one or more orbit, in each intersection,
    // we can fix all but one variable to zero. For now we only test positive
    // literal, and maximize the number of fixing.
    std::vector<int> can_be_fixed_to_false;
    {
      const std::vector<int> orbits = GetOrbits(num_vars, generators);
      std::vector<int> orbit_sizes;
      int max_orbit_size = 0;
      for (const int rep : orbits) {
        if (rep == -1) continue;
        if (rep >= orbit_sizes.size()) orbit_sizes.resize(rep + 1, 0);
        orbit_sizes[rep]++;
        max_orbit_size = std::max(max_orbit_size, orbit_sizes[rep]);
      }
  
      std::vector<int> tmp_to_clear;
      std::vector<int> tmp_sizes(num_vars, 0);
      for (const google::protobuf::RepeatedField<int32>* literals :
           at_most_ones) {
        tmp_to_clear.clear();
  
        // Compute how many variables we can fix with this at most one.
        int num_fixable = 0;
        for (const int literal : *literals) {
          if (!RefIsPositive(literal)) continue;
          if (context->IsFixed(literal)) continue;
  
          const int var = PositiveRef(literal);
          const int rep = orbits[var];
          if (rep == -1) continue;
  
          // We count all but the first one in each orbit.
          if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
          if (tmp_sizes[rep] > 0) ++num_fixable;
          tmp_sizes[rep]++;
        }
  
        // Redo a pass to copy the intersection.
        if (num_fixable > can_be_fixed_to_false.size()) {
          can_be_fixed_to_false.clear();
          for (const int literal : *literals) {
            if (!RefIsPositive(literal)) continue;
            if (context->IsFixed(literal)) continue;
  
            const int var = PositiveRef(literal);
            const int rep = orbits[var];
            if (rep == -1) continue;
  
            // We push all but the first one in each orbit.
            if (tmp_sizes[rep] == 0) can_be_fixed_to_false.push_back(var);
            tmp_sizes[rep] = 0;
          }
        } else {
          // Sparse clean up.
          for (const int rep : tmp_to_clear) tmp_sizes[rep] = 0;
        }
      }
  
      if (log_info) {
        LOG(INFO) << "Num fixable by intersecting at_most_one with orbits: "
                  << can_be_fixed_to_false.size()
                  << " largest_orbit: " << max_orbit_size;
      }
    }
  
    // Orbitope approach.
    //
    // This is basically the same as the generic approach, but because of the
    // extra structure, computing the orbit of any stabilizer subgroup is easy.
    // We look for orbits intersecting at most one constraints, so we can break
    // symmetry by fixing variables.
    //
    // TODO(user): The same effect could be achieved by adding symmetry breaking
    // constraints of the form "a >= b " between Booleans and let the presolve do
    // the reduction. This might be less code, but it is also less efficient.
    // Similarly, when we cannot just fix variables to break symmetries, we could
    // add these constraints, but it is unclear if we should do it all the time or
    // not.
    //
    // TODO(user): code the generic approach with orbits and stabilizer.
    std::vector<std::vector<int>> orbitope = BasicOrbitopeExtraction(generators);
    if (!orbitope.empty() && log_info) {
      LOG(INFO) << "Found orbitope of size " << orbitope.size() << " x "
                << orbitope[0].size();
    }
  
    // Supper simple heuristic to use the orbitope or not.
    //
    // In an orbitope with an at most one on each row, we can fix the upper right
    // triangle. We could use a formula, but the loop is fast enough.
    //
    // TODO(user): Compute the stabilizer under the only non-fixed element and
    // iterate!
    int max_num_fixed_in_orbitope = 0;
    if (!orbitope.empty()) {
      const int num_rows = orbitope[0].size();
      int size_left = num_rows;
      for (int col = 0; size_left > 1 && col < orbitope.size(); ++col) {
        max_num_fixed_in_orbitope += size_left - 1;
        --size_left;
      }
    }
    if (max_num_fixed_in_orbitope < can_be_fixed_to_false.size()) {
      for (int i = 0; i < can_be_fixed_to_false.size(); ++i) {
        const int var = can_be_fixed_to_false[i];
        context->UpdateRuleStats("symmetry: fixed to false in general orbit");
        if (!context->SetLiteralToFalse(var)) return false;
      }
      return true;
    }
    if (orbitope.empty()) return true;
  
    // This will always be kept all zero after usage.
    std::vector<int> tmp_to_clear;
    std::vector<int> tmp_sizes(num_vars, 0);
    std::vector<int> tmp_num_positive(num_vars, 0);
  
    // TODO(user): The code below requires that no variable appears twice in the
    // same at most one. In particular lit and not(lit) cannot appear in the same
    // at most one.
    for (const google::protobuf::RepeatedField<int32>* literals : at_most_ones) {
      for (const int lit : *literals) {
        const int var = PositiveRef(lit);
        CHECK_NE(tmp_sizes[var], 1);
        tmp_sizes[var] = 1;
      }
      for (const int lit : *literals) {
        tmp_sizes[PositiveRef(lit)] = 0;
      }
    }
  
    while (!orbitope.empty() && orbitope[0].size() > 1) {
      const int num_cols = orbitope[0].size();
      const std::vector<int> orbits = GetOrbitopeOrbits(num_vars, orbitope);
  
      // Because in the orbitope case, we have a full symmetry group of the
      // columns, we can infer more than just using the orbits under a general
      // permutation group. If an at most one contains two variables from the
      // orbit, we can infer:
      // 1/ If the two variables appear positively, then there is an at most one
      //    on the full orbit, and we can set n - 1 variables to zero to break the
      //    symmetry.
      // 2/ If the two variables appear negatively, then the opposite situation
      //    arise and there is at most one zero on the orbit, we can set n - 1
      //    variables to one.
      // 3/ If two literals of opposite sign appear, then the only possibility
      //    for the orbit are all at one or all at zero, thus we can mark all
      //    variables as equivalent.
      //
      // These property comes from the fact that when we permute a line of the
      // orbitope in any way, then the position than ends up in the at most one
      // must never be both at one.
      //
      // Note that 1/ can be done without breaking any symmetry, but for 2/ and 3/
      // by choosing which variable is not fixed, we will break some symmetry, and
      // we will need to update the orbitope to stabilize this choice before
      // continuing.
      //
      // TODO(user): for 2/ and 3/ we could add an at most one constraint on the
      // full orbit if it is not already there!
      //
      // Note(user): On the miplib, only 1/ happens currently. Not sure with LNS
      // though.
      std::vector<bool> all_equivalent_rows(orbitope.size(), false);
  
      // The result described above can be generalized if an at most one intersect
      // many of the orbitope rows, each in at leat two positions. We will track
      // the set of best rows on which we have an at most one (or at most one
      // zero) on all their entries.
      bool at_most_one_in_best_rows;  // The alternative is at most one zero.
      int64 best_score = 0;
      std::vector<int> best_rows;
  
      std::vector<int> rows_in_at_most_one;
      for (const google::protobuf::RepeatedField<int32>* literals :
           at_most_ones) {
        tmp_to_clear.clear();
        for (const int literal : *literals) {
          if (context->IsFixed(literal)) continue;
          const int var = PositiveRef(literal);
          const int rep = orbits[var];
          if (rep == -1) continue;
  
          if (tmp_sizes[rep] == 0) tmp_to_clear.push_back(rep);
          tmp_sizes[rep]++;
          if (RefIsPositive(literal)) tmp_num_positive[rep]++;
        }
  
        int num_positive_direction = 0;
        int num_negative_direction = 0;
  
        // An at most one touching two positions in an orbitope row can possibly
        // be extended, depending if it has singleton intersection swith other
        // rows and where.
        bool possible_extension = false;
  
        rows_in_at_most_one.clear();
        for (const int row : tmp_to_clear) {
          const int size = tmp_sizes[row];
          const int num_positive = tmp_num_positive[row];
          const int num_negative = tmp_sizes[row] - tmp_num_positive[row];
          tmp_sizes[row] = 0;
          tmp_num_positive[row] = 0;
  
          if (num_positive > 1 && num_negative == 0) {
            if (size < num_cols) possible_extension = true;
            rows_in_at_most_one.push_back(row);
            ++num_positive_direction;
          } else if (num_positive == 0 && num_negative > 1) {
            if (size < num_cols) possible_extension = true;
            rows_in_at_most_one.push_back(row);
            ++num_negative_direction;
          } else if (num_positive > 0 && num_negative > 0) {
            all_equivalent_rows[row] = true;
          }
        }
  
        if (possible_extension) {
          context->UpdateRuleStats(
              "TODO symmetry: possible at most one extension.");
        }
  
        if (num_positive_direction > 0 && num_negative_direction > 0) {
          return context->NotifyThatModelIsUnsat("Symmetry and at most ones");
        }
        const bool direction = num_positive_direction > 0;
  
        // Because of symmetry, the choice of the column shouldn't matter (they
        // will all appear in the same number of constraints of the same types),
        // however we prefer to fix the variables that seems to touch more
        // constraints.
        //
        // TODO(user): maybe we should simplify the constraint using the variable
        // we fix before choosing the next row to break symmetry on. If there are
        // multiple row involved, we could also take the intersection instead of
        // probably counting the same constraints more than once.
        int64 score = 0;
        for (const int row : rows_in_at_most_one) {
          score +=
              context->VarToConstraints(PositiveRef(orbitope[row][0])).size();
        }
        if (score > best_score) {
          at_most_one_in_best_rows = direction;
          best_score = score;
          best_rows = rows_in_at_most_one;
        }
      }
  
      // Mark all the equivalence.
      // Note that this operation do not change the symmetry group.
      //
      // TODO(user): We could remove these rows from the orbitope. Note that
      // currently this never happen on the miplib (maybe in LNS though).
      for (int i = 0; i < all_equivalent_rows.size(); ++i) {
        if (all_equivalent_rows[i]) {
          for (int j = 1; j < num_cols; ++j) {
            context->StoreBooleanEqualityRelation(orbitope[i][0], orbitope[i][j]);
            context->UpdateRuleStats("symmetry: all equivalent in orbit");
            if (context->ModelIsUnsat()) return false;
          }
        }
      }
  
      // Break the symmetry on our set of best rows by picking one columns
      // and setting all the other entries to zero or one. Note that the at most
      // one applies to all entries in all rows.
      //
      // TODO(user): We don't have any at most one relation on this orbitope,
      // but we could still add symmetry breaking inequality by picking any matrix
      // entry and making it the largest/lowest value on its row. This also work
      // for non-Booleans.
      if (best_score == 0) {
        context->UpdateRuleStats(
            "TODO symmetry: add symmetry breaking inequalities?");
        break;
      }
  
      // If our symmetry group is valid, they cannot be any variable already
      // fixed to one (or zero if !at_most_one_in_best_rows). Otherwise all would
      // be fixed to one and the problem would be unsat.
      for (const int i : best_rows) {
        for (int j = 0; j < num_cols; ++j) {
          const int var = orbitope[i][j];
          if ((at_most_one_in_best_rows && context->LiteralIsTrue(var)) ||
              (!at_most_one_in_best_rows && context->LiteralIsFalse(var))) {
            return context->NotifyThatModelIsUnsat("Symmetry and at most one");
          }
        }
      }
  
      // We have an at most one on a set of rows, we will pick a column, and set
      // all other entries on these rows to zero.
      //
      // TODO(user): All choices should be equivalent, but double check?
      const int best_col = 0;
      for (const int i : best_rows) {
        for (int j = 0; j < num_cols; ++j) {
          if (j == best_col) continue;
          const int var = orbitope[i][j];
          if (at_most_one_in_best_rows) {
            context->UpdateRuleStats("symmetry: fixed to false");
            if (!context->SetLiteralToFalse(var)) return false;
          } else {
            context->UpdateRuleStats("symmetry: fixed to true");
            if (!context->SetLiteralToTrue(var)) return false;
          }
        }
      }
  
      // Remove all best rows.
      for (const int i : best_rows) orbitope[i].clear();
      int new_size = 0;
      for (int i = 0; i < orbitope.size(); ++i) {
        if (!orbitope[i].empty()) orbitope[new_size++] = orbitope[i];
      }
      CHECK_LT(new_size, orbitope.size());
      orbitope.resize(new_size);
  
      // Remove best_col.
      for (int i = 0; i < orbitope.size(); ++i) {
        std::swap(orbitope[i][best_col], orbitope[i].back());
        orbitope[i].pop_back();
      }
    }
  
    // If we are left with a set of variable than can all be permuted, lets
    // break the symmetry by ordering them.
    if (orbitope.size() == 1) {
      const int num_cols = orbitope[0].size();
      for (int i = 0; i + 1 < num_cols; ++i) {
        // Add orbitope[0][i] >= orbitope[0][i+1].
        ConstraintProto* ct = context->working_model->add_constraints();
        ct->mutable_linear()->add_coeffs(1);
        ct->mutable_linear()->add_vars(orbitope[0][i]);
        ct->mutable_linear()->add_coeffs(-1);
        ct->mutable_linear()->add_vars(orbitope[0][i + 1]);
        ct->mutable_linear()->add_domain(0);
        ct->mutable_linear()->add_domain(kint64max);
        context->UpdateRuleStats("symmetry: added symmetry breaking inequality");
      }
      context->UpdateNewConstraintsVariableUsage();
    }
  
    return true;
  }
  
  }  // namespace sat
  }  // namespace operations_research
  