Fix knapsack rational coefficient scaling

cpp/src/cuts/cuts.cpp

@@ -2651,6 +2651,7 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
   f_t objective_constant = 0.0;
   std::vector<i_t> fixed_variables;
   std::vector<f_t> fixed_values;
+  std::vector<f_t> fixed_weights;
   const f_t x_tol = 1e-5;
   for (i_t k = 0; k < knapsack_inequality.size(); k++) {
     const i_t j = knapsack_inequality.index(k);
@@ -2662,6 +2663,7 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
         // if xstar_j is close to 0, then we can fix z to zero
         fixed_variables.push_back(j);
         fixed_values.push_back(0.0);
+        fixed_weights.push_back(knapsack_inequality.vector.x[k]);
         seperation_rhs -= knapsack_inequality.vector.x[k];
         // No need to adjust the objective constant
         continue;
@@ -2670,6 +2672,7 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
         // if xstar_j is close to 1, then we can fix z to 1
         fixed_variables.push_back(j);
         fixed_values.push_back(1.0);
+        fixed_weights.push_back(knapsack_inequality.vector.x[k]);
         // Note seperation rhs is unchanged
         objective_constant += vj;
         continue;
@@ -2711,12 +2714,28 @@ i_t knapsack_generation_t<i_t, f_t>::generate_knapsack_cut(
     return -1;
   }
 
-  i_t cover_size = 0;
+  i_t cover_size   = 0;
+  f_t cover_weight = 0.0;
   for (i_t k = 0; k < solution.size(); k++) {
-    if (solution[k] == 0.0) { cover_size++; }
+    if (solution[k] == 0.0) {
+      cover_size++;
+      cover_weight += weights[k];
+    }
   }
   for (i_t k = 0; k < fixed_values.size(); k++) {
-    if (fixed_values[k] == 1.0) { cover_size++; }
+    if (fixed_values[k] == 1.0) {
+      cover_size++;
+      cover_weight += fixed_weights[k];
+    }
+  }
+
+  // sum_{j in C} a_j > beta is what makes sum_{j in C} x_j <= |C| - 1 valid. The coefficients are
+  // integral here, so demand a full unit rather than letting rounding in the sums decide.
+  const bool is_cover = cover_weight >= knapsack_inequality.rhs + 1.0 - tol;
+  cuopt_assert(is_cover, "knapsack separation produced a set that is not a cover");
+  if (!is_cover) {
+    restore_complemented(complemented_variables);
+    return -1;
   }
 
   cut.reserve(cover_size);
@@ -2887,6 +2906,8 @@ void knapsack_generation_t<i_t, f_t>::minimal_cover_and_partition(
     }
   }
 
