diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index be45ffeecd..ae6d71e78d 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -2406,6 +2406,7 @@ i_t knapsack_generation_t::generate_knapsack_cut( f_t objective_constant = 0.0; std::vector fixed_variables; std::vector fixed_values; + std::vector 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); @@ -2417,6 +2418,7 @@ i_t knapsack_generation_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; @@ -2425,6 +2427,7 @@ i_t knapsack_generation_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; @@ -2466,12 +2469,28 @@ i_t knapsack_generation_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); @@ -2641,6 +2660,9 @@ void knapsack_generation_t::minimal_cover_and_partition( } } + cuopt_assert(cover_sum >= beta + 1.0 - 1e-6, + "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) { @@ -2958,8 +2980,7 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector::max() / 2; + i_t sum_value = std::accumulate(scaled_values.begin(), scaled_values.end(), 0); if (verbose) { settings_.log.printf("sum value %d\n", sum_value); } const i_t max_size = 10000; if (sum_value <= 0.0 || sum_value >= max_size) { @@ -2972,10 +2993,12 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector dp(n + 1, sum_value + 1, INT_INF); + // dp(j, v) = minimum weight using first j items to get value v. + // The weights are carried at full precision: rounding one down would let the DP return a set + // that violates the capacity, and the caller reads the complement of that set as a cover. + dense_matrix_t dp(n + 1, sum_value + 1, inf); dense_matrix_t take(n + 1, sum_value + 1, 0); - dp(0, 0) = 0; + dp(0, 0) = 0.0; // 4. Dynamic programming for (i_t j = 1; j <= n; ++j) { @@ -2985,8 +3008,7 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector= scaled_values[j - 1]) { - i_t candidate = - dp(j - 1, v - scaled_values[j - 1]) + static_cast(std::floor(weights[j - 1])); + f_t candidate = dp(j - 1, v - scaled_values[j - 1]) + weights[j - 1]; if (candidate < dp(j, v)) { dp(j, v) = candidate; take(j, v) = 1; @@ -3012,6 +3034,15 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector& var_types, rational_inequality.scale(scalar); + // The scaled product can land an ulp off the integer it represents. Callers rely on the + // integer-variable coefficients being exact integers: the knapsack cover test + // sum_C a_j > beta is only equivalent to sum_C a_j >= beta + 1 for integral a_j. + constexpr f_t integral_tol = 1e-6; + for (i_t k : indices) { + const f_t scaled = rational_inequality.vector.x[k]; + const f_t rounded = std::round(scaled); + if (std::abs(scaled - rounded) > integral_tol) { return false; } + rational_inequality.vector.x[k] = rounded; + } + return true; } diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index b4fc3e8cc7..88c553bcf1 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -442,6 +442,19 @@ void disable_all_cuts(mip_solver_settings_t& settings) settings.strong_chvatal_gomory_cuts = 0; } +void disable_non_knapsack_cuts(mip_solver_settings_t& 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& cut_vars, const std::vector& incumbent, double tol) @@ -934,6 +947,62 @@ TEST(cuts, test_cuts_2) EXPECT_EQ(solution.get_num_nodes(), 0); } +io::mps_data_model_t 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 settings; + settings.time_limit = 10.; + disable_non_knapsack_cuts(settings); + settings.presolver = presolver_t::None; + + mip_solution_t 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 settings; + settings.time_limit = 10.; + disable_all_cuts(settings); + settings.presolver = presolver_t::None; + + mip_solution_t 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 settings;