diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ba2b63983..1345cb3e7 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -243,9 +243,11 @@ branch_and_bound_t::branch_and_bound_t( const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, f_t start_time, + const probing_implied_bounds_t& probing_implied_bounds, std::shared_ptr> clique_table) : original_problem_(user_problem), settings_(solver_settings), + probing_implied_bounds_(probing_implied_bounds), clique_table_(std::move(clique_table)), original_lp_(user_problem.handle_ptr, 1, 1, 1), Arow_(1, 1, 0), @@ -724,12 +726,6 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& obj, is_maximization ? "Upper" : "Lower", user_bound); - { - const f_t root_lp_obj = root_lp_current_lower_bound_.load(); - if (std::isfinite(root_lp_obj)) { - settings_.log.printf("Root LP dual objective (last): %.16e\n", root_lp_obj); - } - } if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { solver_status_ = mip_status_t::OPTIMAL; @@ -2150,6 +2146,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut new_slacks_, var_types_, original_problem_, + probing_implied_bounds_, clique_table_, &clique_table_future_, &signal_extend_cliques_); @@ -2220,7 +2217,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); if (num_cuts == 0) { break; } cut_info.record_cut_types(cut_types); -#ifdef PRINT_CUT_POOL_TYPES +#if 1 cut_pool.print_cutpool_types(); print_cut_types("In LP ", cut_types, settings_); printf("Cut pool size: %d\n", cut_pool.pool_size()); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 98674b7f9..3c5ee6538 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -77,6 +77,7 @@ class branch_and_bound_t { branch_and_bound_t(const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, f_t start_time, + const probing_implied_bounds_t& probing_implied_bounds, std::shared_ptr> clique_table = nullptr); // Set an initial guess based on the user_problem. This should be called before solve. @@ -148,6 +149,7 @@ class branch_and_bound_t { private: const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; + const probing_implied_bounds_t& probing_implied_bounds_; std::shared_ptr> clique_table_; std::future>> clique_table_future_; std::atomic signal_extend_cliques_{false}; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index c74a12b36..c00341ed3 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -445,6 +445,21 @@ void extend_clique_vertices(std::vector& clique_vertices, } // namespace +template +bool rational_coefficients(const std::vector& var_types, + const inequality_t& inequality, + inequality_t& rational_inequality); + +template +bool rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator); + +int64_t gcd(const std::vector& integers); + +int64_t lcm(const std::vector& integers); + // This function is only used in tests std::vector> find_maximal_cliques_for_test( const std::vector>& adjacency_list, @@ -680,69 +695,76 @@ knapsack_generation_t::knapsack_generation_t( csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types) - : settings_(settings) + : is_slack_(lp.num_cols, 0), + is_complemented_(lp.num_cols, 0), + is_marked_(lp.num_cols, 0), + workspace_(lp.num_cols, 0.0), + complemented_xstar_(lp.num_cols, 0.0), + settings_(settings) { const bool verbose = false; knapsack_constraints_.reserve(lp.num_rows); - is_slack_.resize(lp.num_cols, 0); for (i_t j : new_slacks) { is_slack_[j] = 1; } for (i_t i = 0; i < lp.num_rows; i++) { - const i_t row_start = Arow.row_start[i]; - const i_t row_end = Arow.row_start[i + 1]; - const i_t row_len = row_end - row_start; + inequality_t inequality(Arow, i, lp.rhs[i]); + inequality_t rational_inequality = inequality; + if (!rational_coefficients(var_types, inequality, rational_inequality)) { + continue; + } + inequality = rational_inequality; + + const i_t row_len = rational_inequality.size(); if (row_len < 3) { continue; } bool is_knapsack = true; f_t sum_pos = 0.0; - for (i_t p = row_start; p < row_end; p++) { - const i_t j = Arow.j[p]; + f_t sum_neg = 0.0; + for (i_t p = 0; p < row_len; p++) { + const i_t j = inequality.index(p); if (is_slack_[j]) { continue; } - const f_t aj = Arow.x[p]; - if (std::abs(aj - std::round(aj)) > settings.integer_tol) { - is_knapsack = false; - break; - } + const f_t aj = inequality.coeff(p); if (var_types[j] != variable_type_t::INTEGER || lp.lower[j] != 0.0 || lp.upper[j] != 1.0) { is_knapsack = false; break; } if (aj < 0.0) { - is_knapsack = false; - break; + sum_pos += -aj; + sum_neg += -aj; + } else { + sum_pos += aj; } - sum_pos += aj; } if (is_knapsack) { - const f_t beta = lp.rhs[i]; - if (std::abs(beta - std::round(beta)) <= settings.integer_tol) { - if (beta > 0.0 && beta <= sum_pos && std::abs(sum_pos / (row_len - 1) - beta) > 1e-3) { - if (verbose) { - settings.log.printf( - "Knapsack constraint %d row len %d beta %e sum_pos %e sum_pos / (row_len - 1) %e\n", - i, - row_len, - beta, - sum_pos, - sum_pos / (row_len - 1)); - } - knapsack_constraints_.push_back(i); + const f_t beta = inequality.rhs + sum_neg; + if (beta > 0.0 && beta <= sum_pos && std::abs(sum_pos / (row_len - 1) - beta) > 1e-3) { + if (verbose) { + settings.log.printf( + "Knapsack constraint %d row len %d beta %e sum_neg %e sum_pos %e sum_pos / (row_len - " + "1) %e\n", + i, + row_len, + beta, + sum_neg, + sum_pos, + sum_pos / (row_len - 1)); } + knapsack_constraints_.push_back(i); } } } -#ifdef PRINT_KNAPSACK_INFO +#if 1 i_t num_knapsack_constraints = knapsack_constraints_.size(); settings.log.printf("Number of knapsack constraints %d\n", num_knapsack_constraints); #endif } template -i_t knapsack_generation_t::generate_knapsack_cuts( +i_t knapsack_generation_t::generate_knapsack_cut( const lp_problem_t& lp, const simplex_solver_settings_t& settings, csr_matrix_t& Arow, @@ -754,28 +776,61 @@ i_t knapsack_generation_t::generate_knapsack_cuts( { const bool verbose = false; // Get the row associated with the knapsack constraint - sparse_vector_t knapsack_inequality(Arow, knapsack_row); - f_t knapsack_rhs = lp.rhs[knapsack_row]; + inequality_t knapsack_inequality(Arow, knapsack_row, lp.rhs[knapsack_row]); + inequality_t rational_knapsack_inequality = knapsack_inequality; + if (!rational_coefficients(var_types, knapsack_inequality, rational_knapsack_inequality)) { + return -1; + } + knapsack_inequality = rational_knapsack_inequality; + + // Given the following knapsack constraint: + // sum_j a_j x_j <= beta + // + // We solve the following separation problem: + // minimize sum_j (1 - xstar_j) z_j + // subject to sum_j a_j z_j > beta + // z_j in {0, 1} + // When z_j = 1, then j is in the cover. + // Let phi_star be the optimal objective of this problem. + // We have a violated cover when phi_star < 1.0 + // + // We convert this problem into a 0-1 knapsack problem + // maximize sum_j (1 - xstar_j) zbar_j + // subject to sum_j a_j zbar_j <= sum_j a_j - (beta + 1) + // zbar_j in {0, 1} + // where zbar_j = 1 - z_j + // This problem is in the form of a 0-1 knapsack problem + // which we can solve with dynamic programming or generate + // a heuristic solution with a greedy algorithm. // Remove the slacks from the inequality f_t seperation_rhs = 0.0; if (verbose) { settings.log.printf(" Knapsack : "); } - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; + std::vector complemented_variables; + complemented_variables.reserve(knapsack_inequality.size()); + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); if (is_slack_[j]) { - knapsack_inequality.x[k] = 0.0; + knapsack_inequality.vector.x[k] = 0.0; } else { - if (verbose) { settings.log.printf(" %g x%d +", knapsack_inequality.x[k], j); } - seperation_rhs += knapsack_inequality.x[k]; + const f_t aj = knapsack_inequality.vector.x[k]; + if (aj < 0.0) { + knapsack_inequality.rhs -= aj; + knapsack_inequality.vector.x[k] *= -1.0; + complemented_variables.push_back(j); + is_complemented_[j] = 1; + } + if (verbose) { settings.log.printf(" %g x%d +", knapsack_inequality.vector.x[k], j); } + seperation_rhs += knapsack_inequality.vector.x[k]; } } - if (verbose) { settings.log.printf(" <= %g\n", knapsack_rhs); } - seperation_rhs -= (knapsack_rhs + 1); + if (verbose) { settings.log.printf(" <= %g\n", knapsack_inequality.rhs); } + seperation_rhs -= (knapsack_inequality.rhs + 1); if (verbose) { settings.log.printf("\t"); - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); if (!is_slack_[j]) { if (std::abs(xstar[j]) > 1e-3) { settings.log.printf("x_relax[%d]= %g ", j, xstar[j]); } } @@ -784,71 +839,445 @@ i_t knapsack_generation_t::generate_knapsack_cuts( settings.log.printf("seperation_rhs %g\n", seperation_rhs); } - if (seperation_rhs <= 0.0) { return -1; } + + if (seperation_rhs <= 0.0) { + restore_complemented(complemented_variables); + return -1; + } std::vector values; - values.resize(knapsack_inequality.i.size() - 1); + values.reserve(knapsack_inequality.size() - 1); std::vector weights; - weights.resize(knapsack_inequality.i.size() - 1); - i_t h = 0; + weights.reserve(knapsack_inequality.size() - 1); + std::vector indices; + indices.reserve(knapsack_inequality.size() - 1); f_t objective_constant = 0.0; - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; + std::vector fixed_variables; + std::vector fixed_values; + 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); if (!is_slack_[j]) { - const f_t vj = std::min(1.0, std::max(0.0, 1.0 - xstar[j])); + const f_t xstar_j = is_complemented_[j] ? 1.0 - xstar[j] : xstar[j]; + complemented_xstar_[j] = xstar_j; + const f_t vj = std::min(1.0, std::max(0.0, 1.0 - xstar_j)); + if (xstar_j < x_tol) { + // 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); + seperation_rhs -= knapsack_inequality.vector.x[k]; + // No need to adjust the objective constant + continue; + } + if (xstar_j > 1.0 - x_tol) { + // 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); + // Note seperation rhs is unchanged + objective_constant += vj; + continue; + } objective_constant += vj; - values[h] = vj; - weights[h] = knapsack_inequality.x[k]; - h++; + values.push_back(vj); + weights.push_back(knapsack_inequality.vector.x[k]); + indices.push_back(j); } } std::vector solution; - solution.resize(knapsack_inequality.i.size() - 1); + solution.resize(values.size()); + + if (seperation_rhs <= 0.0) { + restore_complemented(complemented_variables); + return -1; + } + + f_t objective = 0.0; + if (!values.empty()) { + if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } - if (verbose) { settings.log.printf("Calling solve_knapsack_problem\n"); } - f_t objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); - if (std::isnan(objective)) { return -1; } + objective = solve_knapsack_problem(values, weights, seperation_rhs, solution); + } else { + solution.clear(); + } + if (std::isnan(objective)) { restore_complemented(complemented_variables); return -1; } if (verbose) { settings.log.printf("objective %e objective_constant %e\n", objective, objective_constant); } f_t seperation_value = -objective + objective_constant; if (verbose) { settings.log.printf("seperation_value %e\n", seperation_value); } const f_t tol = 1e-6; - if (seperation_value >= 1.0 - tol) { return -1; } + if (seperation_value >= 1.0 - tol) { restore_complemented(complemented_variables); return -1; } i_t cover_size = 0; for (i_t k = 0; k < solution.size(); k++) { if (solution[k] == 0.0) { cover_size++; } } + for (i_t k = 0; k < fixed_values.size(); k++) { + if (fixed_values[k] == 1.0) { cover_size++; } + } cut.reserve(cover_size); cut.clear(); - h = 0; - for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { - const i_t j = knapsack_inequality.i[k]; - if (!is_slack_[j]) { - if (solution[h] == 0.0) { cut.push_back(j, -1.0); } - h++; - } + for (i_t k = 0; k < solution.size(); k++) { + const i_t j = indices[k]; + if (solution[k] == 0.0) { cut.push_back(j, -1.0); } + } + for (i_t k = 0; k < fixed_variables.size(); k++) { + const i_t j = fixed_variables[k]; + if (fixed_values[k] == 1.0) { cut.push_back(j, -1.0); } } cut.rhs = -cover_size + 1; - cut.sort(); // The cut is in the form: - sum_{j in cover} x_j >= -cover_size + 1 // Which is equivalent to: sum_{j in cover} x_j <= cover_size - 1 + + // Compute the minimal cover and partition the variables into C1 and C2 + inequality_t minimal_cover_cut(lp.num_cols); + std::vector c1_partition; + std::vector c2_partition; + minimal_cover_and_partition( + knapsack_inequality, cut, complemented_xstar_, minimal_cover_cut, c1_partition, c2_partition); + + // Lift the cut + inequality_t lifted_cut(lp.num_cols); + lift_knapsack_cut(knapsack_inequality, + minimal_cover_cut, + c1_partition, + c2_partition, + lifted_cut); + lifted_cut.negate(); + + + // The cut is now in the form: + // -\sum_{j in C} x_j - \sum_{j in F} alpha_j x_j >= -cover_size + 1 + for (i_t k = 0; k < lifted_cut.size(); k++) { + const i_t j = lifted_cut.index(k); + // \sum_{k!=j} d_k x_k + d_j xbar_j >= gamma + // xbar_j = 1 - x_j + // \sum_{k!=j} d_k x_k + d_j (1 - x_j) >= gamma + // \sum_{k!=j} d_k x_k + d_j - d_j x_j >= gamma + // \sum_{k!=j} d_k x_k + d_j x_j >= gamma - d_j + if (is_complemented_[j]) { + lifted_cut.rhs -= lifted_cut.vector.x[k]; + lifted_cut.vector.x[k] *= -1.0; + } + } + lifted_cut.sort(); + // Verify the cut is violated - f_t dot = cut.vector.dot(xstar); - f_t violation = dot - cut.rhs; + f_t lifted_dot = lifted_cut.vector.dot(xstar); + f_t lifted_violation = lifted_dot - lifted_cut.rhs; if (verbose) { - settings.log.printf("Knapsack cut %d violation %e < 0\n", knapsack_row, violation); + settings.log.printf("Knapsack cut %d lifted violation %e < 0\n", knapsack_row, lifted_violation); } - if (violation >= -tol) { return -1; } + if (lifted_violation >= -tol) { + restore_complemented(complemented_variables); + return -1; + } + + cut = lifted_cut; + restore_complemented(complemented_variables); return 0; } +template +bool knapsack_generation_t::is_minimal_cover(f_t cover_sum, + f_t beta, + const std::vector& cover_coefficients) +{ + // Check if the cover is minimial + // A set C is a cover if + // sum_{j in C} a_j > beta + // A set C is a minimal cover if + // sum_{k in C \ {j}} a_k <= beta for all j in C + bool minimal = true; + + // cover_sum = sum_{j in C} a_j + + // A cover is minimal if cover_sum - a_j <= beta for all j in C + + for (i_t k = 0; k < cover_coefficients.size(); k++) { + const f_t a_j = cover_coefficients[k]; + if (a_j == 0.0) { continue; } + if (cover_sum - a_j > beta) { + minimal = false; + break; + } + } + return minimal; + +} + +template +void knapsack_generation_t::minimal_cover_and_partition( + const inequality_t& knapsack_inequality, + const inequality_t& negated_base_cut, + const std::vector& xstar, + inequality_t& minimal_cover_cut, + std::vector& c1_partition, + std::vector& c2_partition) +{ + // Compute the minimal cover cut + inequality_t base_cut = negated_base_cut; + base_cut.negate(); + + std::vector cover_indicies; + cover_indicies.reserve(base_cut.size()); + + std::vector cover_coefficients; + cover_coefficients.reserve(base_cut.size()); + + std::vector score; + score.reserve(base_cut.size()); + + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + workspace_[j] = knapsack_inequality.coeff(k); + } + + for (i_t k = 0; k < base_cut.size(); k++) { + const i_t j = base_cut.index(k); + const f_t xstar_j = xstar[j]; + score.push_back((1.0 - xstar_j) / workspace_[j]); + cover_indicies.push_back(j); + cover_coefficients.push_back(workspace_[j]); + } + + f_t cover_sum = std::accumulate(cover_coefficients.begin(), cover_coefficients.end(), 0.0); + + bool is_minimal = is_minimal_cover(cover_sum, knapsack_inequality.rhs, cover_coefficients); + + if (is_minimal) { + minimal_cover_cut = base_cut; + return; + } + + // We don't have a minimal cover. So sort the score from smallest to largest breaking ties by largest to smallest a_j + std::vector permuation(cover_indicies.size()); + std::iota(permuation.begin(), permuation.end(), 0); + std::sort(permuation.begin(), permuation.end(), [&](i_t a, i_t b) { + if (score[a] < score[b]) { + return true; + } else if (score[a] == score[b]) { + return cover_coefficients[a] > cover_coefficients[b]; + } else { + return false; + } + }); + + const f_t beta = knapsack_inequality.rhs; + for (i_t k = 0; k < permuation.size(); k++) { + const i_t h = permuation[k]; + const f_t a_j = cover_coefficients[h]; + if (cover_sum - a_j > beta) { + // sum_{k in C} a_k - a_j > beta + // so sum_{k in C \ {k}} a_k > beta and C \ {k} remains a cover + + cover_sum -= a_j; + // Set the coefficient to 0 to remove it from the cover + cover_coefficients[h] = 0.0; + + is_minimal = is_minimal_cover(cover_sum, beta, cover_coefficients); + if (is_minimal) { + break; + } + } else { + // C \ {j} is no longer a cover. + continue; + } + } + + // Go through and correct cover_indicies and cover_coefficients + for (i_t k = 0; k < cover_coefficients.size(); ) { + if (cover_coefficients[k] == 0.0) { + cover_indicies[k] = cover_indicies.back(); + cover_indicies.pop_back(); + cover_coefficients[k] = cover_coefficients.back(); + cover_coefficients.pop_back(); + } else { + k++; + } + } + + // We now have a minimal cover cut + // sum_{j in C} x_j <= |C| - 1 + minimal_cover_cut.vector.i = cover_indicies; + minimal_cover_cut.vector.x.resize(cover_indicies.size(), 1.0); + minimal_cover_cut.rhs = cover_coefficients.size() - 1; + + // Now we need to partition the variables into C1 and C2 + // C2 = {j in C | x_j = 1} + // C1 = C / C2 + + const f_t x_tol = 1e-5; + for (i_t j : cover_indicies) { + if (xstar[j] > 1.0 - x_tol) { + c2_partition.push_back(j); + } else { + c1_partition.push_back(j); + } + } +} + +template +void knapsack_generation_t::lift_knapsack_cut( + const inequality_t& knapsack_inequality, + const inequality_t& base_cut, + const std::vector& c1_partition, + const std::vector& c2_partition, + inequality_t& lifted_cut) +{ + // The base cut is in the form: sum_{j in cover} x_j <= |cover| - 1 + + // We will attempt to lift the cut by adding a new variable x_k with k not in C to the base cut + // so that the cut becomes + // sum_{j in cover} x_j + alpha_k * x_k <= |cover| - 1 + // + // We can do this for multiple variables so that in the end the cut becomes + // + // sum_{j in cover} x_j + sum_{k in F} alpha_k * x_k <= |cover| - 1 + + // Determine which variables are in the knapsack constraint and not in the cover + std::vector marked_variables; + marked_variables.reserve(knapsack_inequality.size()); + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + if (!is_marked_[j]) { + is_marked_[j] = 1; // is_marked_[j] = 1 for all j in N + marked_variables.push_back(j); + } + } + for (i_t k = 0; k < base_cut.size(); k++) { + const i_t j = base_cut.index(k); + if (is_marked_[j]) { + is_marked_[j] = 0; // is_marked_[j] = 1 for all j in N \ C + // OK to leave marked_variables unchanged as marked_variables will be a superset of all dirty is_marked + } + } + std::vector remaining_variables; + std::vector remaining_indices; + std::vector remaining_coefficients; + remaining_variables.reserve(knapsack_inequality.size()); + remaining_indices.reserve(knapsack_inequality.size()); + remaining_coefficients.reserve(knapsack_inequality.size()); + + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + if (is_marked_[j]) { + remaining_variables.push_back(j); + remaining_indices.push_back(k); + remaining_coefficients.push_back(knapsack_inequality.coeff(k)); + } + } + + // We start with F = {} and lift remaining variables one by one + // For a variable k not in C union F, the inequality + // + // alpha_k * x_k + sum_{j in C} x_j <= |C| - 1 + // is trivially satisfied when x_k = 0. + // If x_k = 1, then the inequality will be valid for all alpha_k + // such that + // alpha_k + maximize sum_{j in C} x_j <= |C| - 1 + // subject to a_k + sum_{j in C} a_j x_j <= beta + // + // where here we require a_k + sum_{j in C} a_j x_j <= beta so that the inequality + // is valid for the set { x_j in {0, 1}^(|C| + 1) | sum_{j in C union k} a_j x_j <= beta} + // + // Let phi^star_k denote the optimal objective value of the problem + // + // maximize sum_{j in C} x_j + // subject to a_k + sum_{j in C} a_j x_j <= beta + // x_j in {0, 1} for all j in C + // Then alpha_k <= |C| - 1 - phi^star_k + // and we can set alpha_k = |C| - 1 - phi^star_k + // + // We can continue this process for each variable k not in C union F + // + // Assume the valid inequality + // sum_{j in C} x_j + sum_{j in F} alpha_j * x_j <= |C| - 1 + // has been obtained so far. We now add the variable x_k with k not in C union F to the inequality. + // So we have + // alpha_k * x_k + sum_{j in C} x_j + sum_{j in F} alpha_j * x_j <= |C| - 1 + // + // Again, this is trivially satisfied when x_k = 0. And we can determine the max value of alpha_k + // by solving the 0-1 knapsack problem: + // + // maximize sum_{j in C} x_j + sum_{j in F} alpha_j * x_j + // subject to sum_{j in C} a_j x_j + sum_{j in F} a_j * x_j <= beta - a_k + // x_j in {0, 1} for all j in C union F + // + // Let phi^star_k denote the optimal objective value of the knapsack problem. + // The lifted coefficient alpha_k = |C| - 1 - phi^star_k + + + // Construct weight and values for C + std::vector values; + values.reserve(knapsack_inequality.size()); + + std::vector weights; + weights.reserve(knapsack_inequality.size()); + + for (i_t k = 0; k < knapsack_inequality.size(); k++) { + const i_t j = knapsack_inequality.index(k); + if (!is_marked_[j]) { + // j is in C + weights.push_back(knapsack_inequality.coeff(k)); + values.push_back(1.0); + } + } + + std::vector F; + std::vector alpha; + + std::vector solution; + + f_t cover_size = base_cut.size(); + + lifted_cut = base_cut; + + // Sort the coefficients such that the largest coefficients are lifted first + // We will pop the largest coefficients from the back of the permutation + std::vector permutation; + best_score_last_permutation(remaining_coefficients, permutation); + + while (permutation.size() > 0) { + const i_t h = permutation.back(); + const i_t k = remaining_variables[h]; + const f_t a_k = remaining_coefficients[h]; + + f_t capacity = knapsack_inequality.rhs - a_k; + + f_t objective = greedy_knapsack_problem(values, weights, capacity, solution); + if (std::isnan(objective)) { settings_.log.printf("lifting knapsack problem failed\n"); break; } + + f_t alpha_k = std::max(0.0, cover_size - 1.0 - objective); + + if (alpha_k > 0.0) { + settings_.log.printf("Lifted variable %d with alpha %g\n", k, alpha_k); + F.push_back(k); + alpha.push_back(alpha_k); + values.push_back(alpha_k); + weights.push_back(a_k); + + lifted_cut.vector.i.push_back(k); + lifted_cut.vector.x.push_back(alpha_k); + } + + // Remove the variable from the permutation + permutation.pop_back(); + } + settings_.log.printf("Lifted %ld variables\n", F.size()); + + // Restore is_marked_ + for (i_t j : marked_variables) { + is_marked_[j] = 0; + } + +} + template f_t knapsack_generation_t::greedy_knapsack_problem(const std::vector& values, const std::vector& weights, @@ -970,6 +1399,8 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector dp(n + 1, sum_value + 1, INT_INF); dense_matrix_t take(n + 1, sum_value + 1, 0); @@ -1014,6 +1445,121 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector +void cut_generation_t::generate_implied_bounds_cuts( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time) +{ + if (probing_implied_bounds_.zero_offsets.empty()) { return; } + + const f_t tol = 1e-4; + i_t num_cuts = 0; + const i_t pib_cols = static_cast(probing_implied_bounds_.zero_offsets.size()) - 1; + const i_t n_cols = std::min(lp.num_cols, pib_cols); + + for (i_t j = 0; j < n_cols; j++) { + if (var_types[j] == variable_type_t::CONTINUOUS) { continue; } + const f_t xstar_j = xstar[j]; + + // x_j = 0 implications + const i_t zero_begin = probing_implied_bounds_.zero_offsets[j]; + const i_t zero_end = probing_implied_bounds_.zero_offsets[j + 1]; + for (i_t p = zero_begin; p < zero_end; p++) { + const i_t i = probing_implied_bounds_.zero_variables[p]; + const f_t l_i = lp.lower[i]; + const f_t u_i = lp.upper[i]; + + // Tightened upper bound: x_j = 0 implies y_i <= b, where b < u_i + // Valid inequality: y_i <= b + (u_i - b)*x_j or -y_i + (u_i - b)*x_j >= -b + const f_t b_ub = probing_implied_bounds_.zero_upper_bound[p]; + if (b_ub < u_i - tol) { + const f_t coeff_j = u_i - b_ub; + const f_t y_i = xstar[i]; + const f_t lhs_val = -y_i + coeff_j * xstar_j; + const f_t rhs_val = -b_ub; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, -1.0); + cut.push_back(j, coeff_j); + cut.rhs = -b_ub; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + + // Tightened lower bound: x_j = 0 implies y_i >= b, where b > l_i + // Valid inequality: y_i >= b - (b - l_i)*x_j or y_i + (b - l_i)*x_j >= b + const f_t b_lb = probing_implied_bounds_.zero_lower_bound[p]; + if (b_lb > l_i + tol) { + const f_t coeff_j = b_lb - l_i; + const f_t y_i = xstar[i]; + const f_t lhs_val = y_i + coeff_j * xstar_j; + const f_t rhs_val = b_lb; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, 1.0); + cut.push_back(j, coeff_j); + cut.rhs = b_lb; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + } + + // x_j = 1 implications + const i_t one_begin = probing_implied_bounds_.one_offsets[j]; + const i_t one_end = probing_implied_bounds_.one_offsets[j + 1]; + for (i_t p = one_begin; p < one_end; p++) { + const i_t i = probing_implied_bounds_.one_variables[p]; + const f_t l_i = lp.lower[i]; + const f_t u_i = lp.upper[i]; + + // Tightened upper bound: x_j = 1 implies y_i <= b, where b < u_i + // Valid inequality: y_i <= u_i - (u_i - b)*x_j or -y_i - (u_i - b)*x_j >= -u_i + const f_t b_ub = probing_implied_bounds_.one_upper_bound[p]; + if (b_ub < u_i - tol) { + const f_t coeff_j = -(u_i - b_ub); + const f_t y_i = xstar[i]; + const f_t lhs_val = -y_i + coeff_j * xstar_j; + const f_t rhs_val = -u_i; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, -1.0); + cut.push_back(j, coeff_j); + cut.rhs = -u_i; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + + // Tightened lower bound: x_j = 1 implies y_i >= b, where b > l_i + // Valid inequality: y_i >= l_i + (b - l_i)*x_j or y_i - (b - l_i)*x_j >= l_i + const f_t b_lb = probing_implied_bounds_.one_lower_bound[p]; + if (b_lb > l_i + tol) { + const f_t coeff_j = -(b_lb - l_i); + const f_t lhs_val = xstar[i] + coeff_j * xstar_j; + const f_t rhs_val = l_i; + if (lhs_val < rhs_val - tol) { + inequality_t cut; + cut.push_back(i, 1.0); + cut.push_back(j, coeff_j); + cut.rhs = rhs_val; + cut_pool_.add_cut(cut_type_t::IMPLIED_BOUNDS, cut); + num_cuts++; + } + } + } + } + + if (num_cuts > 0) { + settings.log.debug("Generated %d implied bounds cuts from probing\n", num_cuts); + } +} + template bool cut_generation_t::generate_cuts(const lp_problem_t& lp, const simplex_solver_settings_t& settings, @@ -1045,8 +1591,8 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, f_t cut_start_time = tic(); generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar); f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings.log.debug("Knapsack cut generation time %.2f seconds\n", cut_generation_time); + if (1 || cut_generation_time > 1.0) { + settings.log.printf("Knapsack cut generation time %.2f seconds\n", cut_generation_time); } } @@ -1073,6 +1619,16 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, settings.log.debug("Clique cut generation time %.2f seconds\n", cut_generation_time); } } + + // Generate implied bounds cuts + if (settings.implied_bounds_cuts != 0) { + f_t cut_start_time = tic(); + generate_implied_bounds_cuts(lp, settings, var_types, xstar, variable_bounds, start_time); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); + } + } return true; } @@ -1088,7 +1644,7 @@ void cut_generation_t::generate_knapsack_cuts( if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { inequality_t cut(lp.num_cols); - i_t knapsack_status = knapsack_generation_.generate_knapsack_cuts( + i_t knapsack_status = knapsack_generation_.generate_knapsack_cut( lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut); if (knapsack_status == 0) { cut_pool_.add_cut(cut_type_t::KNAPSACK, cut); } } @@ -1405,6 +1961,8 @@ void cut_generation_t::generate_mir_cuts( // Transform the relaxation solution std::vector transformed_xstar; complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); + std::vector initial_scores = scores; + std::vector starting_count(lp.num_rows, 0); const i_t max_cuts = std::min(lp.num_rows, 100000); f_t work_estimate = 0.0; @@ -1638,8 +2196,8 @@ void cut_generation_t::generate_mir_cuts( // Clear the aggregated rows aggregated_rows.clear(); - // Set the score of the current row to zero - scores[i] = 0.0; + starting_count[i]++; + scores[i] = initial_scores[i] * std::pow(0.9, starting_count[i]); score_queue.push(std::make_pair(scores[i], i)); work_estimate += std::log2(std::max(1, static_cast(score_queue.size()))); } @@ -1708,7 +2266,7 @@ void cut_generation_t::generate_gomory_cuts( complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_A_float); inequality_t cut_A(lp.num_cols); - if (cut_ok) { cut_ok = gomory_cut.rational_coefficients(var_types, cut_A_float, cut_A); } + if (cut_ok) { cut_ok = rational_coefficients(var_types, cut_A_float, cut_A); } // See if the inequality is violated by the original relaxation solution f_t cut_A_violation = complemented_mir.compute_violation(cut_A, xstar); @@ -1745,7 +2303,7 @@ void cut_generation_t::generate_gomory_cuts( complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_B_float); inequality_t cut_B(lp.num_cols); - if (cut_ok) { cut_ok = gomory_cut.rational_coefficients(var_types, cut_B_float, cut_B); } + if (cut_ok) { cut_ok = rational_coefficients(var_types, cut_B_float, cut_B); } bool B_valid = false; f_t cut_B_distance = 0.0; @@ -1930,11 +2488,11 @@ i_t tableau_equality_t::generate_base_equality( return 0; } -template -bool mixed_integer_gomory_cut_t::rational_approximation(f_t x, - int64_t max_denominator, - int64_t& numerator, - int64_t& denominator) +template +bool rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator) { int64_t a, p0 = 0, q0 = 1, p1 = 1, q1 = 0; f_t val = x; @@ -1970,7 +2528,7 @@ bool mixed_integer_gomory_cut_t::rational_approximation(f_t x, } template -bool mixed_integer_gomory_cut_t::rational_coefficients( +bool rational_coefficients( const std::vector& var_types, const inequality_t& input_inequality, inequality_t& rational_inequality) @@ -2007,8 +2565,7 @@ bool mixed_integer_gomory_cut_t::rational_coefficients( return true; } -template -int64_t mixed_integer_gomory_cut_t::gcd(const std::vector& integers) +int64_t gcd(const std::vector& integers) { if (integers.empty()) { return 0; } @@ -2019,8 +2576,7 @@ int64_t mixed_integer_gomory_cut_t::gcd(const std::vector& in return result; } -template -int64_t mixed_integer_gomory_cut_t::lcm(const std::vector& integers) +int64_t lcm(const std::vector& integers) { if (integers.empty()) { return 0; } int64_t result = diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 91806d81a..2653ee793 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -38,7 +38,8 @@ enum cut_type_t : int8_t { KNAPSACK = 2, CHVATAL_GOMORY = 3, CLIQUE = 4, - MAX_CUT_TYPE = 5 + IMPLIED_BOUNDS = 5, + MAX_CUT_TYPE = 6 }; template @@ -62,6 +63,46 @@ cut_gap_closure_t compute_cut_gap_closure(f_t objective_reference, return {initial_gap, final_gap, gap_closed, gap_closed_ratio}; } +template +struct probing_implied_bounds_t { + // Probing implications stored in CSR format, indexed by binary variable x_j. + // + // "zero" = implications discovered when probing x_j = 0. + // "one" = implications discovered when probing x_j = 1. + // + // For a binary variable x_j, the range + // zero_offsets[j] .. zero_offsets[j+1] + // indexes into the flat arrays zero_variables, zero_lower_bound, zero_upper_bound. + // + // For each position p in that range: + // zero_variables[p] = i if variable y_i bounds were tightened + // when x_j was fixed to 0 and constraints were propagated. + // zero_lower_bound[p] = tightened lower bound on y_i (i.e., x_j = 0 => y_i >= zero_lower_bound[p]). + // zero_upper_bound[p] = tightened upper bound on y_i (i.e., x_j = 0 => y_i <= zero_upper_bound[p]). + // + // The one arrays are analogous for probing x_j = 1. + // + // Non-binary variables have empty ranges (zero_offsets[j] == zero_offsets[j+1]). + // Offsets vectors have size num_cols + 1. + + probing_implied_bounds_t() = default; + + probing_implied_bounds_t(i_t num_cols) + : zero_offsets(num_cols + 1, 0), one_offsets(num_cols + 1, 0) + { + } + + std::vector zero_offsets; + std::vector zero_variables; + std::vector zero_lower_bound; + std::vector zero_upper_bound; + + std::vector one_offsets; + std::vector one_variables; + std::vector one_lower_bound; + std::vector one_upper_bound; +}; + template struct inequality_t { inequality_t() : vector(), rhs(0.0) {} @@ -132,7 +173,7 @@ struct cut_info_t { } } const char* cut_type_names[MAX_CUT_TYPE] = { - "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique "}; + "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique ", "Implied Bounds"}; std::array num_cuts = {0}; }; @@ -288,7 +329,7 @@ class knapsack_generation_t { const std::vector& new_slacks, const std::vector& var_types); - i_t generate_knapsack_cuts(const lp_problem_t& lp, + i_t generate_knapsack_cut(const lp_problem_t& lp, const simplex_solver_settings_t& settings, csr_matrix_t& Arow, const std::vector& new_slacks, @@ -301,6 +342,28 @@ class knapsack_generation_t { const std::vector& get_knapsack_constraints() const { return knapsack_constraints_; } private: + void restore_complemented(const std::vector& complemented_variables) { + for (i_t j : complemented_variables) { + is_complemented_[j] = 0; + } + } + bool is_minimal_cover(f_t cover_sum, + f_t beta, + const std::vector& cover_coefficients); + + void minimal_cover_and_partition(const inequality_t& knapsack_inequality, + const inequality_t& negated_base_cut, + const std::vector& xstar, + inequality_t& minimal_cover_cut, + std::vector& c1_partition, + std::vector& c2_partition); + + void lift_knapsack_cut(const inequality_t& knapsack_inequality, + const inequality_t& base_cut, + const std::vector& c1_partition, + const std::vector& c2_partition, + inequality_t& lifted_cut); + // Generate a heuristic solution to the 0-1 knapsack problem f_t greedy_knapsack_problem(const std::vector& values, const std::vector& weights, @@ -315,6 +378,10 @@ class knapsack_generation_t { std::vector is_slack_; std::vector knapsack_constraints_; + std::vector is_complemented_; + std::vector is_marked_; + std::vector workspace_; + std::vector complemented_xstar_; const simplex_solver_settings_t& settings_; }; @@ -336,12 +403,14 @@ class cut_generation_t { const std::vector& new_slacks, const std::vector& var_types, const user_problem_t& user_problem, + const probing_implied_bounds_t& probing_implied_bounds, std::shared_ptr> clique_table = nullptr, std::future>>* clique_table_future = nullptr, std::atomic* signal_extend = nullptr) : cut_pool_(cut_pool), knapsack_generation_(lp, settings, Arow, new_slacks, var_types), user_problem_(user_problem), + probing_implied_bounds_(probing_implied_bounds), clique_table_(std::move(clique_table)), clique_table_future_(clique_table_future), signal_extend_(signal_extend) @@ -400,9 +469,18 @@ class cut_generation_t { const std::vector& reduced_costs, f_t start_time); + // Generate implied bounds cuts from probing implications + void generate_implied_bounds_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + variable_bounds_t& variable_bounds, + f_t start_time); + cut_pool_t& cut_pool_; knapsack_generation_t knapsack_generation_; const user_problem_t& user_problem_; + const probing_implied_bounds_t& probing_implied_bounds_; std::shared_ptr> clique_table_; std::future>>* clique_table_future_{nullptr}; std::atomic* signal_extend_{nullptr}; @@ -464,19 +542,6 @@ template class mixed_integer_gomory_cut_t { public: mixed_integer_gomory_cut_t() {} - - bool rational_coefficients(const std::vector& var_types, - const inequality_t& input_inequality, - inequality_t& rational_inequality); - - private: - bool rational_approximation(f_t x, - int64_t max_denominator, - int64_t& numerator, - int64_t& denominator); - - int64_t gcd(const std::vector& integers); - int64_t lcm(const std::vector& integers); }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index eadd93040..67e52f9f7 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -100,6 +100,7 @@ struct simplex_solver_settings_t { mir_cuts(-1), mixed_integer_gomory_cuts(-1), knapsack_cuts(-1), + implied_bounds_cuts(-1), clique_cuts(-1), strong_chvatal_gomory_cuts(-1), reduced_cost_strengthening(-1), @@ -180,6 +181,7 @@ struct simplex_solver_settings_t { i_t mixed_integer_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable mixed integer Gomory // cuts i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts + i_t implied_bounds_cuts; // -1 automatic, 0 to disable, >0 to enable implied bounds cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory // cuts diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index d300d6011..c0f919020 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -706,7 +706,8 @@ i_t solve(const user_problem_t& problem, { i_t status; if (is_mip(problem) && !settings.relaxation) { - branch_and_bound_t branch_and_bound(problem, settings, tic()); + probing_implied_bounds_t empty_probing(problem.num_cols); + branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing); mip_solution_t mip_solution(problem.num_cols); mip_status_t mip_status = branch_and_bound.solve(mip_solution); if (mip_status == mip_status_t::OPTIMAL) { @@ -745,7 +746,8 @@ i_t solve_mip_with_guess(const user_problem_t& problem, { i_t status; if (is_mip(problem)) { - branch_and_bound_t branch_and_bound(problem, settings, tic()); + probing_implied_bounds_t empty_probing(problem.num_cols); + branch_and_bound_t branch_and_bound(problem, settings, tic(), empty_probing); branch_and_bound.set_initial_guess(guess); mip_status_t mip_status = branch_and_bound.solve(solution); if (mip_status == mip_status_t::OPTIMAL) { diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index d7d760101..de33d5487 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -22,6 +22,7 @@ #include #include +#include #include namespace cuopt::linear_programming::detail { @@ -274,8 +275,10 @@ void rins_t::run_rins() f_t objective) { rins_solution_queue.push_back(solution); }; + dual_simplex::probing_implied_bounds_t empty_probing( + branch_and_bound_problem.num_cols); dual_simplex::branch_and_bound_t branch_and_bound( - branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic()); + branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing); branch_and_bound.set_initial_guess(cuopt::host_copy(fixed_assignment, rins_handle.get_stream())); branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution); diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index a867141d0..42e33d214 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -117,8 +117,10 @@ class sub_mip_recombiner_t : public recombiner_t { // disable B&B logs, so that it is not interfering with the main B&B thread branch_and_bound_settings.log.log = false; + dual_simplex::probing_implied_bounds_t empty_probing( + branch_and_bound_problem.num_cols); dual_simplex::branch_and_bound_t branch_and_bound( - branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic()); + branch_and_bound_problem, branch_and_bound_settings, dual_simplex::tic(), empty_probing); branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution); if (solution_vector.size() > 0) { cuopt_assert(fixed_assignment.size() == branch_and_bound_solution.x.size(), diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index f25c093af..0e03a529f 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -192,6 +192,8 @@ solution_t mip_solver_t::run_solver() branch_and_bound_solution_helper_t solution_helper(&dm, branch_and_bound_settings); dual_simplex::mip_solution_t branch_and_bound_solution(1); + dual_simplex::probing_implied_bounds_t probing_implied_bounds; + bool run_bb = !context.settings.heuristics_only; if (run_bb) { // Convert the presolved problem to dual_simplex::user_problem_t @@ -199,6 +201,77 @@ solution_t mip_solver_t::run_solver() // Resize the solution now that we know the number of columns/variables branch_and_bound_solution.resize(branch_and_bound_problem.num_cols); + // Extract probing cache into CPU-only CSR struct for implied bounds cuts + { + const i_t num_cols = branch_and_bound_problem.num_cols; + probing_implied_bounds = dual_simplex::probing_implied_bounds_t(num_cols); + auto& pc = dm.ls.constraint_prop.bounds_update.probing_cache.probing_cache; + auto& rev_ids = context.problem_ptr->reverse_original_ids; + + // First pass: count entries per binary variable + for (auto& [var_idx, entries] : pc) { + if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } + i_t j = (var_idx < static_cast(rev_ids.size())) ? rev_ids[var_idx] : -1; + if (j < 0 || j >= num_cols) { continue; } + + for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + probing_implied_bounds.zero_offsets[j + 1]++; + } + for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + probing_implied_bounds.one_offsets[j + 1]++; + } + } + + // Prefix sum + for (i_t j = 0; j < num_cols; j++) { + probing_implied_bounds.zero_offsets[j + 1] += probing_implied_bounds.zero_offsets[j]; + probing_implied_bounds.one_offsets[j + 1] += probing_implied_bounds.one_offsets[j]; + } + + // Allocate flat arrays + i_t zero_nnz = probing_implied_bounds.zero_offsets[num_cols]; + i_t one_nnz = probing_implied_bounds.one_offsets[num_cols]; + probing_implied_bounds.zero_variables.resize(zero_nnz); + probing_implied_bounds.zero_lower_bound.resize(zero_nnz); + probing_implied_bounds.zero_upper_bound.resize(zero_nnz); + probing_implied_bounds.one_variables.resize(one_nnz); + probing_implied_bounds.one_lower_bound.resize(one_nnz); + probing_implied_bounds.one_upper_bound.resize(one_nnz); + + // Second pass: fill flat arrays using write cursors + std::vector zero_cursor(probing_implied_bounds.zero_offsets); + std::vector one_cursor(probing_implied_bounds.one_offsets); + + for (auto& [var_idx, entries] : pc) { + if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } + i_t j = (var_idx < static_cast(rev_ids.size())) ? rev_ids[var_idx] : -1; + if (j < 0 || j >= num_cols) { continue; } + + for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + i_t p = zero_cursor[j]++; + probing_implied_bounds.zero_variables[p] = i; + probing_implied_bounds.zero_lower_bound[p] = bound.lb; + probing_implied_bounds.zero_upper_bound[p] = bound.ub; + } + for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { + i_t i = (imp_var < static_cast(rev_ids.size())) ? rev_ids[imp_var] : -1; + if (i < 0 || i >= num_cols) { continue; } + i_t p = one_cursor[j]++; + probing_implied_bounds.one_variables[p] = i; + probing_implied_bounds.one_lower_bound[p] = bound.lb; + probing_implied_bounds.one_upper_bound[p] = bound.ub; + } + } + + CUOPT_LOG_INFO("Probing implied bounds: %d zero entries, %d one entries", zero_nnz, one_nnz); + } + // Fill in the settings for branch and bound branch_and_bound_settings.time_limit = timer_.get_time_limit(); branch_and_bound_settings.node_limit = context.settings.node_limit; @@ -269,6 +342,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_problem, branch_and_bound_settings, timer_.get_tic_start(), + probing_implied_bounds, context.problem_ptr->clique_table); context.branch_and_bound_ptr = branch_and_bound.get(); auto* stats_ptr = &context.stats;