[CP-SAT] new lns (graph_dec_lns); internal code cleaning

ortools/sat/cp_constraints.cc

@@ -79,13 +79,11 @@ void BooleanXorPropagator::RegisterWith(GenericLiteralWatcher* watcher) {
 }
 
 GreaterThanAtLeastOneOfPropagator::GreaterThanAtLeastOneOfPropagator(
-    IntegerVariable target_var, const absl::Span<const IntegerVariable> vars,
-    const absl::Span<const IntegerValue> offsets,
+    IntegerVariable target_var, const absl::Span<const AffineExpression> exprs,
     const absl::Span<const Literal> selectors,
     const absl::Span<const Literal> enforcements, Model* model)
     : target_var_(target_var),
-      vars_(vars.begin(), vars.end()),
-      offsets_(offsets.begin(), offsets.end()),
+      exprs_(exprs.begin(), exprs.end()),
       selectors_(selectors.begin(), selectors.end()),
       enforcements_(enforcements.begin(), enforcements.end()),
       trail_(model->GetOrCreate<Trail>()),
@@ -103,11 +101,10 @@ bool GreaterThanAtLeastOneOfPropagator::Propagate() {
   // Propagate() calls.
   IntegerValue target_min = kMaxIntegerValue;
   const IntegerValue current_min = integer_trail_->LowerBound(target_var_);
-  for (int i = 0; i < vars_.size(); ++i) {
+  for (int i = 0; i < exprs_.size(); ++i) {
     if (trail_->Assignment().LiteralIsTrue(selectors_[i])) return true;
     if (trail_->Assignment().LiteralIsFalse(selectors_[i])) continue;
-    target_min = std::min(target_min,
-                          integer_trail_->LowerBound(vars_[i]) + offsets_[i]);
+    target_min = std::min(target_min, integer_trail_->LowerBound(exprs_[i]));
 
     // Abort if we can't get a better bound.
     if (target_min <= current_min) return true;
@@ -123,12 +120,13 @@ bool GreaterThanAtLeastOneOfPropagator::Propagate() {
   for (const Literal l : enforcements_) {
     literal_reason_.push_back(l.Negated());
   }
-  for (int i = 0; i < vars_.size(); ++i) {
+  for (int i = 0; i < exprs_.size(); ++i) {
     if (trail_->Assignment().LiteralIsFalse(selectors_[i])) {
       literal_reason_.push_back(selectors_[i]);
     } else {
-      integer_reason_.push_back(
-          IntegerLiteral::GreaterOrEqual(vars_[i], target_min - offsets_[i]));
+      if (!exprs_[i].IsConstant()) {
+        integer_reason_.push_back(exprs_[i].GreaterOrEqual(target_min));
+      }
     }
   }
   return integer_trail_->Enqueue(
@@ -141,7 +139,7 @@ void GreaterThanAtLeastOneOfPropagator::RegisterWith(
   const int id = watcher->Register(this);
   for (const Literal l : selectors_) watcher->WatchLiteral(l.Negated(), id);
   for (const Literal l : enforcements_) watcher->WatchLiteral(l, id);
-  for (const IntegerVariable v : vars_) watcher->WatchLowerBound(v, id);
+  for (const AffineExpression e : exprs_) watcher->WatchLowerBound(e, id);
 }
 
 }  // namespace sat

ortools/sat/cp_constraints.h

@@ -73,8 +73,7 @@ class BooleanXorPropagator : public PropagatorInterface {
 class GreaterThanAtLeastOneOfPropagator : public PropagatorInterface {
  public:
   GreaterThanAtLeastOneOfPropagator(IntegerVariable target_var,
-                                    absl::Span<const IntegerVariable> vars,
-                                    absl::Span<const IntegerValue> offsets,
+                                    absl::Span<const AffineExpression> exprs,
                                     absl::Span<const Literal> selectors,
                                     absl::Span<const Literal> enforcements,
                                     Model* model);
@@ -84,8 +83,7 @@ class GreaterThanAtLeastOneOfPropagator : public PropagatorInterface {
 
  private:
   const IntegerVariable target_var_;
-  const std::vector<IntegerVariable> vars_;
-  const std::vector<IntegerValue> offsets_;
+  const std::vector<AffineExpression> exprs_;
   const std::vector<Literal> selectors_;
   const std::vector<Literal> enforcements_;
 
@@ -124,28 +122,19 @@ inline std::function<void(Model*)> LiteralXorIs(
   };
 }
 
-inline std::function<void(Model*)> GreaterThanAtLeastOneOf(
-    IntegerVariable target_var, const absl::Span<const IntegerVariable> vars,
-    const absl::Span<const IntegerValue> offsets,
-    const absl::Span<const Literal> selectors) {
-  return [=](Model* model) {
-    GreaterThanAtLeastOneOfPropagator* constraint =
-        new GreaterThanAtLeastOneOfPropagator(target_var, vars, offsets,
-                                              selectors, {}, model);
-    constraint->RegisterWith(model->GetOrCreate<GenericLiteralWatcher>());
-    model->TakeOwnership(constraint);
-  };
-}
-
 inline std::function<void(Model*)> GreaterThanAtLeastOneOf(
     IntegerVariable target_var, const absl::Span<const IntegerVariable> vars,
     const absl::Span<const IntegerValue> offsets,
     const absl::Span<const Literal> selectors,
     const absl::Span<const Literal> enforcements) {
   return [=](Model* model) {
+    std::vector<AffineExpression> exprs;
+    for (int i = 0; i < vars.size(); ++i) {
+      exprs.push_back(AffineExpression(vars[i], 1, offsets[i]));
+    }
     GreaterThanAtLeastOneOfPropagator* constraint =
-        new GreaterThanAtLeastOneOfPropagator(target_var, vars, offsets,
-                                              selectors, enforcements, model);
+        new GreaterThanAtLeastOneOfPropagator(target_var, exprs, selectors,
+                                              enforcements, model);
     constraint->RegisterWith(model->GetOrCreate<GenericLiteralWatcher>());
     model->TakeOwnership(constraint);
   };
@@ -168,12 +157,13 @@ inline std::function<void(Model*)> PartialIsOneOfVar(
     const std::vector<IntegerValue> offsets(vars.size(), IntegerValue(0));
     if (vars.size() > 2) {
       // Propagate the min.
-      model->Add(GreaterThanAtLeastOneOf(target_var, vars, offsets, selectors));
+      model->Add(
+          GreaterThanAtLeastOneOf(target_var, vars, offsets, selectors, {}));
     }
     if (vars.size() > 2) {
       // Propagate the max.
-      model->Add(GreaterThanAtLeastOneOf(NegationOf(target_var),
-                                         NegationOf(vars), offsets, selectors));
+      model->Add(GreaterThanAtLeastOneOf(
+          NegationOf(target_var), NegationOf(vars), offsets, selectors, {}));
     }
   };
 }

ortools/sat/cp_model_lns.cc

@@ -51,6 +51,7 @@
 #include "ortools/sat/subsolver.h"
 #include "ortools/sat/synchronization.h"
 #include "ortools/util/adaptative_parameter_value.h"
+#include "ortools/util/integer_pq.h"
 #include "ortools/util/saturated_arithmetic.h"
 #include "ortools/util/sorted_interval_list.h"
 #include "ortools/util/strong_integers.h"
@@ -239,7 +240,9 @@ void NeighborhoodGeneratorHelper::RecomputeHelperData() {
 
     // We replace intervals by their underlying integer variables. Note that
     // this is needed for a correct decomposition into independent part.
+    bool need_sort = false;
     for (const int interval : UsedIntervals(constraints[ct_index])) {
+      need_sort = true;
       for (const int var : UsedVariables(constraints[interval])) {
         if (IsConstant(var)) continue;
         constraint_to_var_[reduced_ct_index].push_back(var);
@@ -254,6 +257,9 @@ void NeighborhoodGeneratorHelper::RecomputeHelperData() {
     }
 
     // Keep this constraint.
+    if (need_sort) {
+      gtl::STLSortAndRemoveDuplicates(&constraint_to_var_[reduced_ct_index]);
+    }
     for (const int var : constraint_to_var_[reduced_ct_index]) {
       var_to_constraint_[var].push_back(reduced_ct_index);
     }
@@ -1401,6 +1407,165 @@ Neighborhood ConstraintGraphNeighborhoodGenerator::Generate(
   return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
 }
 
+Neighborhood DecompositionGraphNeighborhoodGenerator::Generate(
+    const CpSolverResponse& initial_solution, double difficulty,
+    absl::BitGenRef random) {
+  int max_width = 0;
+  int size_at_min_width_after_100;
+  int min_width_after_100 = std::numeric_limits<int>::max();
+  int num_zero_score = 0;
+  std::vector<int> relaxed_variables;
+
+  // Note(user): The algo is slower than the other graph generator, so we
+  // might not want to lock the graph for so long? it is just a reader lock
+  // though.
+  {
+    absl::ReaderMutexLock graph_lock(&helper_.graph_mutex_);
+
+    const int num_active_vars =
+        helper_.ActiveVariablesWhileHoldingLock().size();
+    const int target_size = std::ceil(difficulty * num_active_vars);
+    if (target_size == 0) return helper_.FullNeighborhood();
+
+    const int num_vars = helper_.VarToConstraint().size();
+    const int num_constraints = helper_.ConstraintToVar().size();
+    if (num_constraints == 0 || num_vars == 0) {
+      return helper_.FullNeighborhood();
+    }
+
+    // We will grow this incrementally.
+    // Index in the graph are first variables then constraints.
+    const int num_nodes = num_vars + num_constraints;
+    std::vector<bool> added(num_nodes, false);
+    std::vector<bool> added_or_connected(num_nodes, false);
+
+    // We will process var/constraint node by minimum "score".
+    struct QueueElement {
+      int Index() const { return index; }
+      bool operator<(const QueueElement& o) const {
+        if (score == o.score) return tie_break < o.tie_break;
+        return score < o.score;
+      }
+
+      int index;
+      int score = 0;
+      double tie_break = 0.0;
+    };
+    std::vector<QueueElement> elements(num_nodes);
+    IntegerPriorityQueue<QueueElement> pq(num_nodes);
+
+    // Initialize elements.
+    for (int i = 0; i < num_nodes; ++i) {
+      elements[i].index = i;
+      elements[i].tie_break = absl::Uniform<double>(random, 0.0, 1.0);
+    }
+
+    // We start by a random node.
+    const int first_index = absl::Uniform<int>(random, 0, num_nodes);
+    elements[first_index].score =
+        first_index < num_vars
+            ? helper_.VarToConstraint()[first_index].size()
+            : helper_.ConstraintToVar()[first_index - num_vars].size();
+    pq.Add(elements[first_index]);
+    added_or_connected[first_index] = true;
+
+    // Pop max-degree from queue and update.
+    std::vector<int> to_update;
+    while (!pq.IsEmpty() && relaxed_variables.size() < target_size) {
+      // Just for logging.
+      if (relaxed_variables.size() > 100 && pq.Size() < min_width_after_100) {
+        min_width_after_100 = pq.Size();
+        size_at_min_width_after_100 = relaxed_variables.size();
+      }
+
+      const int index = pq.Top().index;
+      const int score = pq.Top().score;
+      pq.Pop();
+      added[index] = true;
+
+      // When the score is zero, we don't need to update anything since the
+      // frontier does not grow.
+      if (score == 0) {
+        if (index < num_vars) relaxed_variables.push_back(index);
+        ++num_zero_score;
+        continue;
+      }
+
+      // Note that while it might looks bad, the overall complexity of this is
+      // in O(num_edge) since we scan each index once and each newly connected
+      // vertex once.
+      int num_added = 0;
+      to_update.clear();
+      if (index < num_vars) {
+        relaxed_variables.push_back(index);
+        for (const int c : helper_.VarToConstraint()[index]) {
+          const int c_index = num_vars + c;
+          if (added_or_connected[c_index]) continue;
+          ++num_added;
+          added_or_connected[c_index] = true;
+          to_update.push_back(c_index);
+          for (const int v : helper_.ConstraintToVar()[c]) {
+            if (added[v]) continue;
+            if (added_or_connected[v]) {
+              to_update.push_back(v);
+              elements[v].score--;
+            } else {
+              elements[c_index].score++;
+            }
+          }
+        }
+      } else {
+        for (const int v : helper_.ConstraintToVar()[index - num_vars]) {
+          if (added_or_connected[v]) continue;
+          ++num_added;
+          added_or_connected[v] = true;
+          to_update.push_back(v);
+          for (const int c : helper_.VarToConstraint()[v]) {
+            if (added[num_vars + c]) continue;
+            if (added_or_connected[num_vars + c]) {
+              elements[num_vars + c].score--;
+              to_update.push_back(num_vars + c);
+            } else {
+              elements[v].score++;
+            }
+          }
+        }
+      }
+
+      // The score is exactly the frontier increase in size.
+      // This is the same as the min-degree heuristic for the elimination order.
+      // Except we only consider connected nodes.
+      CHECK_EQ(num_added, score);
+
+      gtl::STLSortAndRemoveDuplicates(&to_update);
+      for (const int index : to_update) {
+        DCHECK(!added[index]);
+        if (pq.Contains(index)) {
+          pq.ChangePriority(elements[index]);
+        } else {
+          pq.Add(elements[index]);
+        }
+      }
+
+      max_width = std::max(max_width, pq.Size());
+    }
+
+    // Just for logging.
+    if (pq.Size() < min_width_after_100) {
+      min_width_after_100 = pq.Size();
+      size_at_min_width_after_100 = relaxed_variables.size();
+    }
+
+    VLOG(2) << "#relaxed " << relaxed_variables.size() << " #zero_score "
+            << num_zero_score << " max_width " << max_width
+            << " (size,min_width)_after_100 (" << size_at_min_width_after_100
+            << "," << min_width_after_100 << ") "
+            << " final_width " << pq.Size();
+  }
+
+  return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
+}
+
 Neighborhood RelaxObjectiveVariablesGenerator::Generate(
     const CpSolverResponse& initial_solution, double difficulty,
     absl::BitGenRef random) {

ortools/sat/cp_model_lns.h

@@ -507,6 +507,26 @@ class ConstraintGraphNeighborhoodGenerator : public NeighborhoodGenerator {
                         double difficulty, absl::BitGenRef random) final;
 };
 
+// The idea here is to try to generate a random neighborhood incrementally in
+// such a way that we have at various point a "minimum connection" in term of
+// constraints or variable to the outside world.
+//
+// This is inspired by what would be a good neighborhood if one where to use
+// a tree decomposition of the constraint-variable graph with small treewidth.
+//
+// TODO(user): Doing the full heuristic treewidth decomposition is probably
+// better because when we grow the current neighborhood, just using local
+// connection to the current candidate is probably not enough to orient the
+// search towards a good final neighborhood.
+class DecompositionGraphNeighborhoodGenerator : public NeighborhoodGenerator {
+ public:
+  explicit DecompositionGraphNeighborhoodGenerator(
+      NeighborhoodGeneratorHelper const* helper, const std::string& name)
+      : NeighborhoodGenerator(name, helper) {}
+  Neighborhood Generate(const CpSolverResponse& initial_solution,
+                        double difficulty, absl::BitGenRef random) final;
+};
+
 // Pick a random subset of objective terms.
 class RelaxObjectiveVariablesGenerator : public NeighborhoodGenerator {
  public: