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 505a3b61b..4464a156a 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 4dc6bc67a..f0177991a 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); @@ -1018,7 +1032,7 @@ 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]); return {branch_var, round_dir}; @@ -2100,7 +2114,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); } } @@ -2164,19 +2178,34 @@ bool branch_and_bound_t::launch_rins_worker(const std::vector& so if (!incumbent_.has_incumbent) return false; if (rins_worker_pool_.num_idle() == 0) return false; + bool is_root_heuristic = false; 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 clear 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->recompute_bounds = false; + worker->recompute_basis = true; + 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_, is_root_heuristic); } 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_, is_root_heuristic); } return true; @@ -2185,10 +2214,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_type, 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(); @@ -2203,11 +2234,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; @@ -2220,10 +2252,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), @@ -2241,12 +2273,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, @@ -2257,7 +2296,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_type, settings_, submip_problem); third_party_presolve_t presolver; f_t presolve_time_limit = std::min(0.1 * submip_settings.time_limit, 60.0); @@ -2275,13 +2314,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]; @@ -2293,26 +2331,25 @@ 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; } 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); @@ -2325,24 +2362,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_internal_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; }); } @@ -2380,8 +2417,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(), @@ -2395,7 +2433,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 @@ -2405,9 +2443,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; @@ -2530,11 +2568,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); @@ -2545,41 +2585,27 @@ void branch_and_bound_t::rins(diving_worker_t* rins_worker, const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; 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; @@ -2595,11 +2621,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; } @@ -2624,11 +2650,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; } @@ -2638,7 +2664,7 @@ 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, @@ -2649,21 +2675,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; } @@ -2680,14 +2706,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); @@ -2695,11 +2721,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; @@ -2710,14 +2736,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); - - rins_worker->start_lower = lower; - rins_worker->start_upper = upper; + worker->start_node.packed_vstatus = simplex::compress_vstatus(worker->leaf_vstatus); + worker->start_lower = lower; + 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; @@ -2730,21 +2753,32 @@ 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); } } @@ -2753,7 +2787,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, @@ -2767,7 +2801,89 @@ 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 (!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, + std::list>& heuristics, + omp_atomic_t& worker_count) +{ + if (settings_.deterministic) return; + if (settings_.num_threads < 2) return; + + // 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). + if (worker_count + 4 > settings_.num_threads && !heuristics.empty()) { + heuristics.erase(heuristics.begin()); + } + + bool is_root_heuristic = true; + i_t id = cut_pass; + root_heuristics_t* heuristic = &heuristics.emplace_back(Arow_, var_types_); + + constexpr bool is_cpufj_enabled = true; + if (is_cpufj_enabled) { + fj_cpu_worker_t* worker = &heuristic->fj_cpu_worker_; + + f_t work_limit = 2.0; + f_t time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + + worker->improvement_callback = + [this](f_t obj, const std::vector& assignment, double work_units) { + set_solution_from_cpu_fj(obj, assignment, work_units); + }; + worker->create_worker(lp, var_types_, sol, settings_, "[RootCut CPUFJ] "); + + ++worker_count; +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ + firstprivate(worker, heuristic) shared(heuristics, worker_count) depend(out : *worker -> fj_cpu) + { + worker->run_sync(time_limit, work_limit); + --worker_count; + } + } + + if (settings_.submip_settings.rins != 0 && incumbent_.has_incumbent) { + diving_worker_t* worker = + heuristic->create_submip_worker(id, 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, heuristic->var_types_, is_root_heuristic); + heuristic->Arow_ = csr_matrix_t(1, 1, 1); + heuristic->var_types_ = {}; + heuristic->submip_worker_.reset(); + + } else { + ++worker_count; +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ + shared(heuristics, worker_count) firstprivate(worker, current_incumbent, heuristic) \ + depend(out : *worker) + { + rins(worker, current_incumbent, heuristic->var_types_, is_root_heuristic); + --worker_count; + heuristic->Arow_ = csr_matrix_t(1, 1, 1); + heuristic->var_types_ = {}; + heuristic->submip_worker_.reset(); + } + } + } } template @@ -3316,13 +3432,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) @@ -3330,7 +3447,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(); @@ -3341,7 +3458,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); @@ -3350,7 +3466,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); @@ -3359,7 +3474,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); @@ -3368,7 +3482,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, @@ -3447,12 +3560,9 @@ 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); - }; + omp_atomic_t root_worker_count = 0; + std::list> root_heuristics; + launch_root_heuristics(original_lp_, root_relax_soln_.x, 0, root_heuristics, root_worker_count); f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; @@ -3485,8 +3595,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } cut_pass_result_t cut_pass_result; - root_fj_cpu_worker.run_async(settings_.time_limit - toc(exploration_stats_.start_time)); - cut_pass_result = do_cut_pass(cut_pass, solution, num_fractional, @@ -3505,7 +3613,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) { @@ -3517,15 +3633,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } 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)); - } + launch_root_heuristics( + original_lp_, root_relax_soln_.x, cut_pass + 1, root_heuristics, root_worker_count); } // Publish the post-cuts root LP value. @@ -3541,18 +3650,6 @@ 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("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", @@ -3563,29 +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.clear(); set_uninitialized_steepest_edge_norms(original_lp_, basic_list, 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 96b8a6d8f..75d98f7d6 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,32 @@ 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_type, 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); // 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); + + void launch_root_heuristics(const simplex::lp_problem_t& lp, + const std::vector& sol, + i_t cut_pass, + std::list>& heuristics, + omp_atomic_t& worker_count); // 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 5629360a8..2bc94c8f1 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/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index fd5255f86..f40dcb71b 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -243,6 +243,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/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 7907abd3b..9792ff72d 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_internal_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_internal_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 90c2dbd69..56d032ce9 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_internal_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 57a6a8947..134499bf9 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -1814,13 +1814,20 @@ void fj_cpu_worker_t::create_worker( } template -void fj_cpu_worker_t::run_async(f_t time_limit, double work_unit_limit) +void fj_cpu_worker_t::run_async(f_t time_limit, + double work_unit_limit, + omp_atomic_t* worker_count) { if (!fj_cpu) return; -#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit) \ + if (worker_count) ++(*worker_count); +#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit, worker_count) \ priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : *fj_cpu) - cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); + { + cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); + fj_cpu.reset(); + if (worker_count) --(*worker_count); + } } template 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 ff6022c4c..d594eef8d 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 @@ -47,8 +49,9 @@ struct fj_cpu_worker_t { // Run the worker asynchronously (i.e., launch an openmp task and then continue the // execution). Call `stop()` for stopping the worker - void run_async(f_t time_limit = std::numeric_limits::infinity(), - double work_unit_limit = std::numeric_limits::infinity()); + void run_async(f_t time_limit = std::numeric_limits::infinity(), + double work_unit_limit = std::numeric_limits::infinity(), + omp_atomic_t* worker_count = nullptr); // Run the CPU FJ synchronously (i.e., wait for it to finish before proceeding) void run_sync(f_t time_limit = std::numeric_limits::infinity(), diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 5fd233309..3182839d2 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -1290,9 +1290,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 " @@ -1308,7 +1308,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 60f4f7fc6..79eb28d01 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 000000000..78a62819a --- /dev/null +++ b/cpp/src/mip_heuristics/root_heuristics.hpp @@ -0,0 +1,60 @@ +/* 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 root_heuristics_t { + std::unique_ptr> submip_worker_; + fj_cpu_worker_t fj_cpu_worker_; + std::vector var_types_; + csr_matrix_t Arow_; + + root_heuristics_t(const csr_matrix_t& Arow, + const std::vector& var_types) + : submip_worker_(nullptr), var_types_(var_types), Arow_(Arow) {}; + + ~root_heuristics_t() { stop(); } + + void stop() + { + fj_cpu_worker_.stop(); + if (submip_worker_) { + submip_worker_->halt = true; + diving_worker_t* worker = submip_worker_.get(); +#pragma omp taskwait depend(in : *worker) + } + } + + 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); + 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(); + } +}; +} // namespace cuopt::mathematical_optimization::mip