diff --git a/ortools/sat/2d_orthogonal_packing.cc b/ortools/sat/2d_orthogonal_packing.cc
new file mode 100644
index 0000000000..5a93954b6c
--- /dev/null
+++ b/ortools/sat/2d_orthogonal_packing.cc
@@ -0,0 +1,287 @@
+// Copyright 2010-2024 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/2d_orthogonal_packing.h"
+
+#include <algorithm>
+#include <cstdint>
+#include <optional>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "absl/log/check.h"
+#include "absl/types/span.h"
+#include "ortools/base/vlog_is_on.h"
+#include "ortools/sat/integer.h"
+
+namespace operations_research {
+namespace sat {
+
+OrthogonalPackingInfeasibilityDetector::
+    ~OrthogonalPackingInfeasibilityDetector() {
+  if (!VLOG_IS_ON(1)) return;
+  std::vector<std::pair<std::string, int64_t>> stats;
+  stats.push_back(
+      {"OrthogonalPackingInfeasibilityDetector/called", num_calls_});
+  stats.push_back(
+      {"OrthogonalPackingInfeasibilityDetector/conflicts", num_conflicts_});
+  stats.push_back({"OrthogonalPackingInfeasibilityDetector/trivial_conflicts",
+                   num_trivial_conflicts_});
+  stats.push_back({"OrthogonalPackingInfeasibilityDetector/conflicts_two_items",
+                   num_conflicts_two_items_});
+
+  shared_stats_->AddStats(stats);
+}
+
+namespace {
+std::optional<std::pair<int, int>> FindPairwiseConflict(
+    absl::Span<const IntegerValue> sizes_x,
+    absl::Span<const IntegerValue> sizes_y,
+    std::pair<IntegerValue, IntegerValue> bounding_box_size,
+    const std::vector<int>& index_by_decreasing_x_size,
+    const std::vector<int>& index_by_decreasing_y_size) {
+  // Look for pairwise incompatible pairs by using the logic such conflict can
+  // only happen between a "tall" item a "wide" item.
+  int x_idx = 0;
+  int y_idx = 0;
+  while (x_idx < index_by_decreasing_x_size.size() &&
+         y_idx < index_by_decreasing_y_size.size()) {
+    if (index_by_decreasing_x_size[x_idx] ==
+        index_by_decreasing_y_size[y_idx]) {
+      if (sizes_x[index_by_decreasing_x_size[x_idx]] >
+          sizes_y[index_by_decreasing_x_size[x_idx]]) {
+        y_idx++;
+      } else {
+        x_idx++;
+      }
+      continue;
+    }
+    const bool overlap_on_x = (sizes_x[index_by_decreasing_x_size[x_idx]] +
+                                   sizes_x[index_by_decreasing_y_size[y_idx]] >
+                               bounding_box_size.first);
+    const bool overlap_on_y = (sizes_y[index_by_decreasing_x_size[x_idx]] +
+                                   sizes_y[index_by_decreasing_y_size[y_idx]] >
+                               bounding_box_size.second);
+    if (overlap_on_x && overlap_on_y) {
+      return std::make_pair(index_by_decreasing_x_size[x_idx],
+                            index_by_decreasing_y_size[y_idx]);
+    } else if (overlap_on_x) {
+      x_idx++;
+    } else if (overlap_on_y) {
+      y_idx++;
+    } else {
+      y_idx++;
+    }
+  }
+  return std::nullopt;
+}
+
+// Check for conflict using the f_0 dual feasible function described on Côté,
+// Jean-François, Mohamed Haouari, and Manuel Iori. "A primal decomposition
+// algorithm for the two-dimensional bin packing problem." arXiv preprint
+// arXiv:1909.06835 (2019). This function tries all possible values of the `k`
+// parameter.
+//
+// The `f_0^k(x)` function for a total bin width `C` and a parameter `k` taking
+// values in [0, C/2] is defined as:
+//
+//            / C, if x > C - k,
+// f_0^k(x) = | x, if k <= x <= C - k,
+//            \ 0, if x < k.
+//
+// If a conflict is found, this replaces the conflict in `result` if it detects
+// a conflict using less items or if `result` is empty.
+// The algorithm is the same if we swap the x and y dimension.
+void GetOrReplaceDffConflict(absl::Span<const IntegerValue> sizes_x,
+                             absl::Span<const IntegerValue> sizes_y,
+                             absl::Span<const int> index_by_decreasing_x_size,
+                             IntegerValue x_bb_size, IntegerValue total_energy,
+                             IntegerValue bb_area, std::vector<int>* result) {
+  // If we found a conflict for a k parameter, which is rare, recompute the
+  // total used energy consumed by the items to find the minimal set of
+  // conflicting items.
+  int num_items = sizes_x.size();
+  auto add_result_if_better = [&result, &sizes_x, &sizes_y, num_items,
+                               &x_bb_size, &bb_area](const IntegerValue k) {
+    std::vector<std::pair<int, IntegerValue>> index_to_energy;
+    index_to_energy.reserve(num_items);
+    for (int i = 0; i < num_items; i++) {
+      IntegerValue x_size = sizes_x[i];
+      if (x_size > x_bb_size - k) {
+        x_size = x_bb_size;
+      } else if (x_size < k) {
+        continue;
+      }
+      index_to_energy.push_back({i, x_size * sizes_y[i]});
+    }
+    std::sort(index_to_energy.begin(), index_to_energy.end(),
+              [](const std::pair<int, IntegerValue>& a,
+                 const std::pair<int, IntegerValue>& b) {
+                return a.second > b.second;
+              });
+    IntegerValue recomputed_energy = 0;
+    for (int i = 0; i < index_to_energy.size(); i++) {
+      recomputed_energy += index_to_energy[i].second;
+      if (recomputed_energy > bb_area) {
+        if (result->empty() || i + 1 < result->size()) {
+          result->resize(i + 1);
+          for (int j = 0; j <= i; j++) {
+            (*result)[j] = index_to_energy[j].first;
+          }
+        }
+        break;
+      }
+    }
+  };
+
+  // One thing we use in this implementation is that not all values of k are
+  // interesting: what can cause an energy conflict is increasing the size of
+  // the large items, removing the small ones makes it less constrained and we
+  // do it only to preserve correctness. Thus, it is enough to check the values
+  // of k that are just small enough to enlarge a large item. That means that
+  // large items and small ones are not symmetric with respect to what values of
+  // k are important.
+  IntegerValue current_energy = total_energy;
+  // We keep an index on the largest item yet-to-be enlarged and the smallest
+  // one yet-to-be removed.
+  int removing_item_index = index_by_decreasing_x_size.size() - 1;
+  int enlarging_item_index = 0;
+  while (enlarging_item_index < index_by_decreasing_x_size.size()) {
+    int index = index_by_decreasing_x_size[enlarging_item_index];
+    IntegerValue size = sizes_x[index];
+    if (size <= x_bb_size / 2) {
+      break;
+    }
+    // Note that since `size_x` is decreasing, we test increasingly large
+    // values of k. Also note that a item with size `k` cannot fit alongside
+    // an item with size `size_x`, but smaller ones can.
+    const IntegerValue k = x_bb_size - size + 1;
+    // First, add the area contribution of enlarging the all the items of size
+    // exactly size_x. All larger items were already enlarged in the previous
+    // iterations.
+    do {
+      index = index_by_decreasing_x_size[enlarging_item_index];
+      size = sizes_x[index];
+      current_energy += (x_bb_size - size) * sizes_y[index];
+      enlarging_item_index++;
+    } while (enlarging_item_index < index_by_decreasing_x_size.size() &&
+             sizes_x[index_by_decreasing_x_size[enlarging_item_index]] == size);
+
+    // Now remove the area contribution of removing all the items smaller than
+    // k that were not removed before.
+    while (removing_item_index >= 0 &&
+           sizes_x[index_by_decreasing_x_size[removing_item_index]] < k) {
+      const int remove_idx = index_by_decreasing_x_size[removing_item_index];
+      current_energy -= sizes_x[remove_idx] * sizes_y[remove_idx];
+      removing_item_index--;
+    }
+
+    if (current_energy > bb_area) {
+      add_result_if_better(k);
+    }
+  }
+}
+
+}  // namespace
+
+OrthogonalPackingInfeasibilityDetector::Result
+OrthogonalPackingInfeasibilityDetector::TestFeasibility(
+    absl::Span<const IntegerValue> sizes_x,
+    absl::Span<const IntegerValue> sizes_y,
+    std::pair<IntegerValue, IntegerValue> bounding_box_size) {
+  num_calls_++;
+
+  const int num_items = sizes_x.size();
+  DCHECK_EQ(num_items, sizes_y.size());
+  const IntegerValue bb_area =
+      bounding_box_size.first * bounding_box_size.second;
+  IntegerValue total_energy = 0;
+
+  index_by_decreasing_x_size_.resize(num_items);
+  index_by_decreasing_y_size_.resize(num_items);
+  for (int i = 0; i < num_items; i++) {
+    total_energy += sizes_x[i] * sizes_y[i];
+    index_by_decreasing_x_size_[i] = i;
+    index_by_decreasing_y_size_[i] = i;
+    if (sizes_x[i] > bounding_box_size.first ||
+        sizes_y[i] > bounding_box_size.second) {
+      num_conflicts_++;
+      return Result{.result = Status::INFEASIBLE,
+                    .minimum_unfeasible_subproblem = {i}};
+    }
+  }
+
+  if (num_items <= 1) {
+    return Result{.result = Status::FEASIBLE};
+  }
+
+  std::sort(index_by_decreasing_x_size_.begin(),
+            index_by_decreasing_x_size_.end(),
+            [&sizes_x](int a, int b) { return sizes_x[a] > sizes_x[b]; });
+  std::sort(index_by_decreasing_y_size_.begin(),
+            index_by_decreasing_y_size_.end(),
+            [&sizes_y](int a, int b) { return sizes_y[a] > sizes_y[b]; });
+
+  // First look for pairwise incompatible pairs.
+  if (auto pair = FindPairwiseConflict(sizes_x, sizes_y, bounding_box_size,
+                                       index_by_decreasing_x_size_,
+                                       index_by_decreasing_y_size_);
+      pair.has_value()) {
+    num_conflicts_++;
+    num_conflicts_two_items_++;
+    return Result{.result = Status::INFEASIBLE,
+                  .minimum_unfeasible_subproblem = {pair.value().first,
+                                                    pair.value().second}};
+  }
+
+  if (num_items == 2) {
+    return Result{.result = Status::FEASIBLE};
+  }
+
+  // If there is no pairwise incompatible pairs, this DFF cannot find a conflict
+  // by enlarging a item on both x and y directions: this would create an item
+  // as long as the whole box and another item as high as the whole box, which
+  // is obviously incompatible, and this incompatibility would be present
+  // already before enlarging the items since it is a DFF. So it is enough to
+  // test making items wide or high, but no need to try both.
+
+  std::vector<int> result;
+  GetOrReplaceDffConflict(sizes_x, sizes_y, index_by_decreasing_x_size_,
+                          bounding_box_size.first, total_energy, bb_area,
+                          &result);
+  if (result.size() == 3) {
+    // Since we already checked pairwise conflict, we won't find anything
+    // narrower than a three-item conflict.
+    num_conflicts_++;
+    return Result{.result = Status::INFEASIBLE,
+                  .minimum_unfeasible_subproblem = result};
+  }
+  GetOrReplaceDffConflict(sizes_y, sizes_x, index_by_decreasing_y_size_,
+                          bounding_box_size.second, total_energy, bb_area,
+                          &result);
+
+  // The items have a trivial conflict and we didn't found any tighter one.
+  if (total_energy > bb_area && result.size() == num_items) {
+    num_trivial_conflicts_++;
+  }
+  if (!result.empty()) {
+    num_conflicts_++;
+    return Result{.result = Status::INFEASIBLE,
+                  .minimum_unfeasible_subproblem = result};
+  }
+  return Result{.result = Status::UNKNOWN};
+}
+
+}  // namespace sat
+}  // namespace operations_research
diff --git a/ortools/sat/2d_orthogonal_packing.h b/ortools/sat/2d_orthogonal_packing.h
new file mode 100644
index 0000000000..4a68b7574a
--- /dev/null
+++ b/ortools/sat/2d_orthogonal_packing.h
@@ -0,0 +1,68 @@
+// Copyright 2010-2024 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.
+
+#ifndef OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_H_
+#define OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_H_
+
+#include <cstdint>
+#include <utility>
+#include <vector>
+
+#include "absl/types/span.h"
+#include "ortools/sat/integer.h"
+#include "ortools/sat/synchronization.h"
+
+namespace operations_research {
+namespace sat {
+
+// Class for solving the orthogonal packing problem when it can be done
+// efficiently (ie., not applying any heuristic slower than O(N^2)).
+class OrthogonalPackingInfeasibilityDetector {
+ public:
+  explicit OrthogonalPackingInfeasibilityDetector(
+      SharedStatistics* shared_stats)
+      : shared_stats_(shared_stats) {}
+
+  ~OrthogonalPackingInfeasibilityDetector();
+
+  enum class Status {
+    INFEASIBLE,
+    FEASIBLE,
+    UNKNOWN,
+  };
+  struct Result {
+    Status result;
+    std::vector<int> minimum_unfeasible_subproblem;
+  };
+
+  Result TestFeasibility(
+      absl::Span<const IntegerValue> sizes_x,
+      absl::Span<const IntegerValue> sizes_y,
+      std::pair<IntegerValue, IntegerValue> bounding_box_size);
+
+ private:
+  std::vector<int> index_by_decreasing_x_size_;
+  std::vector<int> index_by_decreasing_y_size_;
+
+  int64_t num_calls_ = 0;
+  int64_t num_conflicts_ = 0;
+  int64_t num_conflicts_two_items_ = 0;
+  int64_t num_trivial_conflicts_ = 0;
+
+  SharedStatistics* shared_stats_;
+};
+
+}  // namespace sat
+}  // namespace operations_research
+
+#endif  // OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_H_
diff --git a/ortools/sat/2d_orthogonal_packing_testing.cc b/ortools/sat/2d_orthogonal_packing_testing.cc
new file mode 100644
index 0000000000..35aceb5c86
--- /dev/null
+++ b/ortools/sat/2d_orthogonal_packing_testing.cc
@@ -0,0 +1,164 @@
+// Copyright 2010-2024 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/2d_orthogonal_packing_testing.h"
+
+#include <algorithm>
+#include <utility>
+#include <vector>
+
+#include "absl/log/check.h"
+#include "absl/random/bit_gen_ref.h"
+#include "absl/random/distributions.h"
+#include "ortools/sat/diffn_util.h"
+#include "ortools/sat/integer.h"
+
+namespace operations_research {
+namespace sat {
+
+std::vector<Rectangle> GenerateNonConflictingRectangles(
+    int num_rectangles, absl::BitGenRef random) {
+  static constexpr int kSizeMax = 1'000'000;
+  std::vector<Rectangle> rectangles;
+  rectangles.reserve(num_rectangles);
+  rectangles.push_back(
+      {.x_min = 0, .x_max = kSizeMax, .y_min = 0, .y_max = kSizeMax});
+  for (int i = 0; i < num_rectangles; ++i) {
+    std::swap(rectangles.back(),
+              rectangles[absl::Uniform(random, 0ull, rectangles.size() - 1)]);
+    const Rectangle& rec = rectangles.back();
+    if (absl::Bernoulli(random, 0.5)) {
+      const IntegerValue cut = IntegerValue(
+          absl::Uniform(random, rec.x_min.value(), rec.x_max.value()));
+      const Rectangle new_range = {.x_min = rec.x_min,
+                                   .x_max = cut,
+                                   .y_min = rec.y_min,
+                                   .y_max = rec.y_max};
+      const Rectangle new_range2 = {.x_min = cut,
+                                    .x_max = rec.x_max,
+                                    .y_min = rec.y_min,
+                                    .y_max = rec.y_max};
+      rectangles.pop_back();
+      rectangles.push_back(new_range);
+      rectangles.push_back(new_range2);
+    } else {
+      const IntegerValue cut = IntegerValue(
+          absl::Uniform(random, rec.y_min.value(), rec.y_max.value()));
+      const Rectangle new_range = {.x_min = rec.x_min,
+                                   .x_max = rec.x_max,
+                                   .y_min = rec.y_min,
+                                   .y_max = cut};
+      const Rectangle new_range2 = {.x_min = rec.x_min,
+                                    .x_max = rec.x_max,
+                                    .y_min = cut,
+                                    .y_max = rec.y_max};
+      rectangles.pop_back();
+      rectangles.push_back(new_range);
+      rectangles.push_back(new_range2);
+    }
+  }
+  return rectangles;
+}
+
+std::vector<RectangleInRange> MakeItemsFromRectangles(
+    const std::vector<Rectangle>& rectangles, double slack_factor,
+    absl::BitGenRef random) {
+  IntegerValue size_max_x = 0;
+  IntegerValue size_max_y = 0;
+  for (const Rectangle& rec : rectangles) {
+    size_max_x = std::max(size_max_x, rec.SizeX());
+    size_max_y = std::max(size_max_y, rec.SizeY());
+  }
+  std::vector<RectangleInRange> ranges;
+  ranges.reserve(rectangles.size());
+  const int max_slack_x = slack_factor * size_max_x.value();
+  const int max_slack_y = slack_factor * size_max_y.value();
+  for (const Rectangle& rec : rectangles) {
+    RectangleInRange range;
+    range.x_size = rec.x_max - rec.x_min;
+    range.y_size = rec.y_max - rec.y_min;
+    range.bounding_area = {
+        .x_min =
+            rec.x_min - IntegerValue(absl::Uniform(random, 0, max_slack_x)),
+        .x_max =
+            rec.x_max + IntegerValue(absl::Uniform(random, 0, max_slack_x)),
+        .y_min =
+            rec.y_min - IntegerValue(absl::Uniform(random, 0, max_slack_y)),
+        .y_max =
+            rec.y_max + IntegerValue(absl::Uniform(random, 0, max_slack_y))};
+    ranges.push_back(range);
+  }
+  return ranges;
+}
+
+std::vector<ItemForPairwiseRestriction>
+GenerateItemsRectanglesWithNoPairwiseConflict(
+    const std::vector<Rectangle>& rectangles, double slack_factor,
+    absl::BitGenRef random) {
+  const std::vector<RectangleInRange> range_items =
+      MakeItemsFromRectangles(rectangles, slack_factor, random);
+  std::vector<ItemForPairwiseRestriction> items;
+  items.reserve(rectangles.size());
+  for (int i = 0; i < range_items.size(); ++i) {
+    const RectangleInRange& rec = range_items[i];
+    items.push_back({.index = i,
+                     .x = {.start_min = rec.bounding_area.x_min,
+                           .start_max = rec.bounding_area.x_max - rec.x_size,
+                           .end_min = rec.bounding_area.x_min + rec.x_size,
+                           .end_max = rec.bounding_area.x_max},
+                     .y = {.start_min = rec.bounding_area.y_min,
+                           .start_max = rec.bounding_area.y_max - rec.y_size,
+                           .end_min = rec.bounding_area.y_min + rec.y_size,
+                           .end_max = rec.bounding_area.y_max}});
+  }
+  return items;
+}
+
+std::vector<ItemForPairwiseRestriction>
+GenerateItemsRectanglesWithNoPairwisePropagation(int num_rectangles,
+                                                 double slack_factor,
+                                                 absl::BitGenRef random) {
+  const std::vector<Rectangle> rectangles =
+      GenerateNonConflictingRectangles(num_rectangles, random);
+  std::vector<ItemForPairwiseRestriction> items =
+      GenerateItemsRectanglesWithNoPairwiseConflict(rectangles, slack_factor,
+                                                    random);
+  bool done = false;
+  // Now run the propagator until there is no more pairwise conditions.
+  do {
+    std::vector<PairwiseRestriction> results;
+    AppendPairwiseRestrictions(items, &results);
+    for (const PairwiseRestriction& restriction : results) {
+      CHECK(restriction.type !=
+            PairwiseRestriction::PairwiseRestrictionType::CONFLICT);
+      // Remove the slack we added
+      for (int index : {restriction.first_index, restriction.second_index}) {
+        const auto& rec = rectangles[index];
+        items[index] = {.index = index,
+                        .x = {.start_min = rec.x_min,
+                              .start_max = rec.x_min,
+                              .end_min = rec.x_max,
+                              .end_max = rec.x_max},
+                        .y = {.start_min = rec.y_min,
+                              .start_max = rec.y_min,
+                              .end_min = rec.y_max,
+                              .end_max = rec.y_max}};
+      }
+    }
+    done = results.empty();
+  } while (!done);
+  return items;
+}
+
+}  // namespace sat
+}  // namespace operations_research
diff --git a/ortools/sat/2d_orthogonal_packing_testing.h b/ortools/sat/2d_orthogonal_packing_testing.h
new file mode 100644
index 0000000000..2f12462456
--- /dev/null
+++ b/ortools/sat/2d_orthogonal_packing_testing.h
@@ -0,0 +1,45 @@
+// Copyright 2010-2024 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.
+
+#ifndef OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_TESTING_H_
+#define OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_TESTING_H_
+
+#include <vector>
+
+#include "absl/random/bit_gen_ref.h"
+#include "ortools/sat/diffn_util.h"
+
+namespace operations_research {
+namespace sat {
+
+std::vector<Rectangle> GenerateNonConflictingRectangles(int num_rectangles,
+                                                        absl::BitGenRef random);
+
+std::vector<RectangleInRange> MakeItemsFromRectangles(
+    const std::vector<Rectangle>& rectangles, double slack_factor,
+    absl::BitGenRef random);
+
+std::vector<ItemForPairwiseRestriction>
+GenerateItemsRectanglesWithNoPairwiseConflict(
+    const std::vector<Rectangle>& rectangles, double slack_factor,
+    absl::BitGenRef random);
+
+std::vector<ItemForPairwiseRestriction>
+GenerateItemsRectanglesWithNoPairwisePropagation(int num_rectangles,
+                                                 double slack_factor,
+                                                 absl::BitGenRef random);
+
+}  // namespace sat
+}  // namespace operations_research
+
+#endif  // OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_TESTING_H_
diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel
index 3fb6d0db1a..69bd5d857b 100644
--- a/ortools/sat/BUILD.bazel
+++ b/ortools/sat/BUILD.bazel
@@ -1014,6 +1014,7 @@ cc_library(
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/strings",
         "@com_google_absl//absl/time",
+        "@com_google_absl//absl/types:span",
     ],
 )
 
