diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index e1097455f..392931d99 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -64,6 +64,7 @@ #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" #define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS "mip_mixed_integer_gomory_cuts" #define CUOPT_MIP_KNAPSACK_CUTS "mip_knapsack_cuts" +#define CUOPT_MIP_IMPLIED_BOUND_CUTS "mip_implied_bound_cuts" #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 8f4a6c1d7..c2cd33df9 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -94,6 +94,7 @@ class mip_solver_settings_t { i_t mixed_integer_gomory_cuts = -1; i_t knapsack_cuts = -1; i_t clique_cuts = -1; + i_t implied_bound_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; f_t cut_change_threshold = -1.0; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index c9b456839..7106d8dce 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -205,7 +205,11 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound) } else { constexpr int BUFFER_LEN = 32; char buffer[BUFFER_LEN]; - snprintf(buffer, BUFFER_LEN - 1, "%5.1f%%", user_mip_gap * 100); + if (user_mip_gap > 1e-3) { + snprintf(buffer, BUFFER_LEN - 1, "%5.1f%%", user_mip_gap * 100); + } else { + snprintf(buffer, BUFFER_LEN - 1, "%5.2f%%", user_mip_gap * 100); + } return std::string(buffer); } } @@ -243,9 +247,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_bound_t& probing_implied_bound, std::shared_ptr> clique_table) : original_problem_(user_problem), settings_(solver_settings), + probing_implied_bound_(probing_implied_bound), clique_table_(std::move(clique_table)), original_lp_(user_problem.handle_ptr, 1, 1, 1), Arow_(1, 1, 0), @@ -556,7 +562,7 @@ template bool branch_and_bound_t::repair_solution(const std::vector& edge_norms, const std::vector& potential_solution, f_t& repaired_obj, - std::vector& repaired_solution) const + std::vector& repaired_solution) { bool feasible = false; repaired_obj = std::numeric_limits::quiet_NaN(); @@ -579,9 +585,10 @@ bool branch_and_bound_t::repair_solution(const std::vector& edge_ i_t iter = 0; f_t lp_start_time = tic(); simplex_solver_settings_t lp_settings = settings_; + lp_settings.concurrent_halt = &node_concurrent_halt_; std::vector vstatus = root_vstatus_; lp_settings.set_log(false); - lp_settings.inside_mip = true; + lp_settings.inside_mip = 2; std::vector leaf_edge_norms = edge_norms; // should probably set the cut off here lp_settings.cut_off dual::status_t lp_status = dual_phase2( @@ -724,12 +731,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; @@ -1336,6 +1337,7 @@ dual::status_t branch_and_bound_t::solve_node_lp( assert(leaf_vstatus.size() == worker->leaf_problem.num_cols); simplex_solver_settings_t lp_settings = settings_; + lp_settings.concurrent_halt = &node_concurrent_halt_; lp_settings.set_log(false); f_t cutoff = get_cutoff(); if (original_lp_.objective_is_integral) { @@ -1432,23 +1434,25 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_trecompute_basis = true; worker->recompute_bounds = true; - while (stack.size() > 0 && solver_status_ == mip_status_t::UNSET) { + f_t lower_bound = get_lower_bound(); + f_t upper_bound = upper_bound_; + f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + f_t abs_gap = upper_bound - lower_bound; + + while (stack.size() > 0 && (solver_status_ == mip_status_t::UNSET && is_running_) && + rel_gap > settings_.relative_mip_gap_tol && abs_gap > settings_.absolute_mip_gap_tol) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); - f_t lower_bound = node_ptr->lower_bound; - f_t upper_bound = upper_bound_; - f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); - // This is based on three assumptions: // - The stack only contains sibling nodes, i.e., the current node and it sibling (if // applicable) // - The current node and its siblings uses the lower bound of the parent before solving the LP // relaxation // - The lower bound of the parent is lower or equal to its children - worker->lower_bound = lower_bound; + worker->lower_bound = node_ptr->lower_bound; - if (lower_bound > get_cutoff()) { + if (node_ptr->lower_bound > get_cutoff()) { search_tree_.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); search_tree_.update(node_ptr, node_status_t::FATHOMED); worker->recompute_basis = true; @@ -1472,6 +1476,9 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_t::plunge_with(branch_and_bound_worker_tget_down_child()); } } + + lower_bound = get_lower_bound(); + upper_bound = upper_bound_; + rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + abs_gap = upper_bound - lower_bound; + } + + lower_bound = get_lower_bound(); + upper_bound = upper_bound_; + rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + abs_gap = upper_bound - lower_bound; + + if (stack.size() > 0 && + (rel_gap <= settings_.relative_mip_gap_tol || abs_gap <= settings_.absolute_mip_gap_tol)) { + // If the solver converged according to the gap rules, but we still have nodes to explore + // in the stack, then we should add all the pending nodes back to the heap so the lower + // bound of the solver is set to the correct value. + while (!stack.empty()) { + auto node = stack.front(); + stack.pop_front(); + node_queue_.push(node); + } } if (settings_.num_threads > 1) { @@ -1549,14 +1578,17 @@ void branch_and_bound_t::dive_with(branch_and_bound_worker_t dive_stats.nodes_explored = 0; dive_stats.nodes_unexplored = 1; - while (stack.size() > 0 && solver_status_ == mip_status_t::UNSET && is_running_) { + f_t lower_bound = get_lower_bound(); + f_t upper_bound = upper_bound_; + f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + f_t abs_gap = upper_bound - lower_bound; + + while (stack.size() > 0 && (solver_status_ == mip_status_t::UNSET && is_running_) && + rel_gap > settings_.relative_mip_gap_tol && abs_gap > settings_.absolute_mip_gap_tol) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); - f_t lower_bound = node_ptr->lower_bound; - f_t upper_bound = upper_bound_; - f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); - worker->lower_bound = lower_bound; + worker->lower_bound = node_ptr->lower_bound; if (node_ptr->lower_bound > get_cutoff()) { worker->recompute_basis = true; @@ -1572,6 +1604,8 @@ void branch_and_bound_t::dive_with(branch_and_bound_worker_t if (lp_status == dual::status_t::TIME_LIMIT) { solver_status_ = mip_status_t::TIME_LIMIT; break; + } else if (lp_status == dual::status_t::CONCURRENT_LIMIT) { + break; } else if (lp_status == dual::status_t::ITERATION_LIMIT) { break; } @@ -1598,6 +1632,11 @@ void branch_and_bound_t::dive_with(branch_and_bound_worker_t if (stack.size() > 1 && stack.front()->depth - stack.back()->depth > diving_backtrack_limit) { stack.pop_back(); } + + lower_bound = get_lower_bound(); + upper_bound = upper_bound_; + rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + abs_gap = upper_bound - lower_bound; } worker_pool_.return_worker_to_pool(worker); @@ -1637,9 +1676,6 @@ void branch_and_bound_t::run_scheduler() rel_gap > settings_.relative_mip_gap_tol && (active_workers_per_strategy_[0] > 0 || node_queue_.best_first_queue_size() > 0)) { bool launched_any_task = false; - lower_bound = get_lower_bound(); - abs_gap = upper_bound_ - lower_bound; - rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound); repair_heuristic_solutions(); @@ -1740,6 +1776,16 @@ void branch_and_bound_t::run_scheduler() } } + lower_bound = get_lower_bound(); + abs_gap = upper_bound_ - lower_bound; + rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound); + + if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) { + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::OPTIMAL; + break; + } + // If no new task was launched in this iteration, suspend temporarily the // execution of the scheduler. As of 8/Jan/2026, GCC does not // implement taskyield, but LLVM does. @@ -1759,10 +1805,6 @@ void branch_and_bound_t::single_threaded_solve() while (solver_status_ == mip_status_t::UNSET && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && node_queue_.best_first_queue_size() > 0) { bool launched_any_task = false; - lower_bound = get_lower_bound(); - abs_gap = upper_bound_ - lower_bound; - rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound); - repair_heuristic_solutions(); f_t now = toc(exploration_stats_.start_time); @@ -1800,6 +1842,15 @@ void branch_and_bound_t::single_threaded_solve() worker.init_best_first(start_node.value(), original_lp_); plunge_with(&worker); + + lower_bound = get_lower_bound(); + abs_gap = upper_bound_ - lower_bound; + rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound); + + if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) { + solver_status_ = mip_status_t::OPTIMAL; + break; + } } } @@ -2151,6 +2202,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut new_slacks_, var_types_, original_problem_, + probing_implied_bound_, clique_table_, &clique_table_future_, &signal_extend_cliques_); @@ -2529,6 +2581,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut node_queue_.push(search_tree_.root.get_up_child()); settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads); + node_concurrent_halt_ = 0; exploration_stats_.nodes_explored = 0; exploration_stats_.nodes_unexplored = 2; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 765f6095b..4faadcc6b 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_bound_t& probing_implied_bound, std::shared_ptr> clique_table = nullptr); // Set an initial guess based on the user_problem. This should be called before solve. @@ -133,7 +134,7 @@ class branch_and_bound_t { bool repair_solution(const std::vector& leaf_edge_norms, const std::vector& potential_solution, f_t& repaired_obj, - std::vector& repaired_solution) const; + std::vector& repaired_solution); f_t get_lower_bound(); bool enable_concurrent_lp_root_solve() const { return enable_concurrent_lp_root_solve_; } @@ -162,6 +163,7 @@ class branch_and_bound_t { private: const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; + const probing_implied_bound_t& probing_implied_bound_; std::shared_ptr> clique_table_; std::future>> clique_table_future_; std::atomic signal_extend_cliques_{false}; @@ -222,6 +224,7 @@ class branch_and_bound_t { omp_atomic_t solving_root_relaxation_{false}; bool enable_concurrent_lp_root_solve_{false}; std::atomic root_concurrent_halt_{0}; + std::atomic node_concurrent_halt_{0}; bool is_root_solution_set{false}; // Pseudocosts diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index c74a12b36..6d7d97ef0 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, @@ -495,7 +510,7 @@ std::vector> find_maximal_cliques_for_test( template void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { - // TODO: Need to deduplicate cuts and only add if the cut is not already in the pool + // TODO: Add fast duplicate check and only add if the cut is not already in the pool for (i_t p = 0; p < cut.size(); p++) { const i_t j = cut.index(p); @@ -572,9 +587,142 @@ f_t cut_pool_t::cut_orthogonality(i_t i, i_t j) return 1.0 - std::abs(dot) / (norm_i * norm_j); } +template +void cut_pool_t::check_for_duplicate_cuts() +{ + // Algorithm from Finding Duplicate Rows in a Linear Programming Model + // by J. A. Tomlin and J.S. Welch + // Operations Research Letters Volume 5, Number 1, June 1986 + std::vector divisors(cut_storage_.m, 0.0); + std::vector sets(cut_storage_.m, 0); + + csc_matrix_t cut_storage_csc(0, 0, 1); + cut_storage_.to_compressed_col(cut_storage_csc); + i_t n = cut_storage_csc.n; + i_t m = cut_storage_csc.m; + + const i_t sentinel = std::numeric_limits::max(); + + i_t new_set = 1; + i_t remaining_potential_duplicates = cut_storage_.m; + for (i_t j = 0; j < n; j++) { + i_t r0 = -1; + i_t new_rows = 0; + i_t new_set_0 = new_set; + new_set++; + const i_t col_start = cut_storage_csc.col_start[j]; + const i_t col_end = cut_storage_csc.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t r = cut_storage_csc.i[p]; + const f_t a_rj = cut_storage_csc.x[p]; + const f_t f_r = divisors[r]; + if (sets[r] == 0) { + r0 = r; // To enable use to find this new set later + sets[r] = new_set_0; + divisors[r] = a_rj; + new_rows++; + } else if (sets[r] < new_set_0) { + // Look over indices a_ij with i > r + for (i_t q = p + 1; q < col_end; q++) { + const i_t i = cut_storage_csc.i[q]; + const f_t a_ij = cut_storage_csc.x[q]; + if (sets[i] == sets[r]) { + // These two rows are currently in the same set + // Check to see if the coefficients still match + const f_t f_i = divisors[i]; + const f_t val = (a_rj / f_r) * (f_i / a_ij); + const f_t epsilon = 1e-10; + if ((val >= 1.0 - epsilon && val <= 1.0 + epsilon)) { + sets[r] = new_set; + sets[i] = new_set; + } + } + } + if (sets[r] >= new_set_0) { // This is only true if a match was found inside the above loop + new_set++; + } else { + sets[r] = sentinel; + remaining_potential_duplicates--; + if (remaining_potential_duplicates == 0) { break; } + } + } + } + if (remaining_potential_duplicates == 0) { break; } + if (new_rows == 1) { + sets[r0] = sentinel; + remaining_potential_duplicates--; + if (remaining_potential_duplicates == 0) { break; } + } + } + + // The cuts are stored in the form: sum_j d_ij x_j >= rhs_i + // We now look for cuts that are duplicates of each other and remove them + std::vector cuts_to_remove(m, 0); + i_t num_cuts_to_remove = 0; + for (i_t r = 0; r < m; r++) { + const i_t set_r = sets[r]; + if (set_r > 0 && set_r < sentinel && cuts_to_remove[r] == 0) { + // This cut has a duplicate + for (i_t i = r + 1; i < m; i++) { + if (sets[i] == set_r) { + const f_t f_r = divisors[r]; + const f_t f_i = divisors[i]; + const f_t theta_r = rhs_storage_[r] / f_r; + const f_t theta_i = rhs_storage_[i] / f_i; + if (f_r > 0 && f_i > 0) { + // We have sum_j d_rj / f_r x_j >= rhs_r / f_r = theta_r + // and sum_j d_ij / f_i x_j >= rhs_i / f_i = theta_i + if (theta_r <= theta_i) { + // Cut i is either the same or stronger than cut r + if (cuts_to_remove[r] == 0) { num_cuts_to_remove++; } + cuts_to_remove[r] = 1; // Remove row r + } else { + // theta_r > theta_i, so cut r is stricly stronger than cut i + if (cuts_to_remove[i] == 0) { num_cuts_to_remove++; } + cuts_to_remove[i] = 1; // Remove row i + } + } else if (f_r < 0 && f_i < 0) { + // We have sum_j d_rj / f_r x_j <= rhs_r / f_r = theta_r + // and sum_j d_ij / f_i x_j <= rhs_i / f_i = theta_i + if (theta_r >= theta_i) { + // Cut i is either the same or stronger than cut r + if (cuts_to_remove[r] == 0) { num_cuts_to_remove++; } + cuts_to_remove[r] = 1; // Remove row r + } else { + // theta_r < theta_i, so cut r is strictly stronger than cut i + if (cuts_to_remove[i] == 0) { num_cuts_to_remove++; } + cuts_to_remove[i] = 1; // Remove row i + } + } + } + } + } + } + + if (num_cuts_to_remove > 0) { + settings_.log.debug("Removing %d duplicate cuts\n", num_cuts_to_remove); + csr_matrix_t new_cut_storage(0, 0, 0); + cut_storage_.remove_rows(cuts_to_remove, new_cut_storage); + cut_storage_ = new_cut_storage; + i_t write = 0; + for (i_t i = 0; i < m; i++) { + if (cuts_to_remove[i] == 0) { + rhs_storage_[write] = rhs_storage_[i]; + cut_type_[write] = cut_type_[i]; + cut_age_[write] = cut_age_[i]; + write++; + } + } + rhs_storage_.resize(write); + cut_type_.resize(write); + cut_age_.resize(write); + } +} + template void cut_pool_t::score_cuts(std::vector& x_relax) { + check_for_duplicate_cuts(); cut_distances_.resize(cut_storage_.m, 0.0); cut_norms_.resize(cut_storage_.m, 0.0); @@ -680,57 +828,68 @@ 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]; - 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; + 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]) { + if (inequality.coeff(p) < 0.0) { + is_knapsack = false; + break; + } + continue; } + 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); } } } @@ -742,7 +901,7 @@ knapsack_generation_t::knapsack_generation_t( } 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 +913,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 +976,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]) { + if (is_slack_[j]) { continue; } + 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); + } + } + + 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 = + exact_knapsack_problem_integer_values_fraction_values(values, weights, capacity, solution); + if (std::isnan(objective)) { + settings_.log.debug("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.debug("Lifted variable %d with alpha %g\n", k, alpha_k); + F.push_back(k); + alpha.push_back(alpha_k); + values.push_back(static_cast(std::round(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(); + } + // 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 +1536,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 +1582,191 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector +f_t knapsack_generation_t::exact_knapsack_problem_integer_values_fraction_values( + const std::vector& values, + const std::vector& weights, + f_t rhs, + std::vector& solution) +{ + // Solve the knapsack problem + // maximize sum_{j=0}^n values[j] * solution[j] + // subject to sum_{j=0}^n weights[j] * solution[j] <= rhs + // values: values of the items + // weights: weights of the items + // return the value of the solution + + const i_t n = weights.size(); + + const bool verbose = false; + i_t sum_value = std::accumulate(values.begin(), 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) { + if (verbose) { settings_.log.printf("sum value %d is negative or too large\n", sum_value); } + return std::numeric_limits::quiet_NaN(); + } + + solution.assign(n, 0.0); + + // dp(j, v) = minimum weight using first j items to get value v + dense_matrix_t dp(n + 1, sum_value + 1, inf); + dense_matrix_t take(n + 1, sum_value + 1, 0); + dp(0, 0) = 0; + + // 4. Dynamic programming + for (i_t j = 1; j <= n; ++j) { + for (i_t v = 0; v <= sum_value; ++v) { + // Do not take item i-1 + dp(j, v) = dp(j - 1, v); + + // Take item j-1 if possible + if (v >= values[j - 1]) { + f_t candidate = dp(j - 1, v - values[j - 1]) + weights[j - 1]; + if (candidate < dp(j, v)) { + dp(j, v) = candidate; + take(j, v) = 1; + } + } + } + } + + // 5. Find best achievable value within capacity + i_t best_value = 0; + for (i_t v = 0; v <= sum_value; ++v) { + if (dp(n, v) <= rhs) { best_value = v; } + } + + // 6. Backtrack to recover solution + i_t v = best_value; + for (i_t j = n; j >= 1; --j) { + if (take(j, v)) { + solution[j - 1] = 1.0; + v -= values[j - 1]; + } else { + solution[j - 1] = 0.0; + } + } + + return f_t(best_value); +} + +template +void cut_generation_t::generate_implied_bound_cuts( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + f_t start_time) +{ + if (probing_implied_bound_.zero_offsets.empty()) { return; } + + const f_t tol = 1e-4; + i_t num_cuts = 0; + const i_t pib_cols = static_cast(probing_implied_bound_.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_bound_.zero_offsets[j]; + const i_t zero_end = probing_implied_bound_.zero_offsets[j + 1]; + for (i_t p = zero_begin; p < zero_end; p++) { + const i_t i = probing_implied_bound_.zero_variables[p]; + if (i == j) { continue; } + 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_bound_.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_BOUND, 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_bound_.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_BOUND, cut); + num_cuts++; + } + } + } + + // x_j = 1 implications + const i_t one_begin = probing_implied_bound_.one_offsets[j]; + const i_t one_end = probing_implied_bound_.one_offsets[j + 1]; + for (i_t p = one_begin; p < one_end; p++) { + const i_t i = probing_implied_bound_.one_variables[p]; + if (i == j) { continue; } + 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_bound_.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_BOUND, 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_bound_.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_BOUND, 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, @@ -1043,7 +1796,7 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, // Generate Knapsack cuts if (settings.knapsack_cuts != 0) { f_t cut_start_time = tic(); - generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar); + generate_knapsack_cuts(lp, settings, Arow, new_slacks, var_types, xstar, start_time); 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); @@ -1073,6 +1826,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 bound cuts + if (settings.implied_bound_cuts != 0) { + f_t cut_start_time = tic(); + generate_implied_bound_cuts(lp, settings, var_types, xstar, 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; } @@ -1083,12 +1846,14 @@ void cut_generation_t::generate_knapsack_cuts( csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types, - const std::vector& xstar) + const std::vector& xstar, + f_t start_time) { if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { + if (toc(start_time) >= settings.time_limit) { return; } 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); } } @@ -1638,7 +2403,6 @@ 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; score_queue.push(std::make_pair(scores[i], i)); work_estimate += std::log2(std::max(1, static_cast(score_queue.size()))); @@ -1708,7 +2472,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 +2509,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 +2694,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,10 +2734,9 @@ bool mixed_integer_gomory_cut_t::rational_approximation(f_t x, } template -bool mixed_integer_gomory_cut_t::rational_coefficients( - const std::vector& var_types, - const inequality_t& input_inequality, - inequality_t& rational_inequality) +bool rational_coefficients(const std::vector& var_types, + const inequality_t& input_inequality, + inequality_t& rational_inequality) { rational_inequality = input_inequality; @@ -2007,8 +2770,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 +2781,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 = @@ -2413,7 +3174,7 @@ bool complemented_mixed_integer_rounding_cut_t::cut_generation_heurist const f_t dist_upper = new_upper_j - x_j; const f_t dist_lower = x_j; const bool between_bounds = x_j > 1e-6 && (new_upper_j == inf || dist_upper > 0.0); - if (between_bounds) { deltas_to_try.push_back(abs_aj); } + if (between_bounds && abs_aj > 1e-6) { deltas_to_try.push_back(abs_aj); } } } if (max_coeff > 1e-6 && max_coeff != 1.0) { diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 91806d81a..2da9760e2 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_BOUND = 5, + MAX_CUT_TYPE = 6 }; template @@ -62,6 +63,47 @@ 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_bound_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_bound_t() = default; + + probing_implied_bound_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) {} @@ -131,9 +173,13 @@ struct cut_info_t { num_cuts[static_cast(cut_type)]++; } } - const char* cut_type_names[MAX_CUT_TYPE] = { - "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique "}; - std::array num_cuts = {0}; + const char* cut_type_names[MAX_CUT_TYPE] = {"Gomory ", + "MIR ", + "Knapsack ", + "Strong CG ", + "Clique ", + "Implied Bounds"}; + std::array num_cuts = {0}; }; template @@ -257,6 +303,8 @@ class cut_pool_t { void print_cutpool_types() { print_cut_types("In cut pool", cut_type_, settings_); } + void check_for_duplicate_cuts(); + private: f_t cut_distance(i_t row, const std::vector& x, f_t& cut_violation, f_t& cut_norm); f_t cut_density(i_t row); @@ -288,19 +336,40 @@ class knapsack_generation_t { const std::vector& new_slacks, const std::vector& var_types); - i_t generate_knapsack_cuts(const lp_problem_t& lp, - const simplex_solver_settings_t& settings, - csr_matrix_t& Arow, - const std::vector& new_slacks, - const std::vector& var_types, - const std::vector& xstar, - i_t knapsack_row, - inequality_t& cut); + 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, + const std::vector& var_types, + const std::vector& xstar, + i_t knapsack_row, + inequality_t& cut); i_t num_knapsack_constraints() const { return knapsack_constraints_.size(); } 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, @@ -313,8 +382,17 @@ class knapsack_generation_t { f_t rhs, std::vector& solution); + f_t exact_knapsack_problem_integer_values_fraction_values(const std::vector& values, + const std::vector& weights, + f_t rhs, + std::vector& solution); + 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 +414,14 @@ class cut_generation_t { const std::vector& new_slacks, const std::vector& var_types, const user_problem_t& user_problem, + const probing_implied_bound_t& probing_implied_bound, 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_bound_(probing_implied_bound), clique_table_(std::move(clique_table)), clique_table_future_(clique_table_future), signal_extend_(signal_extend) @@ -390,7 +470,8 @@ class cut_generation_t { csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types, - const std::vector& xstar); + const std::vector& xstar, + f_t start_time); // Generate clique cuts from conflict graph cliques bool generate_clique_cuts(const lp_problem_t& lp, @@ -400,9 +481,17 @@ class cut_generation_t { const std::vector& reduced_costs, f_t start_time); + // Generate implied bounds cuts from probing implications + void generate_implied_bound_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + f_t start_time); + cut_pool_t& cut_pool_; knapsack_generation_t knapsack_generation_; const user_problem_t& user_problem_; + const probing_implied_bound_t& probing_implied_bound_; std::shared_ptr> clique_table_; std::future>>* clique_table_future_{nullptr}; std::atomic* signal_extend_{nullptr}; @@ -464,19 +553,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/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 43429ba2d..4a9976a81 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -3597,7 +3597,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, 100.0 * dense_delta_z / (sparse_delta_z + dense_delta_z)); ft.print_stats(); } - if (settings.inside_mip && settings.concurrent_halt != nullptr) { + if (settings.inside_mip == 1 && settings.concurrent_halt != nullptr) { settings.log.debug("Setting concurrent halt in Dual Simplex Phase 2\n"); *settings.concurrent_halt = 1; } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index eadd93040..1ba65aad0 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_bound_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_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound 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..5fb2f6c6e 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_bound_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_bound_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/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index a8bc0ea33..183c964a9 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -95,6 +95,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS, &mip_settings.mixed_integer_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_KNAPSACK_CUTS, &mip_settings.knapsack_cuts, -1, 1, -1}, {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, + {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index d7d760101..32ffa778d 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,9 @@ void rins_t::run_rins() f_t objective) { rins_solution_queue.push_back(solution); }; + dual_simplex::probing_implied_bound_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..f6e4f172c 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_bound_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/presolve/conflict_graph/clique_table.cu b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu index 70855267d..82462c11c 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu @@ -146,6 +146,10 @@ void fill_knapsack_constraints(const dual_simplex::user_problem_t& pro std::pair constraint_range = A.get_constraint_range(i); if (constraint_range.second - constraint_range.first < 2) { CUOPT_LOG_DEBUG("Constraint %d has less than 2 variables, skipping", i); + if (problem.row_sense[i] == 'E' && ranged_constraint_counter < problem.num_range_rows && + problem.range_rows[ranged_constraint_counter] == i) { + ranged_constraint_counter++; + } continue; } bool all_binary = true; @@ -158,7 +162,13 @@ void fill_knapsack_constraints(const dual_simplex::user_problem_t& pro } } // if all variables are binary, convert the constraint to a knapsack constraint - if (!all_binary) { continue; } + if (!all_binary) { + if (problem.row_sense[i] == 'E' && ranged_constraint_counter < problem.num_range_rows && + problem.range_rows[ranged_constraint_counter] == i) { + ranged_constraint_counter++; + } + continue; + } knapsack_constraint_t knapsack_constraint; knapsack_constraint.cstr_idx = i; diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 0662584f1..61ab37f6c 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -679,9 +679,18 @@ std::optional> third_party_presolve_t& settings_; }; +// Extract probing cache into CPU-only CSR struct for implied bounds cuts +template +void extract_probing_implied_bounds( + const problem_t& op_problem, + const dual_simplex::user_problem_t& branch_and_bound_problem, + const probing_cache_t& probing_cache, + dual_simplex::probing_implied_bound_t& probing_implied_bound) + +{ + auto& pc = probing_cache.probing_cache; + const i_t num_cols = branch_and_bound_problem.num_cols; + probing_implied_bound = dual_simplex::probing_implied_bound_t(num_cols); + + // First pass: count entries per binary variable + // Probing cache indices are in pre-trivial-presolve space; remap to post-presolve (B&B) space + auto& rev_ids = op_problem.reverse_original_ids; + i_t rev_size = static_cast(rev_ids.size()); + auto remap = [&](i_t raw_idx) -> i_t { + if (rev_size == 0) return raw_idx; + if (raw_idx < 0 || raw_idx >= rev_size) return -1; + return rev_ids[raw_idx]; + }; + auto is_bb_binary = [&](i_t j) { + return branch_and_bound_problem.lower[j] == 0.0 && branch_and_bound_problem.upper[j] == 1.0; + }; + auto bb_bounds_consistent = [&](i_t i, f_t b_lb, f_t b_ub) { + return b_ub >= branch_and_bound_problem.lower[i] - 1e-6 && + b_lb <= branch_and_bound_problem.upper[i] + 1e-6; + }; + for (auto& [var_idx, entries] : pc) { + if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } + i_t j = remap(var_idx); + if (j < 0 || j >= num_cols) { continue; } + if (!is_bb_binary(j)) { continue; } + + for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { + i_t i = remap(imp_var); + if (i < 0 || i >= num_cols) { continue; } + if (!bb_bounds_consistent(i, bound.lb, bound.ub)) { continue; } + probing_implied_bound.zero_offsets[j + 1]++; + } + for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { + i_t i = remap(imp_var); + if (i < 0 || i >= num_cols) { continue; } + if (!bb_bounds_consistent(i, bound.lb, bound.ub)) { continue; } + probing_implied_bound.one_offsets[j + 1]++; + } + } + + // Prefix sum + for (i_t j = 0; j < num_cols; j++) { + probing_implied_bound.zero_offsets[j + 1] += probing_implied_bound.zero_offsets[j]; + probing_implied_bound.one_offsets[j + 1] += probing_implied_bound.one_offsets[j]; + } + + // Allocate flat arrays + i_t zero_nnz = probing_implied_bound.zero_offsets[num_cols]; + i_t one_nnz = probing_implied_bound.one_offsets[num_cols]; + probing_implied_bound.zero_variables.resize(zero_nnz); + probing_implied_bound.zero_lower_bound.resize(zero_nnz); + probing_implied_bound.zero_upper_bound.resize(zero_nnz); + probing_implied_bound.one_variables.resize(one_nnz); + probing_implied_bound.one_lower_bound.resize(one_nnz); + probing_implied_bound.one_upper_bound.resize(one_nnz); + + // Second pass: fill flat arrays using write cursors + std::vector zero_cursor(probing_implied_bound.zero_offsets); + std::vector one_cursor(probing_implied_bound.one_offsets); + + for (auto& [var_idx, entries] : pc) { + if (entries[0].val_interval.interval_type != interval_type_t::EQUALS) { continue; } + i_t j = remap(var_idx); + if (j < 0 || j >= num_cols) { continue; } + if (!is_bb_binary(j)) { continue; } + + for (auto& [imp_var, bound] : entries[0].var_to_cached_bound_map) { + i_t i = remap(imp_var); + if (i < 0 || i >= num_cols) { continue; } + if (!bb_bounds_consistent(i, bound.lb, bound.ub)) { continue; } + i_t p = zero_cursor[j]++; + probing_implied_bound.zero_variables[p] = i; + probing_implied_bound.zero_lower_bound[p] = bound.lb; + probing_implied_bound.zero_upper_bound[p] = bound.ub; + } + for (auto& [imp_var, bound] : entries[1].var_to_cached_bound_map) { + i_t i = remap(imp_var); + if (i < 0 || i >= num_cols) { continue; } + if (!bb_bounds_consistent(i, bound.lb, bound.ub)) { continue; } + i_t p = one_cursor[j]++; + probing_implied_bound.one_variables[p] = i; + probing_implied_bound.one_lower_bound[p] = bound.lb; + probing_implied_bound.one_upper_bound[p] = bound.ub; + } + } + + CUOPT_LOG_INFO("Probing implied bounds: %d zero entries, %d one entries", zero_nnz, one_nnz); +} + template solution_t mip_solver_t::run_solver() { @@ -210,6 +308,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_bound_t probing_implied_bound; + bool run_bb = !context.settings.heuristics_only; if (run_bb) { // Convert the presolved problem to dual_simplex::user_problem_t @@ -217,6 +317,11 @@ 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_implied_bounds(op_problem_, + branch_and_bound_problem, + dm.ls.constraint_prop.bounds_update.probing_cache, + probing_implied_bound); + // 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; @@ -237,8 +342,9 @@ solution_t mip_solver_t::run_solver() } branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; - branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; - branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; + branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; + branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; + branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; branch_and_bound_settings.reduced_cost_strengthening = @@ -287,6 +393,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_problem, branch_and_bound_settings, timer_.get_tic_start(), + probing_implied_bound, context.problem_ptr->clique_table); context.branch_and_bound_ptr = branch_and_bound.get(); diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index c3d671954..33b5457df 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -969,6 +969,55 @@ TEST(cuts, test_cuts_2) EXPECT_EQ(solution.get_num_nodes(), 0); } +TEST(cuts, test_duplicate_cuts_detection) +{ + dual_simplex::simplex_solver_settings_t settings; + dual_simplex::cut_pool_t cut_pool(4, settings); + dual_simplex::inequality_t cut1; + cut1.push_back(0, 1.0); + cut1.push_back(1, 2.0); + cut1.rhs = 1.0; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut1); + dual_simplex::inequality_t cut2; + cut2.push_back(0, 2.0); + cut2.push_back(1, 4.0); + cut2.rhs = 2.0; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut2); + dual_simplex::inequality_t cut3; + cut3.push_back(0, 0.1); + cut3.push_back(2, 0.2); + cut3.rhs = 1.0; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut3); + dual_simplex::inequality_t cut4; + cut4.push_back(0, 0.2); + cut4.push_back(2, 0.4); + cut4.rhs = 1.0; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut4); + dual_simplex::inequality_t cut5; + cut5.push_back(1, 10.0); + cut5.push_back(3, 20.0); + cut5.rhs = 0.1; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut5); + dual_simplex::inequality_t cut6; + cut6.push_back(1, 20.0); + cut6.push_back(3, 40.0); + cut6.rhs = 0.2; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut6); + dual_simplex::inequality_t cut7; + cut7.push_back(0, 1.0); + cut7.push_back(1, 1.0); + cut7.push_back(2, 1.0); + cut7.push_back(3, 1.0); + cut7.rhs = 1.0; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut7); + dual_simplex::inequality_t cut8; + cut8.push_back(1, 3.0); + cut8.rhs = 7.0; + cut_pool.add_cut(dual_simplex::cut_type_t::MIXED_INTEGER_GOMORY, cut8); + + cut_pool.check_for_duplicate_cuts(); +} + TEST(cuts, clique_phase1_smoke_conflict_graph_edges) { const raft::handle_t handle{};