diff --git a/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp b/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp index 505a3b61ba..4464a156a7 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp @@ -41,6 +41,11 @@ struct mip_submip_hyper_params_t { // number of simplex iteration from the parent B&B. f_t iteration_limit_ratio = 0.8; + // If there is not enough variables fixed or we already found an improving solution, + // perform a short DFS to quickly find a feasible solution. This setting controls + // the maximum number of nodes allow for backtracking. + i_t dfs_max_backtrack = 5; + // Run CPU FJ over the sub-MIP bool enable_cpufj = true; }; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index fd7987ae1c..1454884a0e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -41,9 +41,17 @@ #include #include #include +#include #include #include +#define SUBMIP_VERBOSE false +#if SUBMIP_VERBOSE +#define DEBUG_SUBMIP(fmt, ...) settings_.log.print_format(fmt, __VA_ARGS__); +#else +#define DEBUG_SUBMIP(fmt, ...) +#endif + namespace cuopt::mathematical_optimization::mip { using simplex::basis_update_mpf_t; @@ -205,6 +213,7 @@ inline char feasible_solution_symbol(heuristics_origin_t origin) inline char feasible_solution_symbol(search_strategy_t strategy, bool show_diving) { if (strategy == search_strategy_t::BEST_FIRST) return 'B'; + if (strategy == search_strategy_t::RINS) return 'S'; if (!show_diving) return 'D'; switch (strategy) { @@ -215,8 +224,10 @@ inline char feasible_solution_symbol(search_strategy_t strategy, bool show_divin case search_strategy_t::GUIDED_DIVING: return 'G'; case search_strategy_t::FARKAS_DIVING: return 'F'; case search_strategy_t::VECTOR_LENGTH_DIVING: return 'V'; - default: return 'U'; + case search_strategy_t::RINS: return 'S'; } + + return 'U'; } template @@ -654,16 +665,19 @@ void branch_and_bound_t::set_solution_from_cpu_fj(f_t obj, // user space. template void branch_and_bound_t::set_solution_from_submip( + const lp_problem_t& lp, const std::vector& solution, const third_party_presolve_t& presolver, - f_t fixrate, - f_t obj) + f_t fixrate) { + bool check_postsolve = false; std::vector leaf_sol; - presolver.uncrush_primal_solution(solution, leaf_sol); + presolver.uncrush_primal_solution(solution, leaf_sol, check_postsolve); + f_t obj = compute_objective(lp, leaf_sol); + std::vector user_sol; mutex_original_lp_.lock(); - uncrush_primal_solution(original_problem_, original_lp_, leaf_sol, user_sol); + uncrush_primal_solution(original_problem_, lp, leaf_sol, user_sol); mutex_original_lp_.unlock(); settings_.log.debug_format("SubMIP found a feasible solution with obj={:.4g}", obj); bool success = set_solution_from_heuristics(user_sol, heuristics_origin_t::SUBMIP); @@ -991,7 +1005,7 @@ branch_variable_t branch_and_bound_t::variable_selection( branch_var = pc_.variable_selection(fractional, solution); } - round_dir = martin_criteria(solution[branch_var], root_relax_soln_.x[branch_var]); + round_dir = martin_criteria(solution[branch_var], worker->root_solution[branch_var]); return {branch_var, round_dir}; @@ -1000,10 +1014,10 @@ branch_variable_t branch_and_bound_t::variable_selection( original_lp_, fractional, solution, var_up_locks_, var_down_locks_, log); case search_strategy_t::LINE_SEARCH_DIVING: - return line_search_diving(fractional, solution, root_relax_soln_.x, log); + return line_search_diving(fractional, solution, worker->root_solution, log); case search_strategy_t::PSEUDOCOST_DIVING: - return pseudocost_diving(pc_, fractional, solution, root_relax_soln_.x, log); + return pseudocost_diving(pc_, fractional, solution, worker->root_solution, log); case search_strategy_t::GUIDED_DIVING: assert(incumbent_.has_incumbent); @@ -1018,9 +1032,9 @@ branch_variable_t branch_and_bound_t::variable_selection( case search_strategy_t::VECTOR_LENGTH_DIVING: return vector_length_diving(worker->leaf_problem, fractional, solution, log); - case search_strategy_t::SUBMIP: // This is used for solving the DFS of the sub-MIP. + case search_strategy_t::RINS: // This is used for solving the DFS of the sub-MIP. branch_var = pc_.variable_selection(fractional, solution); - round_dir = martin_criteria(solution[branch_var], root_relax_soln_.x[branch_var]); + round_dir = martin_criteria(solution[branch_var], worker->root_solution[branch_var]); return {branch_var, round_dir}; } @@ -1173,7 +1187,7 @@ struct deterministic_bfs_policy_t const std::vector& x) override { i_t var = this->worker.pc_snapshot.variable_selection(fractional, x); - auto dir = martin_criteria(x[var], this->bnb.root_relax_soln_.x[var]); + auto dir = martin_criteria(x[var], this->worker.root_solution[var]); return {var, dir}; } @@ -1255,15 +1269,15 @@ struct deterministic_diving_policy_t switch (this->worker.diving_type) { case search_strategy_t::PSEUDOCOST_DIVING: return pseudocost_diving( - this->worker.pc_snapshot, fractional, x, *this->worker.root_solution, log); + this->worker.pc_snapshot, fractional, x, this->worker.root_solution, log); case search_strategy_t::LINE_SEARCH_DIVING: - return line_search_diving(fractional, x, *this->worker.root_solution, log); + return line_search_diving(fractional, x, this->worker.root_solution, log); case search_strategy_t::GUIDED_DIVING: if (this->worker.incumbent_snapshot.empty()) { return pseudocost_diving( - this->worker.pc_snapshot, fractional, x, *this->worker.root_solution, log); + this->worker.pc_snapshot, fractional, x, this->worker.root_solution, log); } else { return guided_diving( this->worker.pc_snapshot, fractional, x, this->worker.incumbent_snapshot, log); @@ -1595,7 +1609,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( bool feasible = worker->set_lp_variable_bounds(node_ptr, settings_); dual_status_t lp_status = dual_status_t::DUAL_UNBOUNDED; - worker->leaf_edge_norms = edge_norms_; + worker->leaf_edge_norms = worker->root_edge_norm; if (worker->recompute_bounds && worker->orbital_fixing && worker->search_strategy == search_strategy_t::BEST_FIRST) { worker->orbital_fixing->reset(symmetry_, node_ptr); @@ -2113,7 +2127,7 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker, // 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::SUBMIP) { + if (worker->search_strategy != search_strategy_t::RINS) { diving_worker_pool_.return_worker_to_pool(worker); } } @@ -2180,16 +2194,28 @@ bool branch_and_bound_t::launch_rins_worker(const std::vector& so diving_worker_t* worker = rins_worker_pool_.pop_idle_worker(); if (!worker) return false; + std::vector current_incumbent; + mutex_upper_.lock(); + 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) + worker->start_node = mip_node_t(root_objective_, root_vstatus_); + worker->leaf_vstatus = root_vstatus_; + worker->leaf_problem.lower = original_lp_.lower; + worker->leaf_problem.upper = original_lp_.upper; + worker->leaf_solution.x = sol; + worker->search_strategy = search_strategy_t::RINS; worker->set_active(); - worker->search_strategy = search_strategy_t::SUBMIP; if (settings_.inside_submip) { // LLVM libomp's GOMP compatibility path skips GCC's firstprivate copy // function for included tasks. - rins(worker, sol); + rins(worker, current_incumbent, var_types_); } else { -#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) firstprivate(worker, sol) - rins(worker, sol); +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ + firstprivate(worker, current_incumbent) + rins(worker, current_incumbent, var_types_); } return true; @@ -2198,10 +2224,12 @@ bool branch_and_bound_t::launch_rins_worker(const std::vector& so template void branch_and_bound_t::solve_submip(diving_worker_t* worker, const std::vector& current_incumbent, + const std::vector& var_types, i_t num_var_fixed, i_t num_integers, i_t submip_level, - std::string_view log_prefix) + std::string_view log_prefix, + bool is_root_heuristic) { double start_time = tic(); @@ -2216,11 +2244,12 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke if (!feasible) { // This should never happen since we are fixing bounds that are already in the incumbent. rins_stats_.save_infeasible(fixrate); + DEBUG_SUBMIP("{} The problem is infeasible after running bound strengthening!", log_prefix); return; } - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); - f_t user_obj = compute_user_objective(original_lp_, upper_bound_.load()); + f_t user_lower = compute_user_objective(worker->leaf_problem, get_lower_bound()); + f_t user_obj = compute_user_objective(worker->leaf_problem, upper_bound_.load()); f_t rel_gap = user_relative_gap(user_obj, user_lower); i_t explored = exploration_stats_.nodes_explored; @@ -2233,10 +2262,10 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.submip_settings.level = submip_level; - submip_settings.log.log = false; submip_settings.benchmark_info_ptr = nullptr; + submip_settings.log.log = SUBMIP_VERBOSE; -#ifdef DEBUG_SUBMIP +#ifdef SAVE_SUBMIP_TO_FILE submip_settings.log.log_prefix = std::format("{}{}", settings_.log.log_prefix, worker->worker_id); CUOPT_LOG_INFO("Writting submip %s to MPS file", submip_settings.log.log_prefix); worker->leaf_problem.write_mps(std::format("submip-{}.mps", submip_settings.log.log_prefix), @@ -2254,12 +2283,19 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke submip_settings.relative_mip_gap_tol = std::min(settings_.submip_settings.target_mip_gap, rel_gap); - submip_settings.submip_settings.rins = - settings_.submip_settings.rins != 0 && submip_level <= settings_.submip_settings.max_level; + bool max_recursion = submip_level > settings_.submip_settings.max_level; + submip_settings.submip_settings.rins = settings_.submip_settings.rins != 0 && !max_recursion; + + DEBUG_SUBMIP("{}Sub-MIP: num variables fixed={}/{} ({:.2f}%)", + log_prefix, + num_var_fixed, + num_integers, + fixrate * 100); - submip_settings.log.debug_format( - "Sub-MIP solve settings: time_limit={:.2f}, node_limit={}, iter_limit={} (current_iter={}), " + DEBUG_SUBMIP( + "{}Sub-MIP solve settings: time_limit={:.2f}, node_limit={}, iter_limit={} (current_iter={}), " "tol={:g}", + log_prefix, submip_settings.time_limit, submip_settings.node_limit, submip_settings.branch_and_bound_simplex_iteration_limit, @@ -2270,7 +2306,7 @@ 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, 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); @@ -2288,13 +2324,12 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke // Also handle optimal if (submip_problem.num_rows == 0 || submip_problem.num_cols == 0) { - submip_settings.log.debug_format( - "Sub-MIP presolved to a trivial {} x {} problem; solving by bound pushing", - submip_problem.num_rows, - submip_problem.num_cols); + DEBUG_SUBMIP("{}Sub-MIP presolved to a trivial {} x {} problem; solving by bound pushing", + log_prefix, + submip_problem.num_rows, + submip_problem.num_cols); std::vector reduced_sol(submip_problem.num_cols); - f_t obj = 0.0; for (i_t j = 0; j < submip_problem.num_cols; ++j) { const f_t c = submip_problem.objective[j]; @@ -2306,26 +2341,30 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke } else { reduced_sol[j] = std::isfinite(l) ? l : (std::isfinite(u) ? u : 0); } - - obj += reduced_sol[j] * c; } - set_solution_from_submip(reduced_sol, presolver, fixrate, obj); + set_solution_from_submip(worker->leaf_problem, reduced_sol, presolver, fixrate); + return; + } + + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; return; } submip_settings.heuristic_preemption_callback = nullptr; submip_settings.dual_simplex_objective_callback = nullptr; submip_settings.set_simplex_solution_callback = nullptr; - submip_settings.solution_callback = [this, &presolver, fixrate](const std::vector& solution, - f_t obj) { - this->set_solution_from_submip(solution, presolver, fixrate, obj); + submip_settings.solution_callback = [this, &presolver, fixrate, worker]( + const std::vector& solution, f_t obj) { + this->set_solution_from_submip(worker->leaf_problem, solution, presolver, fixrate); }; - submip_settings.log.debug_format("Sub-MIP: {} constraints, {} variables, {} nonzeros\n", - submip_problem.num_rows, - submip_problem.num_cols, - submip_problem.A.nnz()); + DEBUG_SUBMIP("{}Sub-MIP: {} constraints, {} variables, {} nonzeros\n", + log_prefix, + submip_problem.num_rows, + submip_problem.num_cols, + submip_problem.A.nnz()); probing_implied_bound_t empty_probing(submip_problem.num_cols); branch_and_bound_t submip_bnb(submip_problem, submip_settings, tic(), empty_probing); @@ -2338,24 +2377,24 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke presolver.crush_primal_solution(submip_problem, current_incumbent, presolved_incumbent); submip_bnb.set_initial_guess(presolved_incumbent); - const f_t user_upper = compute_user_objective(original_lp_, upper_bound_.load()); - const f_t submip_cutoff = - user_upper / submip_bnb.original_lp_.obj_scale - submip_bnb.original_lp_.obj_constant; + const f_t user_upper = compute_user_objective(worker->leaf_problem, upper_bound_.load()); + const f_t submip_cutoff = compute_presolved_objective(submip_bnb.original_lp_, user_upper); submip_bnb.set_initial_upper_bound(submip_cutoff); - submip_bnb.set_initial_pseudocost(pc_, presolver.get_reduced_to_original_map()); + if (!is_root_heuristic) + submip_bnb.set_initial_pseudocost(pc_, presolver.get_reduced_to_original_map()); if (submip_halt_callback_) { // Copy the halt callback to the deeper level. submip_bnb.set_submip_halt_callback(submip_halt_callback_); } else { // This should only be called by the main solver. - submip_bnb.set_submip_halt_callback([this](f_t, f_t submip_lower_bound) { + submip_bnb.set_submip_halt_callback([this, worker](f_t, f_t submip_lower_bound) { f_t user_upper = compute_user_objective(this->original_lp_, this->upper_bound_.load()); bool is_cutoff = original_lp_.obj_scale > 0 ? submip_lower_bound > user_upper : user_upper > submip_lower_bound; bool is_solver_running = this->solver_status_ == mip_status_t::UNSET && this->is_running_; - return is_cutoff || !is_solver_running; + return is_cutoff || !is_solver_running || worker->halt; }); } @@ -2393,8 +2432,9 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke mip_status_t submip_status = submip_bnb.solve(submip_solution); f_t submip_time = toc(start_time); - submip_settings.log.debug_format( - "Sub-MIP: status={}, iterations={} (total={}), presolve_time={:.2f}, total_time={:.2f} \n", + DEBUG_SUBMIP( + "{}Sub-MIP: status={}, iterations={} (total={}), presolve_time={:.2f}, total_time={:.2f} \n", + log_prefix, mip_status_to_string(submip_status), submip_solution.simplex_iterations, exploration_stats_.total_simplex_iters.load(), @@ -2408,7 +2448,7 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke } if (submip_solution.has_incumbent) { - set_solution_from_submip(submip_solution.x, presolver, fixrate, submip_solution.objective); + set_solution_from_submip(worker->leaf_problem, submip_solution.x, presolver, fixrate); } // Accumulate simplex iterations to determine when to stop exploring the sub-MIP @@ -2418,9 +2458,9 @@ void branch_and_bound_t::solve_submip(diving_worker_t* worke } template -inline f_t submip_get_max_fixrate(const submip_stats_t& stats, - const mip_submip_hyper_params_t& submip_settings, - pcgenerator_t& rng) +f_t submip_get_max_fixrate(const submip_stats_t& stats, + const mip_submip_hyper_params_t& submip_settings, + pcgenerator_t& rng) { // Adaptive fix rate based on previous successes and failures. f_t low = submip_settings.base_target_fixrate; @@ -2543,11 +2583,13 @@ void extend_variable_fixings(const simplex_solver_settings_t& settings } template -void branch_and_bound_t::rins(diving_worker_t* rins_worker, - const std::vector& node_solution) +void branch_and_bound_t::rins(diving_worker_t* worker, + const std::vector& current_incumbent, + const std::vector& var_types, + bool is_root_heuristic) { raft::common::nvtx::range scope("BB::rins_thread"); - if (rins_worker->orbital_fixing) { rins_worker->orbital_fixing->disable(); } + if (worker->orbital_fixing) { worker->orbital_fixing->disable(); } i_t submip_level = settings_.submip_settings.level + 1; std::string log_prefix = std::format("[RINS {}] ", submip_level); @@ -2555,44 +2597,31 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, ++rins_stats_.total_calls; bool has_submip = false; - const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + worker->recompute_bounds = false; + worker->recompute_basis = true; branch_and_bound_stats_t rins_stats; + mip_node_t& node = worker->start_node; + std::vector& lower = worker->leaf_problem.lower; + std::vector& upper = worker->leaf_problem.upper; + std::vector& bounds_changed = worker->bounds_changed; + std::vector& current_sol = worker->leaf_solution.x; - // Note that this node does not have the vstatus (it was clear at the start of B&B exploration) - mip_node_t node = search_tree_.root.detach_copy(); - rins_worker->leaf_vstatus = root_vstatus_; - rins_worker->leaf_problem.lower = original_lp_.lower; - rins_worker->leaf_problem.upper = original_lp_.upper; - rins_worker->leaf_solution.x = node_solution; - rins_worker->recompute_bounds = false; - rins_worker->recompute_basis = true; - - std::vector& lower = rins_worker->leaf_problem.lower; - std::vector& upper = rins_worker->leaf_problem.upper; - std::vector& bounds_changed = rins_worker->bounds_changed; - std::vector& current_sol = rins_worker->leaf_solution.x; std::vector fractional; - i_t num_frac = fractional_variables(settings_, current_sol, var_types_, fractional); - - std::vector current_incumbent; - mutex_upper_.lock(); - current_incumbent = incumbent_.x; - mutex_upper_.unlock(); + i_t num_frac = fractional_variables(settings_, current_sol, 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, var_types, settings_.fixed_tol, integer_list); i_t num_integers = integer_list.size(); - f_t max_fixrate = - submip_get_max_fixrate(rins_stats_, settings_.submip_settings, rins_worker->rng); + f_t max_fixrate = submip_get_max_fixrate(rins_stats_, settings_.submip_settings, worker->rng); f_t min_fixrate = std::min(settings_.submip_settings.min_fixrate, max_fixrate); i_t max_var_fixed = max_fixrate * num_integers; i_t min_var_fixed = min_fixrate * num_integers; i_t num_var_fixed = 0; - while (solver_status_ == mip_status_t::UNSET && is_running_) { + while (solver_status_ == mip_status_t::UNSET && is_running_ && !worker->halt) { // RINS neighbourhood 1: Fix all the integer variables where the starting solution matches the // current incumbent, considering only the fractional values in the current node i_t prev_num_fixed = num_var_fixed; @@ -2608,11 +2637,11 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, // Enough variables has been fixed if (num_var_fixed >= min_var_fixed) { - settings_.log.debug_format("{}Fixed {} variables (max={}, min={})\n", - log_prefix, - num_var_fixed, - max_var_fixed, - min_var_fixed); + DEBUG_SUBMIP("{}Fixed {} variables (max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } @@ -2637,11 +2666,11 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, // Enough variables were fixed if (num_var_fixed >= min_var_fixed) { - settings_.log.debug_format("{}Fixed {} variables (max={}, min={})\n", - log_prefix, - num_var_fixed, - max_var_fixed, - min_var_fixed); + DEBUG_SUBMIP("{}Fixed {} variables (max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } @@ -2651,10 +2680,10 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, // an integer solution first in order to reach the fixing threshold. if (prev_num_fixed == num_var_fixed) { extend_variable_fixings(settings_, - rins_worker->leaf_problem.objective, + worker->leaf_problem.objective, fractional, current_sol, - root_relax_soln_.x, + worker->root_solution, max_var_fixed, lower, upper, @@ -2662,21 +2691,21 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, num_var_fixed); if (num_var_fixed >= min_var_fixed) { - settings_.log.debug_format("{}Fixed {} variables (max={}, min={})\n", - log_prefix, - num_var_fixed, - max_var_fixed, - min_var_fixed); + DEBUG_SUBMIP("{}Fixed {} variables (max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } if (prev_num_fixed == num_var_fixed) { - settings_.log.debug_format("{}Could not fix more variables ({}, max={}, min={})\n", - log_prefix, - num_var_fixed, - max_var_fixed, - min_var_fixed); + DEBUG_SUBMIP("{}Could not fix more variables ({}, max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } @@ -2693,14 +2722,14 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, // We continue to do this until enough variables were fixed or no variable is left to fix. logger_t log; log.log = false; - dual_status_t lp_status = solve_node_lp(&node, rins_worker, rins_stats, log); + dual_status_t lp_status = solve_node_lp(&node, worker, rins_stats, log); if (lp_status != dual_status_t::OPTIMAL) { break; } fractional.clear(); - num_frac = fractional_variables(settings_, current_sol, var_types_, fractional); + num_frac = fractional_variables(settings_, current_sol, var_types, fractional); - f_t leaf_obj = compute_objective(rins_worker->leaf_problem, current_sol); + f_t leaf_obj = compute_objective(worker->leaf_problem, current_sol); node.lower_bound = leaf_obj; snap_to_lattice(&node, leaf_obj); @@ -2708,11 +2737,11 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, if (num_frac == 0) { // We found a feasible solution when fixing the variables in RINS. - add_feasible_solution(leaf_obj, current_sol, -1, search_strategy_t::SUBMIP); + add_feasible_solution(leaf_obj, current_sol, -1, search_strategy_t::RINS); break; } - rins_worker->recompute_basis = false; + worker->recompute_basis = false; } f_t fixrate = (f_t)num_var_fixed / num_integers; @@ -2723,14 +2752,11 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, // levels up to try to find a feasible solution quickly from the neighbourhood. if (fixrate < settings_.submip_settings.min_fixrate_cap || (settings_.inside_submip && rins_stats_.total_success != 0)) { - // We need to re-populate the vstatus of the node since it was previously cleared. - rins_worker->start_node = std::move(node); - rins_worker->start_node.packed_vstatus = simplex::compress_vstatus(rins_worker->leaf_vstatus); + worker->start_node.packed_vstatus = simplex::compress_vstatus(worker->leaf_vstatus); + worker->start_lower = lower; + worker->start_upper = upper; - rins_worker->start_lower = lower; - rins_worker->start_upper = upper; - - bool is_feasible = rins_worker->presolve_start_bounds(settings_); + bool is_feasible = worker->presolve_start_bounds(settings_); if (is_feasible) { fj_cpu_worker_t submip_fj_cpu_worker; @@ -2743,21 +2769,31 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, f_t time_limit = std::max(settings_.time_limit - toc(exploration_stats_.start_time), 0); f_t work_limit = 1.0; - submip_fj_cpu_worker.create_worker(rins_worker->leaf_problem, - var_types_, - rins_worker->leaf_solution.x, + submip_fj_cpu_worker.create_worker(worker->leaf_problem, + var_types, + worker->leaf_solution.x, settings_, std::format("{} [CPU FJ]", log_prefix), - rins_worker->rng.next_i64()); + worker->rng.next_i64()); submip_fj_cpu_worker.run_sync(time_limit, work_limit); } - dive_with(rins_worker, 5); + // We need the pseudocost to do the DFS, which we do not have during the cut passes. + if (!is_root_heuristic) { + DEBUG_SUBMIP("{} Running a quick DFS for the submip!", log_prefix); + dive_with(worker, settings_.submip_settings.dfs_max_backtrack); + } } } else { - solve_submip( - rins_worker, current_incumbent, num_var_fixed, num_integers, submip_level, log_prefix); + solve_submip(worker, + current_incumbent, + var_types, + num_var_fixed, + num_integers, + submip_level, + log_prefix, + is_root_heuristic); } } @@ -2766,7 +2802,7 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, exploration_stats_.total_simplex_iters += rins_stats.total_simplex_iters; } - settings_.log.debug_format( + DEBUG_SUBMIP( "{}success={}, infeasible={}, calls={}, fixrate={:.4g} ({}), max_fixrate={:.4g} ({}), " "min_fixrate={:.4g} ({})\n", log_prefix, @@ -2780,7 +2816,77 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, min_fixrate, min_var_fixed); - rins_worker_pool_.return_worker_to_pool(rins_worker); + // If the pool is uninitialized (i.e., in the root node), then this just inactivate the worker. + if (!is_root_heuristic) { + rins_worker_pool_.return_worker_to_pool(worker); + } else { + worker->set_inactive(); + } +} + +template +void branch_and_bound_t::launch_root_heuristics( + const lp_problem_t& lp, + const std::vector& sol, + i_t cut_pass, + root_heuristics_t& root_heuristics) +{ + if (settings_.deterministic) return; + if (settings_.num_threads < 2) return; + + // 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; + + if (is_cpufj_enabled) { + f_t work_limit = std::numeric_limits::infinity(); + f_t time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + + current_heuristic->fj_cpu_worker_.improvement_callback = + [this](f_t obj, const std::vector& assignment, double work_units) { + set_solution_from_cpu_fj(obj, assignment, work_units); + }; + current_heuristic->fj_cpu_worker_.create_worker( + lp, var_types_, sol, settings_, "[RootCut CPUFJ] "); + ++(*worker_count); + +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) \ + affinity(current_heuristic -> fj_cpu_worker_) firstprivate(current_heuristic, worker_count) \ + depend(out : current_heuristic->fj_cpu_worker_.fj_cpu) + { + current_heuristic->fj_cpu_worker_.run_sync(time_limit, work_limit); + --(*worker_count); + } + } + + if (settings_.submip_settings.rins != 0 && incumbent_.has_incumbent) { + diving_worker_t* worker = current_heuristic->create_submip_worker( + cut_pass, lp, settings_, root_objective_, root_vstatus_, sol); + + std::vector current_incumbent; + mutex_upper_.lock(); + current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + + if (settings_.inside_submip) { + // LLVM libomp's GOMP compatibility path skips GCC's firstprivate copy + // function for included tasks. + rins(worker, current_incumbent, current_heuristic->var_types_, is_root_heuristic); + + } else { + ++(*worker_count); +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ + firstprivate(current_incumbent, current_heuristic, worker_count) depend(out : *worker) + { + rins(worker, current_incumbent, current_heuristic->var_types_, is_root_heuristic); + --(*worker_count); + } + } + } } template @@ -3329,13 +3435,14 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nonbasic_list, edge_norms_); } + settings_.log.printf("\n"); solving_root_relaxation_ = false; f_t root_relax_elapsed_time = toc(root_relax_start_time); exploration_stats_.total_lp_solve_time = root_relax_elapsed_time; if (root_status == lp_status_t::INFEASIBLE) { - settings_.log.printf("\nThe root LP relaxation is infeasible\n", + settings_.log.printf("The root LP relaxation is infeasible\n", lp_status_to_string(root_status).c_str()); signal_extend_cliques_.store(true, std::memory_order_release); #pragma omp taskwait depend(in : *clique_signal) @@ -3343,7 +3450,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::UNBOUNDED) { - settings_.log.printf("\nThe root relaxation is unbounded\n", + settings_.log.printf("The root relaxation is unbounded\n", lp_status_to_string(root_status).c_str()); if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); @@ -3354,7 +3461,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::TIME_LIMIT) { - settings_.log.printf("\n"); solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, -inf); signal_extend_cliques_.store(true, std::memory_order_release); @@ -3363,7 +3469,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::WORK_LIMIT) { - settings_.log.printf("\n"); solver_status_ = mip_status_t::WORK_LIMIT; set_final_solution(solution, -inf); signal_extend_cliques_.store(true, std::memory_order_release); @@ -3372,7 +3477,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::NUMERICAL_ISSUES) { - settings_.log.printf("\n"); solver_status_ = mip_status_t::NUMERICAL; set_final_solution(solution, -inf); signal_extend_cliques_.store(true, std::memory_order_release); @@ -3381,7 +3485,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } assert(root_status == lp_status_t::OPTIMAL); - settings_.log.printf("\n"); settings_.log.print_format("Root relaxation solution found in {} iterations and {:.2f}s by {}\n", root_relax_soln_.iterations, root_relax_elapsed_time, @@ -3460,12 +3563,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut compute_user_objective(original_lp_, root_relax_objective); } - constexpr bool enable_root_cut_cpufj = true; - fj_cpu_worker_t root_fj_cpu_worker; - root_fj_cpu_worker.improvement_callback = - [this](f_t obj, const std::vector& assignment, double work_units) { - set_solution_from_cpu_fj(obj, assignment, work_units); - }; + root_heuristics_t root_heuristics(settings_.num_threads - 1); f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; @@ -3497,9 +3595,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::OPTIMAL; } - cut_pass_result_t cut_pass_result; - root_fj_cpu_worker.run_async(settings_.time_limit - toc(exploration_stats_.start_time)); + launch_root_heuristics(original_lp_, root_relax_soln_.x, cut_pass, root_heuristics); + cut_pass_result_t cut_pass_result; cut_pass_result = do_cut_pass(cut_pass, solution, num_fractional, @@ -3518,7 +3616,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_objective, cut_pool_size, saved_solution); - root_fj_cpu_worker.stop(); + + mutex_upper_.lock(); + if (incumbent_.has_incumbent && incumbent_.x.size() != original_lp_.num_cols) { + std::vector uncrushed_incumbent; + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, uncrushed_incumbent); + crush_primal_solution( + original_problem_, original_lp_, uncrushed_incumbent, new_slacks_, incumbent_.x); + } + mutex_upper_.unlock(); if (cut_pass_result.action == cut_pass_action_t::RETURN) { if (settings_.benchmark_info_ptr != nullptr) { @@ -3529,16 +3635,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return cut_pass_result.status; } if (cut_pass_result.action == cut_pass_action_t::BREAK) { break; } - - if (enable_root_cut_cpufj && !settings_.deterministic && settings_.num_threads >= 2 && - cut_pass + 1 < settings_.max_cut_passes) { - f_t root_cut_cpufj_build_start_time = tic(); - root_fj_cpu_worker.create_worker( - original_lp_, var_types_, root_relax_soln_.x, settings_, "[RootCut CPUFJ] "); - settings_.log.debug("Root cut CPUFJ problem build time after pass %d: %.6f seconds\n", - cut_pass, - toc(root_cut_cpufj_build_start_time)); - } } // Publish the post-cuts root LP value. @@ -3554,19 +3650,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.benchmark_info_ptr->cut_generation_time_sec = cut_generation_time; } if (cut_info.has_cuts()) { - // If the incumbent is set before or during the cut passes, it may not have the correct - // dimensions as cuts add additional constraints/variables to `original_lp_`. - mutex_upper_.lock(); - if (incumbent_.has_incumbent && incumbent_.x.size() != original_lp_.num_cols) { - std::vector uncrushed_incumbent; - uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, uncrushed_incumbent); - crush_primal_solution( - original_problem_, original_lp_, uncrushed_incumbent, new_slacks_, incumbent_.x); - } - - mutex_upper_.unlock(); - - settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time); + settings_.log.printf("Root cut passes time: %.2f seconds\n", cut_generation_time); settings_.log.printf("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", original_lp_.num_rows, @@ -3576,30 +3660,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.printf("\n"); } - if (enable_root_cut_cpufj && cut_info.has_cuts()) { - f_t root_cut_cpufj_build_start_time = tic(); - // In deterministic mode this CPUFJ is built on the B&B task while the LS deterministic - // CPUFJ is being built on the main thread; both would otherwise race on the global - // seed_generator and pick non-reproducible seeds. Pin a stable seed here so this - // climber's behavior depends only on settings_.random_seed. - int64_t root_cut_cpufj_seed = - settings_.deterministic ? static_cast(settings_.random_seed) : -1; - root_fj_cpu_worker.create_worker(original_lp_, - var_types_, - root_relax_soln_.x, - settings_, - "[RootCut CPUFJ] ", - root_cut_cpufj_seed); - settings_.log.debug("Root cut CPUFJ final problem build time: %.6f seconds\n", - toc(root_cut_cpufj_build_start_time)); - f_t remaining_time = f_t(settings_.time_limit - toc(exploration_stats_.start_time)); - // Reserve at least half of the remaining time for B&B exploration; cap absolute spend - // at 1s so generous budgets don't grant CPUFJ more than the historical ceiling. - f_t fj_time_limit = - settings_.deterministic ? remaining_time : std::min(remaining_time * 0.5, 1.0); - root_fj_cpu_worker.run_sync(fj_time_limit, 0.5); - } - + // Stops the root heuristics and clear the associated data + root_heuristics.stop_and_sync(); set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); pc_.resize(original_lp_.num_cols); @@ -3724,9 +3786,23 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut const i_t num_bfs_workers = std::max(num_workers / 2, 1); const i_t num_submip_workers = std::max(num_workers / 8, 1); const i_t num_diving_workers = std::max(num_workers - num_bfs_workers, 1); - bfs_worker_pool_.init(num_bfs_workers, original_lp_, Arow_, var_types_, symmetry_, settings_); - rins_worker_pool_.init( - num_submip_workers, original_lp_, Arow_, var_types_, symmetry_, settings_, num_bfs_workers); + bfs_worker_pool_.init(num_bfs_workers, + original_lp_, + Arow_, + var_types_, + symmetry_, + settings_, + root_relax_soln_.x, + edge_norms_); + rins_worker_pool_.init(num_submip_workers, + original_lp_, + Arow_, + var_types_, + symmetry_, + settings_, + root_relax_soln_.x, + edge_norms_, + num_bfs_workers); if (num_diving_workers > 0) { diving_worker_pool_.init(num_diving_workers, @@ -3735,6 +3811,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut var_types_, symmetry_, settings_, + root_relax_soln_.x, + edge_norms_, num_bfs_workers + num_submip_workers); } @@ -3917,7 +3995,7 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri deterministic_global_termination_status_ = mip_status_t::UNSET; deterministic_workers_ = std::make_unique>( - num_bfs_workers, original_lp_, Arow, var_types_, settings_); + num_bfs_workers, original_lp_, Arow, var_types_, settings_, root_relax_soln_.x, edge_norms_); if (num_diving_workers > 0) { // Extract diving types from search_strategies (skip BEST_FIRST at index 0) @@ -3936,7 +4014,8 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri Arow, var_types_, settings_, - &root_relax_soln_.x); + root_relax_soln_.x, + edge_norms_); } } diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 96b8a6d8fe..fb3f9fce5d 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -34,12 +34,14 @@ #include #include +#include #include #include #include #include +#include #include #include @@ -368,24 +370,31 @@ class branch_and_bound_t { // Launch a new RINS worker bool launch_rins_worker(const std::vector& sol); - void set_solution_from_submip(const std::vector& solution, + void set_solution_from_submip(const simplex::lp_problem_t& lp, + const std::vector& solution, const third_party_presolve_t& presolver, - f_t fixrate, - f_t obj); + f_t fixrate); // Solve the RINS sub-MIP void solve_submip(diving_worker_t* worker, const std::vector& current_incumbent, + const std::vector& var_types, i_t num_var_fixed, i_t num_integers, i_t submip_level, - std::string_view log_prefix); + std::string_view log_prefix, + bool is_root_heuristic = false); // Creates and solves the RINS sub-MIP - void rins(diving_worker_t* rins_worker, const std::vector& node_solution); - - // Get the simplex settings for solving the LP of a single node - simplex::simplex_solver_settings_t get_node_lp_settings(); + void rins(diving_worker_t* worker, + const std::vector& current_incumbent, + const std::vector& var_types, + bool is_root_heuristic = false); + + void launch_root_heuristics(const simplex::lp_problem_t& lp, + const std::vector& sol, + i_t cut_pass, + root_heuristics_t& root_heuristics); // Solve the LP relaxation of a leaf node simplex::dual_status_t solve_node_lp(mip_node_t* node_ptr, diff --git a/cpp/src/branch_and_bound/constants.hpp b/cpp/src/branch_and_bound/constants.hpp index 5629360a86..2bc94c8f1b 100644 --- a/cpp/src/branch_and_bound/constants.hpp +++ b/cpp/src/branch_and_bound/constants.hpp @@ -27,15 +27,31 @@ enum class heuristics_origin_t { // doi: 10.1007/s10107-004-0518-7. enum class search_strategy_t : int { BEST_FIRST = 0, // Best-First + Plunging. - PSEUDOCOST_DIVING = 1, // Pseudocost diving (9.2.5) - LINE_SEARCH_DIVING = 2, // Line search diving (9.2.4) - GUIDED_DIVING = 3, // Guided diving (9.2.3). - COEFFICIENT_DIVING = 4, // Coefficient diving (9.2.1) + PSEUDOCOST_DIVING = 1, // Pseudocost diving [1, Section 9.2.5] + LINE_SEARCH_DIVING = 2, // Line search diving [1, Section 9.2.4] + GUIDED_DIVING = 3, // Guided diving. [1, Section 9.2.3] + COEFFICIENT_DIVING = 4, // Coefficient diving [1, Section 9.2.1] FARKAS_DIVING = 5, // Farkas Diving (see [2]) - VECTOR_LENGTH_DIVING = 6, // Vector Length Diving (9.2.6) - SUBMIP = 7 // RINS (see [3]) + VECTOR_LENGTH_DIVING = 6, // Vector Length Diving [1, Section 9.2.6] + RINS = 7, // RINS (see [3]) }; enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 }; +inline const char* search_strategy_to_string(search_strategy_t search_strategy) +{ + switch (search_strategy) { + case search_strategy_t::BEST_FIRST: return "BEST_FIRST"; + case search_strategy_t::PSEUDOCOST_DIVING: return "PSEUDOCOST_DIVING"; + case search_strategy_t::LINE_SEARCH_DIVING: return "LINE_SEARCH_DIVING"; + case search_strategy_t::GUIDED_DIVING: return "GUIDED_DIVING"; + case search_strategy_t::COEFFICIENT_DIVING: return "COEFFICIENT_DIVING"; + case search_strategy_t::FARKAS_DIVING: return "FARKAS_DIVING"; + case search_strategy_t::VECTOR_LENGTH_DIVING: return "VECTOR_LENGTH_DIVING"; + case search_strategy_t::RINS: return "RINS"; + } + + return "UNKNOWN"; +} + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/deterministic_workers.hpp b/cpp/src/branch_and_bound/deterministic_workers.hpp index 4de3086e61..7c426452a8 100644 --- a/cpp/src/branch_and_bound/deterministic_workers.hpp +++ b/cpp/src/branch_and_bound/deterministic_workers.hpp @@ -89,8 +89,10 @@ 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, + 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), + : base_t(id, original_lp, Arow, var_types, settings, root_solution, root_edge_norm), work_context(context_name), pc_snapshot(1, settings) { @@ -140,8 +142,17 @@ class deterministic_bfs_worker_t const simplex::lp_problem_t& original_lp, const csr_matrix_t& Arow, const std::vector& var_types, - const simplex::simplex_solver_settings_t& settings) - : base_t(id, original_lp, Arow, var_types, settings, "BB_Worker_" + std::to_string(id)) + const simplex::simplex_solver_settings_t& settings, + const std::vector& root_solution, + const std::vector& root_edge_norm) + : base_t(id, + original_lp, + Arow, + var_types, + settings, + root_solution, + root_edge_norm, + "BB_Worker_" + std::to_string(id)) { } @@ -282,9 +293,6 @@ class deterministic_diving_worker_t std::vector dive_lower; std::vector dive_upper; - // Root LP relaxation solution (constant, set once at construction) - const std::vector* root_solution{nullptr}; - // Diving state bool recompute_bounds_and_basis{true}; @@ -300,10 +308,17 @@ class deterministic_diving_worker_t const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, - const std::vector* root_sol) - : base_t(id, original_lp, Arow, var_types, settings, "Diving_Worker_" + std::to_string(id)), - diving_type(type), - root_solution(root_sol) + const std::vector& root_solution, + const std::vector& root_edge_norm) + : base_t(id, + original_lp, + Arow, + var_types, + settings, + root_solution, + root_edge_norm, + "Diving_Worker_" + std::to_string(id)), + diving_type(type) { dive_lower = original_lp.lower; dive_upper = original_lp.upper; @@ -407,11 +422,14 @@ class deterministic_bfs_worker_pool_t const simplex::lp_problem_t& original_lp, const csr_matrix_t& Arow, const std::vector& var_types, - const simplex::simplex_solver_settings_t& settings) + const simplex::simplex_solver_settings_t& settings, + 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); + this->workers_.emplace_back( + i, original_lp, Arow, var_types, settings, root_solution, root_edge_norm); } } @@ -443,12 +461,14 @@ class deterministic_diving_worker_pool_t const csr_matrix_t& Arow, const std::vector& var_types, const simplex::simplex_solver_settings_t& settings, - const std::vector* root_solution) + const std::vector& root_solution, + const std::vector& root_edge_norm) { this->workers_.reserve(num_workers); 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); + this->workers_.emplace_back( + i, type, original_lp, Arow, var_types, settings, root_solution, root_edge_norm); } } diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index fd5255f86f..50e9a4a206 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -78,6 +78,9 @@ class branch_and_bound_worker_t { bool recompute_basis = true; bool recompute_bounds = true; + const std::vector& root_solution; + const std::vector& root_edge_norm; + void ensure_orbital_fixing() { if (orbital_fixing == nullptr && symmetry_ptr != nullptr) { @@ -94,6 +97,8 @@ class branch_and_bound_worker_t { const csr_matrix_t& Arow, const std::vector& var_type, const simplex::simplex_solver_settings_t& settings, + const std::vector& root_solution, + const std::vector& root_edge_norm, uint64_t rng_offset = 0) : worker_id(worker_id), search_strategy(search_strategy_t::BEST_FIRST), @@ -108,7 +113,9 @@ class branch_and_bound_worker_t { node_presolver(leaf_problem, Arow, {}, var_type), bounds_changed(original_lp.num_cols, false), rng(settings.random_seed + pcgenerator_t::default_seed + rng_offset + worker_id, - pcgenerator_t::default_stream ^ (worker_id + rng_offset)) + pcgenerator_t::default_stream ^ (worker_id + rng_offset)), + root_solution(root_solution), + root_edge_norm(root_edge_norm) { } @@ -146,8 +153,11 @@ 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, + 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, rng_offset) + : Base( + worker_id, original_lp, Arow, var_type, settings, root_solution, root_edge_norm, rng_offset) { this->start_lower = original_lp.lower; this->start_upper = original_lp.upper; @@ -243,6 +253,8 @@ class diving_worker_t : public branch_and_bound_worker_t { // The best-first worker that is associated with this diving worker. Used for controlling the // number of active diving workers. bfs_worker_t* bfs_worker{nullptr}; + + std::atomic halt = false; }; struct submip_stats_t { diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp index e9d55bafe3..6977c7882b 100644 --- a/cpp/src/branch_and_bound/worker_pool.hpp +++ b/cpp/src/branch_and_bound/worker_pool.hpp @@ -24,6 +24,8 @@ class worker_pool_t { const std::vector& var_type, mip_symmetry_t* symmetry, const simplex::simplex_solver_settings_t& settings, + const std::vector& root_solution, + const std::vector& root_edge_norm, const uint64_t rng_offset = 0) { assert(!is_initialized_); @@ -33,8 +35,8 @@ class worker_pool_t { num_idle_workers_ = num_workers; 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, rng_offset); + workers_[i] = std::make_unique( + i, original_lp, Arow, var_type, settings, 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(). @@ -61,14 +63,17 @@ class worker_pool_t { void return_worker_to_pool(WorkerType* worker) { - std::lock_guard lock(mutex_); assert(worker != nullptr); + worker->set_inactive(); + assert(!worker->is_active.load()); + + if (!is_initialized_) return; + + std::lock_guard lock(mutex_); assert(workers_[worker->worker_id].get() == worker); assert(static_cast(num_idle_workers_.load()) == idle_workers_.size()); assert(idle_workers_.size() <= workers_.size()); - worker->set_inactive(); - assert(!worker->is_active.load()); idle_workers_.push_back(worker->worker_id); num_idle_workers_++; } diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 7907abd3b9..388bb43b35 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -104,6 +104,12 @@ f_t compute_user_objective(const lp_problem_t& lp, f_t obj) return user_obj; } +template +f_t compute_presolved_objective(const lp_problem_t& lp, f_t user_obj) +{ + return user_obj / lp.obj_scale - lp.obj_constant; +} + template lp_status_t solve_linear_program_advanced(const lp_problem_t& original_lp, const f_t start_time, @@ -813,6 +819,8 @@ template double compute_user_objective(const lp_problem_t& lp, double obj); +template double compute_presolved_objective(const lp_problem_t& lp, double user_obj); + template lp_status_t solve_linear_program_advanced( const lp_problem_t& original_lp, const double start_time, diff --git a/cpp/src/dual_simplex/solve.hpp b/cpp/src/dual_simplex/solve.hpp index 90c2dbd690..308c462de5 100644 --- a/cpp/src/dual_simplex/solve.hpp +++ b/cpp/src/dual_simplex/solve.hpp @@ -63,6 +63,9 @@ f_t compute_user_objective(const lp_problem_t& lp, const std::vector f_t compute_user_objective(const lp_problem_t& lp, f_t obj); +template +f_t compute_presolved_objective(const lp_problem_t& lp, f_t user_obj); + template lp_status_t solve_linear_program_advanced(const lp_problem_t& original_lp, const f_t start_time, diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 5d9b0267b1..b789159953 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -2109,37 +2109,52 @@ void fj_cpu_worker_t::create_worker( fj_cpu.reset(new_climber.release()); fj_cpu->log_prefix = std::move(log_prefix); fj_cpu->improvement_callback = improvement_callback; + fj_cpu->halted = false; + preemption_flag = false; + is_initialized = true; } template void fj_cpu_worker_t::run_async(f_t time_limit, double work_unit_limit) { - if (!fj_cpu) return; + if (!is_initialized) return; -#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit) \ - priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : *fj_cpu) - cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); + auto& fj_ptr = fj_cpu; +#pragma omp task shared(fj_cpu, is_initialized, fj_ptr) firstprivate(time_limit, work_unit_limit) \ + priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : fj_ptr) + { + if (is_initialized) { cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); } + } } template void fj_cpu_worker_t::run_sync(f_t time_limit, double work_unit_limit) { - if (!fj_cpu) return; + if (!is_initialized) return; cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); + is_initialized = false; fj_cpu.reset(); } template void fj_cpu_worker_t::stop() { - if (!fj_cpu) return; + if (!is_initialized) return; + + preemption_flag = true; - fj_cpu->preemption_flag = true; - fj_cpu->halted = true; -#pragma omp taskwait depend(in : *fj_cpu) + auto& fj_ptr = fj_cpu; +#pragma omp taskwait depend(in : fj_ptr) + is_initialized = false; fj_cpu.reset(); } +template +void fj_cpu_worker_t::send_stop_signal() +{ + preemption_flag = true; +} + #if MIP_INSTANTIATE_FLOAT template class fj_t; template struct fj_cpu_worker_t; diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh index ff6022c4c2..bb2c69f81c 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh @@ -10,6 +10,8 @@ #include #include +#include + #include #include #include @@ -28,6 +30,8 @@ struct fj_cpu_worker_t { struct fj_cpu_deleter_t { void operator()(fj_cpu_climber_t* ptr) const; }; + + std::atomic is_initialized{false}; std::atomic preemption_flag{false}; std::unique_ptr, fj_cpu_deleter_t> fj_cpu; std::function&, double)> improvement_callback; @@ -55,6 +59,8 @@ struct fj_cpu_worker_t { double work_unit_limit = std::numeric_limits::infinity()); void stop(); + + void send_stop_signal(); }; } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 3d324d2007..3d5b165abb 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -1293,9 +1293,9 @@ void third_party_presolve_t::undo(std::vector& primal_solution, template void third_party_presolve_t::uncrush_primal_solution( - const std::vector& reduced_primal, std::vector& full_primal) const + const std::vector& reduced_primal, std::vector& full_primal, bool check_postsolve) const { - if (presolver_ == cuopt::mathematical_optimization::presolver_t::PSLP) { + if (presolver_ == PSLP) { cuopt_expects(false, error_type_t::RuntimeError, "This code path should be never called, as this is meant for callbacks and they " @@ -1311,7 +1311,7 @@ void third_party_presolve_t::uncrush_primal_solution( bool is_optimal = false; auto status = post_solver.undo(reduced_sol, full_sol, *papilo_post_solve_storage_, is_optimal); - check_postsolve_status(status); + if (check_postsolve) check_postsolve_status(status); full_primal = std::move(full_sol.primal); } diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index 60f4f7fc6d..79eb28d01e 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -141,7 +141,8 @@ class third_party_presolve_t { bool dual_postsolve); void uncrush_primal_solution(const std::vector& reduced_primal, - std::vector& full_primal) const; + std::vector& full_primal, + bool check_postsolve = true) const; void crush_primal_solution(const optimization_problem_t& reduced_problem, const std::vector& original_primal, diff --git a/cpp/src/mip_heuristics/root_heuristics.hpp b/cpp/src/mip_heuristics/root_heuristics.hpp new file mode 100644 index 0000000000..f54aff7c1b --- /dev/null +++ b/cpp/src/mip_heuristics/root_heuristics.hpp @@ -0,0 +1,134 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include +#include "feasibility_jump/fj_cpu_worker.cuh" + +namespace cuopt::mathematical_optimization::mip { + +template +struct cut_pass_heuristics_t { + std::vector var_types_; + csr_matrix_t Arow_; + std::vector root_solution_; + std::vector root_edge_norm_; + + std::unique_ptr> submip_worker_; + 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) + : var_types_(var_types), + Arow_(Arow), + root_solution_(root_solution), + root_edge_norm_(root_edge_norm), + submip_worker_(nullptr) {}; + + ~cut_pass_heuristics_t() { stop_and_sync(); } + + void send_stop_signal() + { + fj_cpu_worker_.send_stop_signal(); + if (submip_worker_) { submip_worker_->halt = true; } + } + + void stop_and_sync() + { + fj_cpu_worker_.stop(); + if (submip_worker_) { + diving_worker_t* worker = submip_worker_.get(); + worker->halt = true; +#pragma omp taskwait depend(in : *worker) + submip_worker_.reset(); + } + } + + diving_worker_t* create_submip_worker( + i_t id, + const simplex::lp_problem_t& lp, + const simplex::simplex_solver_settings_t& settings, + f_t root_obj, + const std::vector& root_vstatus, + const std::vector& sol) + { + submip_worker_ = std::make_unique>( + id, lp, Arow_, var_types_, settings, 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; + submip_worker_->recompute_bounds = false; + submip_worker_->recompute_basis = true; + submip_worker_->search_strategy = search_strategy_t::RINS; + submip_worker_->set_active(); + + return submip_worker_.get(); + } +}; + +/// \brief Object Representing the heuristics run on the root node. +template +struct root_heuristics_t { + // List of the heuristics that run alongside a single cut pass. + // It holds the workers and all the necessary information. + // + // We use the `shared_ptr` here so the object is only destroyed when the task terminates + // (we declare the `shared_ptr` as firstprivate in the task, so they live until the end the + // task). In this way, we can send the stop signal, destroy the entry in the list and the + // object itself will be destroyed when all related tasks ends. + std::list>> cut_passes_heuristics_; + + // Count the number of active workers. Same reason as above. + std::shared_ptr> worker_count_; + i_t max_workers_; + + root_heuristics_t(i_t max_workers) + : worker_count_(std::make_shared>(0)), max_workers_(max_workers) + { + } + + ~root_heuristics_t() { stop_and_sync(); } + + void stop_and_sync() + { + for (auto& heuristic : cut_passes_heuristics_) { + heuristic->send_stop_signal(); + } + + for (auto& heuristic : cut_passes_heuristics_) { + heuristic->stop_and_sync(); + } + } + + 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) + { + // 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)); + } +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/tests/mip/miplib_test.cu b/cpp/tests/mip/miplib_test.cu index 924763c437..3bdd6f51ff 100644 --- a/cpp/tests/mip/miplib_test.cu +++ b/cpp/tests/mip/miplib_test.cu @@ -97,7 +97,7 @@ TEST(mip_solve, node_limit_test) { mip_solver_settings_t settings; settings.node_limit = 1000; - settings.time_limit = 60; + settings.time_limit = 120; settings.num_cpu_threads = 8; double expect_obj = 3.8151140644999992e+02;