diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 96396e4efc..eca331c661 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -131,7 +131,6 @@ #define CUOPT_MIP_HYPER_DIVING_FARKAS "mip_hyper_diving_farkas" #define CUOPT_MIP_HYPER_DIVING_VECTOR_LENGTH "mip_hyper_diving_vector_length" /* @brief Diving heuristic limits */ -#define CUOPT_MIP_HYPER_DIVING_MIN_NODE_DEPTH "mip_hyper_diving_min_node_depth" #define CUOPT_MIP_HYPER_DIVING_NODE_LIMIT "mip_hyper_diving_node_limit" #define CUOPT_MIP_HYPER_DIVING_ITERATION_LIMIT_FACTOR "mip_hyper_diving_iteration_limit_factor" #define CUOPT_MIP_HYPER_DIVING_BACKTRACK_LIMIT "mip_hyper_diving_backtrack_limit" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp b/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp index 8ae1eda2a2..e61f25cec4 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp @@ -6,6 +6,7 @@ /* clang-format on */ #pragma once +#include namespace cuopt::mathematical_optimization { @@ -22,16 +23,13 @@ struct mip_diving_hyper_params_t { i_t farkas_diving = -1; i_t vector_length_diving = -1; - // The minimum depth to start diving from. - i_t min_node_depth = 10; - // The maximum number of nodes when performing a dive. i_t node_limit = 500; // The maximum number of dual simplex iteration allowed - // in a single dive. This set in terms of the total number of - // iterations in the best-first threads. - f_t iteration_limit_factor = 0.05; + // in a single dive. + f_t iteration_limit_factor = 0.05; + int64_t iteration_limit_offset = 10000; // The maximum backtracking allowed. i_t backtrack_limit = 5; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e2ea33bb6f..6756205e92 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -426,14 +426,19 @@ void branch_and_bound_t::report_heuristic(f_t obj, heuristics_origin_t } template -void branch_and_bound_t::report( - char symbol, f_t obj, f_t lower_bound, i_t node_depth, i_t node_int_infeas, double work_time) +void branch_and_bound_t::report(const lp_problem_t& lp, + char symbol, + f_t obj, + f_t lower_bound, + i_t node_depth, + i_t node_int_infeas, + double work_time) { - update_user_bound(lower_bound); + update_user_bound(lp, lower_bound); const i_t nodes_explored = exploration_stats_.nodes_explored; const i_t nodes_unexplored = exploration_stats_.nodes_unexplored; - const f_t user_obj = compute_user_objective(original_lp_, obj); - const f_t user_lower = compute_user_objective(original_lp_, lower_bound); + const f_t user_obj = compute_user_objective(lp, obj); + const f_t user_lower = compute_user_objective(lp, lower_bound); const f_t iters = static_cast(exploration_stats_.total_simplex_iters); const f_t iter_node = nodes_explored > 0 ? iters / nodes_explored : iters; f_t user_gap = user_relative_gap(user_obj, user_lower); @@ -516,10 +521,11 @@ i_t branch_and_bound_t::find_reduced_cost_fixings(f_t upper_bound, } template -void branch_and_bound_t::update_user_bound(f_t lower_bound) +void branch_and_bound_t::update_user_bound(const lp_problem_t& lp, + f_t lower_bound) { if (user_bound_callback_ == nullptr) { return; } - f_t user_lower = compute_user_objective(original_lp_, lower_bound); + f_t user_lower = compute_user_objective(lp, lower_bound); user_bound_callback_(user_lower); } @@ -932,7 +938,8 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& } template -void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, +void branch_and_bound_t::add_feasible_solution(const lp_problem_t& lp, + f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, search_strategy_t thread_type) @@ -940,7 +947,7 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, bool send_solution = false; settings_.log.debug("%c found a feasible solution with obj=%.10e.\n", feasible_solution_symbol(thread_type, settings_.diving_settings.show_type), - compute_user_objective(original_lp_, leaf_objective)); + compute_user_objective(lp, leaf_objective)); mutex_upper_.lock(); if (!incumbent_.has_incumbent || leaf_objective < incumbent_.objective) { @@ -948,13 +955,13 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, upper_bound_ = std::min(upper_bound_.load(), leaf_objective); char symbol = feasible_solution_symbol(thread_type, settings_.diving_settings.show_type); - report(symbol, leaf_objective, get_lower_bound(), leaf_depth, 0); + report(lp, symbol, leaf_objective, get_lower_bound(), leaf_depth, 0); send_solution = true; } if (send_solution && settings_.solution_callback != nullptr) { std::vector original_x; - uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, original_x); + uncrush_primal_solution(original_problem_, lp, incumbent_.x, original_x); settings_.solution_callback(original_x, leaf_objective); } mutex_upper_.unlock(); @@ -991,24 +998,23 @@ branch_variable_t branch_and_bound_t::variable_selection( log.log = false; i_t branch_var = -1; branch_direction_t round_dir = branch_direction_t::NONE; - std::vector current_incumbent; - std::vector& solution = worker->leaf_solution.x; + std::vector& solution = worker->leaf_solution.x; switch (worker->search_strategy) { case search_strategy_t::BEST_FIRST: if (settings_.reliability_branching != 0) { - branch_var = pc_.reliable_variable_selection(node_ptr, - fractional, - worker, - var_types_, - exploration_stats_, - upper_bound_, - bfs_worker_pool_.num_idle(), - new_slacks_, - original_lp_); + branch_var = worker->pseudo_costs.reliable_variable_selection(node_ptr, + fractional, + worker, + var_types_, + exploration_stats_, + upper_bound_, + bfs_worker_pool_.num_idle(), + new_slacks_, + original_lp_); } else { - branch_var = pc_.variable_selection(fractional, solution); + branch_var = worker->pseudo_costs.variable_selection(fractional, solution); } round_dir = martin_criteria(solution[branch_var], worker->root_solution[branch_var]); @@ -1016,21 +1022,29 @@ branch_variable_t branch_and_bound_t::variable_selection( return {branch_var, round_dir}; case search_strategy_t::COEFFICIENT_DIVING: - return coefficient_diving( - original_lp_, fractional, solution, var_up_locks_, var_down_locks_, log); + if (worker->var_up_locks.empty()) { + calculate_variable_locks( + worker->leaf_problem, worker->var_up_locks, worker->var_down_locks); + } + + return coefficient_diving(worker->leaf_problem, + fractional, + solution, + worker->var_up_locks, + worker->var_down_locks, + log); case search_strategy_t::LINE_SEARCH_DIVING: return line_search_diving(fractional, solution, worker->root_solution, log); case search_strategy_t::PSEUDOCOST_DIVING: - return pseudocost_diving(pc_, fractional, solution, worker->root_solution, log); + return pseudocost_diving( + worker->pseudo_costs, fractional, solution, worker->root_solution, log); case search_strategy_t::GUIDED_DIVING: - assert(incumbent_.has_incumbent); - mutex_upper_.lock(); - current_incumbent = incumbent_.x; - mutex_upper_.unlock(); - return guided_diving(pc_, fractional, solution, current_incumbent, log); + assert(!worker->current_incumbent.empty()); + return guided_diving( + worker->pseudo_costs, fractional, solution, worker->current_incumbent, log); case search_strategy_t::FARKAS_DIVING: return farkas_diving(worker->leaf_problem, fractional, solution, settings_.zero_tol, log); @@ -1040,7 +1054,7 @@ branch_variable_t branch_and_bound_t::variable_selection( case search_strategy_t::RINS: // This is used for solving the DFS of the sub-MIP. case search_strategy_t::RENS: - branch_var = pc_.variable_selection(fractional, solution); + branch_var = worker->pseudo_costs.variable_selection(fractional, solution); round_dir = martin_criteria(solution[branch_var], worker->root_solution[branch_var]); return {branch_var, round_dir}; } @@ -1094,14 +1108,14 @@ struct nondeterministic_policy_t : tree_update_policy_t { void update_pseudo_costs(mip_node_t* node, f_t leaf_obj) override { - bnb.pc_.update_pseudo_costs(node, leaf_obj); + worker->pseudo_costs.update_pseudo_costs(node, leaf_obj); } void handle_integer_solution(mip_node_t* node, f_t obj, const std::vector& x) override { - bnb.add_feasible_solution(obj, x, node->depth, worker->search_strategy); + bnb.add_feasible_solution(worker->leaf_problem, obj, x, node->depth, worker->search_strategy); } branch_variable_t select_branch_variable(mip_node_t* node, @@ -1116,7 +1130,8 @@ struct nondeterministic_policy_t : tree_update_policy_t { const std::vector& x) override { if (worker->search_strategy == search_strategy_t::BEST_FIRST) { - node->objective_estimate = bnb.pc_.obj_estimate(fractional, x, node->lower_bound); + node->objective_estimate = + worker->pseudo_costs.obj_estimate(fractional, x, node->lower_bound); } } @@ -1291,11 +1306,16 @@ struct deterministic_diving_policy_t } case search_strategy_t::COEFFICIENT_DIVING: { + if (this->worker.var_up_locks.empty()) { + calculate_variable_locks( + this->worker.leaf_problem, this->worker.var_up_locks, this->worker.var_down_locks); + } + return coefficient_diving(this->worker.leaf_problem, fractional, x, - this->bnb.var_up_locks_, - this->bnb.var_down_locks_, + this->worker.var_up_locks, + this->worker.var_down_locks, log); } @@ -1345,15 +1365,17 @@ struct deterministic_diving_policy_t // We use the lower bound to decide if we should fathom the // node or branch. template -void branch_and_bound_t::snap_to_lattice(mip_node_t* node_ptr, f_t leaf_obj) +void branch_and_bound_t::snap_to_lattice(const lp_problem_t& lp, + mip_node_t* node_ptr, + f_t leaf_obj) { - if (original_lp_.objective_step.has_step()) { - f_t step = original_lp_.objective_step.step_size; - f_t bias = original_lp_.objective_step.bias; + if (lp.objective_step.has_step()) { + f_t step = lp.objective_step.step_size; + f_t bias = lp.objective_step.bias; // Round up to next value on the lattice: k * step + bias >= leaf_obj f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); node_ptr->lower_bound = k * step + bias; - } else if (original_lp_.objective_is_integral) { + } else if (lp.objective_is_integral) { node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); } } @@ -1392,7 +1414,8 @@ std::pair branch_and_bound_t::updat } else if (lp_status == dual_status_t::OPTIMAL) { std::vector leaf_fractional; - i_t num_frac = fractional_variables(settings_, leaf_solution.x, var_types_, leaf_fractional); + i_t num_frac = + fractional_variables(settings_, leaf_solution.x, worker->var_types, leaf_fractional); #ifdef DEBUG_FRACTIONAL_FIXED for (i_t j : leaf_fractional) { @@ -1415,7 +1438,7 @@ std::pair branch_and_bound_t::updat policy.graphviz(search_tree, node_ptr, "lower bound", leaf_obj); policy.update_pseudo_costs(node_ptr, leaf_obj); node_ptr->lower_bound = leaf_obj; - snap_to_lattice(node_ptr, leaf_obj); + snap_to_lattice(leaf_problem, node_ptr, leaf_obj); if (num_frac == 0) { policy.handle_integer_solution(node_ptr, leaf_obj, leaf_solution.x); @@ -1530,7 +1553,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( branch_and_bound_worker_t* worker, branch_and_bound_stats_t& stats, logger_t& log, - i_t iter_limit) + int64_t iter_limit) { raft::common::nvtx::range scope("BB::solve_node"); #ifdef DEBUG_BRANCHING @@ -1572,13 +1595,13 @@ dual_status_t branch_and_bound_t::solve_node_lp( lp_settings.concurrent_halt = &node_concurrent_halt_; lp_settings.set_log(false); f_t cutoff = upper_bound_.load(); - if (original_lp_.objective_step.has_step()) { - f_t step = original_lp_.objective_step.step_size; - f_t bias = original_lp_.objective_step.bias; + if (worker->leaf_problem.objective_step.has_step()) { + f_t step = worker->leaf_problem.objective_step.step_size; + f_t bias = worker->leaf_problem.objective_step.bias; // Any improving feasible solution must have objective <= cutoff - step. f_t k = std::floor((cutoff - bias) / step + settings_.integer_tol); lp_settings.cut_off = (k - 1) * step + bias + settings_.dual_tol; - } else if (original_lp_.objective_is_integral) { + } else if (worker->leaf_problem.objective_is_integral) { // If the objective is integral, any feasible solution should produce an upper bound that is // (approximately) integral. We add a small tolerance and floor this value to get an integer, // we then subtract 1, to stop simplex on problems that cannot improve the primal objective. @@ -1590,7 +1613,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; } lp_settings.scale_columns = false; - lp_settings.iteration_limit = iter_limit; + lp_settings.iteration_limit = std::min(iter_limit, std::numeric_limits::max()); #ifdef LOG_NODE_SIMPLEX lp_settings.set_log(true); @@ -1689,10 +1712,10 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; - f_t user_obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, lower_bound); + f_t user_obj = compute_user_objective(worker->leaf_problem, upper_bound); + f_t user_lower = compute_user_objective(worker->leaf_problem, lower_bound); f_t rel_gap = user_relative_gap(user_obj, user_lower); - f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); + f_t abs_gap = compute_user_abs_gap(worker->leaf_problem, upper_bound, lower_bound); bool can_launch_new_submip = true; @@ -1742,7 +1765,12 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, if (((nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && time_since_last_log >= 1) || (time_since_last_log > 30) || now > settings_.time_limit) { - report(' ', upper_bound_, lower_bound, node_ptr->depth, node_ptr->integer_infeasible); + report(worker->leaf_problem, + ' ', + upper_bound_, + lower_bound, + node_ptr->depth, + node_ptr->integer_infeasible); exploration_stats_.last_log = tic(); exploration_stats_.nodes_since_last_log = 0; } @@ -1849,10 +1877,10 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, lower_bound = get_lower_bound(); upper_bound = upper_bound_; - user_obj = compute_user_objective(original_lp_, upper_bound); - user_lower = compute_user_objective(original_lp_, lower_bound); + user_obj = compute_user_objective(worker->leaf_problem, upper_bound); + user_lower = compute_user_objective(worker->leaf_problem, lower_bound); rel_gap = user_relative_gap(user_obj, user_lower); - abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); + abs_gap = compute_user_abs_gap(worker->leaf_problem, upper_bound, lower_bound); } if (solver_status_ == mip_status_t::TIME_LIMIT || solver_status_ == mip_status_t::OPTIMAL) { @@ -1927,9 +1955,9 @@ template void branch_and_bound_t::best_first_search_with(bfs_worker_t* worker) { f_t lower_bound = get_lower_bound(); - f_t user_obj = compute_user_objective(original_lp_, upper_bound_.load()); - f_t user_lower = compute_user_objective(original_lp_, lower_bound); - f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound); + f_t user_obj = compute_user_objective(worker->leaf_problem, upper_bound_.load()); + f_t user_lower = compute_user_objective(worker->leaf_problem, lower_bound); + f_t abs_gap = compute_user_abs_gap(worker->leaf_problem, upper_bound_.load(), lower_bound); f_t rel_gap = user_relative_gap(user_obj, user_lower); f_t steal_chance = settings_.bnb_steal_chance >= 0 ? settings_.bnb_steal_chance : MIP_DEFAULT_STEAL_CHANCE; @@ -1942,12 +1970,13 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t if (diving_settings.farkas_diving != 0) { f_t obj_dyn; - if (std::abs(original_lp_.min_abs_obj_coeff) < settings_.zero_tol) { - obj_dyn = std::abs(original_lp_.max_abs_obj_coeff) < settings_.zero_tol + if (std::abs(worker->leaf_problem.min_abs_obj_coeff) < settings_.zero_tol) { + obj_dyn = std::abs(worker->leaf_problem.max_abs_obj_coeff) < settings_.zero_tol ? 0 : std::numeric_limits::infinity(); } else { - obj_dyn = std::log10(original_lp_.max_abs_obj_coeff / original_lp_.min_abs_obj_coeff); + obj_dyn = + std::log10(worker->leaf_problem.max_abs_obj_coeff / worker->leaf_problem.min_abs_obj_coeff); } if (obj_dyn < diving_settings.farkas_obj_dynamism_tol) { diving_settings.farkas_diving = 0; } } @@ -1982,6 +2011,12 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t break; } + if (received_halt_signal()) { + solver_status_ = mip_status_t::HALT; + node_concurrent_halt_ = true; + break; + } + // If the guided diving was disabled previously due to the lack of an incumbent solution, // re-enable as soon as a new incumbent is found. if (diving_worker_pool_.size() > 0 && settings_.diving_settings.guided_diving != 0 && @@ -2015,9 +2050,9 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t plunge_with(worker, start_node); lower_bound = get_lower_bound(); - user_obj = compute_user_objective(original_lp_, upper_bound_.load()); - user_lower = compute_user_objective(original_lp_, lower_bound); - abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound); + user_obj = compute_user_objective(worker->leaf_problem, upper_bound_.load()); + user_lower = compute_user_objective(worker->leaf_problem, lower_bound); + abs_gap = compute_user_abs_gap(worker->leaf_problem, upper_bound_.load(), lower_bound); rel_gap = user_relative_gap(user_obj, user_lower); if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) { @@ -2051,14 +2086,16 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t } template -void branch_and_bound_t::dive_with(diving_worker_t* worker, i_t backtrack_limit) +void branch_and_bound_t::dive_with(diving_worker_t* worker, + const simplex_solver_settings_t& settings) { raft::common::nvtx::range scope("BB::diving_thread"); if (worker->orbital_fixing) { worker->orbital_fixing->disable(); } logger_t log; log.log = false; - const i_t diving_node_limit = settings_.diving_settings.node_limit; + const i_t diving_node_limit = settings.diving_settings.node_limit; + const i_t backtrack_limit = settings.diving_settings.backtrack_limit; worker->recompute_basis = true; worker->recompute_bounds = true; @@ -2073,13 +2110,14 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker, branch_and_bound_stats_t dive_stats; f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; - f_t user_obj = compute_user_objective(original_lp_, upper_bound); - f_t user_lower = compute_user_objective(original_lp_, lower_bound); + f_t user_obj = compute_user_objective(worker->leaf_problem, upper_bound); + f_t user_lower = compute_user_objective(worker->leaf_problem, lower_bound); f_t rel_gap = user_relative_gap(user_obj, user_lower); - f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); + f_t abs_gap = compute_user_abs_gap(worker->leaf_problem, 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) { + rel_gap > settings.relative_mip_gap_tol && abs_gap > settings.absolute_mip_gap_tol && + !(settings.concurrent_halt && settings.concurrent_halt->load(std::memory_order_acquire))) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); @@ -2100,8 +2138,8 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker, int64_t bnb_lp_iters = exploration_stats_.total_simplex_iters; f_t factor = settings_.diving_settings.iteration_limit_factor; - i_t max_iter = std::min(factor * bnb_lp_iters - dive_stats.total_simplex_iters, - std::numeric_limits::max()); + int64_t offset = settings_.diving_settings.iteration_limit_offset; + int64_t max_iter = offset + factor * bnb_lp_iters - dive_stats.total_simplex_iters; if (max_iter <= 0) { break; } decompress_vstatus( @@ -2142,16 +2180,16 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker, lower_bound = get_lower_bound(); upper_bound = upper_bound_; - user_obj = compute_user_objective(original_lp_, upper_bound); - user_lower = compute_user_objective(original_lp_, lower_bound); + user_obj = compute_user_objective(worker->leaf_problem, upper_bound); + user_lower = compute_user_objective(worker->leaf_problem, lower_bound); rel_gap = user_relative_gap(user_obj, user_lower); - abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); + abs_gap = compute_user_abs_gap(worker->leaf_problem, upper_bound, lower_bound); } // This is called from the RINS method which already handle the return to the // pool part. Besides, they do not share the same pool. if (worker->search_strategy != search_strategy_t::RINS && - worker->search_strategy != search_strategy_t::RENS) { + worker->search_strategy != search_strategy_t::RENS && !settings.inside_root_node) { diving_worker_pool_.return_worker_to_pool(worker); } } @@ -2175,8 +2213,7 @@ bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* return false; } - if (upper_bound_.load() < diving_worker->start_node.lower_bound || - diving_worker->start_node.depth < settings_.diving_settings.min_node_depth) { + if (upper_bound_.load() < diving_worker->start_node.lower_bound) { diving_worker_pool_.return_worker_to_pool(diving_worker); return false; } @@ -2196,6 +2233,13 @@ bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* auto strategy = bfs_worker->next_diving_heuristic(); diving_worker->search_strategy = strategy; diving_worker->bfs_worker = bfs_worker; + + if (strategy == search_strategy_t::GUIDED_DIVING) { + mutex_upper_.lock(); + diving_worker->current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + } + diving_worker->set_active(); ++bfs_worker->active_diving_workers; @@ -2203,7 +2247,7 @@ bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* #pragma omp task affinity(*diving_worker) priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) \ firstprivate(diving_worker) - dive_with(diving_worker, settings_.diving_settings.backtrack_limit); + dive_with(diving_worker, settings_); return true; } @@ -2218,10 +2262,9 @@ bool branch_and_bound_t::launch_submip_worker(const std::vector& diving_worker_t* worker = submip_worker_pool_.pop_idle_worker(); if (!worker) return false; - std::vector current_incumbent; mutex_upper_.lock(); bool use_rins = incumbent_.has_incumbent && settings_.submip_settings.rins != 0; - if (use_rins) current_incumbent = incumbent_.x; + if (use_rins) worker->current_incumbent = incumbent_.x; mutex_upper_.unlock(); // Note that this node does not have the vstatus (it was cleared at the start of B&B exploration) @@ -2236,11 +2279,10 @@ bool branch_and_bound_t::launch_submip_worker(const std::vector& if (settings_.inside_submip) { // LLVM libomp's GOMP compatibility path skips GCC's firstprivate copy // function for included tasks. - recursive_submip(worker, current_incumbent, var_types_, settings_); + recursive_submip(worker, settings_); } else { -#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ - firstprivate(worker, current_incumbent) - recursive_submip(worker, current_incumbent, var_types_, settings_); +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) firstprivate(worker) + recursive_submip(worker, settings_); } return true; @@ -2248,8 +2290,6 @@ bool branch_and_bound_t::launch_submip_worker(const std::vector& template void branch_and_bound_t::solve_submip(diving_worker_t* worker, - const std::vector& current_incumbent, - const std::vector& var_types, submip_stats_t& submip_stats, f_t fixrate, i_t simplex_iter_used, @@ -2315,7 +2355,8 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke // there is only equality rows (the range row vector is empty) and it contains // structural + slacks + cuts constraints/variables. user_problem_t submip_problem(original_problem_.handle_ptr); - simplex::convert_lp_to_user_problem(worker->leaf_problem, var_types, settings_, submip_problem); + simplex::convert_lp_to_user_problem( + worker->leaf_problem, worker->var_types, settings_, submip_problem); third_party_presolve_t presolver; f_t presolve_time_limit = std::min(0.1 * submip_settings.time_limit, 60.0); @@ -2375,11 +2416,11 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke std::vector presolved_incumbent; // We do not have an incumbent yet, so skip the initial guess. - if (!current_incumbent.empty()) { + if (!worker->current_incumbent.empty()) { // Crush the incumbent to presolve space. It may not be valid for the sub-MIP since we // may fix integer variables that does not match the current incumbent to reach the target // fix rate. - presolver.crush_primal_solution(submip_problem, current_incumbent, presolved_incumbent); + presolver.crush_primal_solution(submip_problem, worker->current_incumbent, presolved_incumbent); submip_bnb.set_initial_guess(presolved_incumbent); } @@ -2652,10 +2693,7 @@ f_t calculate_fixrate(const std::vector& integer_list, template void branch_and_bound_t::recursive_submip( - diving_worker_t* worker, - const std::vector& current_incumbent, - const std::vector& var_types, - simplex_solver_settings_t submip_settings) + diving_worker_t* worker, simplex_solver_settings_t submip_settings) { raft::common::nvtx::range scope("BB::submip_thread"); if (worker->orbital_fixing) { worker->orbital_fixing->disable(); } @@ -2687,10 +2725,10 @@ void branch_and_bound_t::recursive_submip( std::fill(bounds_changed.begin(), bounds_changed.end(), false); std::vector fractional; - i_t num_frac = fractional_variables(settings_, current_sol, var_types, fractional); + i_t num_frac = fractional_variables(settings_, current_sol, worker->var_types, fractional); std::vector integer_list; - get_unfixed_integer_variables(lower, upper, var_types, settings_.fixed_tol, integer_list); + get_unfixed_integer_variables(lower, upper, worker->var_types, settings_.fixed_tol, integer_list); i_t num_integers = integer_list.size(); f_t max_fixrate = submip_get_max_fixrate(submip_stats, settings_.submip_settings, worker->rng); @@ -2719,7 +2757,7 @@ void branch_and_bound_t::recursive_submip( num_bound_changed = apply_rins_fixings(settings_, current_sol, integer_list, - current_incumbent, + worker->current_incumbent, max_fixrate - prev_fixrate, lower, upper, @@ -2823,8 +2861,7 @@ void branch_and_bound_t::recursive_submip( int64_t simplex_iter = exploration_stats_.total_simplex_iters; f_t iter_ratio = settings_.submip_settings.iteration_limit_ratio; int64_t simplex_iter_limit = iter_offset + simplex_iter * iter_ratio; - i_t max_iter = std::min(simplex_iter_limit - stats.total_simplex_iters, - std::numeric_limits::max()); + int64_t max_iter = simplex_iter_limit - stats.total_simplex_iters; if (max_iter <= 0) { DEBUG_SUBMIP("{}Round {}: max iteration reached! {}/{}", submip_settings.log.log_prefix, @@ -2844,12 +2881,12 @@ void branch_and_bound_t::recursive_submip( } fractional.clear(); - num_frac = fractional_variables(settings_, current_sol, var_types, fractional); + num_frac = fractional_variables(settings_, current_sol, worker->var_types, fractional); f_t leaf_obj = compute_objective(worker->leaf_problem, current_sol); node.lower_bound = leaf_obj; - snap_to_lattice(&node, leaf_obj); + snap_to_lattice(worker->leaf_problem, &node, leaf_obj); if (leaf_obj > upper_bound_.load()) { DEBUG_SUBMIP("{}Round {}: reached cutoff point. obj={:.4g}. upper_bound={:.4g}", submip_settings.log.log_prefix, @@ -2861,7 +2898,8 @@ void branch_and_bound_t::recursive_submip( if (num_frac == 0) { // We found a feasible solution when fixing the variables in RINS/RENS. - add_feasible_solution(leaf_obj, current_sol, -1, worker->search_strategy); + add_feasible_solution( + worker->leaf_problem, leaf_obj, current_sol, -1, worker->search_strategy); DEBUG_SUBMIP("{}Round {}: found a solution with obj={:.4g}. upper_bound={:.4g}", submip_settings.log.log_prefix, round, @@ -2904,7 +2942,7 @@ void branch_and_bound_t::recursive_submip( f_t work_limit = 1.0; submip_fj_cpu_worker.create_worker( worker->leaf_problem, - var_types, + worker->var_types, worker->leaf_solution.x, settings_, std::format("{} [CPU FJ]", submip_settings.log.log_prefix), @@ -2912,25 +2950,19 @@ void branch_and_bound_t::recursive_submip( submip_fj_cpu_worker.run_sync(time_limit, work_limit); } - // We need the pseudocost to do the DFS, which we do not have during the cut passes. - if (!submip_settings.inside_root_node) { - DEBUG_SUBMIP("{}Running a quick DFS. fixrate={:.4g} ({}/{})", - submip_settings.log.log_prefix, - fixrate, - fixrate * num_integers, - num_integers); - dive_with(worker, settings_.submip_settings.dfs_max_backtrack); - } + DEBUG_SUBMIP("{}Running a quick DFS. fixrate={:.4g} ({}/{})", + submip_settings.log.log_prefix, + fixrate, + fixrate * num_integers, + num_integers); + + simplex_solver_settings_t dfs_settings = submip_settings; + dfs_settings.diving_settings.backtrack_limit = settings_.submip_settings.dfs_max_backtrack; + dive_with(worker, dfs_settings); } } else { - solve_submip(worker, - current_incumbent, - var_types, - submip_stats, - fixrate, - stats.total_simplex_iters, - submip_settings); + solve_submip(worker, submip_stats, fixrate, stats.total_simplex_iters, submip_settings); } } @@ -2958,7 +2990,11 @@ void branch_and_bound_t::recursive_submip( template void branch_and_bound_t::launch_root_heuristics( const lp_problem_t& lp, - const std::vector& sol, + const lp_solution_t& lp_solution, + const std::vector& fractional, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factor, i_t cut_pass, root_heuristics_t& root_heuristics) { @@ -2967,13 +3003,17 @@ void branch_and_bound_t::launch_root_heuristics( // Using shared_ptr here, so the lifetime of the object is tied to the related task. This allows // the solver to send the stop signal and immediately continue the execution. - auto current_heuristic = - root_heuristics.create_new_cut_pass_heuristic(cut_pass, Arow_, var_types_, sol, edge_norms_); - auto worker_count = root_heuristics.worker_count_; - constexpr bool is_root_heuristic = true; - constexpr bool is_cpufj_enabled = true; + auto current_heuristic = root_heuristics.create_new_cut_pass_heuristic( + Arow_, var_types_, lp_solution.x, edge_norms_, settings_); + auto worker_count = root_heuristics.worker_count_; + + current_heuristic->initialize_pseudocost( + lp, root_vstatus_, fractional, lp_solution, basic_list, nonbasic_list, basis_factor); + constexpr bool is_cpufj_enabled = true; if (is_cpufj_enabled) { + root_heuristics.stop_old_workers(cut_pass, 1); + f_t work_limit = std::numeric_limits::infinity(); f_t time_limit = settings_.time_limit - toc(exploration_stats_.start_time); @@ -2982,8 +3022,9 @@ void branch_and_bound_t::launch_root_heuristics( set_solution_from_cpu_fj(obj, assignment, work_units); }; current_heuristic->fj_cpu_worker_.create_worker( - lp, var_types_, sol, settings_, "[RootCut CPUFJ] "); + lp, var_types_, lp_solution.x, settings_, "[RootCut CPUFJ] "); ++(*worker_count); + ++current_heuristic->active_workers_; #pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) \ affinity(current_heuristic -> fj_cpu_worker_) firstprivate(current_heuristic, worker_count) \ @@ -2991,19 +3032,26 @@ void branch_and_bound_t::launch_root_heuristics( { current_heuristic->fj_cpu_worker_.run_sync(time_limit, work_limit); --(*worker_count); + --current_heuristic->active_workers_; } } + mutex_upper_.lock(); bool use_rins = settings_.submip_settings.rins != 0 && incumbent_.has_incumbent; + mutex_upper_.unlock(); + if (use_rins || settings_.submip_settings.rens != 0) { + root_heuristics.stop_old_workers(cut_pass, 1); + search_strategy_t strategy = use_rins ? search_strategy_t::RINS : search_strategy_t::RENS; diving_worker_t* worker = current_heuristic->create_submip_worker( - cut_pass, lp, settings_, root_objective_, root_vstatus_, sol, strategy); + cut_pass, lp, settings_, root_objective_, root_vstatus_, lp_solution.x, strategy); - std::vector current_incumbent; - mutex_upper_.lock(); - if (use_rins) current_incumbent = incumbent_.x; - mutex_upper_.unlock(); + if (use_rins) { + mutex_upper_.lock(); + worker->current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + } simplex_solver_settings_t submip_settings = settings_; submip_settings.concurrent_halt = ¤t_heuristic->halt_; @@ -3012,14 +3060,75 @@ void branch_and_bound_t::launch_root_heuristics( if (settings_.inside_submip) { // LLVM libomp's GOMP compatibility path skips GCC's firstprivate copy // function for included tasks. - recursive_submip(worker, current_incumbent, current_heuristic->var_types_, submip_settings); + recursive_submip(worker, submip_settings); } else { + ++current_heuristic->active_workers_; ++(*worker_count); -#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) firstprivate( \ - current_incumbent, current_heuristic, worker_count, submip_settings) depend(out : *worker) +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ + firstprivate(current_heuristic, worker_count, submip_settings) depend(out : *worker) { - recursive_submip(worker, current_incumbent, current_heuristic->var_types_, submip_settings); + recursive_submip(worker, submip_settings); --(*worker_count); + --current_heuristic->active_workers_; + } + } + } + + constexpr bool use_diving = true; + if (use_diving) { + mip_diving_hyper_params_t diving_settings = settings_.diving_settings; + + if (diving_settings.guided_diving != 0 && !has_solver_space_incumbent()) { + diving_settings.guided_diving = 0; + } + + if (diving_settings.farkas_diving != 0) { + f_t obj_dyn; + if (std::abs(lp.min_abs_obj_coeff) < settings_.zero_tol) { + obj_dyn = std::abs(lp.max_abs_obj_coeff) < settings_.zero_tol + ? 0 + : std::numeric_limits::infinity(); + } else { + obj_dyn = std::log10(lp.max_abs_obj_coeff / lp.min_abs_obj_coeff); + } + if (obj_dyn < diving_settings.farkas_obj_dynamism_tol) { diving_settings.farkas_diving = 0; } + } + + std::vector diving_heuristics; + get_diving_heuristic_list(diving_settings, diving_heuristics); + + i_t available = cut_pass == 0 ? settings_.num_threads - 3 : settings_.num_threads - 2; + i_t num_diving_workers = std::min(diving_heuristics.size(), available); + root_heuristics.stop_old_workers(cut_pass, num_diving_workers); + + mip_node_t root_node(root_objective_, root_vstatus_); + + for (i_t k = 0; k < num_diving_workers; ++k) { + i_t j = root_heuristics.next_diving_type_ % diving_heuristics.size(); + root_heuristics.next_diving_type_ = j + 1 % diving_heuristics.size(); + search_strategy_t strategy = diving_heuristics[j]; + + diving_worker_t* worker = + current_heuristic->create_diving_worker(cut_pass, lp, settings_, root_node, strategy); + + if (strategy == search_strategy_t::GUIDED_DIVING) { + mutex_upper_.lock(); + worker->current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + } + + ++current_heuristic->active_workers_; + ++(*worker_count); +#pragma omp task affinity(*worker) priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) \ + firstprivate(worker, current_heuristic, worker_count) depend(out : *worker) + { + simplex_solver_settings_t dive_settings = settings_; + dive_settings.concurrent_halt = ¤t_heuristic->halt_; + dive_settings.inside_root_node = true; + dive_settings.diving_settings.backtrack_limit = 1; + dive_with(worker, dive_settings); + --(*worker_count); + --current_heuristic->active_workers_; } } } @@ -3470,7 +3579,7 @@ auto branch_and_bound_t::do_cut_pass( mutex_upper_.unlock(); } f_t obj = upper_bound_.load(); - report(' ', obj, root_objective_, 0, num_fractional); + report(original_lp_, ' ', obj, root_objective_, 0, num_fractional); f_t user_obj = compute_user_objective(original_lp_, upper_bound_.load()); f_t user_lower = compute_user_objective(original_lp_, root_objective_); @@ -3771,7 +3880,14 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::OPTIMAL; } - launch_root_heuristics(original_lp_, root_relax_soln_.x, cut_pass, root_heuristics); + launch_root_heuristics(original_lp_, + root_relax_soln_, + fractional, + basic_list, + nonbasic_list, + basis_update, + cut_pass, + root_heuristics); cut_pass_action_t cut_pass_action = do_cut_pass(cut_pass, solution, @@ -3952,9 +4068,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.last_log = tic(); min_node_queue_size_ = 20; - if (settings_.diving_settings.coefficient_diving != 0) { - calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); - } print_table_header(); #pragma omp taskgroup @@ -3972,6 +4085,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut var_types_, symmetry_, settings_, + pc_, root_relax_soln_.x, edge_norms_); submip_worker_pool_.init(num_submip_workers, @@ -3980,6 +4094,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut var_types_, symmetry_, settings_, + pc_, root_relax_soln_.x, edge_norms_, num_bfs_workers); @@ -3991,6 +4106,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut var_types_, symmetry_, settings_, + pc_, root_relax_soln_.x, edge_norms_, num_bfs_workers + num_submip_workers); @@ -4186,18 +4302,21 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri deterministic_horizon_number_ = 0; deterministic_global_termination_status_ = mip_status_t::UNSET; - deterministic_workers_ = std::make_unique>( - num_bfs_workers, original_lp_, Arow, var_types_, settings_, root_relax_soln_.x, edge_norms_); + deterministic_workers_ = + std::make_unique>(num_bfs_workers, + original_lp_, + Arow, + var_types_, + settings_, + pc_, + root_relax_soln_.x, + edge_norms_); if (num_diving_workers > 0) { // Extract diving types from search_strategies (skip BEST_FIRST at index 0) std::vector diving_types; get_diving_heuristic_list(settings_.diving_settings, diving_types); - if (settings_.diving_settings.coefficient_diving != 0) { - calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); - } - if (!diving_types.empty()) { deterministic_diving_workers_ = std::make_unique>(num_diving_workers, @@ -4206,6 +4325,7 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri Arow, var_types_, settings_, + pc_, root_relax_soln_.x, edge_norms_); } @@ -4498,7 +4618,7 @@ void branch_and_bound_t::deterministic_sync_callback() f_t time_since_last_log = exploration_stats_.last_log == 0 ? 1.0 : toc(exploration_stats_.last_log); if (time_since_last_log >= 1) { - report(' ', upper_bound, lower_bound, 0, 0, deterministic_current_horizon_); + report(original_lp_, ' ', upper_bound, lower_bound, 0, 0, deterministic_current_horizon_); exploration_stats_.last_log = tic(); } @@ -4664,7 +4784,8 @@ void branch_and_bound_t::deterministic_process_worker_solutions( i_t nodes_unexplored = exploration_stats_.nodes_unexplored.load(); search_strategy_t worker_type = get_worker_type(pool, sol->worker_id); - report(feasible_solution_symbol(worker_type, settings_.diving_settings.show_type), + report(original_lp_, + feasible_solution_symbol(worker_type, settings_.diving_settings.show_type), sol->objective, deterministic_lower, sol->depth, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 17ebdafa33..3ce193a30f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -208,12 +208,6 @@ class branch_and_bound_t { std::vector new_slacks_; std::vector var_types_; - // Variable locks (see definition 3.3 from T. Achterberg, “Constraint Integer Programming,” - // PhD, Technischen Universität Berlin, Berlin, 2007. doi: 10.14279/depositonce-1634). - // Here we assume that the constraints are in the form `Ax = b, l <= x <= u`. - std::vector var_up_locks_; - std::vector var_down_locks_; - // Mutex for the original LP // The heuristics threads look at the original LP. But the main thread modifies the // size of the original LP by adding slacks for cuts. Heuristic threads should lock @@ -296,7 +290,8 @@ class branch_and_bound_t { void print_table_header(); void report_heuristic(f_t obj, heuristics_origin_t origin); - void report(char symbol, + void report(const simplex::lp_problem_t& lp, + char symbol, f_t obj, f_t lower_bound, i_t node_depth, @@ -333,14 +328,15 @@ class branch_and_bound_t { // Set the solution when found at the root node void set_solution_at_root(simplex::mip_solution_t& solution, const cut_info_t& cut_info); - void update_user_bound(f_t lower_bound); + void update_user_bound(const simplex::lp_problem_t& lp, f_t lower_bound); // Set the final solution. void set_final_solution(simplex::mip_solution_t& solution, f_t lower_bound); // Update the incumbent solution with the new feasible solution // found during branch and bound. - void add_feasible_solution(f_t leaf_objective, + void add_feasible_solution(const simplex::lp_problem_t& lp, + f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, search_strategy_t thread_type); @@ -351,7 +347,9 @@ class branch_and_bound_t { // Launch a new diving worker from a given best-first worker. bool launch_diving_worker(bfs_worker_t* bfs_worker); - void snap_to_lattice(mip_node_t* node_ptr, f_t leaf_obj); + void snap_to_lattice(const simplex::lp_problem_t& lp, + mip_node_t* node_ptr, + f_t leaf_obj); // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); @@ -369,7 +367,8 @@ class branch_and_bound_t { // Perform a deep dive in the subtree determined by the `start_node` in order // to find integer feasible solutions. - void dive_with(diving_worker_t* worker, i_t backtrack_limit); + void dive_with(diving_worker_t* worker, + const simplex::simplex_solver_settings_t& settings); // Launch a new RINS worker bool launch_submip_worker(const std::vector& sol); @@ -382,8 +381,6 @@ class branch_and_bound_t { // Solve the RINS sub-MIP. void solve_submip(diving_worker_t* worker, - const std::vector& current_incumbent, - const std::vector& var_types, submip_stats_t& submip_stats, f_t fixrate, i_t simplex_iter_used, @@ -391,12 +388,14 @@ class branch_and_bound_t { // Creates and solves the RINS/RENS sub-MIP. void recursive_submip(diving_worker_t* worker, - const std::vector& current_incumbent, - const std::vector& var_types, simplex::simplex_solver_settings_t submip_settings); void launch_root_heuristics(const simplex::lp_problem_t& lp, - const std::vector& sol, + const simplex::lp_solution_t& lp_solution, + const std::vector& fractional, + const std::vector& basic_list, + const std::vector& nonbasic_list, + simplex::basis_update_mpf_t& basis_factor, i_t cut_pass, root_heuristics_t& root_heuristics); @@ -405,7 +404,7 @@ class branch_and_bound_t { branch_and_bound_worker_t* worker, branch_and_bound_stats_t& stats, simplex::logger_t& log, - i_t iter_limit = std::numeric_limits::max()); + int64_t iter_limit = std::numeric_limits::max()); // Apply symmetry-based bound reductions (orbital fixing and, when // settings_.symmetry == 2, lexical reduction) to the current node. diff --git a/cpp/src/branch_and_bound/deterministic_workers.hpp b/cpp/src/branch_and_bound/deterministic_workers.hpp index 7c426452a8..fae259ac3f 100644 --- a/cpp/src/branch_and_bound/deterministic_workers.hpp +++ b/cpp/src/branch_and_bound/deterministic_workers.hpp @@ -89,10 +89,11 @@ class deterministic_worker_base_t : public branch_and_bound_worker_t { const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm, const std::string& context_name) - : base_t(id, original_lp, Arow, var_types, settings, root_solution, root_edge_norm), + : base_t(id, original_lp, Arow, var_types, settings, pc, root_solution, root_edge_norm), work_context(context_name), pc_snapshot(1, settings) { @@ -143,6 +144,7 @@ class deterministic_bfs_worker_t const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm) : base_t(id, @@ -150,6 +152,7 @@ class deterministic_bfs_worker_t Arow, var_types, settings, + pc, root_solution, root_edge_norm, "BB_Worker_" + std::to_string(id)) @@ -308,6 +311,7 @@ class deterministic_diving_worker_t const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm) : base_t(id, @@ -315,6 +319,7 @@ class deterministic_diving_worker_t Arow, var_types, settings, + pc, root_solution, root_edge_norm, "Diving_Worker_" + std::to_string(id)), @@ -423,13 +428,14 @@ class deterministic_bfs_worker_pool_t const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm) { this->workers_.reserve(num_workers); for (int i = 0; i < num_workers; ++i) { this->workers_.emplace_back( - i, original_lp, Arow, var_types, settings, root_solution, root_edge_norm); + i, original_lp, Arow, var_types, settings, pc, root_solution, root_edge_norm); } } @@ -461,6 +467,7 @@ class deterministic_diving_worker_pool_t const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm) { @@ -468,7 +475,7 @@ class deterministic_diving_worker_pool_t for (int i = 0; i < num_workers; ++i) { search_strategy_t type = diving_types[i % diving_types.size()]; this->workers_.emplace_back( - i, type, original_lp, Arow, var_types, settings, root_solution, root_edge_norm); + i, type, original_lp, Arow, var_types, settings, pc, root_solution, root_edge_norm); } } diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index cdba90f219..c4071bf3b8 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -1369,6 +1369,69 @@ void strong_branching(const lp_problem_t& original_lp, } } +template +void pseudo_costs_t::initialize_with_estimate( + const lp_problem_t& lp, + const std::vector& vstatus, + const std::vector& fractional, + const lp_solution_t& lp_solution, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors) +{ + i_t m = lp.num_rows; + i_t n = lp.num_cols; + + std::vector delta_z(n, 0); + std::vector workspace(n, 0); + + f_t work_estimate = 0; + + std::vector basic_map(n, -1); + for (i_t i = 0; i < m; i++) { + basic_map[basic_list[i]] = i; + } + + // compute_initial_nonbasic_end permutes columns in place; copy so pc.Arow is unchanged + csr_matrix_t local_Arow = Arow; + + std::vector nonbasic_end(m); + compute_initial_nonbasic_end(basic_map, local_Arow, nonbasic_end); + + for (i_t k = 0; k < fractional.size(); k++) { + const i_t j = fractional[k]; + assert(j >= 0); + + if (pseudo_cost_num_up[j] == 0 || pseudo_cost_num_down[j] == 0) { + objective_change_estimate_t estimate = + single_pivot_objective_change_estimate(lp, + settings, + local_Arow, + vstatus, + j, + basic_map[j], + lp_solution, + basic_list, + nonbasic_list, + nonbasic_end, + basis_factors, + workspace, + delta_z, + work_estimate); + + if (pseudo_cost_num_down[j] == 0) { + pseudo_cost_sum_down[j] += estimate.down_obj_change; + ++pseudo_cost_num_down[j]; + } + + if (pseudo_cost_num_up[j] == 0) { + pseudo_cost_sum_up[j] += estimate.up_obj_change; + ++pseudo_cost_num_up[j]; + } + } + } +} + template inline f_t pseudo_costs_t::compute_pseudocost_average_down() { diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 55fc1fff35..a88b8bb3e1 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -216,6 +216,14 @@ class pseudo_costs_t { f_t compute_pseudocost_average_down(); f_t compute_pseudocost_average_up(); + void initialize_with_estimate(const simplex::lp_problem_t& lp, + const std::vector& vstatus, + const std::vector& fractional, + const simplex::lp_solution_t& lp_solution, + const std::vector& basic_list, + const std::vector& nonbasic_list, + simplex::basis_update_mpf_t& basis_factors); + f_t obj_estimate(const std::vector& fractional, const std::vector& solution, f_t lower_bound); diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index aa0ef6bb88..0ec0f74bf9 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -22,6 +22,9 @@ namespace cuopt::mathematical_optimization::mip { +template +class pseudo_costs_t; + template struct branch_and_bound_stats_t { f_t start_time = 0.0; @@ -69,6 +72,17 @@ class branch_and_bound_worker_t { std::vector start_lower; std::vector start_upper; + // The incumbent may change while we are still constructing RINS + // sub-MIP or doing guided diving. Save it here so we always use + // the same value throughout. + std::vector current_incumbent; + + // Variable locks (see definition 3.3 from T. Achterberg, “Constraint Integer Programming,” + // PhD, Technischen Universität Berlin, Berlin, 2007. doi: 10.14279/depositonce-1634). + // Here we assume that the constraints are in the form `Ax = b, l <= x <= u`. + std::vector var_up_locks; + std::vector var_down_locks; + pcgenerator_t rng; std::unique_ptr> orbital_fixing; @@ -78,8 +92,14 @@ class branch_and_bound_worker_t { bool recompute_basis = true; bool recompute_bounds = true; + // During the cut passes, this values can change. So we save a copy on root_heuristics + // and point these attributes to it. During normal exploration, this is set + // to the final values from the root node. const std::vector& root_solution; const std::vector& root_edge_norm; + const std::vector& var_types; + + pseudo_costs_t& pseudo_costs; void ensure_orbital_fixing() { @@ -97,6 +117,7 @@ class branch_and_bound_worker_t { const csr_matrix_t& Arow, const std::vector& var_type, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm, uint64_t rng_offset = 0) @@ -115,7 +136,9 @@ class branch_and_bound_worker_t { rng(settings.random_seed + pcgenerator_t::default_seed + rng_offset + worker_id, pcgenerator_t::default_stream ^ (worker_id + rng_offset)), root_solution(root_solution), - root_edge_norm(root_edge_norm) + root_edge_norm(root_edge_norm), + var_types(var_type), + pseudo_costs(pc) { } @@ -153,11 +176,19 @@ class bfs_worker_t : public branch_and_bound_worker_t { const csr_matrix_t& Arow, const std::vector& var_type, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm, uint64_t rng_offset = 0) - : Base( - worker_id, original_lp, Arow, var_type, settings, root_solution, root_edge_norm, rng_offset) + : Base(worker_id, + original_lp, + Arow, + var_type, + settings, + pc, + root_solution, + root_edge_norm, + rng_offset) { this->start_lower = original_lp.lower; this->start_upper = original_lp.upper; @@ -225,7 +256,30 @@ template class diving_worker_t : public branch_and_bound_worker_t { public: using Base = branch_and_bound_worker_t; - using Base::Base; + + diving_worker_t(i_t worker_id, + const simplex::lp_problem_t& original_lp, + const csr_matrix_t& Arow, + const std::vector& var_type, + const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, + const std::vector& root_solution, + const std::vector& root_edge_norm, + uint64_t rng_offset = 0) + : Base(worker_id, + original_lp, + Arow, + var_type, + settings, + pc, + root_solution, + root_edge_norm, + rng_offset) + { + this->start_lower = original_lp.lower; + this->start_upper = original_lp.upper; + this->search_strategy = search_strategy_t::COEFFICIENT_DIVING; + } // Apply bound strengthening to the starting variable bounds bool presolve_start_bounds(const simplex::simplex_solver_settings_t& settings) diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp index 6977c7882b..c4e54a61f1 100644 --- a/cpp/src/branch_and_bound/worker_pool.hpp +++ b/cpp/src/branch_and_bound/worker_pool.hpp @@ -24,6 +24,7 @@ class worker_pool_t { const std::vector& var_type, mip_symmetry_t* symmetry, const simplex::simplex_solver_settings_t& settings, + pseudo_costs_t& pc, const std::vector& root_solution, const std::vector& root_edge_norm, const uint64_t rng_offset = 0) @@ -36,7 +37,7 @@ class worker_pool_t { idle_workers_.clear_resize(num_workers); for (i_t i = 0; i < num_workers; ++i) { workers_[i] = std::make_unique( - i, original_lp, Arow, var_type, settings, root_solution, root_edge_norm, rng_offset); + i, original_lp, Arow, var_type, settings, pc, root_solution, root_edge_norm, rng_offset); idle_workers_.push_back(i); // Propagate the (possibly null) symmetry pointer; workers lazily build // their orbital_fixing/lexical_reduction state via ensure_orbital_fixing(). diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index b3940df890..f14b9c7b83 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -183,7 +183,6 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_HYPER_DIVING_COEFFICIENT, &mip_settings.diving_params.coefficient_diving, -1, 1, -1, "coefficient diving toggle: -1 automatic, 0 disabled, 1 enabled"}, {CUOPT_MIP_HYPER_DIVING_FARKAS, &mip_settings.diving_params.farkas_diving, -1, 1, -1, "Farkas diving toggle: -1 automatic, 0 disabled, 1 enabled"}, {CUOPT_MIP_HYPER_DIVING_VECTOR_LENGTH, &mip_settings.diving_params.vector_length_diving, -1, 1, -1, "vector-length diving toggle: -1 automatic, 0 disabled, 1 enabled"}, - {CUOPT_MIP_HYPER_DIVING_MIN_NODE_DEPTH, &mip_settings.diving_params.min_node_depth, 0, std::numeric_limits::max(), 10, "minimum depth at which to start diving"}, {CUOPT_MIP_HYPER_DIVING_NODE_LIMIT, &mip_settings.diving_params.node_limit, 0, std::numeric_limits::max(), 500, "maximum nodes explored per dive"}, {CUOPT_MIP_HYPER_DIVING_BACKTRACK_LIMIT, &mip_settings.diving_params.backtrack_limit, 0, std::numeric_limits::max(), 5, "maximum backtracking allowed per dive"}, // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") diff --git a/cpp/src/mip_heuristics/root_heuristics.hpp b/cpp/src/mip_heuristics/root_heuristics.hpp index eac8cb954f..b9645579b5 100644 --- a/cpp/src/mip_heuristics/root_heuristics.hpp +++ b/cpp/src/mip_heuristics/root_heuristics.hpp @@ -19,21 +19,30 @@ struct cut_pass_heuristics_t { csr_matrix_t Arow_; std::vector root_solution_; std::vector root_edge_norm_; + pseudo_costs_t pseudo_costs_; + omp_atomic_t active_workers_; std::atomic halt_; std::unique_ptr> submip_worker_; + std::vector>> diving_workers_; fj_cpu_worker_t fj_cpu_worker_; cut_pass_heuristics_t(const csr_matrix_t& Arow, const std::vector& var_types, const std::vector& root_solution, - const std::vector& root_edge_norm) + const std::vector& root_edge_norm, + const simplex::simplex_solver_settings_t& settings) : var_types_(var_types), Arow_(Arow), root_solution_(root_solution), root_edge_norm_(root_edge_norm), + pseudo_costs_(root_solution.size(), settings), + active_workers_(0), halt_(false), - submip_worker_(nullptr) {}; + submip_worker_(nullptr) + { + pseudo_costs_.Arow = Arow; + }; ~cut_pass_heuristics_t() { stop_and_sync(); } @@ -53,6 +62,14 @@ struct cut_pass_heuristics_t { #pragma omp taskwait depend(in : *worker) submip_worker_.reset(); } + + for (auto& worker : diving_workers_) { + diving_worker_t* w = worker.get(); +#pragma omp taskwait depend(in : *w) + worker.reset(); + } + + diving_workers_.clear(); } diving_worker_t* create_submip_worker( @@ -65,7 +82,7 @@ struct cut_pass_heuristics_t { search_strategy_t type) { submip_worker_ = std::make_unique>( - id, lp, Arow_, var_types_, settings, root_solution_, root_edge_norm_); + id, lp, Arow_, var_types_, settings, pseudo_costs_, root_solution_, root_edge_norm_); submip_worker_->start_node = mip_node_t(root_obj, root_vstatus); submip_worker_->leaf_vstatus = root_vstatus; submip_worker_->leaf_solution.x = sol; @@ -76,6 +93,43 @@ struct cut_pass_heuristics_t { return submip_worker_.get(); } + + void initialize_pseudocost(const simplex::lp_problem_t& lp, + const std::vector& vstatus, + const std::vector& fractional, + const simplex::lp_solution_t& lp_solution, + const std::vector& basic_list, + const std::vector& nonbasic_list, + simplex::basis_update_mpf_t& basis_factors) + { + pseudo_costs_.initialize_with_estimate( + lp, vstatus, fractional, lp_solution, basic_list, nonbasic_list, basis_factors); + } + + diving_worker_t* create_diving_worker( + i_t cut_pass, + const simplex::lp_problem_t& lp, + const simplex::simplex_solver_settings_t& settings, + const mip_node_t& root_node, + search_strategy_t strategy) + { + std::unique_ptr>& worker = diving_workers_.emplace_back( + std::make_unique>(diving_workers_.size(), + lp, + Arow_, + var_types_, + settings, + pseudo_costs_, + root_solution_, + root_edge_norm_)); + worker->start_node = root_node.detach_copy(); + worker->start_lower = lp.lower; + worker->start_upper = lp.upper; + worker->search_strategy = strategy; + worker->set_active(); + + return worker.get(); + } }; /// \brief Object Representing the heuristics run on the root node. @@ -94,8 +148,14 @@ struct root_heuristics_t { std::shared_ptr> worker_count_; i_t max_workers_; + // Keep track of the last diving heuristic used (so we can cycle between them in low thread + // count systems) + i_t next_diving_type_; + root_heuristics_t(i_t max_workers) - : worker_count_(std::make_shared>(0)), max_workers_(max_workers) + : worker_count_(std::make_shared>(0)), + max_workers_(max_workers), + next_diving_type_(0) { } @@ -114,26 +174,37 @@ struct root_heuristics_t { cut_passes_heuristics_.clear(); } + void stop_old_workers(i_t cut_pass, i_t new_workers) + { + if (new_workers <= 0) return; + + // On the first pass we use a thread to generate the clique table + i_t cut_generation = cut_pass == 0 ? 2 : 1; + i_t total_workers = new_workers + worker_count_->load() + cut_generation; + if (total_workers <= max_workers_) { return; } + + for (auto& heuristic : cut_passes_heuristics_) { + // Skip the current heuristic entry + if (&heuristic == &cut_passes_heuristics_.back()) { break; } + + i_t active = heuristic->active_workers_; + if (active > 0 && !heuristic->halt_.load(std::memory_order_acquire)) { + heuristic->send_stop_signal(); + new_workers -= active; + if (new_workers <= 0) return; + } + } + } + std::shared_ptr> create_new_cut_pass_heuristic( - i_t cut_pass, const csr_matrix_t& Arow, const std::vector& var_types, const std::vector& root_solution, - const std::vector& root_edge_norm) + const std::vector& root_edge_norm, + const simplex::simplex_solver_settings_t& settings) { - // If we already exhausted all threads for the root heuristics, stop workers for the - // oldest set of heuristics launched. Leave 2 threads for the cut passes and the clique - // table generation. Add the number of workers that will be launched (1 submip worker + - // 1 CPU FJ worker). - i_t clique_table_generation = cut_pass == 0 ? 1 : 0; - if (*worker_count_ + 3 + clique_table_generation > max_workers_ && - !cut_passes_heuristics_.empty()) { - cut_passes_heuristics_.begin()->get()->send_stop_signal(); - cut_passes_heuristics_.erase(cut_passes_heuristics_.begin()); - } - return cut_passes_heuristics_.emplace_back(std::make_shared>( - Arow, var_types, root_solution, root_edge_norm)); + Arow, var_types, root_solution, root_edge_norm, settings)); } };