@@ -1373,6 +1374,7 @@ cc_library(
         "@com_google_absl//absl/container:flat_hash_set",
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/strings",
+        "@com_google_absl//absl/types:span",
     ],
 )
 
@@ -1695,6 +1697,7 @@ cc_library(
         "@com_google_absl//absl/random:bit_gen_ref",
         "@com_google_absl//absl/strings",
         "@com_google_absl//absl/strings:str_format",
+        "@com_google_absl//absl/types:span",
     ],
 )
 
@@ -1806,11 +1809,38 @@ cc_library(
     ],
 )
 
+cc_library(
+    name = "2d_orthogonal_packing",
+    srcs = ["2d_orthogonal_packing.cc"],
+    hdrs = ["2d_orthogonal_packing.h"],
+    deps = [
+        ":integer",
+        ":synchronization",
+        "@com_google_absl//absl/log:check",
+        "@com_google_absl//absl/types:span",
+    ],
+)
+
+cc_library(
+    name = "2d_orthogonal_packing_testing",
+    testonly = 1,
+    srcs = ["2d_orthogonal_packing_testing.cc"],
+    hdrs = ["2d_orthogonal_packing_testing.h"],
+    deps = [
+        ":diffn_util",
+        ":integer",
+        "@com_google_absl//absl/log:check",
+        "@com_google_absl//absl/random:bit_gen_ref",
+        "@com_google_absl//absl/random:distributions",
+    ],
+)
+
 cc_library(
     name = "diffn",
     srcs = ["diffn.cc"],
     hdrs = ["diffn.h"],
     deps = [
+        ":2d_orthogonal_packing",
         ":cumulative_energy",
         ":diffn_util",
         ":disjunctive",
@@ -1827,6 +1857,7 @@ cc_library(
         "//ortools/util:saturated_arithmetic",
         "//ortools/util:strong_integers",
         "@com_google_absl//absl/container:flat_hash_set",
+        "@com_google_absl//absl/container:inlined_vector",
         "@com_google_absl//absl/log",
         "@com_google_absl//absl/log:check",
         "@com_google_absl//absl/types:span",
@@ -1954,6 +1985,7 @@ cc_library(
     srcs = ["rins.cc"],
     hdrs = ["rins.h"],
     deps = [
+        ":cp_model_mapping",
         ":integer",
         ":linear_programming_constraint",
         ":model",
diff --git a/ortools/sat/cp_model_mapping.h b/ortools/sat/cp_model_mapping.h
index b203c2890e..85acd1ec1e 100644
--- a/ortools/sat/cp_model_mapping.h
+++ b/ortools/sat/cp_model_mapping.h
@@ -112,6 +112,7 @@ class CpModelMapping {
   template <typename List>
   std::vector<IntegerVariable> Integers(const List& list) const {
     std::vector<IntegerVariable> result;
+    result.reserve(list.size());
     for (const auto i : list) result.push_back(Integer(i));
     return result;
   }
@@ -119,6 +120,7 @@ class CpModelMapping {
   template <typename ProtoIndices>
   std::vector<sat::Literal> Literals(const ProtoIndices& indices) const {
     std::vector<sat::Literal> result;
+    result.reserve(indices.size());
     for (const int i : indices) result.push_back(CpModelMapping::Literal(i));
     return result;
   }
@@ -126,6 +128,7 @@ class CpModelMapping {
   template <typename List>
   std::vector<AffineExpression> Affines(const List& list) const {
     std::vector<AffineExpression> result;
+    result.reserve(list.size());
     for (const auto& i : list) result.push_back(Affine(i));
     return result;
   }
@@ -133,6 +136,7 @@ class CpModelMapping {
   template <typename ProtoIndices>
   std::vector<IntervalVariable> Intervals(const ProtoIndices& indices) const {
     std::vector<IntervalVariable> result;
+    result.reserve(indices.size());
     for (const int i : indices) result.push_back(Interval(i));
     return result;
   }
diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc
index 58b92088d2..4540249dc5 100644
--- a/ortools/sat/cp_model_search.cc
+++ b/ortools/sat/cp_model_search.cc
@@ -167,6 +167,7 @@ void AddDualSchedulingHeuristics(SatParameters& new_params) {
   new_params.set_max_pairs_pairwise_reasoning_in_no_overlap_2d(5000);
   new_params.set_use_timetabling_in_no_overlap_2d(true);
   new_params.set_use_energetic_reasoning_in_no_overlap_2d(true);
+  new_params.set_use_area_energetic_reasoning_in_no_overlap_2d(true);
 }
 
 // We want a random tie breaking among variables with equivalent values.
diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc
index c5443456cf..c95e4a1d61 100644
--- a/ortools/sat/cp_model_solver.cc
+++ b/ortools/sat/cp_model_solver.cc
@@ -938,7 +938,8 @@ IntegerVariable GetOrCreateVariableLinkedToSumOf(
 // Adds one LinearProgrammingConstraint per connected component of the model.
 IntegerVariable AddLPConstraints(bool objective_need_to_be_tight,
                                  const CpModelProto& model_proto, Model* m) {
-  const LinearRelaxation relaxation = ComputeLinearRelaxation(model_proto, m);
+  // Non const as we will std::move() stuff out of there.
+  LinearRelaxation relaxation = ComputeLinearRelaxation(model_proto, m);
 
   // The bipartite graph of LP constraints might be disconnected:
   // make a partition of the variables into connected components.
@@ -966,7 +967,8 @@ IntegerVariable AddLPConstraints(bool objective_need_to_be_tight,
            PositiveVariable(var).value();
   };
   for (int i = 0; i < num_lp_constraints; i++) {
-    for (const IntegerVariable var : relaxation.linear_constraints[i].vars) {
+    for (const IntegerVariable var :
+         relaxation.linear_constraints[i].VarsAsSpan()) {
       components.AddEdge(get_constraint_index(i), get_var_index(var));
     }
   }
@@ -1789,28 +1791,6 @@ void LoadCpModel(const CpModelProto& model_proto, Model* model) {
     }
   }
 
-  // Cache the links between model vars, IntegerVariables and lp constraints.
-  // TODO(user): Cache this only if it is actually used.
-  auto* integer_trail = model->GetOrCreate<IntegerTrail>();
-  auto* lp_dispatcher = model->GetOrCreate<LinearProgrammingDispatcher>();
-  auto* lp_vars = model->GetOrCreate<LPVariables>();
-  IntegerVariable size = integer_trail->NumIntegerVariables();
-  for (IntegerVariable positive_var(0); positive_var < size;
-       positive_var += 2) {
-    LPVariable lp_var;
-    lp_var.positive_var = positive_var;
-    lp_var.model_var =
-        mapping->GetProtoVariableFromIntegerVariable(positive_var);
-    const auto& it = lp_dispatcher->find(positive_var);
-    lp_var.lp = it != lp_dispatcher->end() ? it->second : nullptr;
-
-    if (lp_var.model_var >= 0) {
-      lp_vars->vars.push_back(lp_var);
-      lp_vars->model_vars_size =
-          std::max(lp_vars->model_vars_size, lp_var.model_var + 1);
-    }
-  }
-
   // Initialize the search strategies.
   auto* search_heuristics = model->GetOrCreate<SearchHeuristics>();
   search_heuristics->user_search =
diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc
index b9862fb65f..17d86eb7d7 100644
--- a/ortools/sat/cuts.cc
+++ b/ortools/sat/cuts.cc
@@ -33,6 +33,7 @@
 #include "absl/numeric/int128.h"
 #include "absl/strings/str_cat.h"
 #include "absl/strings/string_view.h"
+#include "absl/types/span.h"
 #include "ortools/base/logging.h"
 #include "ortools/base/stl_util.h"
 #include "ortools/base/strong_vector.h"
@@ -137,7 +138,7 @@ bool CutData::FillFromLinearConstraint(
     IntegerTrail* integer_trail) {
   rhs = absl::int128(base_ct.ub.value());
   terms.clear();
-  const int num_terms = base_ct.vars.size();
+  const int num_terms = base_ct.num_terms;
   for (int i = 0; i < num_terms; ++i) {
     const IntegerVariable var = base_ct.vars[i];
     if (!AppendOneTerm(var, base_ct.coeffs[i], lp_values[base_ct.vars[i]],
@@ -150,24 +151,24 @@ bool CutData::FillFromLinearConstraint(
 }
 
 bool CutData::FillFromParallelVectors(
-    const LinearConstraint& base_ct, const std::vector<double>& lp_values,
-    const std::vector<IntegerValue>& lower_bounds,
-    const std::vector<IntegerValue>& upper_bounds) {
-  rhs = absl::int128(base_ct.ub.value());
+    IntegerValue ub, absl::Span<const IntegerVariable> vars,
+    absl::Span<const IntegerValue> coeffs, absl::Span<const double> lp_values,
+    absl::Span<const IntegerValue> lower_bounds,
+    absl::Span<const IntegerValue> upper_bounds) {
+  rhs = absl::int128(ub.value());
   terms.clear();
 
   const int size = lp_values.size();
   if (size == 0) return true;
 
+  CHECK_EQ(vars.size(), size);
+  CHECK_EQ(coeffs.size(), size);
   CHECK_EQ(lower_bounds.size(), size);
   CHECK_EQ(upper_bounds.size(), size);
-  CHECK_EQ(base_ct.vars.size(), size);
-  CHECK_EQ(base_ct.coeffs.size(), size);
-  CHECK_EQ(base_ct.lb, kMinIntegerValue);
 
   for (int i = 0; i < size; ++i) {
-    if (!AppendOneTerm(base_ct.vars[i], base_ct.coeffs[i], lp_values[i],
-                       lower_bounds[i], upper_bounds[i])) {
+    if (!AppendOneTerm(vars[i], coeffs[i], lp_values[i], lower_bounds[i],
+                       upper_bounds[i])) {
       return false;
     }
   }
@@ -304,13 +305,15 @@ bool CutDataBuilder::ConvertToLinearConstraint(const CutData& cut,
 
   output->lb = kMinIntegerValue;
   output->ub = new_rhs;
-  output->vars.clear();
-  output->coeffs.clear();
+  output->resize(tmp_map_.size());
+  int new_size = 0;
   for (const auto [var, coeff] : tmp_map_) {
     if (coeff == 0) continue;
-    output->vars.push_back(var);
-    output->coeffs.push_back(coeff);
+    output->vars[new_size] = var;
+    output->coeffs[new_size] = coeff;
+    ++new_size;
   }
+  output->resize(new_size);
   DivideByGCD(output);
   return true;
 }
diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h
index 8503799f98..f642c8f09d 100644
--- a/ortools/sat/cuts.h
+++ b/ortools/sat/cuts.h
@@ -98,10 +98,12 @@ struct CutData {
       const absl::StrongVector<IntegerVariable, double>& lp_values,
       IntegerTrail* integer_trail);
 
-  bool FillFromParallelVectors(const LinearConstraint& base_ct,
-                               const std::vector<double>& lp_values,
-                               const std::vector<IntegerValue>& lower_bounds,
-                               const std::vector<IntegerValue>& upper_bounds);
+  bool FillFromParallelVectors(IntegerValue ub,
+                               absl::Span<const IntegerVariable> vars,
+                               absl::Span<const IntegerValue> coeffs,
+                               absl::Span<const double> lp_values,
+                               absl::Span<const IntegerValue> lower_bounds,
+                               absl::Span<const IntegerValue> upper_bounds);
 
   bool AppendOneTerm(IntegerVariable var, IntegerValue coeff, double lp_value,
                      IntegerValue lb, IntegerValue ub);
diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc
index 9f8710406d..eca299ab9e 100644
--- a/ortools/sat/diffn.cc
+++ b/ortools/sat/diffn.cc
@@ -18,6 +18,7 @@
 #include <algorithm>
 #include <cstddef>
 #include <cstdint>
+#include <iterator>
 #include <limits>
 #include <optional>
 #include <string>
@@ -25,9 +26,11 @@
 #include <vector>
 
 #include "absl/container/flat_hash_set.h"
+#include "absl/container/inlined_vector.h"
 #include "absl/log/check.h"
 #include "absl/types/span.h"
 #include "ortools/base/logging.h"
+#include "ortools/sat/2d_orthogonal_packing.h"
 #include "ortools/sat/cumulative_energy.h"
 #include "ortools/sat/diffn_util.h"
 #include "ortools/sat/disjunctive.h"
@@ -251,6 +254,9 @@ NonOverlappingRectanglesEnergyPropagator::
        num_conflicts_two_boxes_});
   stats.push_back({"NonOverlappingRectanglesEnergyPropagator/refined",
                    num_refined_conflicts_});
+  stats.push_back(
+      {"NonOverlappingRectanglesEnergyPropagator/orthogonal_packing_conflict",
+       num_orthogonal_packing_conflict_});
 
   shared_stats_->AddStats(stats);
 }
@@ -266,6 +272,8 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
                             .x_max = std::numeric_limits<IntegerValue>::min(),
                             .y_min = std::numeric_limits<IntegerValue>::max(),
                             .y_max = std::numeric_limits<IntegerValue>::min()};
+  IntegerValue max_x_item_size = IntegerValue(0);
+  IntegerValue max_y_item_size = IntegerValue(0);
   std::vector<RectangleInRange> active_box_ranges;
   active_box_ranges.reserve(num_boxes);
   for (int box = 0; box < num_boxes; ++box) {
@@ -285,6 +293,8 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
                           .y_max = y_.StartMax(box) + y_.SizeMin(box)},
         .x_size = x_.SizeMin(box),
         .y_size = y_.SizeMin(box)});
+    max_x_item_size = std::max(max_x_item_size, x_.SizeMin(box));
+    max_y_item_size = std::max(max_y_item_size, y_.SizeMin(box));
   }
 
   if (active_box_ranges.size() < 2) {
@@ -304,11 +314,13 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
     return true;
   }
 
-  std::optional<Conflict> best_conflict =
-      FindConflict(std::move(active_box_ranges));
+  std::optional<Conflict> best_conflict = FindConflict(
+      std::move(active_box_ranges), max_x_item_size, max_y_item_size);
   if (!best_conflict.has_value()) {
     return true;
   }
+  num_conflicts_++;
+
   // We found a conflict, so we can afford to run the propagator again to
   // search for a best explanation. This is specially the case since we only
   // want to re-run it over the items that participate in the conflict, so it is
@@ -316,7 +328,8 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
   IntegerValue best_explanation_size = best_conflict->items.size();
   bool refined = false;
   while (true) {
-    std::optional<Conflict> conflict = FindConflict(best_conflict->items);
+    const std::optional<Conflict> conflict =
+        FindConflict(best_conflict->items, max_x_item_size, max_y_item_size);
     if (!conflict.has_value()) break;
     // We prefer an explanation with the least number of boxes.
     if (conflict->items.size() >= best_explanation_size) {
@@ -329,30 +342,36 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
   }
 
   num_refined_conflicts_ += refined;
-  std::vector<RectangleInRange> generalized_explanation = GeneralizeExplanation(
-      best_conflict->rectangle_with_too_much_energy, best_conflict->items);
+  const std::vector<RectangleInRange> generalized_explanation =
+      GeneralizeExplanation(*best_conflict);
   if (best_explanation_size == 2) {
     num_conflicts_two_boxes_++;
   }
+  if (best_conflict->opp_problem_size > 0) {
+    num_orthogonal_packing_conflict_++;
+  }
   BuildAndReportEnergyTooLarge(generalized_explanation);
   return false;
 }
 
 std::optional<NonOverlappingRectanglesEnergyPropagator::Conflict>
 NonOverlappingRectanglesEnergyPropagator::FindConflict(
-    std::vector<RectangleInRange> active_box_ranges) {
-  const std::vector<Rectangle> rectangles_with_too_much_energy =
-      FindRectanglesWithEnergyConflictMC(active_box_ranges, *random_, 1.0);
+    std::vector<RectangleInRange> active_box_ranges,
+    const IntegerValue& max_x_item_size, const IntegerValue& max_y_item_size) {
+  const auto rectangles_with_too_much_energy =
+      FindRectanglesWithEnergyConflictMC(active_box_ranges, *random_, 1.0, 0.8);
 
-  if (rectangles_with_too_much_energy.empty()) return std::nullopt;
-
-  num_conflicts_++;
-  num_multiple_conflicts_ += rectangles_with_too_much_energy.size() > 1;
+  if (rectangles_with_too_much_energy.conflicts.empty() &&
+      rectangles_with_too_much_energy.candidates.empty()) {
+    return std::nullopt;
+  }
+  num_multiple_conflicts_ +=
+      rectangles_with_too_much_energy.conflicts.size() > 1;
 
   std::vector<RectangleInRange> best_explanation;
   Rectangle best_rectangle;
-  for (const auto& r : rectangles_with_too_much_energy) {
-    std::vector<RectangleInRange> range_for_explanation =
+  for (const Rectangle& r : rectangles_with_too_much_energy.conflicts) {
+    const std::vector<RectangleInRange> range_for_explanation =
         GetEnergyConflictForRectangle(r, active_box_ranges);
     CheckPropagationIsValid(range_for_explanation, r);
     if (best_explanation.empty() ||
@@ -361,26 +380,103 @@ NonOverlappingRectanglesEnergyPropagator::FindConflict(
       best_rectangle = r;
     }
   }
-  return Conflict{best_explanation, best_rectangle};
+
+  // Now try with a DFF on the candidates.
+  int opp_problem_size = 0;
+  // Also try the DFF on known conflicts, hopefully we will get a smaller
+  // explanation.
+  std::vector<Rectangle> filtered_rectangles;
+  filtered_rectangles.reserve(
+      rectangles_with_too_much_energy.conflicts.size() +
+      rectangles_with_too_much_energy.candidates.size());
+  filtered_rectangles = rectangles_with_too_much_energy.conflicts;
+  // Our current DFF does nothing if all elements are smaller than half of the
+  // size of the rectangle. So we filter all rectangles that are twice the
+  // dimensions of the largest item of our problem.
+  for (const Rectangle& r : rectangles_with_too_much_energy.candidates) {
+    if (r.SizeX() <= 2 * max_x_item_size || r.SizeY() <= 2 * max_y_item_size) {
+      filtered_rectangles.push_back(r);
+    }
+  }
+  // Sample 10 rectangles for looking for a conflict with DFF. This avoids
+  // making this heuristic too costly.
+  constexpr int kSampleSize = 10;
+  absl::InlinedVector<Rectangle, kSampleSize> sampled_rectangles;
+  std::sample(filtered_rectangles.begin(), filtered_rectangles.end(),
+              std::back_inserter(sampled_rectangles), kSampleSize, *random_);
+  std::vector<IntegerValue> sizes_x, sizes_y;
+  sizes_x.reserve(active_box_ranges.size());
+  sizes_y.reserve(active_box_ranges.size());
+  std::vector<RectangleInRange> filtered_items;
+  filtered_items.reserve(active_box_ranges.size());
+  for (const Rectangle& r : sampled_rectangles) {
+    sizes_x.clear();
+    sizes_y.clear();
+    filtered_items.clear();
+    for (int i = 0; i < active_box_ranges.size(); ++i) {
+      const RectangleInRange& box = active_box_ranges[i];
+      const Rectangle intersection = box.GetMinimumIntersection(r);
+      if (intersection.SizeX() > 0 && intersection.SizeY() > 0) {
+        sizes_x.push_back(intersection.SizeX());
+        sizes_y.push_back(intersection.SizeY());
+        filtered_items.push_back(box);
+      }
+    }
+    // This check the feasibility of a related orthogonal packing problem where
+    // our rectangle is the bounding box, and we need to fit inside it a set of
+    // items corresponding to the minimum intersection of the original items
+    // with this bounding box.
+    const auto opp_result = orthogonal_packing_checker_.TestFeasibility(
+        sizes_x, sizes_y, {r.SizeX(), r.SizeY()});
+    if (opp_result.result ==
+            OrthogonalPackingInfeasibilityDetector::Status::INFEASIBLE &&
+        (best_explanation.empty() ||
+         best_explanation.size() >
+             opp_result.minimum_unfeasible_subproblem.size())) {
+      best_explanation.clear();
+      for (int idx : opp_result.minimum_unfeasible_subproblem) {
+        best_explanation.push_back(filtered_items[idx]);
+      }
+      best_rectangle = r;
+      opp_problem_size = static_cast<int>(active_box_ranges.size());
+    }
+    // Use the fact that our rectangles are ordered in shrinking order to remove
+    // all items that will never contribute again.
+    filtered_items.swap(active_box_ranges);
+  }
+  if (!best_explanation.empty()) {
+    return Conflict{.items = std::move(best_explanation),
+                    .rectangle_with_too_much_energy = best_rectangle,
+                    .opp_problem_size = opp_problem_size};
+  }
+  return std::nullopt;
 }
 
 std::vector<RectangleInRange>
 NonOverlappingRectanglesEnergyPropagator::GeneralizeExplanation(
-    const Rectangle& rectangle,
-    const std::vector<RectangleInRange>& active_box_ranges) {
+    const Conflict& conflict) {
+  const Rectangle& rectangle = conflict.rectangle_with_too_much_energy;
   IntegerValue available_energy = rectangle.Area();
   IntegerValue used_energy = 0;
-  for (const RectangleInRange& box : active_box_ranges) {
+  for (const RectangleInRange& box : conflict.items) {
     IntegerValue energy = box.GetMinimumIntersectionArea(rectangle);
     used_energy += energy;
   }
-  std::vector<RectangleInRange> ranges_for_explanation(active_box_ranges);
+  std::vector<RectangleInRange> ranges_for_explanation(conflict.items);
 
   IntegerValue slack = used_energy - available_energy - 1;
-  VLOG_EVERY_N_SEC(2, 3) << "Rectangle with energy overflow: " << rectangle
-                         << " (" << used_energy << "/" << available_energy
-                         << "), with " << ranges_for_explanation.size()
-                         << " boxes inside.";
+
+  if (conflict.opp_problem_size > 0) {
+    VLOG_EVERY_N_SEC(2, 3)
+        << "Rectangle with energy overflow after solving OPP: " << rectangle
+        << " using " << ranges_for_explanation.size() << "/"
+        << conflict.opp_problem_size << " items.";
+  } else {
+    VLOG_EVERY_N_SEC(2, 3) << "Rectangle with energy overflow: " << rectangle
+                           << " (" << used_energy << "/" << available_energy
+                           << "), with " << ranges_for_explanation.size()
+                           << " boxes inside.";
+  }
 
   for (auto& range : ranges_for_explanation) {
     Rectangle overlap = range.GetMinimumIntersection(rectangle);
@@ -388,6 +484,7 @@ NonOverlappingRectanglesEnergyPropagator::GeneralizeExplanation(
     DCHECK_GT(overlap.SizeY(), 0);
     range = RectangleInRange::BiggestWithMinIntersection(
         rectangle, range, overlap.SizeX(), overlap.SizeY());
+
     // Now use the energy slack to make the intervals even bigger.
     if (overlap.SizeX() < slack) {
       IntegerValue y_slack = slack / overlap.SizeX();
@@ -434,8 +531,10 @@ NonOverlappingRectanglesEnergyPropagator::GeneralizeExplanation(
       }
     }
   }
-  CHECK_GT(used_energy, available_energy);
-  CHECK_GT(available_energy, 0);
+  if (conflict.opp_problem_size == 0) {
+    CHECK_GT(used_energy, available_energy);
+    CHECK_GT(available_energy, 0);
+  }
 
   return ranges_for_explanation;
 }
diff --git a/ortools/sat/diffn.h b/ortools/sat/diffn.h
index e60e86ee3d..23bfc35738 100644
--- a/ortools/sat/diffn.h
+++ b/ortools/sat/diffn.h
@@ -20,6 +20,7 @@
 
 #include "absl/container/flat_hash_set.h"
 #include "absl/types/span.h"
+#include "ortools/sat/2d_orthogonal_packing.h"
 #include "ortools/sat/diffn_util.h"
 #include "ortools/sat/disjunctive.h"
 #include "ortools/sat/integer.h"
@@ -42,7 +43,8 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
       : x_(*x),
         y_(*y),
         random_(model->GetOrCreate<ModelRandomGenerator>()),
-        shared_stats_(model->GetOrCreate<SharedStatistics>()) {}
+        shared_stats_(model->GetOrCreate<SharedStatistics>()),
+        orthogonal_packing_checker_(shared_stats_) {}
 
   ~NonOverlappingRectanglesEnergyPropagator() override;
 
@@ -53,13 +55,15 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
   struct Conflict {
     std::vector<RectangleInRange> items;
     Rectangle rectangle_with_too_much_energy;
+    // Size of the Orthogonal Packing Problem we used. Can be smaller than
+    // `items.size()` if we found a OPP conflict that didn't need all items.
+    int opp_problem_size;
   };
   std::optional<Conflict> FindConflict(
-      std::vector<RectangleInRange> active_box_ranges);
+      std::vector<RectangleInRange> active_box_ranges,
+      const IntegerValue& max_x_item_size, const IntegerValue& max_y_item_size);
 
-  std::vector<RectangleInRange> GeneralizeExplanation(
-      const Rectangle& rectangle,
-      const std::vector<RectangleInRange>& active_box_ranges);
+  std::vector<RectangleInRange> GeneralizeExplanation(const Conflict& conflict);
 
   bool BuildAndReportEnergyTooLarge(
       const std::vector<RectangleInRange>& ranges);
@@ -73,12 +77,14 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
   SchedulingConstraintHelper& y_;
   ModelRandomGenerator* random_;
   SharedStatistics* shared_stats_;
+  OrthogonalPackingInfeasibilityDetector orthogonal_packing_checker_;
 
   int64_t num_calls_ = 0;
   int64_t num_conflicts_ = 0;
   int64_t num_multiple_conflicts_ = 0;
   int64_t num_conflicts_two_boxes_ = 0;
   int64_t num_refined_conflicts_ = 0;
+  int64_t num_orthogonal_packing_conflict_ = 0;
 
   NonOverlappingRectanglesEnergyPropagator(
       const NonOverlappingRectanglesEnergyPropagator&) = delete;
diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc
index 52f0ddc586..9922ce3b96 100644
--- a/ortools/sat/diffn_util.cc
+++ b/ortools/sat/diffn_util.cc
@@ -21,6 +21,7 @@
 #include <cstddef>
 #include <cstdint>
 #include <limits>
+#include <optional>
 #include <ostream>
 #include <tuple>
 #include <utility>
@@ -1450,10 +1451,10 @@ std::size_t weighted_pick(absl::Span<const double> input,
 
 }  // namespace
 
-std::vector<Rectangle> FindRectanglesWithEnergyConflictMC(
+FindRectanglesResult FindRectanglesWithEnergyConflictMC(
     const std::vector<RectangleInRange>& intervals, absl::BitGenRef random,
-    double temperature) {
-  std::vector<Rectangle> result;
+    double temperature, double candidate_energy_usage_factor) {
+  FindRectanglesResult result;
   ProbingRectangle ranges(intervals);
 
   static const std::vector<double>* cached_probabilities =
@@ -1467,7 +1468,10 @@ std::vector<Rectangle> FindRectanglesWithEnergyConflictMC(
     const IntegerValue rect_area = ranges.GetCurrentRectangleArea();
     const IntegerValue min_energy = ranges.GetMinimumEnergy();
     if (min_energy > rect_area) {
-      result.push_back(ranges.GetCurrentRectangle());
+      result.conflicts.push_back(ranges.GetCurrentRectangle());
+    } else if (min_energy.value() >
+               candidate_energy_usage_factor * rect_area.value()) {
+      result.candidates.push_back(ranges.GetCurrentRectangle());
     }
     if (min_energy == 0) {
       break;
@@ -1502,7 +1506,7 @@ std::vector<Rectangle> FindRectanglesWithEnergyConflictMC(
     ranges.Shrink(candidates[weighted_pick(weights, random)]);
   }
   if (ranges.GetMinimumEnergy() > ranges.GetCurrentRectangleArea()) {
-    result.push_back(ranges.GetCurrentRectangle());
+    result.conflicts.push_back(ranges.GetCurrentRectangle());
   }
   return result;
 }
diff --git a/ortools/sat/diffn_util.h b/ortools/sat/diffn_util.h
index 6e73b71588..227163f761 100644
--- a/ortools/sat/diffn_util.h
+++ b/ortools/sat/diffn_util.h
@@ -572,11 +572,22 @@ class ProbingRectangle {
 // - shrink the rectangle by an edge to the next "interesting" value. Choose
 //   the edge randomly, but biased towards the change that increases the ratio
 //   area_inside / area_rectangle;
-// - collect a result at every conflict;
+// - collect a result at every conflict or every time the ratio
+//   used_energy/available_energy is more than `candidate_energy_usage_factor`;
 // - stop when the rectangle is empty.
-std::vector<Rectangle> FindRectanglesWithEnergyConflictMC(
+struct FindRectanglesResult {
+  // Known conflicts: the minimal energy used inside the rectangle is larger
+  // than the area of the rectangle.
+  std::vector<Rectangle> conflicts;
+  // Rectangles without a conflict but having used_energy/available_energy >
+  // candidate_energy_usage_factor. Those are good candidates for finding
+  // conflicts using more sophisticated heuristics. Those rectangles are
+  // ordered so the n-th rectangle is always fully inside the n-1-th one.
+  std::vector<Rectangle> candidates;
+};
+FindRectanglesResult FindRectanglesWithEnergyConflictMC(
     const std::vector<RectangleInRange>& intervals, absl::BitGenRef random,
-    double temperature);
+    double temperature, double candidate_energy_usage_factor);
 
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/feasibility_pump.cc b/ortools/sat/feasibility_pump.cc
index 1ebfbb1210..e157703049 100644
--- a/ortools/sat/feasibility_pump.cc
+++ b/ortools/sat/feasibility_pump.cc
@@ -84,7 +84,7 @@ FeasibilityPump::~FeasibilityPump() {
 
 void FeasibilityPump::AddLinearConstraint(const LinearConstraint& ct) {
   // We still create the mirror variable right away though.
-  for (const IntegerVariable var : ct.vars) {
+  for (const IntegerVariable var : ct.VarsAsSpan()) {
     GetOrCreateMirrorVariable(PositiveVariable(var));
   }
 
@@ -92,9 +92,8 @@ void FeasibilityPump::AddLinearConstraint(const LinearConstraint& ct) {
   LinearConstraintInternal& new_ct = integer_lp_.back();
   new_ct.lb = ct.lb;
   new_ct.ub = ct.ub;
-  const int size = ct.vars.size();
   CHECK_LE(ct.lb, ct.ub);
-  for (int i = 0; i < size; ++i) {
+  for (int i = 0; i < ct.num_terms; ++i) {
     // We only use positive variable inside this class.
     IntegerVariable var = ct.vars[i];
     IntegerValue coeff = ct.coeffs[i];
diff --git a/ortools/sat/implied_bounds.cc b/ortools/sat/implied_bounds.cc
index 1dcefd4b08..815fe6358d 100644
--- a/ortools/sat/implied_bounds.cc
+++ b/ortools/sat/implied_bounds.cc
@@ -495,7 +495,9 @@ bool ProductDecomposer::TryToLinearize(const AffineExpression& left,
 }
 
 ProductDetector::ProductDetector(Model* model)
-    : enabled_(model->GetOrCreate<SatParameters>()->linearization_level() > 1),
+    : enabled_(
+          model->GetOrCreate<SatParameters>()->detect_linearized_product() &&
+          model->GetOrCreate<SatParameters>()->linearization_level() > 1),
       sat_solver_(model->GetOrCreate<SatSolver>()),
       trail_(model->GetOrCreate<Trail>()),
       integer_trail_(model->GetOrCreate<IntegerTrail>()),
diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc
index 073bf0616a..607a7a6ee0 100644
--- a/ortools/sat/integer_expr.cc
+++ b/ortools/sat/integer_expr.cc
@@ -44,9 +44,9 @@ namespace sat {
 
 template <bool use_int128>
 LinearConstraintPropagator<use_int128>::LinearConstraintPropagator(
-    const std::vector<Literal>& enforcement_literals,
-    const std::vector<IntegerVariable>& vars,
-    const std::vector<IntegerValue>& coeffs, IntegerValue upper, Model* model)
+    absl::Span<const Literal> enforcement_literals,
+    absl::Span<const IntegerVariable> vars,
+    absl::Span<const IntegerValue> coeffs, IntegerValue upper, Model* model)
     : upper_bound_(upper),
       shared_(
           model->GetOrCreate<LinearConstraintPropagator<use_int128>::Shared>()),
@@ -82,6 +82,33 @@ LinearConstraintPropagator<use_int128>::LinearConstraintPropagator(
   rev_lb_fixed_vars_ = IntegerValue(0);
 }
 
+// TODO(user): Avoid duplication with other constructor.
+template <bool use_int128>
+LinearConstraintPropagator<use_int128>::LinearConstraintPropagator(
+    LinearConstraint ct, Model* model)
+    : upper_bound_(ct.ub),
+      shared_(
+          model->GetOrCreate<LinearConstraintPropagator<use_int128>::Shared>()),
+      size_(ct.num_terms),
+      vars_(std::move(ct.vars)),
+      coeffs_(std::move(ct.coeffs)),
+      max_variations_(new IntegerValue[size_]) {
+  // TODO(user): deal with this corner case.
+  CHECK_GT(size_, 0);
+
+  // Handle negative coefficients.
+  for (int i = 0; i < size_; ++i) {
+    if (coeffs_[i] < 0) {
+      vars_[i] = NegationOf(vars_[i]);
+      coeffs_[i] = -coeffs_[i];
+    }
+  }
+
+  // Initialize the reversible numbers.
+  rev_num_fixed_vars_ = 0;
+  rev_lb_fixed_vars_ = IntegerValue(0);
+}
+
 template <bool use_int128>
 void LinearConstraintPropagator<use_int128>::FillIntegerReason() {
   shared_->integer_reason.clear();
diff --git a/ortools/sat/integer_expr.h b/ortools/sat/integer_expr.h
index 65ce7ec3ca..764b543ecb 100644
--- a/ortools/sat/integer_expr.h
+++ b/ortools/sat/integer_expr.h
@@ -70,11 +70,16 @@ class LinearConstraintPropagator : public PropagatorInterface {
   // otherwise we enforce the implication refied_literal => constraint is true.
   // Note that we don't do the reverse implication here, it is usually done by
   // another LinearConstraintPropagator constraint on the negated variables.
-  LinearConstraintPropagator(const std::vector<Literal>& enforcement_literals,
-                             const std::vector<IntegerVariable>& vars,
-                             const std::vector<IntegerValue>& coeffs,
+  LinearConstraintPropagator(absl::Span<const Literal> enforcement_literals,
+                             absl::Span<const IntegerVariable> vars,
+                             absl::Span<const IntegerValue> coeffs,
                              IntegerValue upper_bound, Model* model);
 
+  // This version allow to std::move the memory from the LinearConstraint
+  // directly. It Only uses the upper bound. Id does not support
+  // enforcement_literals.
+  LinearConstraintPropagator(LinearConstraint ct, Model* model);
+
   // We propagate:
   // - If the sum of the individual lower-bound is > upper_bound, we fail.
   // - For all i, upper-bound of i
@@ -626,23 +631,26 @@ inline std::function<void(Model*)> WeightedSumGreaterOrEqualReif(
 
 // LinearConstraint version.
 inline void LoadLinearConstraint(const LinearConstraint& cst, Model* model) {
-  if (cst.vars.empty()) {
+  if (cst.num_terms == 0) {
     if (cst.lb <= 0 && cst.ub >= 0) return;
     model->GetOrCreate<SatSolver>()->NotifyThatModelIsUnsat();
     return;
   }
 
   // TODO(user): Remove the conversion!
-  std::vector<int64_t> converted_coeffs;
+  std::vector<IntegerVariable> vars(cst.num_terms);
+  std::vector<int64_t> converted_coeffs(cst.num_terms);
+  for (int i = 0; i < cst.num_terms; ++i) {
+    vars[i] = cst.vars[i];
+    converted_coeffs[i] = cst.coeffs[i].value();
+  }
 
-  for (const IntegerValue v : cst.coeffs) converted_coeffs.push_back(v.value());
   if (cst.ub < kMaxIntegerValue) {
-    model->Add(
-        WeightedSumLowerOrEqual(cst.vars, converted_coeffs, cst.ub.value()));
+    model->Add(WeightedSumLowerOrEqual(vars, converted_coeffs, cst.ub.value()));
   }
   if (cst.lb > kMinIntegerValue) {
     model->Add(
-        WeightedSumGreaterOrEqual(cst.vars, converted_coeffs, cst.lb.value()));
+        WeightedSumGreaterOrEqual(vars, converted_coeffs, cst.lb.value()));
   }
 }
 
@@ -652,7 +660,7 @@ inline void LoadConditionalLinearConstraint(
   if (enforcement_literals.empty()) {
     return LoadLinearConstraint(cst, model);
   }
-  if (cst.vars.empty()) {
+  if (cst.num_terms == 0) {
     if (cst.lb <= 0 && cst.ub >= 0) return;
 
     // The enforcement literals cannot be all at true.
@@ -666,16 +674,22 @@ inline void LoadConditionalLinearConstraint(
   // TODO(user): Remove the conversion!
   std::vector<Literal> converted_literals(enforcement_literals.begin(),
                                           enforcement_literals.end());
-  std::vector<int64_t> converted_coeffs;
-  for (const IntegerValue v : cst.coeffs) converted_coeffs.push_back(v.value());
+
+  // TODO(user): Remove the conversion!
+  std::vector<IntegerVariable> vars(cst.num_terms);
+  std::vector<int64_t> converted_coeffs(cst.num_terms);
+  for (int i = 0; i < cst.num_terms; ++i) {
+    vars[i] = cst.vars[i];
+    converted_coeffs[i] = cst.coeffs[i].value();
+  }
 
   if (cst.ub < kMaxIntegerValue) {
     model->Add(ConditionalWeightedSumLowerOrEqual(
-        converted_literals, cst.vars, converted_coeffs, cst.ub.value()));
+        converted_literals, vars, converted_coeffs, cst.ub.value()));
   }
   if (cst.lb > kMinIntegerValue) {
     model->Add(ConditionalWeightedSumGreaterOrEqual(
-        converted_literals, cst.vars, converted_coeffs, cst.lb.value()));
+        converted_literals, vars, converted_coeffs, cst.lb.value()));
   }
 }
 
diff --git a/ortools/sat/intervals.cc b/ortools/sat/intervals.cc
index 6fd47d0c08..f582973747 100644
--- a/ortools/sat/intervals.cc
+++ b/ortools/sat/intervals.cc
@@ -840,7 +840,7 @@ IntegerValue SchedulingConstraintHelper::GetMinOverlap(int t,
 IntegerValue ComputeEnergyMinInWindow(
     IntegerValue start_min, IntegerValue start_max, IntegerValue end_min,
     IntegerValue end_max, IntegerValue size_min, IntegerValue demand_min,
-    const std::vector<LiteralValueValue>& filtered_energy,
+    absl::Span<const LiteralValueValue> filtered_energy,
     IntegerValue window_start, IntegerValue window_end) {
   if (window_end <= window_start) return IntegerValue(0);
 
@@ -1075,7 +1075,7 @@ bool SchedulingDemandHelper::AddLinearizedDemand(
 }
 
 void SchedulingDemandHelper::OverrideLinearizedEnergies(
-    const std::vector<LinearExpression>& energies) {
+    absl::Span<const LinearExpression> energies) {
   const int num_tasks = energies.size();
   DCHECK_EQ(num_tasks, helper_->NumTasks());
   linearized_energies_.resize(num_tasks);
diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h
index 1f75b3e0be..cb12e8f333 100644
--- a/ortools/sat/intervals.h
+++ b/ortools/sat/intervals.h
@@ -684,8 +684,7 @@ class SchedulingDemandHelper {
   }
 
   // Visible for testing.
-  void OverrideLinearizedEnergies(
-      const std::vector<LinearExpression>& energies);
+  void OverrideLinearizedEnergies(absl::Span<const LinearExpression> energies);
   void OverrideDecomposedEnergies(
       const std::vector<std::vector<LiteralValueValue>>& energies);
   // Returns the decomposed energy terms compatible with the current literal
@@ -735,7 +734,7 @@ class SchedulingDemandHelper {
 IntegerValue ComputeEnergyMinInWindow(
     IntegerValue start_min, IntegerValue start_max, IntegerValue end_min,
     IntegerValue end_max, IntegerValue size_min, IntegerValue demand_min,
-    const std::vector<LiteralValueValue>& filtered_energy,
+    absl::Span<const LiteralValueValue> filtered_energy,
     IntegerValue window_start, IntegerValue window_end);
 
 // =============================================================================
@@ -978,13 +977,19 @@ inline std::function<IntervalVariable(Model*)> NewIntervalWithVariableSize(
 inline std::function<IntervalVariable(Model*)> NewOptionalInterval(
     int64_t min_start, int64_t max_end, int64_t size, Literal is_present) {
   return [=](Model* model) {
+    CHECK_LE(min_start + size, max_end);
+    const IntegerVariable start =
+        model->Add(NewIntegerVariable(min_start, max_end - size));
     return model->GetOrCreate<IntervalsRepository>()->CreateInterval(
-        model->Add(NewIntegerVariable(min_start, max_end)),
-        model->Add(NewIntegerVariable(min_start, max_end)), kNoIntegerVariable,
-        IntegerValue(size), is_present.Index());
+        AffineExpression(start),
+        AffineExpression(start, IntegerValue(1), IntegerValue(size)),
+        AffineExpression(IntegerValue(size)), is_present.Index(),
+        /*add_linear_relation=*/true);
   };
 }
 
+// TODO(user): Optional variables can be broken with sat_inprocessing, use with
+// care.
 inline std::function<IntervalVariable(Model*)>
 NewOptionalIntervalWithOptionalVariables(int64_t min_start, int64_t max_end,
                                          int64_t size, Literal is_present) {
diff --git a/ortools/sat/lb_tree_search.cc b/ortools/sat/lb_tree_search.cc
index a73fa9af55..13b7c4e503 100644
--- a/ortools/sat/lb_tree_search.cc
+++ b/ortools/sat/lb_tree_search.cc
@@ -25,6 +25,7 @@
 #include "absl/strings/str_cat.h"
 #include "absl/time/clock.h"
 #include "absl/time/time.h"
+#include "absl/types/span.h"
 #include "ortools/base/logging.h"
 #include "ortools/base/strong_vector.h"
 #include "ortools/sat/cp_model_mapping.h"
@@ -622,7 +623,7 @@ SatSolver::Status LbTreeSearch::Search(
 }
 
 std::vector<Literal> LbTreeSearch::ExtractDecisions(
-    int base_level, const std::vector<Literal>& conflict) {
+    int base_level, absl::Span<const Literal> conflict) {
   std::vector<int> num_per_level(sat_solver_->CurrentDecisionLevel() + 1, 0);
   std::vector<bool> is_marked;
   for (const Literal l : conflict) {
diff --git a/ortools/sat/lb_tree_search.h b/ortools/sat/lb_tree_search.h
index 3bdcd9b21f..8a72301aeb 100644
--- a/ortools/sat/lb_tree_search.h
+++ b/ortools/sat/lb_tree_search.h
@@ -24,6 +24,7 @@
 
 #include "absl/strings/string_view.h"
 #include "absl/time/time.h"
+#include "absl/types/span.h"
 #include "ortools/base/strong_vector.h"
 #include "ortools/sat/integer.h"
 #include "ortools/sat/integer_search.h"
@@ -131,7 +132,7 @@ class LbTreeSearch {
   // Returns a small number of decision needed to reach the same conflict.
   // We basically reduce the number of decision at each level to 1.
   std::vector<Literal> ExtractDecisions(int base_level,
-                                        const std::vector<Literal>& conflict);
+                                        absl::Span<const Literal> conflict);
 
   // Used in the solve logs.
   std::string SmallProgressString() const;
diff --git a/ortools/sat/linear_constraint.cc b/ortools/sat/linear_constraint.cc
index ed181f2b6d..589a810e19 100644
--- a/ortools/sat/linear_constraint.cc
+++ b/ortools/sat/linear_constraint.cc
@@ -15,6 +15,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <cstddef>
 #include <cstdint>
 #include <limits>
 #include <string>
@@ -25,6 +26,7 @@
 #include "absl/container/flat_hash_set.h"
 #include "absl/log/check.h"
 #include "absl/strings/str_cat.h"
+#include "absl/types/span.h"
 #include "ortools/base/mathutil.h"
 #include "ortools/base/strong_vector.h"
 #include "ortools/sat/integer.h"
@@ -80,7 +82,7 @@ void LinearConstraintBuilder::AddLinearExpression(const LinearExpression& expr,
 }
 
 ABSL_MUST_USE_RESULT bool LinearConstraintBuilder::AddDecomposedProduct(
-    const std::vector<LiteralValueValue>& product) {
+    absl::Span<const LiteralValueValue> product) {
   if (product.empty()) return true;
 
   IntegerValue product_min = kMaxIntegerValue;
@@ -164,7 +166,7 @@ double ComputeActivity(
     const LinearConstraint& constraint,
     const absl::StrongVector<IntegerVariable, double>& values) {
   int i = 0;
-  const int size = static_cast<int>(constraint.vars.size());
+  const int size = constraint.num_terms;
   const int shifted_size = size - 3;
   double a0 = 0.0;
   double a1 = 0.0;
@@ -196,37 +198,35 @@ double ComputeActivity(
   return activity;
 }
 
-double ComputeL2Norm(const LinearConstraint& constraint) {
+double ComputeL2Norm(const LinearConstraint& ct) {
   double sum = 0.0;
-  for (const IntegerValue coeff : constraint.coeffs) {
-    sum += ToDouble(coeff) * ToDouble(coeff);
+  for (int i = 0; i < ct.num_terms; ++i) {
+    sum += ToDouble(ct.coeffs[i]) * ToDouble(ct.coeffs[i]);
   }
   return std::sqrt(sum);
 }
 
-IntegerValue ComputeInfinityNorm(const LinearConstraint& constraint) {
+IntegerValue ComputeInfinityNorm(const LinearConstraint& ct) {
   IntegerValue result(0);
-  for (const IntegerValue coeff : constraint.coeffs) {
-    result = std::max(result, IntTypeAbs(coeff));
+  for (int i = 0; i < ct.num_terms; ++i) {
+    result = std::max(result, IntTypeAbs(ct.coeffs[i]));
   }
   return result;
 }
 
-double ScalarProduct(const LinearConstraint& constraint1,
-                     const LinearConstraint& constraint2) {
-  DCHECK(std::is_sorted(constraint1.vars.begin(), constraint1.vars.end()));
-  DCHECK(std::is_sorted(constraint2.vars.begin(), constraint2.vars.end()));
+double ScalarProduct(const LinearConstraint& ct1, const LinearConstraint& ct2) {
+  DCHECK(std::is_sorted(ct1.vars.get(), ct1.vars.get() + ct1.num_terms));
+  DCHECK(std::is_sorted(ct2.vars.get(), ct2.vars.get() + ct2.num_terms));
   double scalar_product = 0.0;
   int index_1 = 0;
   int index_2 = 0;
-  while (index_1 < constraint1.vars.size() &&
-         index_2 < constraint2.vars.size()) {
-    if (constraint1.vars[index_1] == constraint2.vars[index_2]) {
-      scalar_product += ToDouble(constraint1.coeffs[index_1]) *
-                        ToDouble(constraint2.coeffs[index_2]);
+  while (index_1 < ct1.num_terms && index_2 < ct2.num_terms) {
+    if (ct1.vars[index_1] == ct2.vars[index_2]) {
+      scalar_product +=
+          ToDouble(ct1.coeffs[index_1]) * ToDouble(ct2.coeffs[index_2]);
       index_1++;
       index_2++;
-    } else if (constraint1.vars[index_1] > constraint2.vars[index_2]) {
+    } else if (ct1.vars[index_1] > ct2.vars[index_2]) {
       index_2++;
     } else {
       index_1++;
@@ -238,7 +238,7 @@ double ScalarProduct(const LinearConstraint& constraint1,
 namespace {
 
 // TODO(user): Template for any integer type and expose this?
-IntegerValue ComputeGcd(const std::vector<IntegerValue>& values) {
+IntegerValue ComputeGcd(absl::Span<const IntegerValue> values) {
   if (values.empty()) return IntegerValue(1);
   int64_t gcd = 0;
   for (const IntegerValue value : values) {
@@ -252,8 +252,9 @@ IntegerValue ComputeGcd(const std::vector<IntegerValue>& values) {
 }  // namespace
 
 void DivideByGCD(LinearConstraint* constraint) {
-  if (constraint->coeffs.empty()) return;
-  const IntegerValue gcd = ComputeGcd(constraint->coeffs);
+  if (constraint->num_terms == 0) return;
+  const IntegerValue gcd = ComputeGcd(
+      {constraint->coeffs.get(), static_cast<size_t>(constraint->num_terms)});
   if (gcd == 1) return;
 
   if (constraint->lb > kMinIntegerValue) {
@@ -262,24 +263,25 @@ void DivideByGCD(LinearConstraint* constraint) {
   if (constraint->ub < kMaxIntegerValue) {
     constraint->ub = FloorRatio(constraint->ub, gcd);
   }
-  for (IntegerValue& coeff : constraint->coeffs) coeff /= gcd;
+  for (int i = 0; i < constraint->num_terms; ++i) {
+    constraint->coeffs[i] /= gcd;
+  }
 }
 
 void RemoveZeroTerms(LinearConstraint* constraint) {
   int new_size = 0;
-  const int size = constraint->vars.size();
+  const int size = constraint->num_terms;
   for (int i = 0; i < size; ++i) {
     if (constraint->coeffs[i] == 0) continue;
     constraint->vars[new_size] = constraint->vars[i];
     constraint->coeffs[new_size] = constraint->coeffs[i];
     ++new_size;
   }
-  constraint->vars.resize(new_size);
-  constraint->coeffs.resize(new_size);
+  constraint->resize(new_size);
 }
 
 void MakeAllCoefficientsPositive(LinearConstraint* constraint) {
-  const int size = constraint->vars.size();
+  const int size = constraint->num_terms;
   for (int i = 0; i < size; ++i) {
     const IntegerValue coeff = constraint->coeffs[i];
     if (coeff < 0) {
@@ -290,7 +292,7 @@ void MakeAllCoefficientsPositive(LinearConstraint* constraint) {
 }
 
 void MakeAllVariablesPositive(LinearConstraint* constraint) {
-  const int size = constraint->vars.size();
+  const int size = constraint->num_terms;
   for (int i = 0; i < size; ++i) {
     const IntegerVariable var = constraint->vars[i];
     if (!VariableIsPositive(var)) {
@@ -355,16 +357,13 @@ std::string LinearExpression::DebugString() const {
   return result;
 }
 
-// TODO(user): it would be better if LinearConstraint natively supported
-// term and not two separated vectors. Fix?
-//
 // TODO(user): This is really similar to CleanTermsAndFillConstraint(), maybe
 // we should just make the later switch negative variable to positive ones to
 // avoid an extra linear scan on each new cuts.
 void CanonicalizeConstraint(LinearConstraint* ct) {
   std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
 
-  const int size = ct->vars.size();
+  const int size = ct->num_terms;
   for (int i = 0; i < size; ++i) {
     if (VariableIsPositive(ct->vars[i])) {
       terms.push_back({ct->vars[i], ct->coeffs[i]});
@@ -374,17 +373,18 @@ void CanonicalizeConstraint(LinearConstraint* ct) {
   }
   std::sort(terms.begin(), terms.end());
 
-  ct->vars.clear();
-  ct->coeffs.clear();
+  int new_size = 0;
   for (const auto& term : terms) {
-    ct->vars.push_back(term.first);
-    ct->coeffs.push_back(term.second);
+    ct->vars[new_size] = term.first;
+    ct->coeffs[new_size] = term.second;
+    ++new_size;
   }
+  ct->resize(new_size);
 }
 
 bool NoDuplicateVariable(const LinearConstraint& ct) {
   absl::flat_hash_set<IntegerVariable> seen_variables;
-  const int size = ct.vars.size();
+  const int size = ct.num_terms;
   for (int i = 0; i < size; ++i) {
     if (VariableIsPositive(ct.vars[i])) {
       if (!seen_variables.insert(ct.vars[i]).second) return false;
@@ -416,7 +416,7 @@ bool ValidateLinearConstraintForOverflow(const LinearConstraint& constraint,
                                          const IntegerTrail& integer_trail) {
   int64_t positive_sum(0);
   int64_t negative_sum(0);
-  for (int i = 0; i < constraint.vars.size(); ++i) {
+  for (int i = 0; i < constraint.num_terms; ++i) {
     const IntegerVariable var = constraint.vars[i];
     const IntegerValue coeff = constraint.coeffs[i];
     const IntegerValue lb = integer_trail.LevelZeroLowerBound(var);
@@ -488,7 +488,7 @@ bool PossibleOverflow(const IntegerTrail& integer_trail,
                       const LinearConstraint& constraint) {
   IntegerValue min_activity(0);
   IntegerValue max_activity(0);
-  const int size = constraint.vars.size();
+  const int size = constraint.num_terms;
   for (int i = 0; i < size; ++i) {
     const IntegerVariable var = constraint.vars[i];
     const IntegerValue coeff = constraint.coeffs[i];
diff --git a/ortools/sat/linear_constraint.h b/ortools/sat/linear_constraint.h
index d483a5e836..bcc671ea01 100644
--- a/ortools/sat/linear_constraint.h
+++ b/ortools/sat/linear_constraint.h
@@ -22,6 +22,7 @@
 
 #include "absl/base/attributes.h"
 #include "absl/strings/str_cat.h"
+#include "absl/types/span.h"
 #include "ortools/base/strong_vector.h"
 #include "ortools/sat/integer.h"
 #include "ortools/sat/model.h"
@@ -42,20 +43,45 @@ struct LinearConstraint {
   IntegerValue lb;
   IntegerValue ub;
 
-  // TODO(user): This class is almost always static, replace by a size
-  // and two [] to save 24 bytes per constraints.
-  std::vector<IntegerVariable> vars;
-  std::vector<IntegerValue> coeffs;
+  // Rather than using two std::vector<> this class is optimized for memory
+  // consumption, given that most of our LinearConstraint are constructed once
+  // and for all.
+  //
+  // It is however up to clients to maintain the invariants that both vars
+  // and coeffs are properly allocated and of size num_terms.
+  //
+  // Also note that we did not add a copy constructor, to make sure that this is
+  // moved as often as possible. This allowed to optimize a few call site and so
+  // far we never copy this.
+  int num_terms = 0;
+  std::unique_ptr<IntegerVariable[]> vars;
+  std::unique_ptr<IntegerValue[]> coeffs;
 
   LinearConstraint() = default;
   LinearConstraint(IntegerValue _lb, IntegerValue _ub) : lb(_lb), ub(_ub) {}
 
+  // Resize the LinearConstraint to have space for num_terms. We always
+  // re-allocate if the size is different to always be tight in memory.
+  void resize(int size) {
+    if (size == num_terms) return;
+    IntegerVariable* tmp_vars = new IntegerVariable[size];
+    IntegerValue* tmp_coeffs = new IntegerValue[size];
+    const int to_copy = std::min(size, num_terms);
+    if (to_copy > 0) {
+      memcpy(tmp_vars, vars.get(), sizeof(IntegerVariable) * to_copy);
+      memcpy(tmp_coeffs, coeffs.get(), sizeof(IntegerValue) * to_copy);
+    }
+    num_terms = size;
+    vars.reset(tmp_vars);
+    coeffs.reset(tmp_coeffs);
+  }
+
   std::string DebugString() const {
     std::string result;
     if (lb.value() > kMinIntegerValue) {
       absl::StrAppend(&result, lb.value(), " <= ");
     }
-    for (int i = 0; i < vars.size(); ++i) {
+    for (int i = 0; i < num_terms; ++i) {
       absl::StrAppend(&result, i > 0 ? " " : "",
                       IntegerTermDebugString(vars[i], coeffs[i]));
     }
@@ -65,12 +91,32 @@ struct LinearConstraint {
     return result;
   }
 
-  bool operator==(const LinearConstraint other) const {
+  bool IsEqualIgnoringBounds(const LinearConstraint& other) const {
+    if (this->num_terms != other.num_terms) return false;
+    if (this->num_terms == 0) return true;
+    if (memcmp(this->vars.get(), other.vars.get(),
+               sizeof(IntegerVariable) * this->num_terms)) {
+      return false;
+    }
+    if (memcmp(this->coeffs.get(), other.coeffs.get(),
+               sizeof(IntegerValue) * this->num_terms)) {
+      return false;
+    }
+    return true;
+  }
+
+  bool operator==(const LinearConstraint& other) const {
     if (this->lb != other.lb) return false;
     if (this->ub != other.ub) return false;
-    if (this->vars != other.vars) return false;
-    if (this->coeffs != other.coeffs) return false;
-    return true;
+    return IsEqualIgnoringBounds(other);
+  }
+
+  absl::Span<const IntegerVariable> VarsAsSpan() const {
+    return absl::MakeSpan(vars.get(), num_terms);
+  }
+
+  absl::Span<const IntegerValue> CoeffsAsSpan() const {
+    return absl::MakeSpan(coeffs.get(), num_terms);
   }
 };
 
@@ -180,7 +226,7 @@ class LinearConstraintBuilder {
   // It returns false if one literal does not have an integer view, as it
   // actually calls AddLiteralTerm().
   ABSL_MUST_USE_RESULT bool AddDecomposedProduct(
-      const std::vector<LiteralValueValue>& product);
+      absl::Span<const LiteralValueValue> product);
 
   // Add literal * coeff to the constaint. Returns false and do nothing if the
   // given literal didn't have an integer view.
@@ -293,12 +339,9 @@ bool NoDuplicateVariable(const LinearConstraint& ct);
 
 // Sorts and merges duplicate IntegerVariable in the given "terms".
 // Fills the given LinearConstraint or LinearExpression with the result.
-//
-// TODO(user): This actually only sort the terms, we don't clean them.
-template <class ClassWithVarsAndCoeffs>
-void CleanTermsAndFillConstraint(
+inline void CleanTermsAndFillConstraint(
     std::vector<std::pair<IntegerVariable, IntegerValue>>* terms,
-    ClassWithVarsAndCoeffs* output) {
+    LinearExpression* output) {
   output->vars.clear();
   output->coeffs.clear();
 
@@ -327,6 +370,39 @@ void CleanTermsAndFillConstraint(
   }
 }
 
+inline void CleanTermsAndFillConstraint(
+    std::vector<std::pair<IntegerVariable, IntegerValue>>* terms,
+    LinearConstraint* output) {
+  // Sort and add coeff of duplicate variables. Note that a variable and
+  // its negation will appear one after another in the natural order.
+  int new_size = 0;
+  output->resize(terms->size());
+  std::sort(terms->begin(), terms->end());
+  IntegerVariable previous_var = kNoIntegerVariable;
+  IntegerValue current_coeff(0);
+  for (const std::pair<IntegerVariable, IntegerValue>& entry : *terms) {
+    if (previous_var == entry.first) {
+      current_coeff += entry.second;
+    } else if (previous_var == NegationOf(entry.first)) {
+      current_coeff -= entry.second;
+    } else {
+      if (current_coeff != 0) {
+        output->vars[new_size] = previous_var;
+        output->coeffs[new_size] = current_coeff;
+        ++new_size;
+      }
+      previous_var = entry.first;
+      current_coeff = entry.second;
+    }
+  }
+  if (current_coeff != 0) {
+    output->vars[new_size] = previous_var;
+    output->coeffs[new_size] = current_coeff;
+    ++new_size;
+  }
+  output->resize(new_size);
+}
+
 }  // namespace sat
 }  // namespace operations_research
 
diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc
index 9881cd0a1e..c0741798ea 100644
--- a/ortools/sat/linear_constraint_manager.cc
+++ b/ortools/sat/linear_constraint_manager.cc
@@ -52,9 +52,9 @@ namespace {
 const LinearConstraintManager::ConstraintIndex kInvalidConstraintIndex(-1);
 
 size_t ComputeHashOfTerms(const LinearConstraint& ct) {
-  DCHECK(std::is_sorted(ct.vars.begin(), ct.vars.end()));
+  DCHECK(std::is_sorted(ct.vars.get(), ct.vars.get() + ct.num_terms));
   size_t hash = 0;
-  const int num_terms = ct.vars.size();
+  const int num_terms = ct.num_terms;
   for (int i = 0; i < num_terms; ++i) {
     hash = util_hash::Hash(ct.vars[i].value(), hash);
     hash = util_hash::Hash(ct.coeffs[i].value(), hash);
@@ -131,7 +131,7 @@ bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
 // we regenerate identical cuts for some reason.
 LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
     LinearConstraint ct, bool* added) {
-  DCHECK(!ct.vars.empty());
+  DCHECK_GT(ct.num_terms, 0);
   DCHECK(!PossibleOverflow(integer_trail_, ct)) << ct.DebugString();
   DCHECK(NoDuplicateVariable(ct));
   SimplifyConstraint(&ct);
@@ -143,8 +143,7 @@ LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
   const size_t key = ComputeHashOfTerms(ct);
   if (equiv_constraints_.contains(key)) {
     const ConstraintIndex ct_index = equiv_constraints_[key];
-    if (constraint_infos_[ct_index].constraint.vars == ct.vars &&
-        constraint_infos_[ct_index].constraint.coeffs == ct.coeffs) {
+    if (constraint_infos_[ct_index].constraint.IsEqualIgnoringBounds(ct)) {
       bool tightened = false;
       if (ct.lb > constraint_infos_[ct_index].constraint.lb) {
         tightened = true;
@@ -218,7 +217,7 @@ void LinearConstraintManager::ComputeObjectiveParallelism(
 
   const LinearConstraint& lc = constraint_infos_[ct_index].constraint;
   double unscaled_objective_parallelism = 0.0;
-  for (int i = 0; i < lc.vars.size(); ++i) {
+  for (int i = 0; i < lc.num_terms; ++i) {
     const IntegerVariable var = lc.vars[i];
     const auto it = objective_map_.find(var);
     if (it == objective_map_.end()) continue;
@@ -233,11 +232,10 @@ void LinearConstraintManager::ComputeObjectiveParallelism(
 
 // Same as Add(), but logs some information about the newly added constraint.
 // Cuts are also handled slightly differently than normal constraints.
-bool LinearConstraintManager::AddCut(const LinearConstraint& ct,
-                                     std::string type_name,
+bool LinearConstraintManager::AddCut(LinearConstraint ct, std::string type_name,
                                      std::string extra_info) {
   ++num_add_cut_calls_;
-  if (ct.vars.empty()) return false;
+  if (ct.num_terms == 0) return false;
 
   const double activity = ComputeActivity(ct, expanded_lp_solution_);
   const double violation =
@@ -247,7 +245,7 @@ bool LinearConstraintManager::AddCut(const LinearConstraint& ct,
   // Only add cut with sufficient efficacy.
   if (violation / l2_norm < 1e-4) {
     VLOG(3) << "BAD Cut '" << type_name << "'"
-            << " size=" << ct.vars.size()
+            << " size=" << ct.num_terms
             << " max_magnitude=" << ComputeInfinityNorm(ct)
             << " norm=" << l2_norm << " violation=" << violation
             << " eff=" << violation / l2_norm << " " << extra_info;
@@ -261,7 +259,7 @@ bool LinearConstraintManager::AddCut(const LinearConstraint& ct,
   if (PossibleOverflow(integer_trail_, ct)) return false;
 
   bool added = false;
-  const ConstraintIndex ct_index = Add(ct, &added);
+  const ConstraintIndex ct_index = Add(std::move(ct), &added);
 
   // We only mark the constraint as a cut if it is not an update of an already
   // existing one.
@@ -272,7 +270,7 @@ bool LinearConstraintManager::AddCut(const LinearConstraint& ct,
   constraint_infos_[ct_index].is_deletable = true;
 
   VLOG(2) << "Cut '" << type_name << "'"
-          << " size=" << constraint_infos_[ct_index].constraint.vars.size()
+          << " size=" << constraint_infos_[ct_index].constraint.num_terms
           << " max_magnitude="
           << ComputeInfinityNorm(constraint_infos_[ct_index].constraint)
           << " norm=" << l2_norm << " violation=" << violation
@@ -364,7 +362,7 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
   IntegerValue max_magnitude(0);
   IntegerValue min_magnitude = kMaxIntegerValue;
   int new_size = 0;
-  const int num_terms = ct->vars.size();
+  const int num_terms = ct->num_terms;
   for (int i = 0; i < num_terms; ++i) {
     const IntegerVariable var = ct->vars[i];
     const IntegerValue coeff = ct->coeffs[i];
@@ -408,15 +406,13 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
       ct->coeffs[new_size] = coeff;
       ++new_size;
     }
-    ct->vars.resize(new_size);
-    ct->coeffs.resize(new_size);
+    ct->resize(new_size);
   }
 
   // Clear constraints that are always true.
   // We rely on the deletion code to remove them eventually.
   if (min_sum >= ct->lb && max_sum <= ct->ub) {
-    ct->vars.clear();
-    ct->coeffs.clear();
+    ct->resize(0);
     ct->lb = 0;
     ct->ub = 0;
     return true;
@@ -453,7 +449,7 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
   if (max_magnitude > second_threshold) {
     term_changed = true;
     ++num_coeff_strenghtening_;
-    const int num_terms = ct->vars.size();
+    const int num_terms = ct->num_terms;
     for (int i = 0; i < num_terms; ++i) {
       // In all cases, we reason on a transformed constraint where the term
       // is max_value - |coeff| * positive_X. If we change coeff, and
@@ -489,7 +485,7 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
 void LinearConstraintManager::FillDerivedFields(ConstraintInfo* info) {
   IntegerValue min_sum(0);
   IntegerValue max_sum(0);
-  const int num_terms = info->constraint.vars.size();
+  const int num_terms = info->constraint.num_terms;
   for (int i = 0; i < num_terms; ++i) {
     const IntegerVariable var = info->constraint.vars[i];
     const IntegerValue coeff = info->constraint.coeffs[i];
@@ -516,7 +512,7 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
   VLOG(3) << "Enter ChangeLP, scan " << constraint_infos_.size()
           << " constraints";
   const double saved_dtime = dtime_;
-  std::vector<ConstraintIndex> new_constraints;
+  std::vector<std::pair<ConstraintIndex, double>> new_constraints_by_score;
   std::vector<double> new_constraints_efficacies;
   std::vector<double> new_constraints_orthogonalities;
 
@@ -560,8 +556,8 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
 
     // ComputeActivity() often represent the bulk of the time spent in
     // ChangeLP().
-    dtime_ += 1.7e-9 *
-              static_cast<double>(constraint_infos_[i].constraint.vars.size());
+    dtime_ +=
+        1.7e-9 * static_cast<double>(constraint_infos_[i].constraint.num_terms);
     const double activity =
         ComputeActivity(constraint_infos_[i].constraint, expanded_lp_solution_);
     const double lb_violation =
@@ -571,11 +567,6 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
     const double violation = std::max(lb_violation, ub_violation);
     if (violation >= tolerance) {
       constraint_infos_[i].inactive_count = 0;
-      new_constraints.push_back(i);
-      new_constraints_efficacies.push_back(violation /
-                                           constraint_infos_[i].l2_norm);
-      new_constraints_orthogonalities.push_back(1.0);
-
       if (objective_is_defined_ &&
           !constraint_infos_[i].objective_parallelism_computed) {
         ComputeObjectiveParallelism(i);
@@ -583,9 +574,9 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
         constraint_infos_[i].objective_parallelism = 0.0;
       }
 
-      constraint_infos_[i].current_score =
-          new_constraints_efficacies.back() +
-          constraint_infos_[i].objective_parallelism;
+      const double score = violation / constraint_infos_[i].l2_norm +
+                           constraint_infos_[i].objective_parallelism;
+      new_constraints_by_score.push_back({i, score});
 
       if (constraint_infos_[i].is_deletable) {
         constraint_infos_[i].active_count += constraint_active_count_increase_;
@@ -646,53 +637,58 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
   // TODO(user): This blowup factor could be adaptative w.r.t. the constraint
   // limit.
   const int kBlowupFactor = 4;
-  int constraint_limit = std::min(sat_parameters_.new_constraints_batch_size(),
-                                  static_cast<int>(new_constraints.size()));
+  const int current_size = static_cast<int>(new_constraints_by_score.size());
+  int constraint_limit =
+      std::min(sat_parameters_.new_constraints_batch_size(), current_size);
   if (lp_constraints_.empty()) {
-    constraint_limit = std::min(1000, static_cast<int>(new_constraints.size()));
+    constraint_limit = std::min(1000, current_size);
   }
-  VLOG(3) << "   - size = " << new_constraints.size()
-          << ", limit = " << constraint_limit;
+  VLOG(3) << "   - size = " << current_size << ", limit = " << constraint_limit;
 
-  std::stable_sort(new_constraints.begin(), new_constraints.end(),
-                   [&](ConstraintIndex a, ConstraintIndex b) {
-                     return constraint_infos_[a].current_score >
-                            constraint_infos_[b].current_score;
+  std::stable_sort(new_constraints_by_score.begin(),
+                   new_constraints_by_score.end(),
+                   [&](std::pair<ConstraintIndex, double> a,
+                       std::pair<ConstraintIndex, double> b) {
+                     return a.second > b.second;
                    });
-  if (new_constraints.size() > kBlowupFactor * constraint_limit) {
-    VLOG(3) << "Resize candidate constraints from " << new_constraints.size()
-            << " down to " << kBlowupFactor * constraint_limit;
-    new_constraints.resize(kBlowupFactor * constraint_limit);
+  if (new_constraints_by_score.size() > kBlowupFactor * constraint_limit) {
+    VLOG(3) << "Resize candidate constraints from "
+            << new_constraints_by_score.size() << " down to "
+            << kBlowupFactor * constraint_limit;
+    new_constraints_by_score.resize(kBlowupFactor * constraint_limit);
   }
 
   int num_added = 0;
   int num_skipped_checks = 0;
   const int kCheckFrequency = 100;
   ConstraintIndex last_added_candidate = kInvalidConstraintIndex;
+  std::vector<double> orthogonality_score(new_constraints_by_score.size(), 1.0);
   for (int i = 0; i < constraint_limit; ++i) {
     // Iterate through all new constraints and select the one with the best
     // score.
+    //
+    // TODO(user): find better algo, this does 1000 * 4000 scalar product!
     double best_score = 0.0;
     ConstraintIndex best_candidate = kInvalidConstraintIndex;
-    for (int j = 0; j < new_constraints.size(); ++j) {
+    for (int j = 0; j < new_constraints_by_score.size(); ++j) {
       // Checks the time limit, and returns if the lp has changed.
       if (++num_skipped_checks >= kCheckFrequency) {
         if (time_limit_->LimitReached()) return current_lp_is_changed_;
         num_skipped_checks = 0;
       }
 
-      const ConstraintIndex new_constraint = new_constraints[j];
-      if (constraint_infos_[new_constraint].is_in_lp) continue;
+      const ConstraintIndex new_index = new_constraints_by_score[j].first;
+      if (constraint_infos_[new_index].is_in_lp) continue;
 
       if (last_added_candidate != kInvalidConstraintIndex) {
         const double current_orthogonality =
             1.0 - (std::abs(ScalarProduct(
                        constraint_infos_[last_added_candidate].constraint,
-                       constraint_infos_[new_constraint].constraint)) /
+                       constraint_infos_[new_index].constraint)) /
                    (constraint_infos_[last_added_candidate].l2_norm *
-                    constraint_infos_[new_constraint].l2_norm));
-        new_constraints_orthogonalities[j] =
-            std::min(new_constraints_orthogonalities[j], current_orthogonality);
+                    constraint_infos_[new_index].l2_norm));
+        orthogonality_score[j] =
+            std::min(orthogonality_score[j], current_orthogonality);
       }
 
       // NOTE(user): It is safe to not add this constraint as the constraint
@@ -700,19 +696,19 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
       // inactive for a long time and is removed from the LP. In either case,
       // this constraint is not adding significant value and is only making the
       // LP larger.
-      if (new_constraints_orthogonalities[j] <
+      if (orthogonality_score[j] <
           sat_parameters_.min_orthogonality_for_lp_constraints()) {
         continue;
       }
 
       // TODO(user): Experiment with different weights or different
       // functions for computing score.
-      const double score = new_constraints_orthogonalities[j] +
-                           constraint_infos_[new_constraint].current_score;
+      const double score =
+          orthogonality_score[j] + new_constraints_by_score[j].second;
       CHECK_GE(score, 0.0);
       if (score > best_score || best_candidate == kInvalidConstraintIndex) {
         best_score = score;
-        best_candidate = new_constraint;
+        best_candidate = new_index;
       }
     }
 
@@ -726,6 +722,9 @@ bool LinearConstraintManager::ChangeLp(glop::BasisState* solution_state,
       current_lp_is_changed_ = true;
       lp_constraints_.push_back(best_candidate);
       last_added_candidate = best_candidate;
+    } else {
+      // Abort the addition loop.
+      break;
     }
   }
 
@@ -770,7 +769,7 @@ bool LinearConstraintManager::DebugCheckConstraint(
   const auto& debug_solution = *(model_->Get<DebugSolution>());
 
   IntegerValue activity(0);
-  for (int i = 0; i < cut.vars.size(); ++i) {
+  for (int i = 0; i < cut.num_terms; ++i) {
     const IntegerVariable var = cut.vars[i];
     const IntegerValue coeff = cut.coeffs[i];
     CHECK(debug_solution.ivar_has_value[var]);
@@ -788,17 +787,17 @@ bool LinearConstraintManager::DebugCheckConstraint(
 void TopNCuts::AddCut(
     LinearConstraint ct, absl::string_view name,
     const absl::StrongVector<IntegerVariable, double>& lp_solution) {
-  if (ct.vars.empty()) return;
+  if (ct.num_terms == 0) return;
   const double activity = ComputeActivity(ct, lp_solution);
   const double violation =
       std::max(activity - ToDouble(ct.ub), ToDouble(ct.lb) - activity);
   const double l2_norm = ComputeL2Norm(ct);
-  cuts_.Add({std::string(name), ct}, violation / l2_norm);
+  cuts_.Add({std::string(name), std::move(ct)}, violation / l2_norm);
 }
 
 void TopNCuts::TransferToManager(LinearConstraintManager* manager) {
-  for (const CutCandidate& candidate : cuts_.UnorderedElements()) {
-    manager->AddCut(candidate.cut, candidate.name);
+  for (CutCandidate& candidate : *cuts_.MutableUnorderedElements()) {
+    manager->AddCut(std::move(candidate.cut), candidate.name);
   }
   cuts_.Clear();
 }
diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h
index f32c808734..ebf799238f 100644
--- a/ortools/sat/linear_constraint_manager.h
+++ b/ortools/sat/linear_constraint_manager.h
@@ -69,20 +69,25 @@ class LinearConstraintManager {
     LinearConstraint constraint;
 
     double l2_norm = 0.0;
-    int64_t inactive_count = 0;
     double objective_parallelism = 0.0;
-    bool objective_parallelism_computed = false;
-    bool is_in_lp = false;
-    bool ub_is_trivial = false;
-    bool lb_is_trivial = false;
     size_t hash;
-    double current_score = 0.0;
 
     // Updated only for deletable constraints. This is incremented every time
     // ChangeLp() is called and the constraint is active in the LP or not in the
     // LP and violated.
     double active_count = 0.0;
 
+    // TODO(user): This is the number of time the constraint was consecutively
+    // inactive, and go up to 100 with the default param, so we could optimize
+    // the space used more.
+    uint16_t inactive_count = 0;
+
+    // TODO(user): Pack bool and in general optimize the memory of this class.
+    bool objective_parallelism_computed = false;
+    bool is_in_lp = false;
+    bool ub_is_trivial = false;
+    bool lb_is_trivial = false;
+
     // For now, we mark all the generated cuts as deletable and the problem
     // constraints as undeletable.
     // TODO(user): We can have a better heuristics. Some generated good cuts
@@ -113,7 +118,7 @@ class LinearConstraintManager {
   //
   // Returns true if a new cut was added and false if this cut is not
   // efficacious or if it is a duplicate of an already existing one.
-  bool AddCut(const LinearConstraint& ct, std::string type_name,
+  bool AddCut(LinearConstraint ct, std::string type_name,
               std::string extra_info = "");
 
   // These must be level zero bounds.
diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc
index f2f0e18173..bde918e983 100644
--- a/ortools/sat/linear_programming_constraint.cc
+++ b/ortools/sat/linear_programming_constraint.cc
@@ -22,6 +22,7 @@
 #include <limits>
 #include <memory>
 #include <numeric>
+#include <optional>
 #include <random>
 #include <string>
 #include <utility>
@@ -140,30 +141,63 @@ bool ScatteredIntegerVector::AddLinearExpressionMultiple(
   return true;
 }
 
-void ScatteredIntegerVector::ConvertToLinearConstraint(
-    const std::vector<IntegerVariable>& integer_variables,
-    IntegerValue upper_bound, LinearConstraint* result) {
-  result->vars.clear();
-  result->coeffs.clear();
+LinearConstraint ScatteredIntegerVector::ConvertToLinearConstraint(
+    absl::Span<const IntegerVariable> integer_variables,
+    IntegerValue upper_bound,
+    std::optional<std::pair<IntegerVariable, IntegerValue>> extra_term) {
+  // We first do one pass to compute the exact size and not overallocate.
+  int final_size = 0;
+  if (is_sparse_) {
+    for (const glop::ColIndex col : non_zeros_) {
+      const IntegerValue coeff = dense_vector_[col];
+      if (coeff == 0) continue;
+      ++final_size;
+    }
+  } else {
+    for (const IntegerValue coeff : dense_vector_) {
+      if (coeff != 0) ++final_size;
+    }
+  }
+  if (extra_term != std::nullopt) ++final_size;
+
+  // Allocate once.
+  LinearConstraint result;
+  result.resize(final_size);
+
+  // Copy terms.
+  int new_size = 0;
   if (is_sparse_) {
     std::sort(non_zeros_.begin(), non_zeros_.end());
     for (const glop::ColIndex col : non_zeros_) {
       const IntegerValue coeff = dense_vector_[col];
       if (coeff == 0) continue;
-      result->vars.push_back(integer_variables[col.value()]);
-      result->coeffs.push_back(coeff);
+      result.vars[new_size] = integer_variables[col.value()];
+      result.coeffs[new_size] = coeff;
+      ++new_size;
     }
   } else {
     const int size = dense_vector_.size();
     for (glop::ColIndex col(0); col < size; ++col) {
       const IntegerValue coeff = dense_vector_[col];
       if (coeff == 0) continue;
-      result->vars.push_back(integer_variables[col.value()]);
-      result->coeffs.push_back(coeff);
+      result.vars[new_size] = integer_variables[col.value()];
+      result.coeffs[new_size] = coeff;
+      ++new_size;
     }
   }
-  result->lb = kMinIntegerValue;
-  result->ub = upper_bound;
+
+  result.lb = kMinIntegerValue;
+  result.ub = upper_bound;
+
+  if (extra_term != std::nullopt) {
+    result.vars[new_size] += extra_term->first;
+    result.coeffs[new_size] += extra_term->second;
+    ++new_size;
+  }
+
+  CHECK_EQ(new_size, final_size);
+  DivideByGCD(&result);
+  return result;
 }
 
 void ScatteredIntegerVector::ConvertToCutData(
@@ -336,7 +370,7 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
     new_ct.ub = ct.ub;
     new_ct.lb_is_trivial = all_constraints[index].lb_is_trivial;
     new_ct.ub_is_trivial = all_constraints[index].ub_is_trivial;
-    const int size = ct.vars.size();
+    const int size = ct.num_terms;
     if (ct.lb > ct.ub) {
       VLOG(1) << "Trivial infeasible bound in an LP constraint";
       return false;
@@ -604,7 +638,12 @@ void LinearProgrammingConstraint::RegisterWith(Model* model) {
   std::sort(integer_objective_.begin(), integer_objective_.end());
 
   // Set the LP to its initial content.
-  if (!parameters_.add_lp_constraints_lazily()) {
+  //
+  // Note that we always add LP constraint lazily if we have A LOT of them.
+  // This is because currently on large problem with millions of constraints,
+  // our LP is usually not fast enough anyway.
+  if (!parameters_.add_lp_constraints_lazily() &&
+      constraint_manager_.num_constraints() < 1e6) {
     constraint_manager_.AddAllConstraintsToLp();
   }
   if (!CreateLpFromConstraintManager()) {
@@ -1247,17 +1286,20 @@ bool LinearProgrammingConstraint::PostprocessAndAddCut(
     }
   }
 
-  tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, cut_ub,
-                                                  &cut_);
-  DivideByGCD(&cut_);
+  // TODO(user): avoid allocating memory if it turns out this is a duplicate of
+  // something we already added. This tends to happen if the duplicate was
+  // already a generated cut which is currently not part of the LP.
+  LinearConstraint converted_cut =
+      tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_,
+                                                      cut_ub);
 
   // TODO(user): We probably generate too many cuts, keep best one only?
   // Note that we do need duplicate removal and maybe some orthogonality here?
   if (/*DISABLES CODE*/ (false)) {
-    top_n_cuts_.AddCut(cut_, name, expanded_lp_solution_);
+    top_n_cuts_.AddCut(std::move(converted_cut), name, expanded_lp_solution_);
     return true;
   }
-  return constraint_manager_.AddCut(cut_, name, info);
+  return constraint_manager_.AddCut(std::move(converted_cut), name, info);
 }
 
 // TODO(user): This can be still too slow on some problems like
@@ -1364,23 +1406,26 @@ void LinearProgrammingConstraint::AddObjectiveCut() {
   objective_ct.ub = integer_objective_offset_ -
                     integer_trail_->LevelZeroLowerBound(objective_cp_);
   IntegerValue obj_coeff_magnitude(0);
+  objective_ct.resize(integer_objective_.size());
+  int i = 0;
   for (const auto& [col, coeff] : integer_objective_) {
     const IntegerVariable var = integer_variables_[col.value()];
-    objective_ct.vars.push_back(var);
-    objective_ct.coeffs.push_back(-coeff);
+    objective_ct.vars[i] = var;
+    objective_ct.coeffs[i] = -coeff;
     obj_coeff_magnitude = std::max(obj_coeff_magnitude, IntTypeAbs(coeff));
+    ++i;
+  }
+
+  if (!base_ct_.FillFromLinearConstraint(objective_ct, expanded_lp_solution_,
+                                         integer_trail_)) {
+    return;
   }
 
   // If the magnitude is small enough, just try to add the full objective. Other
   // cuts will be derived in subsequent passes. Otherwise, try normal cut
   // heuristic that should result in a cut with reasonable coefficients.
   if (obj_coeff_magnitude < 1e9 &&
-      constraint_manager_.AddCut(objective_ct, "Objective")) {
-    return;
-  }
-
-  if (!base_ct_.FillFromLinearConstraint(objective_ct, expanded_lp_solution_,
-                                         integer_trail_)) {
+      constraint_manager_.AddCut(std::move(objective_ct), "Objective")) {
     return;
   }
 
@@ -1911,7 +1956,7 @@ bool LinearProgrammingConstraint::Propagate() {
 absl::int128 LinearProgrammingConstraint::GetImpliedLowerBound(
     const LinearConstraint& terms) const {
   absl::int128 lower_bound(0);
-  const int size = terms.vars.size();
+  const int size = terms.num_terms;
   for (int i = 0; i < size; ++i) {
     const IntegerVariable var = terms.vars[i];
     const IntegerValue coeff = terms.coeffs[i];
@@ -2032,7 +2077,7 @@ LinearProgrammingConstraint::ScaleLpMultiplier(
 
 template <bool check_overflow>
 bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
-    const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers,
+    absl::Span<const std::pair<RowIndex, IntegerValue>> integer_multipliers,
     ScatteredIntegerVector* scattered_vector, IntegerValue* upper_bound) const {
   // Initialize the new constraint.
   *upper_bound = 0;
@@ -2225,16 +2270,18 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint(
   if (adjusted) ++num_adjusts_;
 }
 
-bool LinearProgrammingConstraint::PropagateLpConstraint(
-    const LinearConstraint& ct) {
+bool LinearProgrammingConstraint::PropagateLpConstraint(LinearConstraint ct) {
   DCHECK(constraint_manager_.DebugCheckConstraint(ct));
-  if (ct.vars.empty()) {
+
+  // We need to cache this before we std::move() the constraint!
+  const int num_terms = ct.num_terms;
+  if (num_terms == 0) {
     if (ct.ub >= 0) return true;
     return integer_trail_->ReportConflict({});  // Unsat.
   }
 
   std::unique_ptr<IntegerSumLE128> cp_constraint(
-      new IntegerSumLE128({}, ct.vars, ct.coeffs, ct.ub, model_));
+      new IntegerSumLE128(std::move(ct), model_));
 
   // We always propagate level zero bounds first.
   // If we are at level zero, there is nothing else to do.
@@ -2252,7 +2299,7 @@ bool LinearProgrammingConstraint::PropagateLpConstraint(
           ? 0
           : cumulative_optimal_constraint_sizes_.back();
   optimal_constraints_.push_back(std::move(cp_constraint));
-  cumulative_optimal_constraint_sizes_.push_back(current_size + ct.vars.size());
+  cumulative_optimal_constraint_sizes_.push_back(current_size + num_terms);
   rev_optimal_constraints_size_ = optimal_constraints_.size();
   return no_conflict;
 }
@@ -2327,19 +2374,18 @@ bool LinearProgrammingConstraint::PropagateExactLpReason() {
 
   // Create the IntegerSumLE that will allow to propagate the objective and more
   // generally do the reduced cost fixing.
-  tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, rc_ub,
-                                                  &tmp_constraint_);
-  tmp_constraint_.vars.push_back(objective_cp_);
-  tmp_constraint_.coeffs.push_back(-obj_scale);
-  DivideByGCD(&tmp_constraint_);
+  LinearConstraint explanation =
+      tmp_scattered_vector_.ConvertToLinearConstraint(
+          integer_variables_, rc_ub,
+          /*extra_term=*/{{objective_cp_, -obj_scale}});
 
   // Corner case where prevent overflow removed all terms.
-  if (tmp_constraint_.vars.empty()) {
+  if (explanation.num_terms == 0) {
     trail_->MutableConflict()->clear();
-    return tmp_constraint_.ub >= 0;
+    return explanation.ub >= 0;
   }
 
-  return PropagateLpConstraint(tmp_constraint_);
+  return PropagateLpConstraint(std::move(explanation));
 }
 
 bool LinearProgrammingConstraint::PropagateExactDualRay() {
@@ -2367,17 +2413,22 @@ bool LinearProgrammingConstraint::PropagateExactDualRay() {
   AdjustNewLinearConstraint(&tmp_integer_multipliers_, &tmp_scattered_vector_,
                             &new_constraint_ub);
 
-  tmp_scattered_vector_.ConvertToLinearConstraint(
-      integer_variables_, new_constraint_ub, &tmp_constraint_);
-  DivideByGCD(&tmp_constraint_);
+  LinearConstraint explanation =
+      tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_,
+                                                      new_constraint_ub);
+
+  std::string message;
+  if (VLOG_IS_ON(1)) {
+    // Unfortunately, we need to set this up before we std::move() it.
+    message = absl::StrCat("LP exact dual ray not infeasible,",
+                           " implied_lb: ", GetImpliedLowerBound(explanation),
+                           " ub: ", explanation.ub.value());
+  }
 
   // This should result in a conflict if the precision is good enough.
-  if (!PropagateLpConstraint(tmp_constraint_)) return false;
+  if (!PropagateLpConstraint(std::move(explanation))) return false;
 
-  const absl::int128 implied_lb = GetImpliedLowerBound(tmp_constraint_);
-  VLOG(1) << "LP exact dual ray not infeasible,"
-          << " implied_lb: " << implied_lb
-          << " ub: " << tmp_constraint_.ub.value();
+  VLOG(1) << message;
   return true;
 }
 
diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h
index f8fe01672a..c7118de432 100644
--- a/ortools/sat/linear_programming_constraint.h
+++ b/ortools/sat/linear_programming_constraint.h
@@ -86,9 +86,11 @@ class ScatteredIntegerVector {
   //
   // TODO(user): Ideally we should convert to IntegerVariable as late as
   // possible. Prefer to use GetTerms().
-  void ConvertToLinearConstraint(
-      const std::vector<IntegerVariable>& integer_variables,
-      IntegerValue upper_bound, LinearConstraint* result);
+  LinearConstraint ConvertToLinearConstraint(
+      absl::Span<const IntegerVariable> integer_variables,
+      IntegerValue upper_bound,
+      std::optional<std::pair<IntegerVariable, IntegerValue>> extra_term =
+          std::nullopt);
 
   void ConvertToCutData(absl::int128 rhs,
                         const std::vector<IntegerVariable>& integer_variables,
@@ -332,7 +334,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
 
   // Called by PropagateExactLpReason() and PropagateExactDualRay() to finish
   // propagation.
-  bool PropagateLpConstraint(const LinearConstraint& ct);
+  bool PropagateLpConstraint(LinearConstraint ct);
 
   // Returns number of non basic variables with zero reduced costs.
   int64_t CalculateDegeneracy();
@@ -364,7 +366,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // is false, we do not check for a bit of extra speed.
   template <bool check_overflow = true>
   bool ComputeNewLinearConstraint(
-      const std::vector<std::pair<glop::RowIndex, IntegerValue>>&
+      absl::Span<const std::pair<glop::RowIndex, IntegerValue>>
           integer_multipliers,
       ScatteredIntegerVector* scattered_vector,
       IntegerValue* upper_bound) const;
@@ -484,9 +486,6 @@ class LinearProgrammingConstraint : public PropagatorInterface,
 
   bool problem_proven_infeasible_by_cuts_ = false;
   CutData base_ct_;
-  LinearConstraint cut_;
-  LinearConstraint saved_cut_;
-  LinearConstraint tmp_constraint_;
 
   ScatteredIntegerVector tmp_scattered_vector_;
 
diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc
index e1c139577b..5fdf18a9db 100644
--- a/ortools/sat/linear_relaxation.cc
+++ b/ortools/sat/linear_relaxation.cc
@@ -209,7 +209,7 @@ void AppendRelaxationForEqualityEncoding(IntegerVariable var,
     // precondition but not the Var = sum bool * value one. In this case, we
     // just don't encode it this way. Hopefully, most normal model will not run
     // into this.
-    const LinearConstraint lc = encoding_ct.Build();
+    LinearConstraint lc = encoding_ct.Build();
     if (!PossibleOverflow(*integer_trail, lc)) {
       relaxation->linear_constraints.push_back(at_least_one.Build());
       relaxation->linear_constraints.push_back(std::move(lc));
@@ -230,7 +230,7 @@ void AppendRelaxationForEqualityEncoding(IntegerVariable var,
                                        rhs - value_literal.value));
     }
 
-    const LinearConstraint lc = encoding_ct.Build();
+    LinearConstraint lc = encoding_ct.Build();
     if (!PossibleOverflow(*integer_trail, lc)) {
       relaxation->at_most_ones.push_back(at_most_one_ct);
       relaxation->linear_constraints.push_back(std::move(lc));
@@ -248,7 +248,7 @@ void AppendRelaxationForEqualityEncoding(IntegerVariable var,
     CHECK(lower_bound_ct.AddLiteralTerm(value_literal.literal,
                                         d_min - value_literal.value));
   }
-  const LinearConstraint built_lb_ct = lower_bound_ct.Build();
+  LinearConstraint built_lb_ct = lower_bound_ct.Build();
   if (!PossibleOverflow(*integer_trail, built_lb_ct)) {
     relaxation->linear_constraints.push_back(std::move(built_lb_ct));
     ++*num_loose;
@@ -263,7 +263,7 @@ void AppendRelaxationForEqualityEncoding(IntegerVariable var,
     CHECK(upper_bound_ct.AddLiteralTerm(value_literal.literal,
                                         d_max - value_literal.value));
   }
-  const LinearConstraint built_ub_ct = upper_bound_ct.Build();
+  LinearConstraint built_ub_ct = upper_bound_ct.Build();
   if (!PossibleOverflow(*integer_trail, built_ub_ct)) {
     relaxation->linear_constraints.push_back(std::move(built_ub_ct));
     ++*num_loose;
@@ -335,7 +335,7 @@ void AppendPartialGreaterThanEncodingRelaxation(IntegerVariable var,
 namespace {
 
 bool AllLiteralsHaveViews(const IntegerEncoder& encoder,
-                          const std::vector<Literal>& literals) {
+                          absl::Span<const Literal> literals) {
   for (const Literal lit : literals) {
     if (!encoder.LiteralOrNegationHasView(lit)) return false;
   }
@@ -1718,7 +1718,7 @@ void AppendElementEncodingRelaxation(Model* m, LinearRelaxation* relaxation) {
         }
       }
       ++num_exactly_one_elements;
-      const LinearConstraint lc = builder.Build();
+      LinearConstraint lc = builder.Build();
       if (!PossibleOverflow(*integer_trail, lc)) {
         relaxation->linear_constraints.push_back(std::move(lc));
       }
@@ -1844,7 +1844,7 @@ LinearRelaxation ComputeLinearRelaxation(const CpModelProto& model_proto,
       std::remove_if(
           relaxation.linear_constraints.begin(),
           relaxation.linear_constraints.end(),
-          [](const LinearConstraint& lc) { return lc.vars.size() <= 1; }),
+          [](const LinearConstraint& lc) { return lc.num_terms <= 1; }),
       relaxation.linear_constraints.end());
 
   // We add a clique cut generation over all Booleans of the problem.
diff --git a/ortools/sat/max_hs.cc b/ortools/sat/max_hs.cc
index 2c3e65810b..9ef9e66f39 100644
--- a/ortools/sat/max_hs.cc
+++ b/ortools/sat/max_hs.cc
@@ -29,6 +29,7 @@
 #include "absl/meta/type_traits.h"
 #include "absl/random/random.h"
 #include "absl/strings/string_view.h"
+#include "absl/types/span.h"
 #include "ortools/base/logging.h"
 #include "ortools/base/strong_vector.h"
 #include "ortools/sat/presolve_util.h"
@@ -208,7 +209,7 @@ void HittingSetOptimizer::ExtractObjectiveVariables() {
 }
 
 void HittingSetOptimizer::ExtractAdditionalVariables(
-    const std::vector<IntegerVariable>& to_extract) {
+    absl::Span<const IntegerVariable> to_extract) {
   MPModelProto* hs_model = request_.mutable_model();
 
   VLOG(2) << "Extract " << to_extract.size() << " additional variables";
@@ -268,14 +269,14 @@ HittingSetOptimizer::ComputeAdditionalVariablesToExtract() {
 
   for (const LinearConstraint& linear : relaxation_.linear_constraints) {
     bool found_at_least_one = extract_all;
-    for (const IntegerVariable var : linear.vars) {
+    for (const IntegerVariable var : linear.VarsAsSpan()) {
       if (GetExtractedIndex(var) != kUnextracted) {
         found_at_least_one = true;
       }
       if (found_at_least_one) break;
     }
     if (!found_at_least_one) continue;
-    for (const IntegerVariable var : linear.vars) {
+    for (const IntegerVariable var : linear.VarsAsSpan()) {
       if (GetExtractedIndex(var) == kUnextracted) {
         result_set.insert(PositiveVariable(var));
       }
@@ -305,7 +306,7 @@ void HittingSetOptimizer::ProjectAndAddAtMostOne(
 MPConstraintProto* HittingSetOptimizer::ProjectAndAddLinear(
     const LinearConstraint& linear) {
   int num_extracted_variables = 0;
-  for (int i = 0; i < linear.vars.size(); ++i) {
+  for (int i = 0; i < linear.num_terms; ++i) {
     if (GetExtractedIndex(PositiveVariable(linear.vars[i])) != kUnextracted) {
       num_extracted_variables++;
     }
@@ -322,7 +323,7 @@ void HittingSetOptimizer::ProjectLinear(const LinearConstraint& linear,
   IntegerValue lb = linear.lb;
   IntegerValue ub = linear.ub;
 
-  for (int i = 0; i < linear.vars.size(); ++i) {
+  for (int i = 0; i < linear.num_terms; ++i) {
     const IntegerVariable var = linear.vars[i];
     const IntegerValue coeff = linear.coeffs[i];
     const int index = GetExtractedIndex(PositiveVariable(var));
diff --git a/ortools/sat/max_hs.h b/ortools/sat/max_hs.h
index 179b8bd266..fa5e1546aa 100644
--- a/ortools/sat/max_hs.h
+++ b/ortools/sat/max_hs.h
@@ -22,6 +22,7 @@
 #include <vector>
 
 #include "absl/container/flat_hash_map.h"
+#include "absl/types/span.h"
 #include "ortools/base/strong_vector.h"
 #include "ortools/linear_solver/linear_solver.pb.h"
 #include "ortools/sat/cp_model.pb.h"
@@ -76,8 +77,7 @@ class HittingSetOptimizer {
 
   // Calls ComputeAdditionalVariablesToExtract() and extract all new variables.
   // This must be called after the linear relaxation has been filled.
-  void ExtractAdditionalVariables(
-      const std::vector<IntegerVariable>& to_extract);
+  void ExtractAdditionalVariables(absl::Span<const IntegerVariable> to_extract);
 
   // Heuristic to decide which variables (in addition to the objective
   // variables) to extract.
diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc
index 99bcb49da3..78932127d8 100644
--- a/ortools/sat/optimization.cc
+++ b/ortools/sat/optimization.cc
@@ -434,9 +434,8 @@ SatSolver::Status FindCores(std::vector<Literal> assumptions,
 }  // namespace
 
 CoreBasedOptimizer::CoreBasedOptimizer(
-    IntegerVariable objective_var,
-    const std::vector<IntegerVariable>& variables,
-    const std::vector<IntegerValue>& coefficients,
+    IntegerVariable objective_var, absl::Span<const IntegerVariable> variables,
+    absl::Span<const IntegerValue> coefficients,
     std::function<void()> feasible_solution_observer, Model* model)
     : parameters_(model->GetOrCreate<SatParameters>()),
       sat_solver_(model->GetOrCreate<SatSolver>()),
diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h
index 61522b34e9..5820783906 100644
--- a/ortools/sat/optimization.h
+++ b/ortools/sat/optimization.h
@@ -17,6 +17,7 @@
 #include <functional>
 #include <vector>
 
+#include "absl/types/span.h"
 #include "ortools/sat/clause.h"
 #include "ortools/sat/cp_model_mapping.h"
 #include "ortools/sat/integer.h"
@@ -92,8 +93,8 @@ void PresolveBooleanLinearExpression(std::vector<Literal>* literals,
 class CoreBasedOptimizer {
  public:
   CoreBasedOptimizer(IntegerVariable objective_var,
-                     const std::vector<IntegerVariable>& variables,
-                     const std::vector<IntegerValue>& coefficients,
+                     absl::Span<const IntegerVariable> variables,
+                     absl::Span<const IntegerValue> coefficients,
                      std::function<void()> feasible_solution_observer,
                      Model* model);
 
diff --git a/ortools/sat/pb_constraint.cc b/ortools/sat/pb_constraint.cc
index 76a79d35cd..6a4553c82e 100644
--- a/ortools/sat/pb_constraint.cc
+++ b/ortools/sat/pb_constraint.cc
@@ -239,7 +239,7 @@ bool CanonicalBooleanLinearProblem::AddLinearConstraint(
 }
 
 bool CanonicalBooleanLinearProblem::AddConstraint(
-    const std::vector<LiteralWithCoeff>& cst, Coefficient max_value,
+    absl::Span<const LiteralWithCoeff> cst, Coefficient max_value,
     Coefficient rhs) {
   if (rhs < 0) return false;          // Trivially unsatisfiable.
   if (rhs >= max_value) return true;  // Trivially satisfiable.
diff --git a/ortools/sat/pb_constraint.h b/ortools/sat/pb_constraint.h
index 576275c780..65f803d100 100644
--- a/ortools/sat/pb_constraint.h
+++ b/ortools/sat/pb_constraint.h
@@ -174,7 +174,7 @@ class CanonicalBooleanLinearProblem {
   }
 
  private:
-  bool AddConstraint(const std::vector<LiteralWithCoeff>& cst,
+  bool AddConstraint(absl::Span<const LiteralWithCoeff> cst,
                      Coefficient max_value, Coefficient rhs);
 
   std::vector<Coefficient> rhs_;
diff --git a/ortools/sat/rins.cc b/ortools/sat/rins.cc
index d6602a884d..c5bccd2736 100644
--- a/ortools/sat/rins.cc
+++ b/ortools/sat/rins.cc
@@ -25,6 +25,7 @@
 #include "absl/log/check.h"
 #include "absl/random/bit_gen_ref.h"
 #include "absl/random/distributions.h"
+#include "ortools/sat/cp_model_mapping.h"
 #include "ortools/sat/integer.h"
 #include "ortools/sat/linear_programming_constraint.h"
 #include "ortools/sat/model.h"
@@ -37,20 +38,26 @@ void RecordLPRelaxationValues(Model* model) {
   auto* lp_solutions = model->Mutable<SharedLPSolutionRepository>();
   if (lp_solutions == nullptr) return;
 
-  const LPVariables& lp_vars = *model->GetOrCreate<LPVariables>();
+  auto* mapping = model->GetOrCreate<CpModelMapping>();
+  auto* lp_values = model->GetOrCreate<ModelLpValues>();
+
+  // TODO(user): The default of ::infinity() for variable for which we do not
+  // have any LP solution is weird and inconsistent with ModelLpValues default
+  // which is zero. Fix. Note that in practice, at linearization level 2, all
+  // variable will eventually have an lp relaxation value, so it shoulnd't
+  // matter much to just use zero in RINS/RENS.
   std::vector<double> relaxation_values(
-      lp_vars.model_vars_size, std::numeric_limits<double>::infinity());
+      mapping->NumProtoVariables(), std::numeric_limits<double>::infinity());
 
-  auto* integer_trail = model->GetOrCreate<IntegerTrail>();
-  for (const LPVariable& lp_var : lp_vars.vars) {
-    const IntegerVariable positive_var = lp_var.positive_var;
-    if (integer_trail->IsCurrentlyIgnored(positive_var)) continue;
-
-    LinearProgrammingConstraint* lp = lp_var.lp;
-    if (lp == nullptr || !lp->HasSolution()) continue;
-
-    relaxation_values[lp_var.model_var] = lp->GetSolutionValue(positive_var);
+  // We only loop over the positive variables.
+  const int size = lp_values->size();
+  for (IntegerVariable var(0); var < size; var += 2) {
+    const int proto_var = mapping->GetProtoVariableFromIntegerVariable(var);
+    if (proto_var != -1) {
+      relaxation_values[proto_var] = (*lp_values)[var];
+    }
   }
+
   lp_solutions->NewLPSolution(std::move(relaxation_values));
 }
 
diff --git a/ortools/sat/rins.h b/ortools/sat/rins.h
index e7969ef502..29513bfec5 100644
--- a/ortools/sat/rins.h
+++ b/ortools/sat/rins.h
@@ -29,24 +29,6 @@
 namespace operations_research {
 namespace sat {
 
-// Links IntegerVariable with model variable and its lp constraint if any.
-struct LPVariable {
-  IntegerVariable positive_var = kNoIntegerVariable;
-  LinearProgrammingConstraint* lp = nullptr;
-  int model_var;
-
-  bool operator==(const LPVariable other) const {
-    return (positive_var == other.positive_var && lp == other.lp &&
-            model_var == other.model_var);
-  }
-};
-
-// This is used to "cache" in the model the set of relevant LPVariable.
-struct LPVariables {
-  std::vector<LPVariable> vars;
-  int model_vars_size = 0;
-};
-
 // A RINS Neighborhood is actually just a generic neighborhood where the domain
 // of some variable have been reduced (fixed or restricted in [lb, ub]).
 //
diff --git a/ortools/sat/sat_decision.cc b/ortools/sat/sat_decision.cc
index 5a6cebeb75..2691265350 100644
--- a/ortools/sat/sat_decision.cc
+++ b/ortools/sat/sat_decision.cc
@@ -48,8 +48,6 @@ void SatDecisionPolicy::IncreaseNumVariables(int num_variables) {
   num_bumps_.clear();
   pq_need_update_for_var_at_trail_index_.IncreaseSize(num_variables);
 
-  weighted_sign_.resize(num_variables, 0.0);
-
   has_forced_polarity_.resize(num_variables, false);
   forced_polarity_.resize(num_variables);
   has_target_polarity_.resize(num_variables, false);
@@ -164,12 +162,6 @@ void SatDecisionPolicy::ResetDecisionHeuristic() {
 
 void SatDecisionPolicy::ResetInitialPolarity(int from, bool inverted) {
   // Sets the initial polarity.
-  //
-  // TODO(user): The WEIGHTED_SIGN one are currently slightly broken because the
-  // weighted_sign_ is updated after this has been called. It requires a call
-  // to ResetDecisionHeuristic() after all the constraint have been added. Fix.
-  // On another hand, this is only used with SolveWithRandomParameters() that
-  // does call this function.
   const int num_variables = activities_.size();
   for (BooleanVariable var(from); var < num_variables; ++var) {
     switch (parameters_.initial_polarity()) {
@@ -182,12 +174,6 @@ void SatDecisionPolicy::ResetInitialPolarity(int from, bool inverted) {
       case SatParameters::POLARITY_RANDOM:
         var_polarity_[var] = std::uniform_int_distribution<int>(0, 1)(*random_);
         break;
-      case SatParameters::POLARITY_WEIGHTED_SIGN:
-        var_polarity_[var] = weighted_sign_[var] > 0;
-        break;
-      case SatParameters::POLARITY_REVERSE_WEIGHTED_SIGN:
-        var_polarity_[var] = weighted_sign_[var] < 0;
-        break;
     }
   }
 }
@@ -226,8 +212,7 @@ void SatDecisionPolicy::InitializeVariableOrdering() {
   for (BooleanVariable var(0); var < num_variables; ++var) {
     if (!trail_.Assignment().VariableIsAssigned(var)) {
       if (activities_[var] > 0.0) {
-        var_ordering_.Add(
-            {var, static_cast<float>(tie_breakers_[var]), activities_[var]});
+        var_ordering_.Add({var, tie_breakers_[var], activities_[var]});
       } else {
         tmp_variables_.push_back(var);
       }
@@ -253,7 +238,7 @@ void SatDecisionPolicy::InitializeVariableOrdering() {
 
   // Add the variables without activity to the queue (in the default order)
   for (const BooleanVariable var : tmp_variables_) {
-    var_ordering_.Add({var, static_cast<float>(tie_breakers_[var]), 0.0});
+    var_ordering_.Add({var, tie_breakers_[var], 0.0});
   }
 
   // Finish the queue initialization.
@@ -262,8 +247,7 @@ void SatDecisionPolicy::InitializeVariableOrdering() {
   var_ordering_is_initialized_ = true;
 }
 
-void SatDecisionPolicy::SetAssignmentPreference(Literal literal,
-                                                double weight) {
+void SatDecisionPolicy::SetAssignmentPreference(Literal literal, float weight) {
   if (!parameters_.use_optimization_hints()) return;
   DCHECK_GE(weight, 0.0);
   DCHECK_LE(weight, 1.0);
@@ -277,13 +261,13 @@ void SatDecisionPolicy::SetAssignmentPreference(Literal literal,
   var_ordering_is_initialized_ = false;
 }
 
-std::vector<std::pair<Literal, double>> SatDecisionPolicy::AllPreferences()
+std::vector<std::pair<Literal, float>> SatDecisionPolicy::AllPreferences()
     const {
-  std::vector<std::pair<Literal, double>> prefs;
+  std::vector<std::pair<Literal, float>> prefs;
   for (BooleanVariable var(0); var < var_polarity_.size(); ++var) {
     // TODO(user): we currently assume that if the tie_breaker is zero then
     // no preference was set (which is not 100% correct). Fix that.
-    const double value = var_ordering_.GetElement(var.value()).tie_breaker;
+    const float value = var_ordering_.GetElement(var.value()).tie_breaker;
     if (value > 0.0) {
       prefs.push_back(std::make_pair(Literal(var, var_polarity_[var]), value));
     }
@@ -291,16 +275,6 @@ std::vector<std::pair<Literal, double>> SatDecisionPolicy::AllPreferences()
   return prefs;
 }
 
-void SatDecisionPolicy::UpdateWeightedSign(
-    const std::vector<LiteralWithCoeff>& terms, Coefficient rhs) {
-  for (const LiteralWithCoeff& term : terms) {
-    const double weight = static_cast<double>(term.coefficient.value()) /
-                          static_cast<double>(rhs.value());
-    weighted_sign_[term.literal.Variable()] +=
-        term.literal.IsPositive() ? -weight : weight;
-  }
-}
-
 void SatDecisionPolicy::BumpVariableActivities(
     absl::Span<const Literal> literals) {
   if (parameters_.use_erwa_heuristic()) {
@@ -402,8 +376,8 @@ Literal SatDecisionPolicy::NextBranch() {
 }
 
 void SatDecisionPolicy::PqInsertOrUpdate(BooleanVariable var) {
-  const WeightedVarQueueElement element{
-      var, static_cast<float>(tie_breakers_[var]), activities_[var]};
+  const WeightedVarQueueElement element{var, tie_breakers_[var],
+                                        activities_[var]};
   if (var_ordering_.Contains(var.value())) {
     // Note that the new weight should always be higher than the old one.
     var_ordering_.IncreasePriority(element);
diff --git a/ortools/sat/sat_decision.h b/ortools/sat/sat_decision.h
index 7f8f1f756f..eed1d2f703 100644
--- a/ortools/sat/sat_decision.h
+++ b/ortools/sat/sat_decision.h
@@ -54,11 +54,6 @@ class SatDecisionPolicy {
   // variables are assigned.
   Literal NextBranch();
 
-  // Updates statistics about literal occurrences in constraints.
-  // Input is a canonical linear constraint of the form (terms <= rhs).
-  void UpdateWeightedSign(const std::vector<LiteralWithCoeff>& terms,
-                          Coefficient rhs);
-
   // Bumps the activity of all variables appearing in the conflict. All literals
   // must be currently assigned. See VSIDS decision heuristic: Chaff:
   // Engineering an Efficient SAT Solver. M.W. Moskewicz et al. ANNUAL ACM IEEE
@@ -100,10 +95,10 @@ class SatDecisionPolicy {
   //
   // Note(user): Having a lot of different weights may slow down the priority
   // queue operations if there is millions of variables.
-  void SetAssignmentPreference(Literal literal, double weight);
+  void SetAssignmentPreference(Literal literal, float weight);
 
   // Returns the vector of the current assignment preferences.
-  std::vector<std::pair<Literal, double>> AllPreferences() const;
+  std::vector<std::pair<Literal, float>> AllPreferences() const;
 
   // Returns the current activity of a BooleanVariable.
   double Activity(Literal l) const {
@@ -210,7 +205,7 @@ class SatDecisionPolicy {
   // Stores variable activity and the number of time each variable was "bumped".
   // The later is only used with the ERWA heuristic.
   absl::StrongVector<BooleanVariable, double> activities_;
-  absl::StrongVector<BooleanVariable, double> tie_breakers_;
+  absl::StrongVector<BooleanVariable, float> tie_breakers_;
   absl::StrongVector<BooleanVariable, int64_t> num_bumps_;
 
   // If the polarity if forced (externally) we always use this first.
@@ -234,9 +229,6 @@ class SatDecisionPolicy {
   // The longest partial assignment since the last reset.
   std::vector<Literal> best_partial_assignment_;
 
-  // Used in initial polarity computation.
-  absl::StrongVector<BooleanVariable, double> weighted_sign_;
-
   // Used in InitializeVariableOrdering().
   std::vector<BooleanVariable> tmp_variables_;
 };
diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto
index 46b4d5edb8..fdc44b1c51 100644
--- a/ortools/sat/sat_parameters.proto
+++ b/ortools/sat/sat_parameters.proto
@@ -23,7 +23,7 @@ option csharp_namespace = "Google.OrTools.Sat";
 // Contains the definitions for all the sat algorithm parameters and their
 // default values.
 //
-// NEXT TAG: 277
+// NEXT TAG: 278
 message SatParameters {
   // In some context, like in a portfolio of search, it makes sense to name a
   // given parameters set for logging purpose.
@@ -52,15 +52,6 @@ message SatParameters {
     POLARITY_TRUE = 0;
     POLARITY_FALSE = 1;
     POLARITY_RANDOM = 2;
-
-    // Choose the sign that tends to satisfy the most constraints. This is
-    // computed using a weighted sum: if a literal l appears in a constraint of
-    // the form: ... + coeff * l +... <= rhs with positive coefficients and
-    // rhs, then -sign(l) * coeff / rhs is added to the weight of l.variable().
-    POLARITY_WEIGHTED_SIGN = 3;
-
-    // The opposite choice of POLARITY_WEIGHTED_SIGN.
-    POLARITY_REVERSE_WEIGHTED_SIGN = 4;
   }
   optional Polarity initial_polarity = 2 [default = POLARITY_FALSE];
 
@@ -535,7 +526,7 @@ message SatParameters {
 
   // Enable or disable "inprocessing" which is some SAT presolving done at
   // each restart to the root level.
-  optional bool use_sat_inprocessing = 163 [default = false];
+  optional bool use_sat_inprocessing = 163 [default = true];
 
   // Proportion of deterministic time we should spend on inprocessing.
   // At each "restart", if the proportion is below this ratio, we will do some
@@ -837,123 +828,6 @@ message SatParameters {
   // as it modifies the reduced_cost, lb_tree_search, and probing workers.
   optional bool use_dual_scheduling_heuristics = 214 [default = true];
 
-  // A non-negative level indicating the type of constraints we consider in the
-  // LP relaxation. At level zero, no LP relaxation is used. At level 1, only
-  // the linear constraint and full encoding are added. At level 2, we also add
-  // all the Boolean constraints.
-  optional int32 linearization_level = 90 [default = 1];
-
-  // A non-negative level indicating how much we should try to fully encode
-  // Integer variables as Boolean.
-  optional int32 boolean_encoding_level = 107 [default = 1];
-
-  // When loading a*x + b*y ==/!= c when x and y are both fully encoded.
-  // The solver may decide to replace the linear equation by a set of clauses.
-  // This is triggered if the sizes of the domains of x and y are below the
-  // threshold.
-  optional int32 max_domain_size_when_encoding_eq_neq_constraints = 191
-      [default = 16];
-
-  // The limit on the number of cuts in our cut pool. When this is reached we do
-  // not generate cuts anymore.
-  //
-  // TODO(user): We should probably remove this parameters, and just always
-  // generate cuts but only keep the best n or something.
-  optional int32 max_num_cuts = 91 [default = 10000];
-
-  // Control the global cut effort. Zero will turn off all cut. For now we just
-  // have one level. Note also that most cuts are only used at linearization
-  // level >= 2.
-  optional int32 cut_level = 196 [default = 1];
-
-  // For the cut that can be generated at any level, this control if we only
-  // try to generate them at the root node.
-  optional bool only_add_cuts_at_level_zero = 92 [default = false];
-
-  // When the LP objective is fractional, do we add the cut that forces the
-  // linear objective expression to be greater or equal to this fractional value
-  // rounded up? We can always do that since our objective is integer, and
-  // combined with MIR heuristic to reduce the coefficient of such cut, it can
-  // help.
-  optional bool add_objective_cut = 197 [default = false];
-
-  // Whether we generate and add Chvatal-Gomory cuts to the LP at root node.
-  // Note that for now, this is not heavily tuned.
-  optional bool add_cg_cuts = 117 [default = true];
-
-  // Whether we generate MIR cuts at root node.
-  // Note that for now, this is not heavily tuned.
-  optional bool add_mir_cuts = 120 [default = true];
-
-  // Whether we generate Zero-Half cuts at root node.
-  // Note that for now, this is not heavily tuned.
-  optional bool add_zero_half_cuts = 169 [default = true];
-
-  // Whether we generate clique cuts from the binary implication graph. Note
-  // that as the search goes on, this graph will contains new binary clauses
-  // learned by the SAT engine.
-  optional bool add_clique_cuts = 172 [default = true];
-
-  // Cut generator for all diffs can add too many cuts for large all_diff
-  // constraints. This parameter restricts the large all_diff constraints to
-  // have a cut generator.
-  optional int32 max_all_diff_cut_size = 148 [default = 64];
-
-  // For the lin max constraints, generates the cuts described in "Strong
-  // mixed-integer programming formulations for trained neural networks" by Ross
-  // Anderson et. (https://arxiv.org/pdf/1811.01988.pdf)
-  optional bool add_lin_max_cuts = 152 [default = true];
-
-  // In the integer rounding procedure used for MIR and Gomory cut, the maximum
-  // "scaling" we use (must be positive). The lower this is, the lower the
-  // integer coefficients of the cut will be. Note that cut generated by lower
-  // values are not necessarily worse than cut generated by larger value. There
-  // is no strict dominance relationship.
-  //
-  // Setting this to 2 result in the "strong fractional rouding" of Letchford
-  // and Lodi.
-  optional int32 max_integer_rounding_scaling = 119 [default = 600];
-
-  // If true, we start by an empty LP, and only add constraints not satisfied
-  // by the current LP solution batch by batch. A constraint that is only added
-  // like this is known as a "lazy" constraint in the literature, except that we
-  // currently consider all constraints as lazy here.
-  optional bool add_lp_constraints_lazily = 112 [default = true];
-
-  // Even at the root node, we do not want to spend too much time on the LP if
-  // it is "difficult". So we solve it in "chunks" of that many iterations. The
-  // solve will be continued down in the tree or the next time we go back to the
-  // root node.
-  optional int32 root_lp_iterations = 227 [default = 2000];
-
-  // While adding constraints, skip the constraints which have orthogonality
-  // less than 'min_orthogonality_for_lp_constraints' with already added
-  // constraints during current call. Orthogonality is defined as 1 -
-  // cosine(vector angle between constraints). A value of zero disable this
-  // feature.
-  optional double min_orthogonality_for_lp_constraints = 115 [default = 0.05];
-
-  // Max number of time we perform cut generation and resolve the LP at level 0.
-  optional int32 max_cut_rounds_at_level_zero = 154 [default = 1];
-
-  // If a constraint/cut in LP is not active for that many consecutive OPTIMAL
-  // solves, remove it from the LP. Note that it might be added again later if
-  // it become violated by the current LP solution.
-  optional int32 max_consecutive_inactive_count = 121 [default = 100];
-
-  // These parameters are similar to sat clause management activity parameters.
-  // They are effective only if the number of generated cuts exceed the storage
-  // limit. Default values are based on a few experiments on miplib instances.
-  optional double cut_max_active_count_value = 155 [default = 1e10];
-  optional double cut_active_count_decay = 156 [default = 0.8];
-
-  // Target number of constraints to remove during cleanup.
-  optional int32 cut_cleanup_target = 157 [default = 1000];
-
-  // Add that many lazy constraints (or cuts) at once in the LP. Note that at
-  // the beginning of the solve, we do add more than this.
-  optional int32 new_constraints_batch_size = 122 [default = 50];
-
   // The search branching will be used to decide how to branch on unfixed nodes.
   enum SearchBranching {
     // Try to fix all literals using the underlying SAT solver's heuristics,
@@ -1013,37 +887,6 @@ message SatParameters {
   // hinted value.
   optional bool fix_variables_to_their_hinted_value = 192 [default = false];
 
-  // All the "exploit_*" parameters below work in the same way: when branching
-  // on an IntegerVariable, these parameters affect the value the variable is
-  // branched on. Currently the first heuristic that triggers win in the order
-  // in which they appear below.
-  //
-  // TODO(user): Maybe do like for the restart algorithm, introduce an enum
-  // and a repeated field that control the order on which these are applied?
-
-  // If true and the Lp relaxation of the problem has an integer optimal
-  // solution, try to exploit it. Note that since the LP relaxation may not
-  // contain all the constraints, such a solution is not necessarily a solution
-  // of the full problem.
-  optional bool exploit_integer_lp_solution = 94 [default = true];
-
-  // If true and the Lp relaxation of the problem has a solution, try to exploit
-  // it. This is same as above except in this case the lp solution might not be
-  // an integer solution.
-  optional bool exploit_all_lp_solution = 116 [default = true];
-
-  // When branching on a variable, follow the last best solution value.
-  optional bool exploit_best_solution = 130 [default = false];
-
-  // When branching on a variable, follow the last best relaxation solution
-  // value. We use the relaxation with the tightest bound on the objective as
-  // the best relaxation solution.
-  optional bool exploit_relaxation_solution = 161 [default = false];
-
-  // When branching an a variable that directly affect the objective,
-  // branch on the value that lead to the best objective first.
-  optional bool exploit_objective = 131 [default = true];
-
   // If true, search will continuously probe Boolean variables, and integer
   // variable bounds. This parameter is set to true in parallel on the probing
   // worker.
@@ -1341,7 +1184,7 @@ message SatParameters {
 
   // All at_most_one constraints with a size <= param will be replaced by a
   // quadratic number of binary implications.
-  optional int32 at_most_one_max_expansion_size = 270 [default = 9];
+  optional int32 at_most_one_max_expansion_size = 270 [default = 3];
 
   // Indicates if the CP-SAT layer should catch Control-C (SIGINT) signals
   // when calling solve. If set, catching the SIGINT signal will terminate the
@@ -1386,6 +1229,163 @@ message SatParameters {
   // faster propation in the CP engine.
   optional int32 linear_split_size = 256 [default = 100];
 
+  // ==========================================================================
+  // Linear programming relaxation
+  // ==========================================================================
+
+  // A non-negative level indicating the type of constraints we consider in the
+  // LP relaxation. At level zero, no LP relaxation is used. At level 1, only
+  // the linear constraint and full encoding are added. At level 2, we also add
+  // all the Boolean constraints.
+  optional int32 linearization_level = 90 [default = 1];
+
+  // A non-negative level indicating how much we should try to fully encode
+  // Integer variables as Boolean.
+  optional int32 boolean_encoding_level = 107 [default = 1];
+
+  // When loading a*x + b*y ==/!= c when x and y are both fully encoded.
+  // The solver may decide to replace the linear equation by a set of clauses.
+  // This is triggered if the sizes of the domains of x and y are below the
+  // threshold.
+  optional int32 max_domain_size_when_encoding_eq_neq_constraints = 191
+      [default = 16];
+
+  // The limit on the number of cuts in our cut pool. When this is reached we do
+  // not generate cuts anymore.
+  //
+  // TODO(user): We should probably remove this parameters, and just always
+  // generate cuts but only keep the best n or something.
+  optional int32 max_num_cuts = 91 [default = 10000];
+
+  // Control the global cut effort. Zero will turn off all cut. For now we just
+  // have one level. Note also that most cuts are only used at linearization
+  // level >= 2.
+  optional int32 cut_level = 196 [default = 1];
+
+  // For the cut that can be generated at any level, this control if we only
+  // try to generate them at the root node.
+  optional bool only_add_cuts_at_level_zero = 92 [default = false];
+
+  // When the LP objective is fractional, do we add the cut that forces the
+  // linear objective expression to be greater or equal to this fractional value
+  // rounded up? We can always do that since our objective is integer, and
+  // combined with MIR heuristic to reduce the coefficient of such cut, it can
+  // help.
+  optional bool add_objective_cut = 197 [default = false];
+
+  // Whether we generate and add Chvatal-Gomory cuts to the LP at root node.
+  // Note that for now, this is not heavily tuned.
+  optional bool add_cg_cuts = 117 [default = true];
+
+  // Whether we generate MIR cuts at root node.
+  // Note that for now, this is not heavily tuned.
+  optional bool add_mir_cuts = 120 [default = true];
+
+  // Whether we generate Zero-Half cuts at root node.
+  // Note that for now, this is not heavily tuned.
+  optional bool add_zero_half_cuts = 169 [default = true];
+
+  // Whether we generate clique cuts from the binary implication graph. Note
+  // that as the search goes on, this graph will contains new binary clauses
+  // learned by the SAT engine.
+  optional bool add_clique_cuts = 172 [default = true];
+
+  // Cut generator for all diffs can add too many cuts for large all_diff
+  // constraints. This parameter restricts the large all_diff constraints to
+  // have a cut generator.
+  optional int32 max_all_diff_cut_size = 148 [default = 64];
+
+  // For the lin max constraints, generates the cuts described in "Strong
+  // mixed-integer programming formulations for trained neural networks" by Ross
+  // Anderson et. (https://arxiv.org/pdf/1811.01988.pdf)
+  optional bool add_lin_max_cuts = 152 [default = true];
+
+  // In the integer rounding procedure used for MIR and Gomory cut, the maximum
+  // "scaling" we use (must be positive). The lower this is, the lower the
+  // integer coefficients of the cut will be. Note that cut generated by lower
+  // values are not necessarily worse than cut generated by larger value. There
+  // is no strict dominance relationship.
+  //
+  // Setting this to 2 result in the "strong fractional rouding" of Letchford
+  // and Lodi.
+  optional int32 max_integer_rounding_scaling = 119 [default = 600];
+
+  // If true, we start by an empty LP, and only add constraints not satisfied
+  // by the current LP solution batch by batch. A constraint that is only added
+  // like this is known as a "lazy" constraint in the literature, except that we
+  // currently consider all constraints as lazy here.
+  optional bool add_lp_constraints_lazily = 112 [default = true];
+
+  // Even at the root node, we do not want to spend too much time on the LP if
+  // it is "difficult". So we solve it in "chunks" of that many iterations. The
+  // solve will be continued down in the tree or the next time we go back to the
+  // root node.
+  optional int32 root_lp_iterations = 227 [default = 2000];
+
+  // While adding constraints, skip the constraints which have orthogonality
+  // less than 'min_orthogonality_for_lp_constraints' with already added
+  // constraints during current call. Orthogonality is defined as 1 -
+  // cosine(vector angle between constraints). A value of zero disable this
+  // feature.
+  optional double min_orthogonality_for_lp_constraints = 115 [default = 0.05];
+
+  // Max number of time we perform cut generation and resolve the LP at level 0.
+  optional int32 max_cut_rounds_at_level_zero = 154 [default = 1];
+
+  // If a constraint/cut in LP is not active for that many consecutive OPTIMAL
+  // solves, remove it from the LP. Note that it might be added again later if
+  // it become violated by the current LP solution.
+  optional int32 max_consecutive_inactive_count = 121 [default = 100];
+
+  // These parameters are similar to sat clause management activity parameters.
+  // They are effective only if the number of generated cuts exceed the storage
+  // limit. Default values are based on a few experiments on miplib instances.
+  optional double cut_max_active_count_value = 155 [default = 1e10];
+  optional double cut_active_count_decay = 156 [default = 0.8];
+
+  // Target number of constraints to remove during cleanup.
+  optional int32 cut_cleanup_target = 157 [default = 1000];
+
+  // Add that many lazy constraints (or cuts) at once in the LP. Note that at
+  // the beginning of the solve, we do add more than this.
+  optional int32 new_constraints_batch_size = 122 [default = 50];
+
+  // All the "exploit_*" parameters below work in the same way: when branching
+  // on an IntegerVariable, these parameters affect the value the variable is
+  // branched on. Currently the first heuristic that triggers win in the order
+  // in which they appear below.
+  //
+  // TODO(user): Maybe do like for the restart algorithm, introduce an enum
+  // and a repeated field that control the order on which these are applied?
+
+  // If true and the Lp relaxation of the problem has an integer optimal
+  // solution, try to exploit it. Note that since the LP relaxation may not
+  // contain all the constraints, such a solution is not necessarily a solution
+  // of the full problem.
+  optional bool exploit_integer_lp_solution = 94 [default = true];
+
+  // If true and the Lp relaxation of the problem has a solution, try to exploit
+  // it. This is same as above except in this case the lp solution might not be
+  // an integer solution.
+  optional bool exploit_all_lp_solution = 116 [default = true];
+
+  // When branching on a variable, follow the last best solution value.
+  optional bool exploit_best_solution = 130 [default = false];
+
+  // When branching on a variable, follow the last best relaxation solution
+  // value. We use the relaxation with the tightest bound on the objective as
+  // the best relaxation solution.
+  optional bool exploit_relaxation_solution = 161 [default = false];
+
+  // When branching an a variable that directly affect the objective,
+  // branch on the value that lead to the best objective first.
+  optional bool exploit_objective = 131 [default = true];
+
+  // Infer products of Boolean or of Boolean time IntegerVariable from the
+  // linear constrainst in the problem. This can be used in some cuts, altough
+  // for now we don't really exploit it.
+  optional bool detect_linearized_product = 277 [default = false];
+
   // ==========================================================================
   // MIP -> CP-SAT (i.e. IP with integer coeff) conversion parameters that are
   // used by our automatic "scaling" algorithm.
diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc
index f0ecc0fec2..db42e94425 100644
--- a/ortools/sat/sat_solver.cc
+++ b/ortools/sat/sat_solver.cc
@@ -289,10 +289,6 @@ bool SatSolver::AddLinearConstraintInternal(
   if (rhs < 0) return SetModelUnsat();  // Unsatisfiable constraint.
   if (rhs >= max_value) return true;    // Always satisfied constraint.
 
-  // The case "rhs = 0" will just fix variables, so there is no need to
-  // updates the weighted sign.
-  if (rhs > 0) decision_policy_->UpdateWeightedSign(cst, rhs);
-
   // Since the constraint is in canonical form, the coefficients are sorted.
   const Coefficient min_coeff = cst.front().coefficient;
   const Coefficient max_coeff = cst.back().coefficient;
diff --git a/ortools/sat/sat_solver.h b/ortools/sat/sat_solver.h
index 923edee5a1..b15d0b481c 100644
--- a/ortools/sat/sat_solver.h
+++ b/ortools/sat/sat_solver.h
@@ -161,22 +161,15 @@ class SatSolver {
   //
   // TODO(user): Clean this up by making clients directly talk to
   // SatDecisionPolicy.
-  void SetAssignmentPreference(Literal literal, double weight) {
+  void SetAssignmentPreference(Literal literal, float weight) {
     decision_policy_->SetAssignmentPreference(literal, weight);
   }
-  std::vector<std::pair<Literal, double>> AllPreferences() const {
+  std::vector<std::pair<Literal, float>> AllPreferences() const {
     return decision_policy_->AllPreferences();
   }
   void ResetDecisionHeuristic() {
     return decision_policy_->ResetDecisionHeuristic();
   }
-  void ResetDecisionHeuristicAndSetAllPreferences(
-      const std::vector<std::pair<Literal, double>>& prefs) {
-    decision_policy_->ResetDecisionHeuristic();
-    for (const std::pair<Literal, double>& p : prefs) {
-      decision_policy_->SetAssignmentPreference(p.first, p.second);
-    }
-  }
 
   // Solves the problem and returns its status.
   // An empty problem is considered to be SAT.
diff --git a/ortools/sat/util.h b/ortools/sat/util.h
index 4c3f51dcc2..b722b1dc90 100644
--- a/ortools/sat/util.h
+++ b/ortools/sat/util.h
@@ -573,6 +573,7 @@ class TopN {
   bool empty() const { return elements_.empty(); }
 
   const std::vector<Element>& UnorderedElements() const { return elements_; }
+  std::vector<Element>* MutableUnorderedElements() { return &elements_; }
 
  private:
   const int n_;