fix: incorrect GF2 presolve constraint addressing

cpp/src/mip_heuristics/presolve/gf2_presolve.cpp

@@ -203,11 +203,12 @@ papilo::PresolveStatus GF2Presolve<f_t>::execute(const papilo::Problem<f_t>& pro
   std::vector<std::vector<int>> A(gf2_constraints.size(),
                                   std::vector<int>(gf2_constraints.size(), 0));
   std::vector<int> b(gf2_constraints.size());
-  for (const auto& cons : gf2_constraints) {
+  for (size_t gf2_cstr_idx = 0; gf2_cstr_idx < gf2_constraints.size(); ++gf2_cstr_idx) {
+    const auto& cons = gf2_constraints[gf2_cstr_idx];
     for (auto [bin_var, _] : cons.bin_vars) {
-      A[cons.cstr_idx][gf2_bin_vars[bin_var]] = 1;
+      A[gf2_cstr_idx][gf2_bin_vars[bin_var]] = 1;
     }
-    b[cons.cstr_idx] = cons.rhs;
+    b[gf2_cstr_idx] = cons.rhs;
   }
 
   std::vector<int> solution(gf2_constraints.size());

cpp/tests/mip/presolve_test.cu

@@ -22,6 +22,8 @@
 
 #include <gtest/gtest.h>
 
+#include <algorithm>
+#include <array>
 #include <cstdint>
 #include <sstream>
 #include <string>
@@ -63,4 +65,89 @@ TEST(problem, find_implied_integers)
             ((int)mip::problem_t<int, double>::var_flags_t::VAR_IMPLIED_INTEGER));
 }
 
+TEST(gf2_presolve, uses_compact_constraint_indices)
+{
+  constexpr int num_packing_vars = 128;
+  constexpr int num_gf2_vars     = 6;
+  constexpr int num_key_vars     = 6;
+  constexpr int num_packing_rows = 128;
+  constexpr int num_key_rows     = 2 * num_key_vars;
+  constexpr int num_gf2_rows     = 6;
+  constexpr int num_vars         = num_packing_vars + num_gf2_vars + num_key_vars;
+  constexpr int num_rows         = num_packing_rows + num_key_rows + num_gf2_rows;
+  constexpr int x_offset         = num_packing_vars;
+  constexpr int y_offset         = x_offset + num_gf2_vars;
+
+  std::vector<double> values;
+  std::vector<int> indices;
+  std::vector<int> offsets{0};
+  std::vector<double> constraint_lb(num_rows, 1.0);
+  std::vector<double> constraint_ub(num_rows, 2.0);
+
+  auto add_entry = [&](int column, double value) {
+    indices.push_back(column);
+    values.push_back(value);
+  };
+  auto finish_row = [&] { offsets.push_back(static_cast<int>(values.size())); };
+
+  // A normal binary MIP block keeps the GF2 rows at high raw row indices.
+  for (int row = 0; row < num_packing_rows; ++row) {
+    std::array columns{row, (row + 1) % num_packing_vars, (row + 2) % num_packing_vars};
+    std::sort(columns.begin(), columns.end());
+    for (int column : columns) {
+      add_entry(column, 1.0);
+    }
+    finish_row();
+  }
+
+  // Keep every GF2 key column non-singleton without forcing it.
+  for (int key = 0; key < num_key_vars; ++key) {
+    add_entry(3 * key, 1.0);
+    add_entry(3 * key + 1, 1.0);
+    add_entry(y_offset + key, 1.0);
+    finish_row();
+
+    add_entry(3 * key + 1, 1.0);
+    add_entry(3 * key + 2, 1.0);
+    add_entry(y_offset + key, 1.0);
+    finish_row();
+  }
+
+  // Over GF(2), this is J-I for even dimension 6, hence nonsingular. Three positive and two
+  // negative coefficients per row prevent ordinary bound propagation from fixing the key.
+  for (int row = 0; row < num_gf2_rows; ++row) {
+    int term = 0;
+    for (int col = 0; col < num_gf2_vars; ++col) {
+      if (col == row) { continue; }
+      add_entry(x_offset + col, term < 3 ? 1.0 : -1.0);
+      ++term;
+    }
+    add_entry(y_offset + row, 2.0);
+    finish_row();
+    constraint_lb[num_packing_rows + num_key_rows + row] = 1.0;
+    constraint_ub[num_packing_rows + num_key_rows + row] = 1.0;
+  }
+
+  const raft::handle_t handle_{};
+  optimization_problem_t<int, double> problem(&handle_);
+  std::vector<double> objective(num_vars, 1.0);
+  std::vector<double> variable_lb(num_vars, 0.0);
+  std::vector<double> variable_ub(num_vars, 1.0);
+  std::vector<var_t> variable_types(num_vars, var_t::INTEGER);
+  problem.set_csr_constraint_matrix(
+    values.data(), values.size(), indices.data(), indices.size(), offsets.data(), offsets.size());
+  problem.set_objective_coefficients(objective.data(), objective.size());
+  problem.set_variable_lower_bounds(variable_lb.data(), variable_lb.size());
+  problem.set_variable_upper_bounds(variable_ub.data(), variable_ub.size());
+  problem.set_variable_types(variable_types.data(), variable_types.size());
+  problem.set_constraint_lower_bounds(constraint_lb.data(), constraint_lb.size());
+  problem.set_constraint_upper_bounds(constraint_ub.data(), constraint_ub.size());
+
+  auto presolver = std::make_unique<mip::third_party_presolve_t<int, double>>();
+  auto result    = presolver->apply(
+    problem, problem_category_t::MIP, presolver_t::Papilo, false, 1e-6, 1e-12, 20, 1);
+
+  EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED);
+}
+
 }  // namespace cuopt::mathematical_optimization::test