+  cuopt_assert(cover_sum > beta, "minimal cover reduction dropped an item the cover needed");
+
   // Go through and correct cover_indicies and cover_coefficients
   for (i_t k = 0; k < cover_coefficients.size();) {
     if (cover_coefficients[k] == 0.0) {
@@ -3220,7 +3241,10 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_
 
   solution.assign(n, 0.0);
 
-  // dp(j, v) = minimum weight using first j items to get value v
+  // dp(j, v) = minimum weight using first j items to get value v.
+  // The weights are rounded up: rounding one down would let the dynamic programming algorithm
+  // return a set that violates the capacity, and the caller reads the complement of that set
+  // as a cover.
   dense_matrix_t<i_t, i_t> dp(n + 1, sum_value + 1, INT_INF);
   dense_matrix_t<i_t, uint8_t> take(n + 1, sum_value + 1, 0);
   dp(0, 0) = 0;
@@ -3233,8 +3257,7 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_
 
       // Take item j-1 if possible
       if (v >= scaled_values[j - 1]) {
-        i_t candidate =
-          dp(j - 1, v - scaled_values[j - 1]) + static_cast<i_t>(std::floor(weights[j - 1]));
+        i_t candidate = dp(j - 1, v - scaled_values[j - 1]) + std::ceil(weights[j - 1]);
         if (candidate < dp(j, v)) {
           dp(j, v)   = candidate;
           take(j, v) = 1;
@@ -3260,6 +3283,15 @@ f_t knapsack_generation_t<i_t, f_t>::solve_knapsack_problem(const std::vector<f_
     }
   }
 
+#ifdef ASSERT_MODE
+  f_t selected_weight = 0.0;
+  for (i_t j = 0; j < n; ++j) {
+    selected_weight += solution[j] * weights[j];
+  }
+  cuopt_assert(selected_weight <= rhs + settings_.primal_tol,
+               "knapsack dynamic program returned a solution over capacity");
+#endif
+
   objective = best_value * scale;
   return objective;
 }
@@ -4918,6 +4950,7 @@ bool rational_coefficients(const std::vector<variable_type_t>& var_types,
   if (scalar < 0) { return false; }
   if (std::abs(scalar) > 1000) { return false; }
 
+  // The scaled product can land an ulp off the integer it represents.
   rational_inequality.scale(scalar);
 
   return true;

cpp/tests/mip/cuts_test.cu

@@ -444,6 +444,19 @@ void disable_all_cuts(mip_solver_settings_t<int, double>& settings)
   settings.strong_chvatal_gomory_cuts = 0;
 }
 
+void disable_non_knapsack_cuts(mip_solver_settings_t<int, double>& settings)
+{
+  settings.max_cut_passes             = 10;
+  settings.knapsack_cuts              = 1;
+  settings.clique_cuts                = 0;
+  settings.zero_half_cuts             = 0;
+  settings.mixed_integer_gomory_cuts  = 0;
+  settings.mir_cuts                   = 0;
+  settings.strong_chvatal_gomory_cuts = 0;
+  settings.flow_cover_cuts            = 0;
+  settings.implied_bound_cuts         = 0;
+}
+
 bool cut_is_invalid_for_incumbent(const std::vector<int>& cut_vars,
                                   const std::vector<double>& incumbent,
                                   double tol)
@@ -936,6 +949,62 @@ TEST(cuts, test_cuts_2)
   EXPECT_EQ(solution.get_num_nodes(), 0);
 }
 
+io::mps_data_model_t<int, double> create_knapsack_cover_floor_problem()
+{
+  // The odd cycle over z1, z2, z3 makes z = (0.5, 0.5, 0.5), w = 0 the unique LP optimum, which
+  // puts y at 0.412078. Integrality moves the optimum to w = 1 with y = 0 and objective 3.
+  //
+  // The capacity coefficients are load bearing: scaling that row to integers multiplies by 100,
+  // and 135.45 and 135.42 do not land on integers when scaled, which is what the knapsack
+  // separator needs them to do.
+  return cuopt::test::parse_inline_lp(R"LP(
+Minimize
+  obj: 2 y + z1 + z2 + z3 + 3 w
+Subject To
+  capacity: -450 y + 135.45 z1 + 135.42 z2 + 100 z3 <= 0
+  tri12: z1 + z2 + w >= 1
+  tri13: z1 + z3 + w >= 1
+  tri23: z2 + z3 + w >= 1
+Binaries
+  y
+  z1
+  z2
+  z3
+  w
+End
+)LP");
+}
+
+TEST(cuts, knapsack_cover_floor_regression)
+{
+  const raft::handle_t handle_{};
+  auto problem = create_knapsack_cover_floor_problem();
+
+  mip_solver_settings_t<int, double> settings;
+  settings.time_limit = 10.;
+  disable_non_knapsack_cuts(settings);
+  settings.presolver = presolver_t::None;
+
+  mip_solution_t<int, double> solution = solve_mip(&handle_, problem, settings);
+  EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal);
+  EXPECT_NEAR(3.0, solution.get_objective_value(), 1e-6);
+}
+
+TEST(cuts, knapsack_cover_floor_regression_reference)
+{
+  const raft::handle_t handle_{};
+  auto problem = create_knapsack_cover_floor_problem();
+
+  mip_solver_settings_t<int, double> settings;
+  settings.time_limit = 10.;
+  disable_all_cuts(settings);
+  settings.presolver = presolver_t::None;
+
+  mip_solution_t<int, double> solution = solve_mip(&handle_, problem, settings);
+  EXPECT_EQ(solution.get_termination_status(), mip_termination_status_t::Optimal);
+  EXPECT_NEAR(3.0, solution.get_objective_value(), 1e-6);
+}
+
 TEST(cuts, test_duplicate_cuts_detection)
 {
   simplex::simplex_solver_settings_t<int, double> settings;