diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index ba80244b3b..1acc16af54 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include @@ -27,6 +28,7 @@ #include #include +#include #include @@ -2023,6 +2025,283 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( return root_status; } +template +auto branch_and_bound_t::do_cut_pass( + [[maybe_unused]] i_t cut_pass, + mip_solution_t& solution, + i_t& num_fractional, + std::vector& fractional, + cut_generation_t& cut_generation, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + variable_bounds_t& variable_bounds, + cut_pool_t& cut_pool, + cut_info_t& cut_info, + simplex_solver_settings_t& lp_settings, + i_t original_rows, + f_t& last_upper_bound, + f_t& last_objective, + f_t root_relax_objective, + i_t& cut_pool_size, + [[maybe_unused]] const std::vector& saved_solution) -> cut_pass_result_t +{ +#ifdef PRINT_FRACTIONAL_INFO + settings_.log.printf("Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); + for (i_t j : fractional) { + settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", + j, + original_lp_.lower[j], + root_relax_soln_.x[j], + original_lp_.upper[j]); + } +#endif + + f_t cut_start_time = tic(); + bool problem_feasible = cut_generation.generate_cuts(original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + basis_update, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + variable_bounds, + exploration_stats_.start_time); + if (!problem_feasible) { + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; + } + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); + } + // Score the cuts + f_t score_start_time = tic(); + cut_pool.score_cuts(root_relax_soln_.x); + f_t score_time = toc(score_start_time); + if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } + // Get the best cuts from the cut pool + csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); + std::vector cut_rhs; + std::vector cut_types; + i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); + if (num_cuts == 0) { return {cut_pass_action_t::BREAK, mip_status_t::UNSET}; } + cut_info.record_cut_types(cut_types); +#ifdef PRINT_CUT_POOL_TYPES + cut_pool.print_cutpool_types(); + print_cut_types("In LP ", cut_types, settings_); + printf("Cut pool size: %d\n", cut_pool.pool_size()); +#endif + +#ifdef CHECK_CUT_MATRIX + if (cuts_to_add.check_matrix() != 0) { + settings_.log.printf("Bad cuts matrix\n"); + for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { + settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); + } + return {cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } +#endif +#ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION + verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); +#endif + cut_pool_size = cut_pool.pool_size(); + + // Resolve the LP with the new cuts + settings_.log.debug( + "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", + num_cuts, + cuts_to_add.row_start[cuts_to_add.m], + cut_pool.pool_size(), + cuts_to_add.m + original_lp_.num_rows); + lp_settings.log.log = false; + + f_t add_cuts_start_time = tic(); + mutex_original_lp_.lock(); + i_t add_cuts_status = add_cuts(settings_, + cuts_to_add, + cut_rhs, + original_lp_, + new_slacks_, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + variable_bounds.resize(original_lp_.num_cols); + mutex_original_lp_.unlock(); + f_t add_cuts_time = toc(add_cuts_start_time); + if (add_cuts_time > 1.0) { settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); } + if (add_cuts_status != 0) { + settings_.log.printf("Failed to add cuts\n"); + return {cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } + + if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { + mutex_upper_.lock(); + last_upper_bound = upper_bound_.load(); + std::vector lower_bounds; + std::vector upper_bounds; + find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + mutex_upper_.unlock(); + mutex_original_lp_.lock(); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; + mutex_original_lp_.unlock(); + } + + // Try to do bound strengthening + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; +#ifdef CHECK_MATRICES + settings_.log.printf("Before A check\n"); + original_lp_.A.check_matrix(); +#endif + original_lp_.A.to_compressed_row(Arow_); + + f_t node_presolve_start_time = tic(); + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + std::vector new_lower = original_lp_.lower; + std::vector new_upper = original_lp_.upper; + bool feasible = + node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + f_t node_presolve_time = toc(node_presolve_start_time); + if (node_presolve_time > 1.0) { + settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); + } + if (!feasible) { + settings_.log.printf("Bound strengthening detected infeasibility\n"); +#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS + original_lp_.write_mps("bound_strengthening_infeasible.mps"); +#endif + return {cut_pass_action_t::RETURN, mip_status_t::INFEASIBLE}; + } + + i_t iter = 0; + bool initialize_basis = false; + lp_settings.concurrent_halt = NULL; + f_t dual_phase2_start_time = tic(); + dual::status_t cut_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + exploration_stats_.total_lp_iters += iter; + f_t dual_phase2_time = toc(dual_phase2_start_time); + if (dual_phase2_time > 1.0) { + settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); + } + if (cut_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, solver_status_}; + } + + if (cut_status != dual::status_t::OPTIMAL) { + settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); + lp_status_t scratch_status = + solve_linear_program_with_advanced_basis(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + if (scratch_status == lp_status_t::OPTIMAL) { + // We recovered + cut_status = convert_lp_status_to_dual_status(scratch_status); + exploration_stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else { + settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); +#ifdef WRITE_CUT_INFEASIBLE_MPS + original_lp_.write_mps("cut_infeasible.mps"); +#endif + return {cut_pass_action_t::RETURN, mip_status_t::NUMERICAL}; + } + } + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + + f_t remove_cuts_start_time = tic(); + mutex_original_lp_.lock(); + remove_cuts(original_lp_, + settings_, + exploration_stats_.start_time, + Arow_, + new_slacks_, + original_rows, + var_types_, + root_vstatus_, + edge_norms_, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + basis_update); + variable_bounds.resize(original_lp_.num_cols); + mutex_original_lp_.unlock(); + f_t remove_cuts_time = toc(remove_cuts_start_time); + if (remove_cuts_time > 1.0) { + settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); + } + fractional.clear(); + num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + + if (num_fractional == 0) { + upper_bound_ = root_objective_; + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + mutex_upper_.unlock(); + } + f_t obj = upper_bound_.load(); + report(' ', obj, root_objective_, 0, num_fractional); + + f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); + f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_); + if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { + if (num_fractional == 0) { set_solution_at_root(solution, cut_info); } + set_final_solution(solution, root_objective_); + return {cut_pass_action_t::RETURN, mip_status_t::OPTIMAL}; + } + + f_t change_in_objective = root_objective_ - last_objective; + const f_t factor = settings_.cut_change_threshold; + const f_t min_objective = 1e-3; + if (factor > 0.0 && + change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { + settings_.log.printf( + "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", + change_in_objective, + root_relax_objective); + return {cut_pass_action_t::BREAK, mip_status_t::UNSET}; + } + last_objective = root_objective_; + return {cut_pass_action_t::CONTINUE, mip_status_t::UNSET}; +} + template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { @@ -2228,273 +2507,92 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t last_objective = root_objective_; f_t root_relax_objective = root_objective_; + constexpr bool enable_root_cut_cpufj = true; + std::unique_ptr> root_cut_cpufj_task; + auto root_cut_cpufj_improvement_callback = + [this](f_t obj, const std::vector& assignment, double work_units) { + std::vector user_assignment; + mutex_original_lp_.lock(); + uncrush_primal_solution(original_problem_, original_lp_, assignment, user_assignment); + mutex_original_lp_.unlock(); + settings_.log.debug("Root cut CPUFJ found solution with objective %.16e\n", obj); + // In deterministic mode the solution must be ordered by its work-unit timestamp so + // B&B sees incumbents in a reproducible sequence; otherwise apply it immediately. + if (settings_.deterministic) { + queue_external_solution_deterministic(user_assignment, work_units); + } else { + set_new_solution(user_assignment); + } + }; + auto stop_root_cut_cpufj = [&]() { + if (!root_cut_cpufj_task) { return; } + detail::stop_fj_cpu_task(*root_cut_cpufj_task); + root_cut_cpufj_task.reset(); + }; + cuopt::scope_guard root_cut_cpufj_guard([&]() { stop_root_cut_cpufj(); }); + f_t cut_generation_start_time = tic(); i_t cut_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { if (num_fractional == 0) { set_solution_at_root(solution, cut_info); + signal_extend_cliques_.store(true, std::memory_order_release); +#pragma omp taskwait depend(in : *clique_signal) return mip_status_t::OPTIMAL; - } else { -#ifdef PRINT_FRACTIONAL_INFO - settings_.log.printf( - "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); - for (i_t j : fractional) { - settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", - j, - original_lp_.lower[j], - root_relax_soln_.x[j], - original_lp_.upper[j]); - } -#endif + } - // Generate cuts and add them to the cut pool - f_t cut_start_time = tic(); - bool problem_feasible = cut_generation.generate_cuts(original_lp_, - settings_, - Arow_, - new_slacks_, - var_types_, - basis_update, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - variable_bounds, - exploration_stats_.start_time); - if (!problem_feasible) { - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); - } - signal_extend_cliques_.store(true, std::memory_order_release); + cut_pass_result_t cut_pass_result; + if (root_cut_cpufj_task) { +#pragma omp task shared(root_cut_cpufj_task) default(none) depend(out : *root_cut_cpufj_task) + detail::run_fj_cpu_task(*root_cut_cpufj_task, + std::numeric_limits::infinity(), + std::numeric_limits::infinity()); + } + + cut_pass_result = do_cut_pass(cut_pass, + solution, + num_fractional, + fractional, + cut_generation, + basis_update, + basic_list, + nonbasic_list, + variable_bounds, + cut_pool, + cut_info, + lp_settings, + original_rows, + last_upper_bound, + last_objective, + root_relax_objective, + cut_pool_size, + saved_solution); + + if (root_cut_cpufj_task) { + detail::stop_fj_cpu_task(*root_cut_cpufj_task); +#pragma omp taskwait depend(in : *root_cut_cpufj_task) + } + + if (cut_pass_result.action == cut_pass_action_t::RETURN) { + signal_extend_cliques_.store(true, std::memory_order_release); #pragma omp taskwait depend(in : *clique_signal) - return mip_status_t::INFEASIBLE; - } - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); - } - // Score the cuts - f_t score_start_time = tic(); - cut_pool.score_cuts(root_relax_soln_.x); - f_t score_time = toc(score_start_time); - if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } - // Get the best cuts from the cut pool - csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); - std::vector cut_rhs; - std::vector cut_types; - i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); - if (num_cuts == 0) { break; } - cut_info.record_cut_types(cut_types); -#ifdef PRINT_CUT_POOL_TYPES - cut_pool.print_cutpool_types(); - print_cut_types("In LP ", cut_types, settings_); - printf("Cut pool size: %d\n", cut_pool.pool_size()); -#endif - -#ifdef CHECK_CUT_MATRIX - if (cuts_to_add.check_matrix() != 0) { - settings_.log.printf("Bad cuts matrix\n"); - for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { - settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); - } - return mip_status_t::NUMERICAL; - } -#endif - // Check against saved solution -#ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); -#endif - cut_pool_size = cut_pool.pool_size(); - - // Resolve the LP with the new cuts - settings_.log.debug( - "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", - num_cuts, - cuts_to_add.row_start[cuts_to_add.m], - cut_pool.pool_size(), - cuts_to_add.m + original_lp_.num_rows); - lp_settings.log.log = false; - - f_t add_cuts_start_time = tic(); - mutex_original_lp_.lock(); - i_t add_cuts_status = add_cuts(settings_, - cuts_to_add, - cut_rhs, - original_lp_, - new_slacks_, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); - variable_bounds.resize(original_lp_.num_cols); - mutex_original_lp_.unlock(); - f_t add_cuts_time = toc(add_cuts_start_time); - if (add_cuts_time > 1.0) { - settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); - } - if (add_cuts_status != 0) { - settings_.log.printf("Failed to add cuts\n"); - return mip_status_t::NUMERICAL; - } - - if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { - mutex_upper_.lock(); - last_upper_bound = upper_bound_.load(); - std::vector lower_bounds; - std::vector upper_bounds; - find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - mutex_upper_.unlock(); - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - mutex_original_lp_.unlock(); - } - - // Try to do bound strengthening - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; -#ifdef CHECK_MATRICES - settings_.log.printf("Before A check\n"); - original_lp_.A.check_matrix(); -#endif - original_lp_.A.to_compressed_row(Arow_); - - f_t node_presolve_start_time = tic(); - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - std::vector new_lower = original_lp_.lower; - std::vector new_upper = original_lp_.upper; - bool feasible = - node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); - mutex_original_lp_.lock(); - original_lp_.lower = new_lower; - original_lp_.upper = new_upper; - mutex_original_lp_.unlock(); - f_t node_presolve_time = toc(node_presolve_start_time); - if (node_presolve_time > 1.0) { - settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); - } - if (!feasible) { - settings_.log.printf("Bound strengthening detected infeasibility\n"); -#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS - original_lp_.write_mps("bound_strengthening_infeasible.mps"); -#endif - return mip_status_t::INFEASIBLE; - } - - i_t iter = 0; - bool initialize_basis = false; - lp_settings.concurrent_halt = NULL; - f_t dual_phase2_start_time = tic(); - dual::status_t cut_status = dual_phase2_with_advanced_basis(2, - 0, - initialize_basis, - exploration_stats_.start_time, - original_lp_, - lp_settings, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - root_relax_soln_, - iter, - edge_norms_); - exploration_stats_.total_lp_iters += iter; - f_t dual_phase2_time = toc(dual_phase2_start_time); - if (dual_phase2_time > 1.0) { - settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); - } - if (cut_status == dual::status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return solver_status_; - } - - if (cut_status != dual::status_t::OPTIMAL) { - settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); - lp_status_t scratch_status = - solve_linear_program_with_advanced_basis(original_lp_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - if (scratch_status == lp_status_t::OPTIMAL) { - // We recovered - cut_status = convert_lp_status_to_dual_status(scratch_status); - exploration_stats_.total_lp_iters += root_relax_soln_.iterations; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - } else { - settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); -#ifdef WRITE_CUT_INFEASIBLE_MPS - original_lp_.write_mps("cut_infeasible.mps"); -#endif - return mip_status_t::NUMERICAL; - } - } - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - - f_t remove_cuts_start_time = tic(); - mutex_original_lp_.lock(); - remove_cuts(original_lp_, - settings_, - exploration_stats_.start_time, - Arow_, - new_slacks_, - original_rows, - var_types_, - root_vstatus_, - edge_norms_, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - basis_update); - variable_bounds.resize(original_lp_.num_cols); - mutex_original_lp_.unlock(); - f_t remove_cuts_time = toc(remove_cuts_start_time); - if (remove_cuts_time > 1.0) { - settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); - } - fractional.clear(); - num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - - if (num_fractional == 0) { - upper_bound_ = root_objective_; - mutex_upper_.lock(); - incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); - mutex_upper_.unlock(); - } - f_t obj = upper_bound_.load(); - report(' ', obj, root_objective_, 0, num_fractional); - - f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); - f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_); - if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { - if (num_fractional == 0) { set_solution_at_root(solution, cut_info); } - set_final_solution(solution, root_objective_); - return mip_status_t::OPTIMAL; - } + return cut_pass_result.status; + } + if (cut_pass_result.action == cut_pass_action_t::BREAK) { break; } - f_t change_in_objective = root_objective_ - last_objective; - const f_t factor = settings_.cut_change_threshold; - const f_t min_objective = 1e-3; - if (factor > 0.0 && - change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { - settings_.log.printf( - "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", - change_in_objective, - root_relax_objective); - break; - } - last_objective = root_objective_; + 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_cut_cpufj_task = + detail::make_fj_cpu_task_from_host_lp(original_lp_, + var_types_, + root_relax_soln_.x, + settings_, + root_cut_cpufj_improvement_callback, + "[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)); } } @@ -2509,6 +2607,33 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_.A.col_start[original_lp_.A.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_cut_cpufj_task = + detail::make_fj_cpu_task_from_host_lp(original_lp_, + var_types_, + root_relax_soln_.x, + settings_, + root_cut_cpufj_improvement_callback, + "[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 * f_t{0.5}, f_t{1}); + detail::run_fj_cpu_task(*root_cut_cpufj_task, fj_time_limit, 0.5); + root_cut_cpufj_task.reset(); + } + set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); pc_.resize(original_lp_.num_cols); @@ -3084,6 +3209,19 @@ void branch_and_bound_t::deterministic_sync_callback() f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); + // Apply limit-based statuses first so a definitive answer (gap closure or tree exhaustion) + // detected in the same callback can override them. Otherwise a long producer wait that + // pushes the wall clock past time_limit would clobber a true INFEASIBLE/OPTIMAL conclusion + // and the solver would report TIME_LIMIT for an already-solved instance. + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + deterministic_global_termination_status_ = mip_status_t::TIME_LIMIT; + } + + // Stop early if next horizon exceeds work limit + if (deterministic_current_horizon_ > settings_.work_limit) { + deterministic_global_termination_status_ = mip_status_t::WORK_LIMIT; + } + if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) { deterministic_global_termination_status_ = mip_status_t::OPTIMAL; } @@ -3097,15 +3235,6 @@ void branch_and_bound_t::deterministic_sync_callback() } } - if (toc(exploration_stats_.start_time) > settings_.time_limit) { - deterministic_global_termination_status_ = mip_status_t::TIME_LIMIT; - } - - // Stop early if next horizon exceeds work limit - if (deterministic_current_horizon_ > settings_.work_limit) { - deterministic_global_termination_status_ = mip_status_t::WORK_LIMIT; - } - // Signal shutdown to prevent threads from entering barriers after termination if (deterministic_global_termination_status_ != mip_status_t::UNSET) { deterministic_scheduler_->signal_shutdown(); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index ae1a225e9a..bb4e7a1040 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -269,6 +269,31 @@ class branch_and_bound_t { i_t node_int_infeas, double work_time = -1); + enum class cut_pass_action_t { CONTINUE, BREAK, RETURN }; + struct cut_pass_result_t { + cut_pass_action_t action{cut_pass_action_t::CONTINUE}; + mip_status_t status{mip_status_t::UNSET}; + }; + + cut_pass_result_t do_cut_pass(i_t cut_pass, + mip_solution_t& solution, + i_t& num_fractional, + std::vector& fractional, + cut_generation_t& cut_generation, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + variable_bounds_t& variable_bounds, + cut_pool_t& cut_pool, + cut_info_t& cut_info, + simplex_solver_settings_t& lp_settings, + i_t original_rows, + f_t& last_upper_bound, + f_t& last_objective, + f_t root_relax_objective, + i_t& cut_pool_size, + const std::vector& saved_solution); + // Set the solution when found at the root node void set_solution_at_root(mip_solution_t& solution, const cut_info_t& cut_info); diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh new file mode 100644 index 0000000000..040674e47a --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh @@ -0,0 +1,56 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace cuopt::linear_programming::detail { + +template +struct fj_cpu_climber_t; + +template +struct fj_cpu_task_t { + struct fj_cpu_deleter_t { + void operator()(fj_cpu_climber_t* ptr) const; + }; + std::atomic preemption_flag{false}; + std::unique_ptr, fj_cpu_deleter_t> fj_cpu; +}; + +// `seed` selects the FJ RNG seed: pass a non-negative value for a deterministic seed, +// or -1 to draw from the global cuopt::seed_generator (the historical behavior). +// In deterministic mode the caller MUST pass an explicit seed, otherwise the underlying +// seed_generator::get_seed() racing with concurrent callers breaks reproducibility. +template +std::unique_ptr> make_fj_cpu_task_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::function&, double)> improvement_callback, + std::string log_prefix, + int64_t seed = -1); + +template +void run_fj_cpu_task(fj_cpu_task_t& task, + f_t time_limit = std::numeric_limits::infinity(), + double work_unit_limit = std::numeric_limits::infinity()); + +template +void stop_fj_cpu_task(fj_cpu_task_t& task); + +} // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 60cf271a55..575228895b 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -7,6 +7,10 @@ #include +#include +#include + +#include "cpu_fj_thread.cuh" #include "feasibility_jump.cuh" #include "feasibility_jump_impl_common.cuh" #include "fj_cpu.cuh" @@ -18,7 +22,9 @@ #include #include +#include #include +#include #include #include #include @@ -41,6 +47,15 @@ namespace cuopt::linear_programming::detail { +template +void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + i_t n_variables, + i_t n_constraints, + i_t n_integer_vars, + i_t nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances); + template thrust::tuple get_mtm_for_bound(const typename fj_t::climber_data_t::view_t& fj, i_t var_idx, @@ -792,9 +807,8 @@ static void apply_move(fj_cpu_climber_t& fj_cpu, fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; fj_cpu.h_best_assignment = fj_cpu.h_assignment; fj_cpu.iterations_since_best = 0; - CUOPT_LOG_TRACE("%sCPUFJ: new best objective: %g", - fj_cpu.log_prefix.c_str(), - fj_cpu.pb_ptr->get_user_obj_from_solver_obj(fj_cpu.h_incumbent_objective)); + CUOPT_LOG_TRACE( + "%sCPUFJ: new best objective: %g", fj_cpu.log_prefix.c_str(), fj_cpu.h_incumbent_objective); if (fj_cpu.improvement_callback) { double current_work_units = fj_cpu.work_units_elapsed.load(std::memory_order_acquire); fj_cpu.improvement_callback( @@ -829,7 +843,6 @@ static thrust::tuple find_mtm_move( fj_cpu_climber_t& fj_cpu, const std::vector& target_cstrs, bool localmin = false) { CPUFJ_NVTX_RANGE("CPUFJ::find_mtm_move"); - auto& problem = *fj_cpu.pb_ptr; raft::random::PCGenerator rng(fj_cpu.settings.seed + fj_cpu.iterations, 0, 0); @@ -1258,33 +1271,29 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, fj_cpu.h_tabu_lastinc.resize(fj_cpu.pb_ptr->n_variables, 0); fj_cpu.iterations = 0; - // set pointers to host copies - // technically not 'device_span's but raft doesn't have a universal span. - // cuda::std::span? - fj_cpu.view.cstr_left_weights = - raft::device_span(fj_cpu.h_cstr_left_weights.data(), fj_cpu.h_cstr_left_weights.size()); - fj_cpu.view.cstr_right_weights = - raft::device_span(fj_cpu.h_cstr_right_weights.data(), fj_cpu.h_cstr_right_weights.size()); - fj_cpu.view.objective_weight = &fj_cpu.h_objective_weight; - fj_cpu.view.incumbent_assignment = - raft::device_span(fj_cpu.h_assignment.data(), fj_cpu.h_assignment.size()); - fj_cpu.view.incumbent_lhs = raft::device_span(fj_cpu.h_lhs.data(), fj_cpu.h_lhs.size()); - fj_cpu.view.incumbent_lhs_sumcomp = - raft::device_span(fj_cpu.h_lhs_sumcomp.data(), fj_cpu.h_lhs_sumcomp.size()); - fj_cpu.view.tabu_nodec_until = - raft::device_span(fj_cpu.h_tabu_nodec_until.data(), fj_cpu.h_tabu_nodec_until.size()); - fj_cpu.view.tabu_noinc_until = - raft::device_span(fj_cpu.h_tabu_noinc_until.data(), fj_cpu.h_tabu_noinc_until.size()); - fj_cpu.view.tabu_lastdec = - raft::device_span(fj_cpu.h_tabu_lastdec.data(), fj_cpu.h_tabu_lastdec.size()); - fj_cpu.view.tabu_lastinc = - raft::device_span(fj_cpu.h_tabu_lastinc.data(), fj_cpu.h_tabu_lastinc.size()); - fj_cpu.view.objective_vars = - raft::device_span(fj_cpu.h_objective_vars.data(), fj_cpu.h_objective_vars.size()); - fj_cpu.view.incumbent_objective = &fj_cpu.h_incumbent_objective; - fj_cpu.view.best_objective = &fj_cpu.h_best_objective; + finalize_fj_cpu_host_initialization(fj_cpu, + problem.n_variables, + problem.n_constraints, + problem.n_integer_vars, + problem.nnz, + problem.tolerances); +} + +template +static void set_host_data_view( + fj_cpu_climber_t& fj_cpu, + i_t n_variables, + i_t n_constraints, + i_t n_integer_vars, + i_t nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + fj_cpu.view.pb.tolerances = tolerances; + fj_cpu.view.pb.n_variables = n_variables; + fj_cpu.view.pb.n_integer_vars = n_integer_vars; + fj_cpu.view.pb.n_constraints = n_constraints; + fj_cpu.view.pb.nnz = nnz; - fj_cpu.view.settings = &fj_cpu.settings; fj_cpu.view.pb.constraint_lower_bounds = raft::device_span(fj_cpu.h_cstr_lb.data(), fj_cpu.h_cstr_lb.size()); fj_cpu.view.pb.constraint_upper_bounds = @@ -1295,6 +1304,8 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, raft::device_span(fj_cpu.h_var_types.data(), fj_cpu.h_var_types.size()); fj_cpu.view.pb.is_binary_variable = raft::device_span(fj_cpu.h_is_binary_variable.data(), fj_cpu.h_is_binary_variable.size()); + fj_cpu.view.pb.binary_indices = + raft::device_span(fj_cpu.h_binary_indices.data(), fj_cpu.h_binary_indices.size()); fj_cpu.view.pb.coefficients = raft::device_span(fj_cpu.h_coefficients.data(), fj_cpu.h_coefficients.size()); fj_cpu.view.pb.offsets = raft::device_span(fj_cpu.h_offsets.data(), fj_cpu.h_offsets.size()); @@ -1308,13 +1319,61 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, raft::device_span(fj_cpu.h_reverse_offsets.data(), fj_cpu.h_reverse_offsets.size()); fj_cpu.view.pb.objective_coefficients = raft::device_span(fj_cpu.h_obj_coeffs.data(), fj_cpu.h_obj_coeffs.size()); - fj_cpu.h_objective_vars.resize(problem.n_variables); +} + +template +void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + i_t n_variables, + i_t n_constraints, + i_t n_integer_vars, + i_t nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + raft::common::nvtx::range scope("finalize_fj_cpu_host_initialization"); + + cuopt_assert(n_variables >= 0, "invalid variable count"); + cuopt_assert(n_constraints >= 0, "invalid constraint count"); + cuopt_assert(fj_cpu.h_offsets.size() == static_cast(n_constraints + 1), + "invalid CSR offsets"); + cuopt_assert(fj_cpu.h_reverse_offsets.size() == static_cast(n_variables + 1), + "invalid reverse offsets"); + cuopt_assert(fj_cpu.h_assignment.size() == static_cast(n_variables), + "seed assignment size mismatch"); + + set_host_data_view(fj_cpu, n_variables, n_constraints, n_integer_vars, nnz, tolerances); + + fj_cpu.view.cstr_left_weights = + raft::device_span(fj_cpu.h_cstr_left_weights.data(), fj_cpu.h_cstr_left_weights.size()); + fj_cpu.view.cstr_right_weights = + raft::device_span(fj_cpu.h_cstr_right_weights.data(), fj_cpu.h_cstr_right_weights.size()); + fj_cpu.view.objective_weight = &fj_cpu.h_objective_weight; + fj_cpu.view.incumbent_assignment = + raft::device_span(fj_cpu.h_assignment.data(), fj_cpu.h_assignment.size()); + fj_cpu.view.incumbent_lhs = raft::device_span(fj_cpu.h_lhs.data(), fj_cpu.h_lhs.size()); + fj_cpu.view.incumbent_lhs_sumcomp = + raft::device_span(fj_cpu.h_lhs_sumcomp.data(), fj_cpu.h_lhs_sumcomp.size()); + fj_cpu.view.tabu_nodec_until = + raft::device_span(fj_cpu.h_tabu_nodec_until.data(), fj_cpu.h_tabu_nodec_until.size()); + fj_cpu.view.tabu_noinc_until = + raft::device_span(fj_cpu.h_tabu_noinc_until.data(), fj_cpu.h_tabu_noinc_until.size()); + fj_cpu.view.tabu_lastdec = + raft::device_span(fj_cpu.h_tabu_lastdec.data(), fj_cpu.h_tabu_lastdec.size()); + fj_cpu.view.tabu_lastinc = + raft::device_span(fj_cpu.h_tabu_lastinc.data(), fj_cpu.h_tabu_lastinc.size()); + fj_cpu.view.incumbent_objective = &fj_cpu.h_incumbent_objective; + fj_cpu.view.best_objective = &fj_cpu.h_best_objective; + fj_cpu.view.settings = &fj_cpu.settings; + + fj_cpu.h_objective_vars.resize(n_variables); auto end = std::copy_if( thrust::counting_iterator(0), - thrust::counting_iterator(problem.n_variables), + thrust::counting_iterator(n_variables), fj_cpu.h_objective_vars.begin(), [&fj_cpu](i_t idx) { return !fj_cpu.view.pb.integer_equal(fj_cpu.h_obj_coeffs[idx], (f_t)0); }); fj_cpu.h_objective_vars.resize(end - fj_cpu.h_objective_vars.begin()); + fj_cpu.view.objective_vars = + raft::device_span(fj_cpu.h_objective_vars.data(), fj_cpu.h_objective_vars.size()); fj_cpu.h_best_objective = +std::numeric_limits::infinity(); @@ -1323,7 +1382,7 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, std::make_pair(0, fj_staged_score_t::zero())); fj_cpu.cached_cstr_bounds.resize(fj_cpu.h_reverse_coefficients.size()); - for (i_t var_idx = 0; var_idx < (i_t)fj_cpu.view.pb.n_variables; ++var_idx) { + for (i_t var_idx = 0; var_idx < n_variables; ++var_idx) { auto [offset_begin, offset_end] = reverse_range_for_var(fj_cpu, var_idx); for (i_t i = offset_begin; i < offset_end; ++i) { fj_cpu.cached_cstr_bounds[i] = @@ -1332,9 +1391,9 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, } } - fj_cpu.flip_move_computed.resize(fj_cpu.view.pb.n_variables, false); - fj_cpu.var_bitmap.resize(fj_cpu.view.pb.n_variables, false); - fj_cpu.iter_mtm_vars.reserve(fj_cpu.view.pb.n_variables); + fj_cpu.flip_move_computed.resize(n_variables, false); + fj_cpu.var_bitmap.resize(n_variables, false); + fj_cpu.iter_mtm_vars.reserve(n_variables); recompute_lhs(fj_cpu); @@ -1342,6 +1401,119 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, precompute_problem_features(fj_cpu); } +template +static std::unique_ptr> init_fj_cpu_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::atomic& preemption_flag, + int64_t seed) +{ + using f_t2 = typename type_2::type; + + cuopt_assert(variable_types.size() >= static_cast(problem.num_cols), + "variable type size mismatch"); + + typename mip_solver_settings_t::tolerances_t tolerances{}; + tolerances.absolute_tolerance = settings.primal_tol; + tolerances.relative_tolerance = settings.zero_tol; + tolerances.integrality_tolerance = settings.integer_tol; + tolerances.absolute_mip_gap = settings.absolute_mip_gap_tol; + tolerances.relative_mip_gap = settings.relative_mip_gap_tol; + + const i_t n_variables = problem.num_cols; + const i_t n_constraints = problem.num_rows; + + dual_simplex::csr_matrix_t csr_A(problem.num_rows, problem.num_cols, problem.A.nnz()); + problem.A.to_compressed_row(csr_A); + std::vector coefficients = csr_A.x; + std::vector variables = csr_A.j; + std::vector offsets = csr_A.row_start; + std::vector constraint_lower_bounds = problem.rhs; + std::vector constraint_upper_bounds = problem.rhs; + std::vector variable_bounds(n_variables); + std::vector cpufj_variable_types(n_variables); + std::vector is_binary_variable(n_variables, 0); + i_t n_integer_vars = 0; + + for (i_t j = 0; j < n_variables; ++j) { + variable_bounds[j] = f_t2{problem.lower[j], problem.upper[j]}; + const auto var_type = variable_types[j]; + cpufj_variable_types[j] = + var_type == dual_simplex::variable_type_t::CONTINUOUS ? var_t::CONTINUOUS : var_t::INTEGER; + + const bool is_integer = cpufj_variable_types[j] == var_t::INTEGER; + const bool is_binary = is_integer && + integer_equal(problem.lower[j], f_t{0}, settings.integer_tol) && + integer_equal(problem.upper[j], f_t{1}, settings.integer_tol); + if (is_integer) { ++n_integer_vars; } + if (is_binary) { is_binary_variable[j] = 1; } + } + + const i_t nnz = static_cast(variables.size()); + dual_simplex::csc_matrix_t reverse_csc(n_constraints, n_variables, nnz); + csr_A.to_compressed_col(reverse_csc); + std::vector reverse_coefficients = std::move(reverse_csc.x); + std::vector reverse_constraints = std::move(reverse_csc.i); + std::vector reverse_offsets = std::move(reverse_csc.col_start); + + std::vector projected_seed(n_variables, f_t{0}); + for (i_t j = 0; j < n_variables; ++j) { + f_t value = j < static_cast(seed_assignment.size()) ? seed_assignment[j] : f_t{0}; + value = std::clamp(value, problem.lower[j], problem.upper[j]); + if (variable_types[j] != dual_simplex::variable_type_t::CONTINUOUS) { + value = std::clamp(std::round(value), problem.lower[j], problem.upper[j]); + } + projected_seed[j] = value; + } + + fj_settings_t fj_settings; + fj_settings.mode = fj_mode_t::EXIT_NON_IMPROVING; + fj_settings.n_of_minimums_for_exit = std::numeric_limits::max(); + fj_settings.time_limit = std::numeric_limits::infinity(); + fj_settings.iteration_limit = std::numeric_limits::max(); + fj_settings.update_weights = true; + fj_settings.feasibility_run = false; + fj_settings.seed = seed >= 0 ? seed : cuopt::seed_generator::get_seed(); + + auto fj_cpu = std::make_unique>(preemption_flag); + fj_cpu->view = typename fj_t::climber_data_t::view_t{}; + fj_cpu->pb_ptr = nullptr; + fj_cpu->settings = fj_settings; + + fj_cpu->h_reverse_coefficients = std::move(reverse_coefficients); + fj_cpu->h_reverse_constraints = std::move(reverse_constraints); + fj_cpu->h_reverse_offsets = std::move(reverse_offsets); + fj_cpu->h_coefficients = std::move(coefficients); + fj_cpu->h_offsets = std::move(offsets); + fj_cpu->h_variables = std::move(variables); + fj_cpu->h_obj_coeffs = problem.objective; + fj_cpu->h_var_bounds = std::move(variable_bounds); + fj_cpu->h_cstr_lb = std::move(constraint_lower_bounds); + fj_cpu->h_cstr_ub = std::move(constraint_upper_bounds); + fj_cpu->h_var_types = std::move(cpufj_variable_types); + fj_cpu->h_is_binary_variable = std::move(is_binary_variable); + + fj_cpu->h_cstr_left_weights.resize(n_constraints, 1.0); + fj_cpu->h_cstr_right_weights.resize(n_constraints, 1.0); + fj_cpu->max_weight = 1.0; + fj_cpu->h_objective_weight = 0.0; + fj_cpu->h_assignment = projected_seed; + fj_cpu->h_best_assignment = std::move(projected_seed); + fj_cpu->h_lhs.resize(n_constraints); + fj_cpu->h_lhs_sumcomp.resize(n_constraints, 0); + fj_cpu->h_tabu_nodec_until.resize(n_variables, 0); + fj_cpu->h_tabu_noinc_until.resize(n_variables, 0); + fj_cpu->h_tabu_lastdec.resize(n_variables, 0); + fj_cpu->h_tabu_lastinc.resize(n_variables, 0); + fj_cpu->iterations = 0; + + finalize_fj_cpu_host_initialization( + *fj_cpu, n_variables, n_constraints, n_integer_vars, nnz, tolerances); + return fj_cpu; +} + template static void sanity_checks(fj_cpu_climber_t& fj_cpu) { @@ -1417,7 +1589,7 @@ std::unique_ptr> fj_t::create_cpu_climber( } template -void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit) +void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit, double work_unit_limit) { i_t local_mins = 0; auto loop_start = std::chrono::high_resolution_clock::now(); @@ -1518,7 +1690,7 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit) fj_cpu->total_violations += fj_cpu->view.excess_score(cstr_idx, fj_cpu->h_lhs[cstr_idx]); } if (fj_cpu->iterations % fj_cpu->log_interval == 0) { - CUOPT_LOG_TRACE( + CUOPT_LOG_DEBUG( "%sCPUFJ iteration: %d/%d, local mins: %d, best_objective: %g, viol: %zu, obj weight %g, " "maxw %g", fj_cpu->log_prefix.c_str(), @@ -1527,7 +1699,7 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit) ? fj_cpu->settings.iteration_limit : -1, local_mins, - fj_cpu->pb_ptr->get_user_obj_from_solver_obj(fj_cpu->h_best_objective), + fj_cpu->h_best_objective, fj_cpu->violated_constraints.size(), fj_cpu->h_objective_weight, fj_cpu->max_weight); @@ -1547,12 +1719,21 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit) #endif if (fj_cpu->iterations % 100 == 0 && fj_cpu->iterations > 0) { - // Collect memory statistics + // Use cumulative byte counts (collect() without flush). Each window's contribution to + // work_units_elapsed therefore grows roughly with the running total of bytes touched, + // i.e. quadratically in iterations rather than linearly. This is intentional: the + // memory_aggregator is calibrated for medium/large MIPs, and a strictly-linear scheme + // forces tiny instances (few KB per iteration) to run for tens of seconds before the + // accumulated bytes cross a 0.5 horizon, causing the deterministic producer_sync to + // stall and B&B to time out on instances that should solve in milliseconds. The + // accumulation is still deterministic across runs of the same problem, which is what + // the producer_sync contract actually requires. auto [loads, stores] = fj_cpu->memory_aggregator.collect(); double biased_work = (loads + stores) * fj_cpu->work_unit_bias / 1e10; fj_cpu->work_units_elapsed += biased_work; if (fj_cpu->producer_sync != nullptr) { fj_cpu->producer_sync->notify_progress(); } + if (fj_cpu->work_units_elapsed >= work_unit_limit) { break; } } cuopt_func_call(sanity_checks(*fj_cpu)); @@ -1593,24 +1774,110 @@ std::unique_ptr> init_fj_cpu_standalone( return fj_cpu; } +template +void fj_cpu_task_t::fj_cpu_deleter_t::operator()(fj_cpu_climber_t* ptr) const +{ + delete ptr; +} + +template +std::unique_ptr> make_fj_cpu_task_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::function&, double)> improvement_callback, + std::string log_prefix, + int64_t seed) +{ + auto task = std::make_unique>(); + auto fj_cpu = init_fj_cpu_from_host_lp( + problem, variable_types, seed_assignment, settings, task->preemption_flag, seed); + fj_cpu->log_prefix = std::move(log_prefix); + fj_cpu->improvement_callback = std::move(improvement_callback); + task->fj_cpu.reset(fj_cpu.release()); + return task; +} + +template +void run_fj_cpu_task(fj_cpu_task_t& task, f_t time_limit, double work_unit_limit) +{ + cuopt_assert(task.fj_cpu != nullptr, "CPUFJ task has no climber"); + cpufj_solve(task.fj_cpu.get(), time_limit, work_unit_limit); +} + +template +void stop_fj_cpu_task(fj_cpu_task_t& task) +{ + if (task.fj_cpu) { + auto& fj_cpu = *task.fj_cpu; + fj_cpu.preemption_flag = true; + fj_cpu.halted = true; + } +} + #if MIP_INSTANTIATE_FLOAT template class fj_t; -template void cpufj_solve(fj_cpu_climber_t* fj_cpu, float in_time_limit); +template struct fj_cpu_task_t; +template void cpufj_solve(fj_cpu_climber_t* fj_cpu, + float in_time_limit, + double work_unit_limit); template std::unique_ptr> init_fj_cpu_standalone( problem_t& problem, solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); +template std::unique_ptr> make_fj_cpu_task_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::function&, double)> improvement_callback, + std::string log_prefix, + int64_t seed); +template void run_fj_cpu_task(fj_cpu_task_t& task, + float time_limit, + double work_unit_limit); +template void stop_fj_cpu_task(fj_cpu_task_t& task); +template void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + int n_variables, + int n_constraints, + int n_integer_vars, + int nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances); #endif #if MIP_INSTANTIATE_DOUBLE template class fj_t; -template void cpufj_solve(fj_cpu_climber_t* fj_cpu, double in_time_limit); +template struct fj_cpu_task_t; +template void cpufj_solve(fj_cpu_climber_t* fj_cpu, + double in_time_limit, + double work_unit_limit); template std::unique_ptr> init_fj_cpu_standalone( problem_t& problem, solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); +template std::unique_ptr> make_fj_cpu_task_from_host_lp( + const dual_simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const dual_simplex::simplex_solver_settings_t& settings, + std::function&, double)> improvement_callback, + std::string log_prefix, + int64_t seed); +template void run_fj_cpu_task(fj_cpu_task_t& task, + double time_limit, + double work_unit_limit); +template void stop_fj_cpu_task(fj_cpu_task_t& task); +template void finalize_fj_cpu_host_initialization( + fj_cpu_climber_t& fj_cpu, + int n_variables, + int n_constraints, + int n_integer_vars, + int nnz, + const typename mip_solver_settings_t::tolerances_t& tolerances); #endif } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index 76bf158f9e..cdf3a2f58a 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -191,7 +192,8 @@ struct fj_cpu_climber_t { template void cpufj_solve(fj_cpu_climber_t* fj_cpu, - f_t in_time_limit = std::numeric_limits::infinity()); + f_t in_time_limit = std::numeric_limits::infinity(), + double work_unit_limit = std::numeric_limits::infinity()); // Standalone CPUFJ init for running without full fj_t infrastructure (avoids GPU allocations). // Used for early CPUFJ during presolve.