diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 67412162c2..58a42d079c 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -252,14 +252,14 @@ endif () FetchContent_Declare( papilo - GIT_REPOSITORY "https://github.com/scipopt/papilo.git" + GIT_REPOSITORY "https://github.com/nguidotti/papilo.git" # We would want to get the main branch. However, the main branch # does not have some of the presolvers and settings that we need # Mainly, probing and clique merging. # This is the reason we are using the development branch # from Oct 12, 2025. Once these changes are merged into the main branch, #we can switch to the main branch. - GIT_TAG "741a2b9c8155b249d6df574d758b4d97d4417520" + GIT_TAG "0cfea20c5655249174dbd04f0fb7bd6b1f6e9a0c" GIT_PROGRESS TRUE EXCLUDE_FROM_ALL SYSTEM @@ -293,13 +293,13 @@ set(BUILD_SHARED_LIBS ${BUILD_SHARED_LIBS_SAVED}) # dejavu - header-only graph automorphism library for MIP symmetry detection # https://github.com/markusa4/dejavu (header-only, skip its CMakeLists.txt) FetchContent_Declare( - dejavu - GIT_REPOSITORY "https://github.com/markusa4/dejavu.git" - GIT_TAG "v2.1" - GIT_PROGRESS TRUE - EXCLUDE_FROM_ALL - SYSTEM - SOURCE_SUBDIR _nonexistent + dejavu + GIT_REPOSITORY "https://github.com/markusa4/dejavu.git" + GIT_TAG "v2.1" + GIT_PROGRESS TRUE + EXCLUDE_FROM_ALL + SYSTEM + SOURCE_SUBDIR _nonexistent ) FetchContent_MakeAvailable(dejavu) message(STATUS "dejavu (graph automorphism): ${dejavu_SOURCE_DIR}") @@ -369,7 +369,7 @@ if (NOT SKIP_GRPC_BUILD) # Proto search paths: manual protos in src/grpc, generated data proto in src/grpc/codegen/generated set(PROTO_PATH_MANUAL "${CMAKE_CURRENT_SOURCE_DIR}/src/grpc") - set(PROTO_PATH_GEN "${CMAKE_CURRENT_SOURCE_DIR}/src/grpc/codegen/generated") + set(PROTO_PATH_GEN "${CMAKE_CURRENT_SOURCE_DIR}/src/grpc/codegen/generated") # Generate C++ code from cuopt_remote_data.proto (auto-generated data definitions) set(DATA_PROTO_FILE "${PROTO_PATH_GEN}/cuopt_remote_data.proto") @@ -441,8 +441,8 @@ add_subdirectory(src) # nvcc 13.0.3 ICE (signal 11) compiling sliding_window.cu with 7 GPU architectures; # --split-compile breaks the codegen into per-arch sub-jobs to avoid the crash set_source_files_properties( - ${CMAKE_CURRENT_SOURCE_DIR}/src/routing/local_search/sliding_window.cu - PROPERTIES COMPILE_OPTIONS "--split-compile=0") + ${CMAKE_CURRENT_SOURCE_DIR}/src/routing/local_search/sliding_window.cu + PROPERTIES COMPILE_OPTIONS "--split-compile=0") if (HOST_LINEINFO) set_source_files_properties(${CUOPT_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g1") @@ -495,8 +495,8 @@ set_target_properties(cuopt ) target_compile_definitions(cuopt - PUBLIC "CUOPT_LOG_ACTIVE_LEVEL=RAPIDS_LOGGER_LOG_LEVEL_${LIBCUOPT_LOGGING_LEVEL}" - PUBLIC CUSPARSE_ENABLE_EXPERIMENTAL_API + PUBLIC "CUOPT_LOG_ACTIVE_LEVEL=RAPIDS_LOGGER_LOG_LEVEL_${LIBCUOPT_LOGGING_LEVEL}" + PUBLIC CUSPARSE_ENABLE_EXPERIMENTAL_API ) target_compile_options(cuopt @@ -528,8 +528,8 @@ target_include_directories(cuopt PRIVATE ) target_include_directories(cuopt SYSTEM PRIVATE - "${pslp_SOURCE_DIR}/include" - "${dejavu_SOURCE_DIR}" + "${pslp_SOURCE_DIR}/include" + "${dejavu_SOURCE_DIR}" ) target_include_directories(cuopt diff --git a/cpp/src/branch_and_bound/CMakeLists.txt b/cpp/src/branch_and_bound/CMakeLists.txt index 1e40c1bbf1..47ca6b98f5 100644 --- a/cpp/src/branch_and_bound/CMakeLists.txt +++ b/cpp/src/branch_and_bound/CMakeLists.txt @@ -4,11 +4,11 @@ # cmake-format: on set(BRANCH_AND_BOUND_SRC_FILES - ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp - ) + ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp +) set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} - ${BRANCH_AND_BOUND_SRC_FILES} PARENT_SCOPE) + ${BRANCH_AND_BOUND_SRC_FILES} PARENT_SCOPE) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e4ce4dfd7b..094913baee 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -14,7 +14,7 @@ #include // benchmark_info_t #include -#include +#include #include #include @@ -41,6 +41,7 @@ #include #include #include +#include #include #include @@ -191,6 +192,25 @@ dual_status_t convert_lp_status_to_dual_status(lp_status_t status) } } +// When `log_diving_type` is true, each diving strategy gets its own letter; +// otherwise every dive collapses to 'D'. +inline char feasible_solution_symbol(worker_type_t type, bool log_diving_type) +{ + if (type == HEURISTICS) return 'H'; + if (type == BEST_FIRST) return 'B'; + if (type == SUBMIP) return 'S'; + if (!log_diving_type) { return 'D'; } + switch (type) { + case COEFFICIENT_DIVING: return 'C'; + case LINE_SEARCH_DIVING: return 'L'; + case PSEUDOCOST_DIVING: return 'P'; + case GUIDED_DIVING: return 'G'; + case FARKAS_DIVING: return 'F'; + case VECTOR_LENGTH_DIVING: return 'V'; + default: return 'U'; + } +} + template f_t sgn(f_t x) { @@ -313,6 +333,15 @@ void branch_and_bound_t::set_initial_upper_bound(f_t bound) upper_bound_ = bound; } +template +void branch_and_bound_t::warm_start(const pseudo_costs_t& parent_pc, + const std::vector& reduced_to_original) +{ + pc_.resize(original_lp_.num_cols); + pc_.warm_start_from(parent_pc, reduced_to_original); + root_warm_start_ = true; +} + template void branch_and_bound_t::print_table_header() { @@ -332,7 +361,7 @@ void branch_and_bound_t::print_table_header() } template -void branch_and_bound_t::report_heuristic(f_t obj) +void branch_and_bound_t::report_heuristic(f_t obj, worker_type_t type) { if (is_running_) { f_t lower_bound = get_lower_bound(); @@ -342,7 +371,8 @@ void branch_and_bound_t::report_heuristic(f_t obj) std::string user_gap_text = to_percentage(user_gap); std::string log_line = - std::format("H {:>12} {:>12} {:^+19.6e} {:^+15.6e} {:>8} {:>7} {:^11} {:^11}", + std::format("{} {:>12} {:>12} {:^+19.6e} {:^+15.6e} {:>8} {:>7} {:^11} {:^11}", + feasible_solution_symbol(type, false), "", // nodes explored "", // nodes unexplored user_obj, @@ -383,7 +413,7 @@ void branch_and_bound_t::report( const i_t nodes_unexplored = exploration_stats_.nodes_unexplored; const f_t user_obj = compute_user_objective(original_lp_, obj); const f_t user_lower = compute_user_objective(original_lp_, lower_bound); - const f_t iters = static_cast(exploration_stats_.total_lp_iters); + const f_t iters = static_cast(exploration_stats_.total_simplex_iters); const f_t iter_node = nodes_explored > 0 ? iters / nodes_explored : iters; f_t user_gap = user_relative_gap(user_obj, user_lower); std::string user_gap_text = to_percentage(user_gap); @@ -473,7 +503,8 @@ void branch_and_bound_t::update_user_bound(f_t lower_bound) } template -void branch_and_bound_t::set_solution_from_heuristics(const std::vector& solution) +bool branch_and_bound_t::set_solution_from_heuristics(const std::vector& solution, + worker_type_t heuristic_type) { mutex_original_lp_.lock(); if (solution.size() != original_problem_.num_cols) { @@ -484,9 +515,15 @@ void branch_and_bound_t::set_solution_from_heuristics(const std::vecto crush_primal_solution( original_problem_, original_lp_, solution, new_slacks_, crushed_solution); f_t obj = compute_objective(original_lp_, crushed_solution); + mutex_original_lp_.unlock(); bool is_feasible = false; bool attempt_repair = false; + bool success = false; + + settings_.log.debug_format("{} found solution with obj={:.4g}", + feasible_solution_symbol(heuristic_type, false), + compute_user_objective(original_lp_, obj)); if (improves_incumbent(obj)) { f_t primal_err; @@ -507,7 +544,11 @@ void branch_and_bound_t::set_solution_from_heuristics(const std::vecto f_t current_upper_bound = upper_bound_.load(); upper_bound_ = std::min(current_upper_bound, obj); incumbent_.set_incumbent_solution(obj, crushed_solution); - if (current_upper_bound > upper_bound_.load()) { report_heuristic(obj); } + if (current_upper_bound > upper_bound_.load()) { + report_heuristic(obj, heuristic_type); + success = true; + } + } else { attempt_repair = true; constexpr bool verbose = false; @@ -530,6 +571,8 @@ void branch_and_bound_t::set_solution_from_heuristics(const std::vecto repair_queue_.push_back(solution); mutex_repair_.unlock(); } + + return success; } template @@ -575,6 +618,28 @@ void branch_and_bound_t::queue_external_solution_deterministic( mutex_heuristic_queue_.unlock(); } +template +void branch_and_bound_t::set_solution_from_cpu_fj(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_format("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 { + if (settings_.solution_callback != nullptr) { + settings_.solution_callback(user_assignment, obj); + } + set_solution_from_heuristics(user_assignment, HEURISTICS); + } +} + template bool branch_and_bound_t::repair_solution(const std::vector& edge_norms, const std::vector& potential_solution, @@ -668,7 +733,7 @@ void branch_and_bound_t::repair_heuristic_solutions() if (improves_incumbent(repaired_obj)) { upper_bound_ = std::min(upper_bound_.load(), repaired_obj); incumbent_.set_incumbent_solution(repaired_obj, repaired_solution); - report_heuristic(repaired_obj); + report_heuristic(repaired_obj, HEURISTICS); if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -716,6 +781,8 @@ template void branch_and_bound_t::set_final_solution(mip_solution_t& solution, f_t lower_bound) { + if (solver_status_ == mip_status_t::HALT) { settings_.log.printf("Stopping submip solve...\n"); } + if (solver_status_ == mip_status_t::NUMERICAL) { settings_.log.printf("Numerical issue encountered. Stopping the solver...\n"); } @@ -732,6 +799,10 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& settings_.log.printf("Node limit reached. Stopping the solver...\n"); } + if (solver_status_ == mip_status_t::ITERATION_LIMIT) { + settings_.log.printf("Simplex iteration limit reached. Stopping the solver...\n"); + } + if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } @@ -744,7 +815,7 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& settings_.log.print_format("Explored {} nodes ({} simplex iterations) in {:.2f}s.", exploration_stats_.nodes_explored.load(), - exploration_stats_.total_lp_iters.load(), + exploration_stats_.total_simplex_iters.load(), toc(exploration_stats_.start_time)); if (exploration_stats_.orbital_fixing_nodes.load() > 0 || @@ -767,7 +838,7 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { solver_status_ = mip_status_t::OPTIMAL; #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - if (settings_.sub_mip == 0 && has_solver_space_incumbent()) { + if (settings_.inside_submip == 0 && has_solver_space_incumbent()) { write_solution_for_cut_verification(original_lp_, incumbent_.x); } #endif @@ -797,18 +868,19 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& if (has_solver_space_incumbent()) { uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); - solution.objective = incumbent_.objective; + solution.objective = incumbent_.objective; + solution.has_incumbent = true; } solution.lower_bound = lower_bound; solution.nodes_explored = exploration_stats_.nodes_explored; - solution.simplex_iterations = exploration_stats_.total_lp_iters; + solution.simplex_iterations = exploration_stats_.total_simplex_iters; } template void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, - search_strategy_t thread_type) + worker_type_t thread_type) { bool send_solution = false; @@ -915,10 +987,20 @@ branch_variable_t branch_and_bound_t::variable_selection( case VECTOR_LENGTH_DIVING: return vector_length_diving(worker->leaf_problem, fractional, solution, log); - default: + case SUBMIP: // 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}; + + // These should never be called for selecting which variable to branch + case HEURISTICS: + case NUM_WORKER_TYPES: log.debug("Unknown variable selection method: %d\n", worker->search_strategy); return {-1, branch_direction_t::NONE}; } + + log.debug("Unknown variable selection method: %d\n", worker->search_strategy); + return {-1, branch_direction_t::NONE}; } // ============================================================================ @@ -947,7 +1029,6 @@ struct tree_update_policy_t { branch_direction_t dir) = 0; virtual void on_numerical_issue(mip_node_t*) = 0; virtual void graphviz(search_tree_t&, mip_node_t*, const char*, f_t) = 0; - virtual void on_optimal_callback(const std::vector&, f_t) = 0; }; template @@ -988,14 +1069,14 @@ struct nondeterministic_policy_t : tree_update_policy_t { const std::vector& fractional, const std::vector& x) override { - if (worker->search_strategy == search_strategy_t::BEST_FIRST) { + if (worker->search_strategy == worker_type_t::BEST_FIRST) { node->objective_estimate = bnb.pc_.obj_estimate(fractional, x, node->lower_bound); } } void on_numerical_issue(mip_node_t* node) override { - if (worker->search_strategy == search_strategy_t::BEST_FIRST) { + if (worker->search_strategy == worker_type_t::BEST_FIRST) { fetch_min(bnb.lower_bound_numerical_, node->lower_bound); log.printf("LP returned numerical issue on node %d. Best bound set to %+10.6e.\n", node->node_id, @@ -1011,16 +1092,6 @@ struct nondeterministic_policy_t : tree_update_policy_t { tree.graphviz_node(log, node, label, value); } - void on_optimal_callback(const std::vector& x, f_t objective) override - { - if (worker->search_strategy == search_strategy_t::BEST_FIRST && - bnb.settings_.node_processed_callback != nullptr) { - std::vector original_x; - uncrush_primal_solution(bnb.original_problem_, bnb.original_lp_, x, original_x); - bnb.settings_.node_processed_callback(original_x, objective); - } - } - void on_node_completed(mip_node_t*, node_status_t, branch_direction_t) override {} }; @@ -1051,7 +1122,6 @@ struct deterministic_policy_base_t : tree_update_policy_t { void on_numerical_issue(mip_node_t*) override {} void graphviz(search_tree_t&, mip_node_t*, const char*, f_t) override {} - void on_optimal_callback(const std::vector&, f_t) override {} }; template @@ -1156,14 +1226,14 @@ struct deterministic_diving_policy_t log.log = false; switch (this->worker.diving_type) { - case search_strategy_t::PSEUDOCOST_DIVING: + case worker_type_t::PSEUDOCOST_DIVING: return pseudocost_diving( this->worker.pc_snapshot, fractional, x, *this->worker.root_solution, log); - case search_strategy_t::LINE_SEARCH_DIVING: + case worker_type_t::LINE_SEARCH_DIVING: return line_search_diving(fractional, x, *this->worker.root_solution, log); - case search_strategy_t::GUIDED_DIVING: + case worker_type_t::GUIDED_DIVING: if (this->worker.incumbent_snapshot.empty()) { return pseudocost_diving( this->worker.pc_snapshot, fractional, x, *this->worker.root_solution, log); @@ -1172,7 +1242,7 @@ struct deterministic_diving_policy_t this->worker.pc_snapshot, fractional, x, this->worker.incumbent_snapshot, log); } - case search_strategy_t::COEFFICIENT_DIVING: { + case worker_type_t::COEFFICIENT_DIVING: { return coefficient_diving(this->worker.leaf_problem, fractional, x, @@ -1220,6 +1290,21 @@ struct deterministic_diving_policy_t } }; +template +void branch_and_bound_t::apply_objective_step(mip_node_t* node_ptr, + f_t leaf_obj) +{ + if (original_lp_.objective_step.has_step()) { + f_t step = original_lp_.objective_step.step_size; + f_t bias = original_lp_.objective_step.bias; + // Round up to next value on the lattice: k * step + bias >= leaf_obj + f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); + node_ptr->lower_bound = k * step + bias; + } else if (original_lp_.objective_is_integral) { + node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); + } +} + template template std::pair branch_and_bound_t::update_tree_impl( @@ -1277,22 +1362,7 @@ std::pair branch_and_bound_t::updat policy.graphviz(search_tree, node_ptr, "lower bound", leaf_obj); policy.update_pseudo_costs(node_ptr, leaf_obj); node_ptr->lower_bound = leaf_obj; - // If the objective is integral or must move in steps than - // the lower bound will be different from the leaf objective. - // We use the leaf objective for RINS (on_optimal_callback) - // and if we are integer feasible (handle_integer_solution). - // We use the lower bound to decide if we should fathom the - // node or branch. - if (original_lp_.objective_step.has_step()) { - f_t step = original_lp_.objective_step.step_size; - f_t bias = original_lp_.objective_step.bias; - // Round up to next value on the lattice: k * step + bias >= leaf_obj - f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); - node_ptr->lower_bound = k * step + bias; - } else if (original_lp_.objective_is_integral) { - node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); - } - policy.on_optimal_callback(leaf_solution.x, leaf_obj); + apply_objective_step(node_ptr, leaf_obj); if (num_frac == 0) { policy.handle_integer_solution(node_ptr, leaf_obj, leaf_solution.x); @@ -1469,11 +1539,11 @@ dual_status_t branch_and_bound_t::solve_node_lp( lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); lp_settings.scale_columns = false; - if (worker->search_strategy != search_strategy_t::BEST_FIRST) { - int64_t bnb_lp_iters = exploration_stats_.total_lp_iters; + if (worker->search_strategy != BEST_FIRST && worker->search_strategy != SUBMIP) { + int64_t bnb_lp_iters = exploration_stats_.total_simplex_iters; f_t factor = settings_.diving_settings.iteration_limit_factor; int64_t max_iter = factor * bnb_lp_iters; - lp_settings.iteration_limit = max_iter - stats.total_lp_iters; + lp_settings.iteration_limit = max_iter - stats.total_simplex_iters; if (lp_settings.iteration_limit <= 0) { return dual_status_t::ITERATION_LIMIT; } } @@ -1544,7 +1614,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( } stats.total_lp_solve_time += toc(lp_start_time); - stats.total_lp_iters += node_iter; + stats.total_simplex_iters += node_iter; } } @@ -1578,11 +1648,13 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, f_t rel_gap = user_relative_gap(user_obj, user_lower); f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); + bool can_launch_rins = true; + while (stack.size() > 0 && (solver_status_ == mip_status_t::UNSET && is_running_) && rel_gap > settings_.relative_mip_gap_tol && abs_gap > settings_.absolute_mip_gap_tol) { if (worker->worker_id == 0) { repair_heuristic_solutions(); } - if (worker->total_active_diving_workers < worker->total_max_diving_workers && + if (worker->active_diving_workers < worker->max_diving_workers && worker->node_queue.diving_queue_size() > 0) { launch_diving_worker(worker); } @@ -1618,10 +1690,10 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, if (worker->worker_id == 0) { f_t time_since_last_log = - exploration_stats_.last_log == 0 ? 1.0 : toc(exploration_stats_.last_log); + exploration_stats_.last_log == 0 ? 50.0 : toc(exploration_stats_.last_log); i_t nodes_since_last_log = exploration_stats_.nodes_since_last_log; - if (((nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + if (((nodes_since_last_log >= 100 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && time_since_last_log >= 1) || (time_since_last_log > 30) || now > settings_.time_limit) { report(' ', upper_bound_, lower_bound, node_ptr->depth, node_ptr->integer_infeasible); @@ -1645,6 +1717,13 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, break; } + if (exploration_stats_.total_simplex_iters > settings_.bnb_iteration_limit) { + solver_status_ = mip_status_t::ITERATION_LIMIT; + stack.push_front(node_ptr); + --exploration_stats_.nodes_being_solved; + break; + } + dual_status_t lp_status = solve_node_lp(node_ptr, worker, exploration_stats_, settings_.log); ++exploration_stats_.nodes_since_last_log; ++exploration_stats_.nodes_explored; @@ -1674,6 +1753,10 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, worker->recompute_bounds = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::HAS_CHILDREN) { + if (can_launch_rins) { + can_launch_rins = !launch_local_branching_worker(worker->leaf_solution.x); + } + // The stack should only contain the children of the current parent. // If the stack size is greater than 0, // we pop the current node from the stack and place it in the global heap, @@ -1803,19 +1886,31 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t if (obj_dyn < diving_settings.farkas_obj_dynamism_tol) { diving_settings.farkas_diving = 0; } } - worker->calculate_num_diving_workers( - bfs_worker_pool_.size(), diving_worker_pool_.size(), diving_settings); + worker->calculate_max_diving_workers(bfs_worker_pool_.size(), diving_worker_pool_.size()); + worker->update_diving_heuristic_list(diving_settings); while (solver_status_ == mip_status_t::UNSET && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && node_queue.best_first_queue_size() > 0) { + if (halt_callback_) { + bool stop = halt_callback_(user_obj, user_lower); + if (stop) { + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::HALT; + settings_.log.print_format( + "Received halt signal. Current best obj={:.6e} and best bound={:.6e}\n", + user_obj, + user_lower); + break; + } + } + // If the guided diving was disabled previously due to the lack of an incumbent solution, // re-enable as soon as a new incumbent is found. if (diving_worker_pool_.size() > 0 && settings_.diving_settings.guided_diving != 0 && diving_settings.guided_diving == 0) { if (has_solver_space_incumbent()) { diving_settings.guided_diving = 1; - worker->calculate_num_diving_workers( - bfs_worker_pool_.size(), diving_worker_pool_.size(), diving_settings); + worker->update_diving_heuristic_list(diving_settings); } } @@ -1865,37 +1960,36 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t if (worker->node_queue.best_first_queue_size() == 0) { bfs_worker_pool_.return_worker_to_pool(worker); } + + // We explore the entire tree and none worker is running. Set is_running_ to false to stop + // the submip. + if (exploration_stats_.nodes_unexplored == 0 && + bfs_worker_pool_.num_idle() == bfs_worker_pool_.size()) { + is_running_ = false; + } } template -void branch_and_bound_t::dive_with(diving_worker_t* worker) +void branch_and_bound_t::dive_with(diving_worker_t* worker, i_t backtrack_limit) { raft::common::nvtx::range scope("BB::diving_thread"); if (worker->orbital_fixing) { worker->orbital_fixing->disable(); } logger_t log; log.log = false; - search_strategy_t search_strategy = worker->search_strategy; - const i_t diving_node_limit = settings_.diving_settings.node_limit; - const i_t diving_backtrack_limit = settings_.diving_settings.backtrack_limit; - - worker->recompute_basis = true; - worker->recompute_bounds = true; + const i_t diving_node_limit = settings_.diving_settings.node_limit; + worker->recompute_basis = true; + worker->recompute_bounds = true; search_tree_t dive_tree(std::move(worker->start_node)); // Since we are perform a DFS with a limit amount of backtracking, the - // stack can hold at most `diving_backtrack_limit` + 2 siblings nodes of the + // stack can hold at most `backtrack_limit` + 2 siblings nodes of the // current level - circular_deque_t*> stack(diving_backtrack_limit + 4); + circular_deque_t*> stack(backtrack_limit + 4); stack.push_front(&dive_tree.root); branch_and_bound_stats_t dive_stats; - dive_stats.total_lp_iters = 0; - dive_stats.total_lp_solve_time = 0; - dive_stats.nodes_explored = 0; - dive_stats.nodes_unexplored = 1; - f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; f_t user_obj = compute_user_objective(original_lp_, upper_bound); @@ -1949,8 +2043,7 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker) // Remove nodes that we can no longer backtrack to (i.e., from the current node, we can only // backtrack to a node that is has a depth of at most 5 levels lower than the current node). - while (stack.size() > 1 && - stack.front()->depth - stack.back()->depth > diving_backtrack_limit) { + while (stack.size() > 1 && stack.front()->depth - stack.back()->depth > backtrack_limit) { stack.pop_back(); } @@ -1962,12 +2055,16 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker) abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); } - diving_worker_pool_.return_worker_to_pool(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 != SUBMIP) { diving_worker_pool_.return_worker_to_pool(worker); } } template bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* bfs_worker) { + if (!bfs_worker->is_diving_enabled()) return false; + // Get an idle worker. diving_worker_t* diving_worker = diving_worker_pool_.pop_idle_worker(); if (diving_worker == nullptr) { return false; } @@ -2000,31 +2097,558 @@ bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* return false; } - for (int i = 1; i < num_search_strategies; ++i) { - auto strategy = search_strategies[i]; - - if (bfs_worker->active_diving_workers[strategy] < bfs_worker->max_diving_workers[strategy]) { - diving_worker->search_strategy = strategy; - diving_worker->bfs_worker = bfs_worker; - diving_worker->set_active(); - bfs_worker->active_diving_workers[strategy]++; - bfs_worker->total_active_diving_workers++; + auto strategy = bfs_worker->next_diving_heuristic(); + diving_worker->search_strategy = strategy; + diving_worker->bfs_worker = bfs_worker; + diving_worker->set_active(); + ++bfs_worker->active_diving_workers; - assert(bfs_worker->active_diving_workers[strategy].load() <= - bfs_worker->max_diving_workers[strategy]); - assert(bfs_worker->total_active_diving_workers.load() <= - bfs_worker->total_max_diving_workers); + assert(bfs_worker->active_diving_workers.load() <= bfs_worker->max_diving_workers); #pragma omp task affinity(*diving_worker) priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) \ firstprivate(diving_worker) - dive_with(diving_worker); + dive_with(diving_worker, settings_.diving_settings.backtrack_limit); + + return true; +} + +template +bool branch_and_bound_t::launch_local_branching_worker(const std::vector& sol) +{ + if (settings_.submip_settings.enable_rins == 0) return false; + if (!incumbent_.has_incumbent) return false; + if (submip_worker_pool_.num_idle() == 0) return false; + + diving_worker_t* worker = submip_worker_pool_.pop_idle_worker(); + if (!worker) return false; + + worker->set_active(); + worker->search_strategy = SUBMIP; + +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) affinity(worker) \ + firstprivate(worker, sol) if (!settings_.inside_submip) + rins(worker, sol); + + return true; +} + +template +void branch_and_bound_t::solve_submip(diving_worker_t* worker, + const std::vector& current_incumbent, + i_t num_var_fixed, + i_t num_integers, + i_t submip_level, + std::string_view log_prefix) +{ + double start_time = tic(); + + std::vector& lower = worker->leaf_problem.lower; + std::vector& upper = worker->leaf_problem.upper; + std::vector& bounds_changed = worker->bounds_changed; + f_t fixrate = (f_t)num_var_fixed / num_integers; + + bool feasible = + worker->node_presolver.bounds_strengthening(settings_, bounds_changed, lower, upper); + + if (!feasible) { + // This should never happen since we are fixing bounds that are already in the incumbent. + submip_stats_.save_infeasible(fixrate); + 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 rel_gap = user_relative_gap(user_obj, user_lower); + i_t explored = exploration_stats_.nodes_explored; + + simplex_solver_settings_t submip_settings = settings_; + submip_settings.print_presolve_stats = false; + submip_settings.num_threads = 1; + submip_settings.reliability_branching = 0; + submip_settings.clique_cuts = 0; + submip_settings.zero_half_cuts = 0; + 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; + +#ifdef DEBUG_SUBMIP + 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), + var_types_); +#else + submip_settings.log.log_prefix = log_prefix; +#endif + + submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; + submip_settings.bnb_iteration_limit = + exploration_stats_.total_simplex_iters * settings_.submip_settings.iteration_limit_ratio; + submip_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + if (submip_settings.time_limit < 0) { return; } + + submip_settings.relative_mip_gap_tol = + std::min(settings_.submip_settings.target_mip_gap, rel_gap); + + submip_settings.submip_settings.enable_rins = settings_.submip_settings.enable_rins != 0 && + submip_level <= settings_.submip_settings.max_level; + + submip_settings.log.print_format( + "Sub-MIP solve settings: time_limit={:.2f}, node_limit={}, iter_limit={} (current_iter={}), " + "tol={:g}", + submip_settings.time_limit, + submip_settings.node_limit, + submip_settings.bnb_iteration_limit, + exploration_stats_.total_simplex_iters.load(), + submip_settings.relative_mip_gap_tol); + + user_problem_t submip_problem(original_problem_.handle_ptr); + simplex::convert_simplex_problem( + worker->leaf_problem, var_types_, settings_, new_slacks_, submip_problem); + + third_party_presolve_t presolver; + f_t presolve_time_limit = std::min(0.1 * submip_settings.time_limit, 60.0); + third_party_presolve_status_t presolver_status = + presolver.apply(submip_problem, submip_settings, presolve_time_limit, 1); + + double presolve_time = toc(start_time); + + if (presolver_status == third_party_presolve_status_t::INFEASIBLE || + presolver_status == third_party_presolve_status_t::UNBNDORINFEAS || + presolver_status == third_party_presolve_status_t::UNBOUNDED) { + submip_stats_.save_infeasible(fixrate); + return; + } + + // Also handle optimal + if (submip_problem.num_rows == 0 || submip_problem.num_cols == 0) { + submip_settings.log.print_format( + "Sub-MIP presolved to a trivial {} x {} problem; solving by bound pushing", + submip_problem.num_rows, + submip_problem.num_cols); + + std::vector reduced_sol(submip_problem.num_cols); + + for (i_t j = 0; j < submip_problem.num_cols; ++j) { + const f_t c = submip_problem.objective[j]; + const f_t l = submip_problem.lower[j]; + const f_t u = submip_problem.upper[j]; + // Minimize c_j x_j over [l, u]; fall back to any finite bound (0 if both are infinite). + if (c < -settings_.zero_tol) { + reduced_sol[j] = std::isfinite(u) ? u : (std::isfinite(l) ? l : 0); + } else { + reduced_sol[j] = std::isfinite(l) ? l : (std::isfinite(u) ? u : 0); + } + } + + std::vector full_sol(original_problem_.num_cols); + presolver.uncrush_primal_solution(reduced_sol, full_sol); + bool success = set_solution_from_heuristics(full_sol, SUBMIP); + if (success) submip_stats_.save_success(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, fixrate, &presolver](std::vector& solution, + f_t obj) { + std::vector full_sol(original_problem_.num_cols); + presolver.uncrush_primal_solution(solution, full_sol); + settings_.log.debug_format("SubMIP found a feasible solution with obj={:.4g}", obj); + bool success = set_solution_from_heuristics(full_sol, SUBMIP); + if (success) submip_stats_.save_success(fixrate); + }; + + submip_settings.log.print_format("Sub-MIP: {} constraints, {} variables, {} nonzeros\n", + 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); + mip_solution_t submip_solution(submip_problem.num_cols); + + // Map the current incumbent to the user space, removing the variables related to the slacks/cuts + // added. + std::vector uncrushed_incumbent; + mutex_original_lp_.lock(); + uncrush_primal_solution(original_problem_, original_lp_, current_incumbent, uncrushed_incumbent); + mutex_original_lp_.unlock(); + + // Crush the incumbent to presolve space. It may not be valid for the sub-MIP since we + // may fix integer variables that does not match the current incumbent to reach the target + // fix rate. + std::vector presolved_incumbent; + presolver.crush_primal_solution(uncrushed_incumbent, presolved_incumbent); + submip_bnb.set_initial_guess(presolved_incumbent); + submip_bnb.set_initial_upper_bound(upper_bound_.load()); + + submip_bnb.warm_start(pc_, presolver.get_reduced_to_original_map()); + + if (halt_callback_) { + // Copy the halt callback to the deeper level. + submip_bnb.set_halt_callback(halt_callback_); + } else { + // This should only be called by the main solver. + submip_bnb.set_halt_callback([this](f_t, f_t submip_lower_bound) { + f_t user_upper = compute_user_objective(this->original_lp_, this->upper_bound_.load()); + bool is_suboptimal = 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_suboptimal || !is_solver_running; + }); + } - return true; + fj_cpu_worker_t submip_fj_cpu_worker; + scope_guard cpufj_guard([&]() { submip_fj_cpu_worker.stop(); }); + + if (settings_.submip_settings.enable_cpufj) { + std::vector initial_guess; + crush_primal_solution(submip_problem, + submip_bnb.original_lp_, + presolved_incumbent, + submip_bnb.new_slacks_, + initial_guess); + + submip_fj_cpu_worker.improvement_callback = + [&submip_bnb](f_t obj, const std::vector& assignment, double work_units) { + submip_bnb.set_solution_from_cpu_fj(obj, assignment, work_units); + }; + + f_t time_limit = submip_settings.time_limit; + f_t work_limit = 1.0; + submip_fj_cpu_worker.from_simplex_lp(submip_bnb.original_lp_, + submip_bnb.var_types_, + initial_guess, + submip_bnb.settings_, + std::format("{} [CPU FJ]", log_prefix), + worker->rng.next_i64()); + submip_fj_cpu_worker.run_async(time_limit, work_limit); + } + + mip_status_t submip_status = submip_bnb.solve(submip_solution); + double submip_time = toc(start_time); + + submip_settings.log.print_format( + "Sub-MIP: status={}, iterations={} (total={}), presolve_time={:.2f}, total_time={:.2f} \n", + mip_status_to_string(submip_status), + submip_solution.simplex_iterations, + exploration_stats_.total_simplex_iters.load(), + presolve_time, + submip_time); + + if (submip_status == mip_status_t::NUMERICAL) { return; } + if (submip_status == mip_status_t::INFEASIBLE || submip_status == mip_status_t::UNBOUNDED) { + submip_stats_.save_infeasible(fixrate); + return; + } + + if (submip_solution.has_incumbent) { + std::vector full_sol(original_problem_.num_cols); + presolver.uncrush_primal_solution(submip_solution.x, full_sol); + bool success = set_solution_from_heuristics(full_sol, SUBMIP); + if (success) submip_stats_.save_success(fixrate); + } + + // Accumulate simplex iterations to determine when to stop exploring the sub-MIP + if (settings_.inside_submip) { + exploration_stats_.total_simplex_iters += submip_solution.simplex_iterations; + } +} + +inline double submip_get_max_fixrate(const submip_stats_t& stats, + const simplex::submip_settings_t& submip_settings, + pcgenerator_t& rng) +{ + // Adaptive fix rate based on previous successes and failures. + double low = submip_settings.base_target_fixrate; + double high = submip_settings.base_target_fixrate; + + if (stats.total_infeasible > 0) { + double infeasible_avg_fixrate = stats.average_infeasible_fixrate(); + high = 0.9 * infeasible_avg_fixrate; + low = std::min(low, high); + } + + if (stats.total_success > 0) { + double success_avg_fixrate = stats.average_success_fixrate(); + low = std::min(low, 0.9 * success_avg_fixrate); + high = std::max(high, 1.1 * success_avg_fixrate); + } + + double fixrate = high > low ? rng.uniform(low, high) : low; + return fixrate; +} + +template +void branch_and_bound_t::rins(diving_worker_t* rins_worker, + const std::vector& node_solution) +{ + raft::common::nvtx::range scope("BB::rins_thread"); + if (rins_worker->orbital_fixing) { rins_worker->orbital_fixing->disable(); } + logger_t log; + log.log = false; + + i_t submip_level = settings_.submip_settings.level + 1; + std::string log_prefix = std::format("[RINS {}] ", submip_level); + + ++submip_stats_.total_calls; + + bool has_submip = false; + const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + + branch_and_bound_stats_t rins_stats; + + mip_node_t node = search_tree_.root.detach_copy(); + node.packed_vstatus = compress_vstatus(root_vstatus_); + 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->skip_set_bounds = 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(); + + std::vector integer_list; + for (i_t j = 0; j < var_types_.size(); ++j) { + if (var_types_[j] == variable_type_t::CONTINUOUS) { continue; } + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + integer_list.push_back(j); + } + i_t num_integers = integer_list.size(); + + f_t max_fixrate = + submip_get_max_fixrate(submip_stats_, settings_.submip_settings, rins_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_) { + // 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; + for (i_t j : fractional) { + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + if (std::abs(node_solution[j] - current_incumbent[j]) <= settings_.integer_tol) { + f_t fixed_val = std::round(node_solution[j]); + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + } + } + + // 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); + has_submip = true; + break; + } + + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; + break; + } + + if (prev_num_fixed == num_var_fixed) { + // RINS neighbourhood 2: Search into the entire list of integer variables where the current + // LP solution matches the current incumbent. + for (i_t j : integer_list) { + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + if (std::abs(current_sol[j] - current_incumbent[j]) <= settings_.integer_tol) { + f_t fixed_val = std::round(current_sol[j]); + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + } + } + + // 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); + has_submip = true; + break; + } + + // Even considering the entire integer list, we were unable to fix a single variable in this + // iteration. Iterate over the fractional variables again and fixing those that closest to + // an integer solution first in order to reach the fixing threshold. + if (prev_num_fixed == num_var_fixed) { + std::vector> candidates; + for (i_t j : fractional) { + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + + f_t root_change = current_sol[j] - root_relax_soln_.x[j]; + f_t obj_coeff = rins_worker->leaf_problem.objective[j]; + f_t fixed_val = 0; + + if (root_change >= 0.4) { + fixed_val = std::ceil(current_sol[j]); + } else if (root_change <= -0.4) { + fixed_val = std::floor(current_sol[j]); + } else if (obj_coeff > settings_.zero_tol) { + fixed_val = std::ceil(current_sol[j]); + } else if (obj_coeff < -settings_.zero_tol) { + fixed_val = std::floor(current_sol[j]); + } else { + fixed_val = std::round(current_sol[j]); + } + + candidates.push_back(std::make_tuple(std::abs(fixed_val - current_sol[j]), j, fixed_val)); + } + + std::sort(candidates.begin(), candidates.end(), [](auto a, auto b) { + return std::get<0>(a) < std::get<0>(b); + }); + + f_t change = 0; + for (auto [dist, j, fixed_val] : candidates) { + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + + // Limit the amount of fixing to the current LP. + change += dist; + if (change >= 0.5) { break; } + } + + 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); + 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); + has_submip = true; + break; + } + } + } + + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; + break; + } + + dual_status_t lp_status = solve_node_lp(&node, rins_worker, rins_stats, log); + if (lp_status != dual_status_t::OPTIMAL) { break; } + + fractional.clear(); + num_frac = fractional_variables(settings_, current_sol, var_types_, fractional); + + f_t leaf_obj = compute_objective(rins_worker->leaf_problem, current_sol); + node.lower_bound = leaf_obj; + + apply_objective_step(&node, leaf_obj); + if (leaf_obj > upper_bound_.load()) { break; } + + if (num_frac == 0) { + // We found a feasible solution when fixing the variables in RINS. + add_feasible_solution(leaf_obj, current_sol, -1, SUBMIP); + break; } } - diving_worker_pool_.return_worker_to_pool(diving_worker); - return false; + f_t fixrate = (f_t)num_var_fixed / num_integers; + + if (has_submip) { + // If no enough variables was fixed (the neighbourhood is too loose) or the sub-MIP already + // found a solution that improved the incumbent, then do a DFS with a backtrack_limit of 5 + // levels up to try to find a feasible solution quickly from the neighbourhood. + if (fixrate < settings_.submip_settings.min_fixrate_cap || + (settings_.inside_submip && submip_stats_.total_success != 0)) { + rins_worker->skip_set_bounds = false; + rins_worker->start_node = std::move(node); + rins_worker->start_lower = lower; + rins_worker->start_upper = upper; + bool is_feasible = rins_worker->presolve_start_bounds(settings_); + if (is_feasible) { + fj_cpu_worker_t submip_fj_cpu_worker; + + if (settings_.submip_settings.enable_cpufj) { + submip_fj_cpu_worker.improvement_callback = + [this](f_t obj, const std::vector& assignment, double work_units) { + this->set_solution_from_cpu_fj(obj, assignment, work_units); + }; + + 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.from_simplex_lp(rins_worker->leaf_problem, + var_types_, + rins_worker->leaf_solution.x, + settings_, + std::format("{} [CPU FJ]", log_prefix), + rins_worker->rng.next_i64()); + submip_fj_cpu_worker.run_sync(time_limit, work_limit); + } + + dive_with(rins_worker, 5); + } + + } else { + solve_submip( + rins_worker, current_incumbent, num_var_fixed, num_integers, submip_level, log_prefix); + } + } + + // Accumulate the iterations for sub-MIP so it stops when it reaches the allocated budget. + if (settings_.inside_submip) { + exploration_stats_.total_simplex_iters += rins_stats.total_simplex_iters; + } + + settings_.log.debug_format( + "{}success={}, infeasible={}, calls={}, fixrate={:.4g} ({}), max_fixrate={:.4g} ({}), " + "min_fixrate={:.4g} ({})\n", + log_prefix, + submip_stats_.total_success.load(), + submip_stats_.total_infeasible.load(), + submip_stats_.total_calls.load(), + fixrate, + num_var_fixed, + max_fixrate, + max_var_fixed, + min_fixrate, + min_var_fixed); + + submip_worker_pool_.return_worker_to_pool(rins_worker); } template @@ -2099,6 +2723,12 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( #pragma omp taskwait depend(in : root_status) set_root_concurrent_halt(0); // Clear the concurrent halt flag + + // Since Barrier/PDLP iterations are not comparable with the simplex iterations + // used in the remaining of the B&B, use the iterations of dual simplex before it + // being stopped as an approximation. + exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; + // Override the root relaxation solution with the crossover solution root_relax_soln = root_crossover_soln_; root_vstatus = crossover_vstatus_; @@ -2147,12 +2777,14 @@ lp_status_t branch_and_bound_t::solve_root_relaxation( } else { // Wait for the dual simplex to finish (after telling PDLP/Barrier to stop) #pragma omp taskwait depend(in : root_status) - root_relax_solved_by = DualSimplex; + root_relax_solved_by = DualSimplex; + exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; } } else { // Wait for the dual simplex to finish (crossover do not produced a solution) #pragma omp taskwait depend(in : root_status) - root_relax_solved_by = DualSimplex; + root_relax_solved_by = DualSimplex; + exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; } is_root_solution_set = true; @@ -2352,7 +2984,7 @@ auto branch_and_bound_t::do_cut_pass( root_relax_soln_, iter, edge_norms_); - exploration_stats_.total_lp_iters += iter; + exploration_stats_.total_simplex_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); @@ -2378,7 +3010,7 @@ auto branch_and_bound_t::do_cut_pass( 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; + exploration_stats_.total_simplex_iters += root_relax_soln_.iterations; root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); } else { settings_.log.printf("Cut status %s\n", simplex::dual_status_to_string(cut_status).c_str()); @@ -2491,6 +3123,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut incumbent_.set_incumbent_solution(computed_obj, crushed_guess); upper_bound_ = computed_obj; mutex_upper_.unlock(); + + settings_.log.print_format("Warm starting B&B with an initial guess. Obj={:.4g}", + compute_user_objective(original_lp_, computed_obj)); } } @@ -2530,14 +3165,16 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut std::vector basic_list(original_lp_.num_rows); std::vector nonbasic_list; basis_update_mpf_t basis_update(original_lp_.num_rows, settings_.refactor_frequency); - lp_status_t root_status; + lp_status_t root_status = lp_status_t::UNSET; solving_root_relaxation_ = true; f_t root_relax_start_time = tic(); + if (!enable_concurrent_lp_root_solve()) { // RINS/SUBMIP path - settings_.log.printf("\nSolving LP root relaxation with dual simplex\n"); - root_status = solve_linear_program_with_advanced_basis(original_lp_, + settings_.log.printf("\n"); + settings_.log.printf("Solving LP root relaxation with dual simplex\n"); + root_status = solve_linear_program_with_advanced_basis(original_lp_, exploration_stats_.start_time, lp_settings, root_relax_soln_, @@ -2546,9 +3183,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nonbasic_list, root_vstatus_, edge_norms_); + root_relax_solved_by = DualSimplex; + exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; } else { - settings_.log.printf("\nSolving LP root relaxation in concurrent mode\n"); + settings_.log.printf("\n"); + settings_.log.printf("Solving LP root relaxation in concurrent mode\n"); root_status = solve_root_relaxation(lp_settings, root_relax_soln_, root_vstatus_, @@ -2557,8 +3197,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nonbasic_list, edge_norms_); } + solving_root_relaxation_ = false; - exploration_stats_.total_lp_iters = root_relax_soln_.iterations; f_t root_relax_elapsed_time = toc(root_relax_start_time); exploration_stats_.total_lp_solve_time = root_relax_elapsed_time; @@ -2609,11 +3249,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } assert(root_status == lp_status_t::OPTIMAL); - settings_.log.print_format( - "\nRoot relaxation solution found in {} iterations and {:.2f}s by {}\n", - root_relax_soln_.iterations, - root_relax_elapsed_time, - method_to_string(root_relax_solved_by)); + 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, + method_to_string(root_relax_solved_by)); settings_.log.printf("Root relaxation objective %+.8e\n\n", root_relax_soln_.user_objective); assert(root_vstatus_.size() == original_lp_.num_cols); @@ -2689,31 +3329,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } constexpr bool enable_root_cut_cpufj = true; - std::unique_ptr> root_cut_cpufj_task; - auto root_cut_cpufj_improvement_callback = + fj_cpu_worker_t root_fj_cpu_worker; + scope_guard cpufj_guard([&]() { root_fj_cpu_worker.stop(); }); + root_fj_cpu_worker.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 { - if (settings_.solution_callback != nullptr) { - settings_.solution_callback(user_assignment, obj); - } - set_solution_from_heuristics(user_assignment); - } + set_solution_from_cpu_fj(obj, assignment, work_units); }; - auto stop_root_cut_cpufj = [&]() { - if (!root_cut_cpufj_task) { return; } - mip::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; @@ -2746,13 +3367,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } cut_pass_result_t cut_pass_result; - if (root_cut_cpufj_task) { -#pragma omp task shared(root_cut_cpufj_task) priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) \ - depend(out : *root_cut_cpufj_task) - mip::run_fj_cpu_task(*root_cut_cpufj_task, - std::numeric_limits::infinity(), - std::numeric_limits::infinity()); - } + root_fj_cpu_worker.run_async(settings_.time_limit - toc(exploration_stats_.start_time)); cut_pass_result = do_cut_pass(cut_pass, solution, @@ -2772,11 +3387,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_objective, cut_pool_size, saved_solution); - - if (root_cut_cpufj_task) { - mip::stop_fj_cpu_task(*root_cut_cpufj_task); -#pragma omp taskwait depend(in : *root_cut_cpufj_task) - } + root_fj_cpu_worker.stop(); if (cut_pass_result.action == cut_pass_action_t::RETURN) { if (settings_.benchmark_info_ptr != nullptr) { @@ -2791,13 +3402,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut 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 = - mip::make_fj_cpu_task_from_host_lp(original_lp_, - var_types_, - root_relax_soln_.x, - settings_, - root_cut_cpufj_improvement_callback, - "[RootCut CPUFJ] "); + root_fj_cpu_worker.from_simplex_lp( + 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)); @@ -2835,30 +3441,28 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // 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 = - mip::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); + root_fj_cpu_worker.from_simplex_lp(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 * f_t{0.5}, f_t{1}); - mip::run_fj_cpu_task(*root_cut_cpufj_task, fj_time_limit, 0.5); - root_cut_cpufj_task.reset(); + settings_.deterministic ? remaining_time : std::min(remaining_time * 0.5, 1.0); + root_fj_cpu_worker.run_sync(fj_time_limit, 0.5); } set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); pc_.resize(original_lp_.num_cols); pc_.Arow = Arow_; - { + + if (!root_warm_start_) { raft::common::nvtx::range scope_sb("BB::strong_branching"); strong_branching(original_lp_, settings_, @@ -2975,8 +3579,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } else { const i_t num_workers = settings_.num_threads; const i_t num_bfs_workers = std::max(settings_.num_threads / 2, 1); - const i_t num_diving_workers = num_workers - num_bfs_workers; + 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 - num_submip_workers, 1); bfs_worker_pool_.init(num_bfs_workers, original_lp_, Arow_, var_types_, symmetry_, settings_); + submip_worker_pool_.init( + num_submip_workers, original_lp_, Arow_, var_types_, symmetry_, settings_, num_bfs_workers); if (num_diving_workers > 0) { diving_worker_pool_.init(num_diving_workers, @@ -2985,7 +3593,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut var_types_, symmetry_, settings_, - num_bfs_workers); + num_bfs_workers + num_submip_workers); } bfs_worker_t* initial_worker = bfs_worker_pool_.pop_idle_worker(); @@ -3090,8 +3698,8 @@ Work Units: 0 0.5 1. ──────────────────────────────────────────────────────────────────────────────────────────► Work Unit Time -Legend: ▓▓▓ = actively working ░░░ = waiting at barrier [hash] = state hash for verification - wut = work unit timestamp PC = pseudo-costs snap = snapshot (local copy) +Legend: ▓▓▓ = actively working ░░░ = waiting at barrier [hash] = state hash for +verification wut = work unit timestamp PC = pseudo-costs snap = snapshot (local copy) */ @@ -3129,23 +3737,22 @@ Producer Sync: Producing solutions in the past would break determinism, therefore this unidirectional sync ensures no such thing can occur. Instrumentation Aggregator: Collects multiple instrument vectors into a single aggregation point for estimating work from memory operations. Worker Context: Object -representing the "context" (e.g.: the worker) that should register the amount of work recorded There -is a 1context:1worker mapping. The Work Unit Scheduler registers such contexts and ensure they -remained synchronized together. Queued Integer Solutions: New integer solutions found within -horizons are queued with a work unit timestamp, in order to be sorted and played in order during the -sync callback. Creation Sequence: In nondeterministic mode, a single global atomic integer is used -to generate sequential IDs for the nodes. Since this is a global atomic, it is inherently +representing the "context" (e.g.: the worker) that should register the amount of work recorded +There is a 1context:1worker mapping. The Work Unit Scheduler registers such contexts and ensure +they remained synchronized together. Queued Integer Solutions: New integer solutions found within +horizons are queued with a work unit timestamp, in order to be sorted and played in order during +the sync callback. Creation Sequence: In nondeterministic mode, a single global atomic integer is +used to generate sequential IDs for the nodes. Since this is a global atomic, it is inherently nondeterministic. To fix this, in deterministic mode, nodes are addressed by a tuple - where "worker_id" is the ID of the worker that created this node, and "seq_id" is a sequential ID -local to the worker.\ This sequential ID is similar in principle to the global atomic ID sequence of -the nondeterminsitic mode but since it is local to each worker, it is updated serially and thus is -deterministic. worker IDs are unique, and sequence IDs are unique to their workers, therefor - is a globally unique node identifier. -Pseudocost Update: - Each worker updates its local pseudocosts when branching. These updates are queued within -horizons. During the horizon sync, these updates are all played in order, and the newly updated -global pseudocosts are broadcast to the worker's pseudocost snapshots for the coming horizon. + where "worker_id" is the ID of the worker that created this node, and "seq_id" is a sequential +ID local to the worker.\ This sequential ID is similar in principle to the global atomic ID +sequence of the nondeterminsitic mode but since it is local to each worker, it is updated serially +and thus is deterministic. worker IDs are unique, and sequence IDs are unique to their workers, +therefor is a globally unique node identifier. Pseudocost Update: Each worker +updates its local pseudocosts when branching. These updates are queued within horizons. During the +horizon sync, these updates are all played in order, and the newly updated global pseudocosts are +broadcast to the worker's pseudocost snapshots for the coming horizon. */ @@ -3171,8 +3778,8 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri if (num_diving_workers > 0) { // Extract diving types from search_strategies (skip BEST_FIRST at index 0) - std::vector diving_types(search_strategies + 1, - search_strategies + num_search_strategies); + std::vector diving_types; + get_diving_heuristic_list(settings_.diving_settings, diving_types); if (settings_.diving_settings.coefficient_diving != 0) { calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); @@ -3255,7 +3862,8 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri "Sync%% | NoWork\n"); settings_.log.printf( " " - "-------+---------+----------+--------+---------+--------+----------+----------+-------+-------" + "-------+---------+----------+--------+---------+--------+----------+----------+-------+-----" + "--" "\n"); for (const auto& worker : *deterministic_workers_) { double sync_time = worker.work_context.total_sync_time; @@ -3283,10 +3891,10 @@ void branch_and_bound_t::run_deterministic_coordinator(const csr_matri for (const auto& worker : *deterministic_diving_workers_) { const char* type_str = "???"; switch (worker.diving_type) { - case search_strategy_t::PSEUDOCOST_DIVING: type_str = "PC"; break; - case search_strategy_t::LINE_SEARCH_DIVING: type_str = "LS"; break; - case search_strategy_t::GUIDED_DIVING: type_str = "GD"; break; - case search_strategy_t::COEFFICIENT_DIVING: type_str = "CD"; break; + case worker_type_t::PSEUDOCOST_DIVING: type_str = "PC"; break; + case worker_type_t::LINE_SEARCH_DIVING: type_str = "LS"; break; + case worker_type_t::GUIDED_DIVING: type_str = "GD"; break; + case worker_type_t::COEFFICIENT_DIVING: type_str = "CD"; break; default: break; } settings_.log.printf(" %6d | %6s | %7d | %6d | %6d | %7.3fs | %5.2fs\n", @@ -3603,7 +4211,7 @@ node_status_t branch_and_bound_t::solve_node_deterministic( worker.clock += work_performed; exploration_stats_.total_lp_solve_time += toc(lp_start_time); - exploration_stats_.total_lp_iters += node_iter; + exploration_stats_.total_simplex_iters += node_iter; ++exploration_stats_.nodes_explored; --exploration_stats_.nodes_unexplored; @@ -3642,7 +4250,7 @@ void branch_and_bound_t::deterministic_process_worker_solutions( i_t nodes_explored = exploration_stats_.nodes_explored.load(); i_t nodes_unexplored = exploration_stats_.nodes_unexplored.load(); - search_strategy_t worker_type = get_worker_type(pool, sol->worker_id); + worker_type_t worker_type = get_worker_type(pool, sol->worker_id); report(feasible_solution_symbol(worker_type, settings_.diving_settings.show_type), sol->objective, deterministic_lower, @@ -3690,10 +4298,10 @@ void branch_and_bound_t::deterministic_broadcast_snapshots( PoolT& pool, const std::vector& incumbent_snapshot) { deterministic_snapshot_t snap{ - .upper_bound = upper_bound_, - .pc_snapshot = pc_, - .incumbent = incumbent_snapshot, - .total_lp_iters = exploration_stats_.total_lp_iters, + .upper_bound = upper_bound_, + .pc_snapshot = pc_, + .incumbent = incumbent_snapshot, + .total_simplex_iters = exploration_stats_.total_simplex_iters, }; for (auto& worker : pool) { @@ -3822,7 +4430,7 @@ void branch_and_bound_t::deterministic_sort_replay_events( } if (new_upper < std::numeric_limits::infinity()) { - report_heuristic(new_upper); + report_heuristic(new_upper, HEURISTICS); if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -3836,7 +4444,7 @@ void branch_and_bound_t::deterministic_sort_replay_events( // Merge integer solutions from BFS workers and update global incumbent deterministic_process_worker_solutions(*deterministic_workers_, [](const deterministic_bfs_worker_pool_t&, int) { - return search_strategy_t::BEST_FIRST; + return worker_type_t::BEST_FIRST; }); // Merge and apply pseudo-cost updates from BFS workers diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 12c93fcd91..e3225e2e57 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -24,6 +24,8 @@ #include #include +#include + #include #include #include @@ -41,29 +43,13 @@ #include namespace cuopt::mathematical_optimization::mip { + template struct clique_table_t; -} - -namespace cuopt::mathematical_optimization::mip { template struct mip_symmetry_t; -enum class mip_status_t { - OPTIMAL = 0, // The optimal integer solution was found - UNBOUNDED = 1, // The problem is unbounded - INFEASIBLE = 2, // The problem is infeasible - TIME_LIMIT = 3, // The solver reached a time limit - NODE_LIMIT = 4, // The maximum number of nodes was reached (not implemented) - NUMERICAL = 5, // The solver encountered a numerical error - UNSET = 6, // The status is not set - WORK_LIMIT = 7, // The solver reached a deterministic work limit -}; - -template -void upper_bound_callback(f_t upper_bound); - template struct nondeterministic_policy_t; template @@ -86,6 +72,11 @@ class branch_and_bound_t { // Set an initial guess based on the user_problem. This should be called before solve. void set_initial_guess(const std::vector& user_guess) { guess_ = user_guess; } + void set_halt_callback(std::function callback) + { + halt_callback_ = std::move(callback); + } + // Set the root solution found by PDLP void set_root_relaxation_solution(const std::vector& primal, const std::vector& dual, @@ -109,7 +100,10 @@ class branch_and_bound_t { } // Set a solution based on the user problem during the course of the solve - void set_solution_from_heuristics(const std::vector& solution); + bool set_solution_from_heuristics(const std::vector& solution, worker_type_t heuristic_type); + + // Apply a solution found by a CPU FJ worker. + void set_solution_from_cpu_fj(f_t obj, const std::vector& assignment, double work_units); // This queues the solution to be processed at the correct work unit timestamp void queue_external_solution_deterministic(const std::vector& solution, double work_unit_ts); @@ -125,6 +119,9 @@ class branch_and_bound_t { // `bound` must be in B&B's internal objective space. void set_initial_upper_bound(f_t bound); + void warm_start(const pseudo_costs_t& parent_pc, + const std::vector& reduced_to_original); + f_t get_upper_bound() const { return upper_bound_.load(); } bool has_solver_space_incumbent() const { return incumbent_.has_incumbent; } @@ -151,6 +148,15 @@ class branch_and_bound_t { std::vector& lower_bounds, std::vector& upper_bounds); + // Whether obj should replace the stored incumbent. Must be called under mutex_upper_. + // Compares against the stored incumbent's objective, NOT against upper_bound_, because + // set_initial_upper_bound can set a tighter bound from an OG-space solution that has no + // corresponding solver-space incumbent (e.g. papilo can't crush it back). + bool improves_incumbent(f_t obj) const + { + return !incumbent_.has_incumbent || obj < incumbent_.objective; + } + // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(simplex::mip_solution_t& solution); @@ -199,18 +205,13 @@ class branch_and_bound_t { // original-space in the mip_solver_context_t), but does NOT imply incumbent_.has_incumbent. omp_atomic_t upper_bound_; + // Callback for halting the solver. This passes the current upper and lower bound of the solver + // in user space. + std::function halt_callback_; + // Solver-space incumbent tracked directly by B&B. simplex::mip_solution_t incumbent_; - // Whether obj should replace the stored incumbent. Must be called under mutex_upper_. - // Compares against the stored incumbent's objective, NOT against upper_bound_, because - // set_initial_upper_bound can set a tighter bound from an OG-space solution that has no - // corresponding solver-space incumbent (e.g. papilo can't crush it back). - bool improves_incumbent(f_t obj) const - { - return !incumbent_.has_incumbent || obj < incumbent_.objective; - } - // Structure with the general info of the solver. branch_and_bound_stats_t exploration_stats_; @@ -233,6 +234,7 @@ class branch_and_bound_t { std::atomic root_concurrent_halt_{0}; std::atomic node_concurrent_halt_{0}; bool is_root_solution_set{false}; + bool root_warm_start_{false}; // Pseudocosts pseudo_costs_t pc_; @@ -246,6 +248,9 @@ class branch_and_bound_t { // Worker pool dedicated to diving diving_worker_pool_t diving_worker_pool_; + diving_worker_pool_t submip_worker_pool_; + submip_stats_t submip_stats_; + // Global status of the solver. omp_atomic_t solver_status_; omp_atomic_t is_running_{false}; @@ -262,7 +267,7 @@ class branch_and_bound_t { std::function user_bound_callback_; void print_table_header(); - void report_heuristic(f_t obj); + void report_heuristic(f_t obj, worker_type_t type); void report(char symbol, f_t obj, f_t lower_bound, @@ -308,7 +313,7 @@ class branch_and_bound_t { void add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, - search_strategy_t thread_type); + worker_type_t thread_type); // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); @@ -316,6 +321,14 @@ class branch_and_bound_t { // Launch a new diving worker from a given best-first worker. bool launch_diving_worker(bfs_worker_t* bfs_worker); + // If the objective is integral or must move in steps than + // the lower bound will be different from the leaf objective. + // We use the leaf objective for RINS (on_optimal_callback) + // and if we are integer feasible (handle_integer_solution). + // We use the lower bound to decide if we should fathom the + // node or branch. + void apply_objective_step(mip_node_t* node_ptr, f_t leaf_obj); + // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); @@ -332,7 +345,21 @@ class branch_and_bound_t { // Perform a deep dive in the subtree determined by the `start_node` in order // to find integer feasible solutions. - void dive_with(diving_worker_t* worker); + void dive_with(diving_worker_t* worker, i_t backtrack_limit); + + // Launch a new RINS/RENS worker + bool launch_local_branching_worker(const std::vector& sol); + + // Solve the RINS/RENS sub-MIP + void solve_submip(diving_worker_t* worker, + const std::vector& current_incumbent, + i_t num_var_fixed, + i_t num_integers, + i_t submip_level, + std::string_view log_prefix); + + // Creates and solves the RINS sub-MIP + void rins(diving_worker_t* rins_worker, const std::vector& node_solution); // 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 0194e23da4..2146109961 100644 --- a/cpp/src/branch_and_bound/constants.hpp +++ b/cpp/src/branch_and_bound/constants.hpp @@ -7,9 +7,9 @@ #pragma once -namespace cuopt::mathematical_optimization::mip { +#include -constexpr int num_search_strategies = 7; +namespace cuopt::mathematical_optimization::mip { // Indicate the search and variable selection algorithms used by each thread // in B&B (See [1]). @@ -19,26 +19,51 @@ constexpr int num_search_strategies = 7; // [2] J. Witzig and A. Gleixner, “Conflict-Driven Heuristics for Mixed Integer Programming,” // Feb. 07, 2019, _arXiv_: arXiv:1902.02615. doi: // [10.48550/arXiv.1902.02615](https://doi.org/10.48550/arXiv.1902.02615). -enum search_strategy_t : int { +enum worker_type_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) FARKAS_DIVING = 5, // Farkas Diving (see [2]) - VECTOR_LENGTH_DIVING = 6 // Vector Length Diving (9.2.6) + VECTOR_LENGTH_DIVING = 6, // Vector Length Diving (9.2.6) + SUBMIP = 7, // RINS/RENS (akin to a guided diving, see HiGHS) + HEURISTICS = 8, // Other heuristics + NUM_WORKER_TYPES }; -constexpr search_strategy_t search_strategies[] = {BEST_FIRST, - PSEUDOCOST_DIVING, - LINE_SEARCH_DIVING, - GUIDED_DIVING, - COEFFICIENT_DIVING, - FARKAS_DIVING, - VECTOR_LENGTH_DIVING}; - enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 }; enum class branch_and_bound_mode_t { PARALLEL = 0, DETERMINISTIC = 1 }; +enum class mip_status_t { + OPTIMAL = 0, // The optimal integer solution was found + UNBOUNDED = 1, // The problem is unbounded + INFEASIBLE = 2, // The problem is infeasible + TIME_LIMIT = 3, // The solver reached a time limit + NODE_LIMIT = 4, // The maximum number of nodes was reached + ITERATION_LIMIT = 5, // The maximum number of simplex iterations was reached + NUMERICAL = 6, // The solver encountered a numerical error + UNSET = 7, // The status is not set + WORK_LIMIT = 8, // The solver reached a deterministic work limit + HALT = 9 // Halt the solver +}; + +inline std::string mip_status_to_string(mip_status_t status) +{ + switch (status) { + case mip_status_t::OPTIMAL: return "OPTIMAL"; + case mip_status_t::UNBOUNDED: return "UNBOUNDED"; + case mip_status_t::INFEASIBLE: return "INFEASIBLE"; + case mip_status_t::TIME_LIMIT: return "TIME_LIMIT"; + case mip_status_t::NODE_LIMIT: return "NODE_LIMIT"; + case mip_status_t::ITERATION_LIMIT: return "ITERATION_LIMIT"; + case mip_status_t::NUMERICAL: return "NUMERICAL"; + case mip_status_t::UNSET: return "UNSET"; + case mip_status_t::WORK_LIMIT: return "WORK_LIMIT"; + case mip_status_t::HALT: return "SUBMIP_HALT"; + } + 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 60d3436401..3e09d6e606 100644 --- a/cpp/src/branch_and_bound/deterministic_workers.hpp +++ b/cpp/src/branch_and_bound/deterministic_workers.hpp @@ -58,7 +58,7 @@ struct deterministic_snapshot_t { f_t upper_bound; pseudo_cost_snapshot_t pc_snapshot; std::vector incumbent; - int64_t total_lp_iters; + int64_t total_simplex_iters; }; template @@ -102,7 +102,7 @@ class deterministic_worker_base_t : public branch_and_bound_worker_t { local_upper_bound = snap.upper_bound; pc_snapshot = snap.pc_snapshot; incumbent_snapshot = snap.incumbent; - total_lp_iters_snapshot = snap.total_lp_iters; + total_lp_iters_snapshot = snap.total_simplex_iters; } bool has_work() const { return static_cast(this)->has_work_impl(); } @@ -275,7 +275,7 @@ class deterministic_diving_worker_t using base_t = deterministic_worker_base_t>; public: - search_strategy_t diving_type{search_strategy_t::PSEUDOCOST_DIVING}; + worker_type_t diving_type{worker_type_t::PSEUDOCOST_DIVING}; // Diving-specific node management std::deque> dive_queue; @@ -295,7 +295,7 @@ class deterministic_diving_worker_t explicit deterministic_diving_worker_t( int id, - search_strategy_t type, + worker_type_t type, const simplex::lp_problem_t& original_lp, const csr_matrix_t& Arow, const std::vector& var_types, @@ -438,7 +438,7 @@ class deterministic_diving_worker_pool_t public: deterministic_diving_worker_pool_t(int num_workers, - const std::vector& diving_types, + const std::vector& diving_types, const simplex::lp_problem_t& original_lp, const csr_matrix_t& Arow, const std::vector& var_types, @@ -447,7 +447,7 @@ class deterministic_diving_worker_pool_t { this->workers_.reserve(num_workers); for (int i = 0; i < num_workers; ++i) { - search_strategy_t type = diving_types[i % diving_types.size()]; + worker_type_t type = diving_types[i % diving_types.size()]; this->workers_.emplace_back(i, type, original_lp, Arow, var_types, settings, root_solution); } } diff --git a/cpp/src/branch_and_bound/diving_heuristics.hpp b/cpp/src/branch_and_bound/diving_heuristics.hpp index f778ae38dc..a3c48808e4 100644 --- a/cpp/src/branch_and_bound/diving_heuristics.hpp +++ b/cpp/src/branch_and_bound/diving_heuristics.hpp @@ -18,38 +18,17 @@ namespace cuopt::mathematical_optimization::mip { -// When `log_diving_type` is true, each diving strategy gets its own letter; -// otherwise every dive collapses to 'D'. -inline char feasible_solution_symbol(search_strategy_t strategy, bool log_diving_type) -{ - if (strategy == BEST_FIRST) return 'B'; - if (!log_diving_type) { return 'D'; } - switch (strategy) { - case COEFFICIENT_DIVING: return 'C'; - case LINE_SEARCH_DIVING: return 'L'; - case PSEUDOCOST_DIVING: return 'P'; - case GUIDED_DIVING: return 'G'; - case FARKAS_DIVING: return 'F'; - case VECTOR_LENGTH_DIVING: return 'V'; - default: return 'U'; - } -} - template -bool is_search_strategy_enabled(search_strategy_t strategy, - const mip_diving_hyper_params_t& settings) +void get_diving_heuristic_list(const mip_diving_hyper_params_t& settings, + std::vector& heuristic_list) { - switch (strategy) { - case BEST_FIRST: return true; - case PSEUDOCOST_DIVING: return settings.pseudocost_diving != 0; - case LINE_SEARCH_DIVING: return settings.line_search_diving != 0; - case GUIDED_DIVING: return settings.guided_diving != 0; - case COEFFICIENT_DIVING: return settings.coefficient_diving != 0; - case FARKAS_DIVING: return settings.farkas_diving != 0; - case VECTOR_LENGTH_DIVING: return settings.vector_length_diving != 0; - } - - return false; + heuristic_list.clear(); + if (settings.pseudocost_diving != 0) heuristic_list.push_back(PSEUDOCOST_DIVING); + if (settings.line_search_diving != 0) heuristic_list.push_back(LINE_SEARCH_DIVING); + if (settings.guided_diving != 0) heuristic_list.push_back(GUIDED_DIVING); + if (settings.coefficient_diving != 0) heuristic_list.push_back(COEFFICIENT_DIVING); + if (settings.farkas_diving != 0) heuristic_list.push_back(FARKAS_DIVING); + if (settings.vector_length_diving != 0) heuristic_list.push_back(VECTOR_LENGTH_DIVING); } template diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp index c763e5ada5..3e291124c3 100644 --- a/cpp/src/branch_and_bound/mip_node.hpp +++ b/cpp/src/branch_and_bound/mip_node.hpp @@ -147,7 +147,7 @@ class mip_node_t { { update_branched_variable_bounds(lower, upper, bounds_changed); - mip_node_t* parent_ptr = parent; + mip_node_t* parent_ptr = parent; while (parent_ptr != nullptr && parent_ptr->node_id != 0) { parent_ptr->update_branched_variable_bounds(lower, upper, bounds_changed); parent_ptr = parent_ptr->parent; @@ -160,6 +160,9 @@ class mip_node_t { std::vector& upper, std::vector& bounds_changed) const { + // We in the root node and the lower/upper are already set to their starting value + if (parent == nullptr) return; + assert(branch_var >= 0); assert(lower.size() > branch_var); assert(upper.size() > branch_var); diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 82aca6c458..cdba90f219 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -1040,7 +1040,7 @@ void strong_branching_reduced(const lp_problem_t& original_lp, i_t effective_batch_pdlp = settings.mip_batch_pdlp_strong_branching; // Disable for sub MIP - if (settings.sub_mip) { effective_batch_pdlp = 0; } + if (settings.inside_submip) { effective_batch_pdlp = 0; } // Disable if running in deterministic mode if (settings.deterministic && settings.mip_batch_pdlp_strong_branching == 1) { @@ -1053,7 +1053,7 @@ void strong_branching_reduced(const lp_problem_t& original_lp, } if (settings.mip_batch_pdlp_strong_branching != 0 && - (settings.sub_mip || settings.deterministic)) { + (settings.inside_submip || settings.deterministic)) { settings.log.printf( "Batch PDLP strong branching is disabled because sub-MIP or deterministic mode is enabled\n"); } @@ -1423,6 +1423,10 @@ template void pseudo_costs_t::update_pseudo_costs(mip_node_t* node_ptr, f_t leaf_objective) { + // The root node carries no branching decision (branch_var < 0), so there is no variable to + // attribute the objective change to. + if (node_ptr->branch_var < 0) return; + const f_t change_in_obj = std::max(leaf_objective - node_ptr->lower_bound, 0.0); const f_t frac = node_ptr->branch_dir == branch_direction_t::DOWN ? node_ptr->fractional_val - std::floor(node_ptr->fractional_val) @@ -1489,7 +1493,7 @@ i_t pseudo_costs_t::reliable_variable_selection( f_t avg_up{0}; lp_solution_t& leaf_solution = worker->leaf_solution; - const int64_t branch_and_bound_lp_iters = bnb_stats.total_lp_iters; + const int64_t branch_and_bound_lp_iters = bnb_stats.total_simplex_iters; const i_t branch_and_bound_lp_iter_per_node = bnb_stats.nodes_explored > 0 ? branch_and_bound_lp_iters / bnb_stats.nodes_explored : 0; const i_t iter_limit_per_trial = std::clamp(2 * branch_and_bound_lp_iter_per_node, @@ -1567,7 +1571,7 @@ i_t pseudo_costs_t::reliable_variable_selection( // Use the heuristic to decide if it should be used (in case it is set to automatic) if (!use_pdlp && rb_mode != 0) { // Check if it is a sub MIP or the determinism mode is on. - use_pdlp = !settings.sub_mip; + use_pdlp = !settings.inside_submip; use_pdlp &= !settings.deterministic; // Check if the warm cache was filled at the root @@ -1587,7 +1591,7 @@ i_t pseudo_costs_t::reliable_variable_selection( // Use the heuristic to decide if it should be used (in case it is set to automatic) if (!use_pdlp && rb_mode != 0) { // Check if it is a sub MIP or the determinism mode is on. - use_pdlp = !settings.sub_mip; + use_pdlp = !settings.inside_submip; use_pdlp &= !settings.deterministic; // Check if the warm cache was filled at the root @@ -1603,7 +1607,7 @@ i_t pseudo_costs_t::reliable_variable_selection( if (rb_mode != 0 && !pdlp_warm_cache->populated) { settings.log.debug("PDLP warm start data not populated, using DS only\n"); - } else if (rb_mode != 0 && settings.sub_mip) { + } else if (rb_mode != 0 && settings.inside_submip) { settings.log.debug("Batch PDLP reliability branching is disabled because sub-MIP is enabled\n"); } else if (rb_mode != 0 && settings.deterministic) { settings.log.debug( diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index eddf2cf28f..6cc34957fd 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -156,12 +157,45 @@ class pseudo_costs_t { void resize(i_t num_variables) { - pseudo_cost_sum_down.assign(num_variables, 0); - pseudo_cost_sum_up.assign(num_variables, 0); - pseudo_cost_num_down.assign(num_variables, 0); - pseudo_cost_num_up.assign(num_variables, 0); pseudo_cost_mutex_up.resize(num_variables); pseudo_cost_mutex_down.resize(num_variables); + + if (!warm_start) { + std::fill(pseudo_cost_sum_down.begin(), pseudo_cost_sum_down.end(), 0); + std::fill(pseudo_cost_sum_up.begin(), pseudo_cost_sum_up.end(), 0); + std::fill(pseudo_cost_num_down.begin(), pseudo_cost_num_down.end(), 0); + std::fill(pseudo_cost_num_up.begin(), pseudo_cost_num_up.end(), 0); + } + + pseudo_cost_sum_down.resize(num_variables, 0); + pseudo_cost_sum_up.resize(num_variables, 0); + pseudo_cost_num_down.resize(num_variables, 0); + pseudo_cost_num_up.resize(num_variables, 0); + } + + void warm_start_from(const pseudo_costs_t& parent, const std::vector& reduced_to_original) + { + warm_start = true; + + for (i_t k = 0; k < reduced_to_original.size(); ++k) { + const i_t orig = reduced_to_original[k]; + assert(orig >= 0); + assert(orig < parent.pseudo_cost_num_up.size()); + + const i_t parent_num_up = parent.pseudo_cost_num_up[orig]; + if (parent_num_up > 0) { + const f_t value = parent.pseudo_cost_sum_up[orig] / parent_num_up; + pseudo_cost_num_up[k] = 1; + pseudo_cost_sum_up[k] = value; + } + + const i_t parent_num_down = parent.pseudo_cost_num_down[orig]; + if (parent_num_down > 0) { + const f_t value = parent.pseudo_cost_sum_down[orig] / parent_num_down; + pseudo_cost_num_down[k] = 1; + pseudo_cost_sum_down[k] = value; + } + } } f_t get_pseudocost_down(i_t j, f_t avg) const @@ -229,6 +263,7 @@ class pseudo_costs_t { std::vector pseudo_cost_mutex_down; omp_atomic_t strong_branching_lp_iter = 0; + bool warm_start = false; }; template diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index 4fd278417e..b238dee523 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -31,9 +31,9 @@ struct branch_and_bound_stats_t { // Tracks the number of nodes being solved by the workers at a given time omp_atomic_t nodes_being_solved = 0; - omp_atomic_t total_lp_iters = 0; - omp_atomic_t nodes_since_last_log = 0; - omp_atomic_t last_log = 0.0; + omp_atomic_t total_simplex_iters = 0; + omp_atomic_t nodes_since_last_log = 0; + omp_atomic_t last_log = 0.0; omp_atomic_t orbital_fixing_nodes = 0; omp_atomic_t orbital_fixings_applied = 0; @@ -50,7 +50,7 @@ class branch_and_bound_worker_t { using int_type = i_t; i_t worker_id; - omp_atomic_t search_strategy; + omp_atomic_t search_strategy; omp_atomic_t is_active; omp_atomic_t lower_bound; @@ -75,6 +75,10 @@ class branch_and_bound_worker_t { std::unique_ptr> lexical_reduction; mip_symmetry_t* symmetry_ptr = nullptr; + bool skip_set_bounds = false; + bool recompute_basis = true; + bool recompute_bounds = true; + void ensure_orbital_fixing() { if (orbital_fixing == nullptr && symmetry_ptr != nullptr) { @@ -86,9 +90,6 @@ class branch_and_bound_worker_t { } } - bool recompute_basis = true; - bool recompute_bounds = true; - branch_and_bound_worker_t(i_t worker_id, const simplex::lp_problem_t& original_lp, const csr_matrix_t& Arow, @@ -116,6 +117,9 @@ class branch_and_bound_worker_t { bool set_lp_variable_bounds(mip_node_t* node_ptr, const simplex::simplex_solver_settings_t& settings) { + // In RINS/RENS, the bounds are already set in the parent method. + if (skip_set_bounds) return true; + // Reset the bound_changed markers std::fill(bounds_changed.begin(), bounds_changed.end(), false); @@ -152,10 +156,9 @@ class bfs_worker_t : public branch_and_bound_worker_t { this->start_lower = original_lp.lower; this->start_upper = original_lp.upper; this->search_strategy = BEST_FIRST; - - max_diving_workers.fill(0); - active_diving_workers.fill(0); - total_active_diving_workers = 0; + max_diving_workers = 0; + active_diving_workers = 0; + next_heuristic = worker_id; } void set_inactive() { this->is_active = false; } @@ -172,53 +175,44 @@ class bfs_worker_t : public branch_and_bound_worker_t { return node_queue.steal_from(victim->node_queue, nodes_to_steal); } - // Calculate the number of diving workers that this worker can launch. Having a fixed number - // of workers allows the solver to be more deterministic. - void calculate_num_diving_workers(i_t num_bfs_workers, - i_t total_diving_workers, - const mip_diving_hyper_params_t& settings) + void calculate_max_diving_workers(i_t num_bfs_workers, i_t total_diving_workers) { - i_t num_active = 0; - for (i_t i = 1; i < num_search_strategies; ++i) { - num_active += is_search_strategy_enabled(search_strategies[i], settings); - } + auto [start, end] = + calculate_index_range(this->worker_id, total_diving_workers, num_bfs_workers); + max_diving_workers = end - start; + } - total_max_diving_workers = 0; - max_diving_workers.fill(0); - if (num_active == 0) { return; } - - for (size_t i = 1, k = 0; i < num_search_strategies; ++i) { - if (is_search_strategy_enabled(search_strategies[i], settings)) { - // Calculate the number of workers for a given diving heuristic - auto [type_start, type_end] = calculate_index_range(k, total_diving_workers, num_active); - i_t workers_per_type = type_end - type_start; - - // Calculate the number of diving workers allocated to this (best-first) worker - auto [start, end] = - calculate_index_range(this->worker_id, workers_per_type, num_bfs_workers); - max_diving_workers[i] = end - start; - total_max_diving_workers += max_diving_workers[i]; - ++k; - } - } + // Calculate the number of diving workers that this worker can launch. And update the list + // with all active diving heuristics. + void update_diving_heuristic_list(const mip_diving_hyper_params_t& settings) + { + get_diving_heuristic_list(settings, diving_heuristics); } - // The worker-local node heap. - node_queue_t node_queue; + // Get the next diving heuristic from the list + worker_type_t next_diving_heuristic() + { + assert(is_diving_enabled()); + next_heuristic = next_heuristic % diving_heuristics.size(); + return diving_heuristics[next_heuristic++]; + } - // The number of diving workers of each type that this (best-first) worker can launch. - std::array max_diving_workers; + bool is_diving_enabled() { return !diving_heuristics.empty(); } - // The number of active diving workers of each type associated with this (best-first) worker. - std::array, num_search_strategies> active_diving_workers; + // The worker-local node heap. + node_queue_t node_queue; - // Keep track of the total number of active diving worker that are associated with this + // Keep track of the number of active diving worker that are associated with this // (best-first) worker - omp_atomic_t total_active_diving_workers{0}; + omp_atomic_t active_diving_workers; // The maximum number of diving worker that are associated with this // (best-first) worker - i_t total_max_diving_workers{0}; + i_t max_diving_workers; + + private: + std::vector diving_heuristics; + i_t next_heuristic; }; template @@ -238,13 +232,12 @@ class diving_worker_t : public branch_and_bound_worker_t { void set_inactive() { if (!this->is_active.load()) { return; } - assert(bfs_worker != nullptr); - assert(bfs_worker->active_diving_workers[this->search_strategy].load() > 0); - assert(bfs_worker->total_active_diving_workers.load() > 0); - this->is_active = false; - --bfs_worker->active_diving_workers[this->search_strategy]; - --bfs_worker->total_active_diving_workers; + + if (bfs_worker) { + assert(bfs_worker->active_diving_workers.load() > 0); + --bfs_worker->active_diving_workers; + } } f_t get_lower_bound() { return this->lower_bound; } @@ -256,4 +249,27 @@ class diving_worker_t : public branch_and_bound_worker_t { bfs_worker_t* bfs_worker{nullptr}; }; +struct submip_stats_t { + omp_atomic_t total_success = 0; + omp_atomic_t success_fixrate_sum = 0; + omp_atomic_t total_infeasible = 0; + omp_atomic_t infeasible_fixrate_sum = 0; + omp_atomic_t total_calls = 0; + + void save_success(double fixrate) + { + ++total_success; + success_fixrate_sum += fixrate; + } + + void save_infeasible(double fixrate) + { + ++total_infeasible; + infeasible_fixrate_sum += fixrate; + } + + double average_infeasible_fixrate() const { return infeasible_fixrate_sum / total_infeasible; } + double average_success_fixrate() const { return success_fixrate_sum / total_success; } +}; + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index ef82c0a940..de7042806c 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -1485,11 +1485,19 @@ bool flow_cover_is_zero_one_integer_variable(const flow_cover_context_t -f_t flow_cover_arc_tol(const flow_cover_context_t& context, - const single_node_flow_arc_t& arc) +f_t flow_cover_arc_lower_tol(const flow_cover_context_t& context) { - return flow_cover_scaled_primal_tol(context, arc.u); + return std::max(10 * context.settings.primal_tol, static_cast(1e-5)); +} + +template +f_t flow_cover_arc_upper_tol(const flow_cover_context_t& context, f_t u) +{ + return std::max(100 * context.settings.primal_tol, static_cast(1e-4)) * + std::max(1.0, u); } template @@ -1878,11 +1886,8 @@ bool flow_cover_generation_t::build_single_node_flow_relaxation( scratch.candidates.end(), [](const single_node_flow_candidate_t& a, const single_node_flow_candidate_t& b) { return a.distance < b.distance; }); - const f_t arc_lower_tol = - std::max(10 * context.settings.primal_tol, static_cast(1e-5)); - const f_t arc_upper_tol = - std::max(100 * context.settings.primal_tol, static_cast(1e-4)) * - std::max(1.0, best->arc.u); + const f_t arc_lower_tol = flow_cover_arc_lower_tol(context); + const f_t arc_upper_tol = flow_cover_arc_upper_tol(context, best->arc.u); if (best->arc.y_value < -arc_lower_tol || best->arc.y_value > best->arc.u * best->arc.x_value + arc_upper_tol) { return false; @@ -1899,7 +1904,15 @@ bool flow_cover_generation_t::build_single_node_flow_relaxation( const f_t coeff = scratch.binary_coefficients[j]; if (std::abs(coeff) <= coefficient_tol) { continue; } const f_t u = std::abs(coeff); - scratch.arcs.push_back(flow_cover_build_arc(context, u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u)); + single_node_flow_arc_t arc = + flow_cover_build_arc(context, u, coeff < 0.0, j, 0.0, 0.0, -1, 0.0, u); + // y_value is built from the raw LP value (u * xstar[j]) while x_value is clamped to [0, 1]. + // When xstar[j] sits marginally below its bound (e.g. -5e-7, normal LP feasibility slop), + // y_value picks up a tiny negative (u * xstar) that trips the arc's lower gate even though the + // flow is really 0. Continuous-term arcs are already clamped this way above (see the + // `best->arc.y_value < 0.0` clamp); the binary arcs bypassed both the gate and the clamp. + if (arc.y_value < 0.0) { arc.y_value = 0.0; } + scratch.arcs.push_back(arc); } if (scratch.arcs.empty()) { return false; } @@ -1909,16 +1922,25 @@ bool flow_cover_generation_t::build_single_node_flow_relaxation( [&]() { f_t single_node_flow_activity = 0.0; f_t single_node_flow_scale = std::max(1.0, std::abs(single_node_flow_b)); + // Each accepted arc may sit up to arc_{lower,upper}_tol outside its capacity — the same + // slack the acceptance gate above tolerates — and the y_value<0 -> 0 clamp can shift it by + // another arc_lower_tol. Those per-arc slacks accumulate into the aggregate activity, so the + // relaxation tolerance must sum them; primal_tol * scale alone (which only covers FP rounding + // in the summation) is far too tight and rejects LP points that are genuinely inside. + f_t single_node_flow_slack = 0.0; for (const auto& arc : scratch.arcs) { - const f_t arc_tol = flow_cover_arc_tol(context, arc); - if (arc.y_value < -arc_tol) { return false; } - if (arc.y_value > arc.u * arc.x_value + arc_tol) { return false; } + const f_t arc_lower_tol = flow_cover_arc_lower_tol(context); + const f_t arc_upper_tol = flow_cover_arc_upper_tol(context, arc.u); + if (arc.y_value < -arc_lower_tol) { return false; } + if (arc.y_value > arc.u * arc.x_value + arc_upper_tol) { return false; } const f_t signed_y = arc.in_n2 ? -arc.y_value : arc.y_value; single_node_flow_activity += signed_y; single_node_flow_scale += std::abs(signed_y); + single_node_flow_slack += arc_lower_tol + arc_upper_tol; } - return single_node_flow_activity <= - single_node_flow_b + context.settings.primal_tol * single_node_flow_scale; + return single_node_flow_activity <= single_node_flow_b + + context.settings.primal_tol * single_node_flow_scale + + single_node_flow_slack; }(), "Flow cover single-node-flow relaxation excludes LP solution"); @@ -4157,7 +4179,7 @@ void cut_generation_t::generate_gomory_cuts( mixed_integer_gomory_cut_t gomory_cut; complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); simplex_solver_settings_t variable_settings = settings; - variable_settings.sub_mip = 1; + variable_settings.inside_submip = 1; variable_bounds_t variable_bounds(lp, variable_settings, var_types, Arow, new_slacks); strong_cg_cut_t cg(lp, var_types, xstar); std::vector transformed_xstar; @@ -4492,7 +4514,7 @@ variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, num_pos_inf_(lp.num_rows, 0), num_neg_inf_(lp.num_rows, 0) { - if (settings.sub_mip) { + if (settings.inside_submip) { return; // Don't compute the variable upper/lower bounds inside sub-MIP } f_t start_time = tic(); diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 9adf5bf0e9..2048df00fa 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -676,6 +676,138 @@ i_t add_artifical_variables(lp_problem_t& problem, return 0; } +// Inverse of convert_user_problem (LP/MIP only -- no quadratic or SOC handling): +// take a problem in simplex standard form +// minimize c^T x +// subject to A x = b +// l <= x <= u +// and express it as a user_problem in range (row-bounded) form +// minimize c^T x +// subject to l_row <= A x <= u_row +// l <= x <= u +template +void convert_simplex_problem(const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const std::vector& new_slacks, + user_problem_t& user_problem) +{ + constexpr bool verbose = false; + if (verbose) { + settings.log.printf("Converting simplex problem with %d rows and %d columns and %d nonzeros\n", + simplex_problem.num_rows, + simplex_problem.num_cols, + simplex_problem.A.col_start[simplex_problem.num_cols]); + } + + const i_t m = simplex_problem.num_rows; + const i_t n = simplex_problem.num_cols; + const csc_matrix_t& simplex_A = simplex_problem.A; + assert((var_types.empty() || var_types.size() == n) && + "var_types must span the full simplex problem (structural + slack columns)"); + + const i_t num_slacks = new_slacks.size(); + const i_t new_n = n - num_slacks; + const i_t new_nnz = simplex_A.col_start[n] - num_slacks; + + // Reset the destination so we can populate user_problem in place. + user_problem.num_rows = m; + user_problem.num_cols = new_n; + user_problem.range_rows.clear(); + user_problem.range_value.clear(); + user_problem.objective.clear(); + user_problem.lower.clear(); + user_problem.upper.clear(); + user_problem.var_types.clear(); + user_problem.objective.reserve(new_n); + user_problem.lower.reserve(new_n); + user_problem.upper.reserve(new_n); + if (!var_types.empty()) { user_problem.var_types.reserve(new_n); } + + // Rows with no associated slack stay equalities at their rhs; the slack loop + // below overrides the rows that had one. + user_problem.row_sense.assign(m, 'E'); + user_problem.rhs = simplex_problem.rhs; + + // Recover each row's constraint from its slack. The forward equality + // a_i^T x_struct + c * s = b_i (with l_s <= s <= u_s) gives + // a_i^T x_struct in [min(b - c l_s, b - c u_s), max(b - c l_s, b - c u_s)], + // which maps back to row_sense / rhs / range. Flag the slack columns so they + // can be dropped from the matrix below. + std::vector is_slack(n, false); + for (i_t s : new_slacks) { + assert(s >= 0 && s < n && "slack index out of range"); + is_slack[s] = true; + + const i_t col_start = simplex_A.col_start[s]; + const i_t col_end = simplex_A.col_start[s + 1]; + assert(col_end - col_start == 1 && "slack/artificial column must have a single nonzero"); + const i_t row = simplex_A.i[col_start]; + const f_t c = simplex_A.x[col_start]; + const f_t b = simplex_problem.rhs[row]; + const f_t t1 = b - c * simplex_problem.lower[s]; + const f_t t2 = b - c * simplex_problem.upper[s]; + const f_t lo = std::min(t1, t2); + const f_t hi = std::max(t1, t2); + if (lo == hi) { + user_problem.row_sense[row] = 'E'; + user_problem.rhs[row] = lo; + } else if (lo == -inf) { + user_problem.row_sense[row] = 'L'; + user_problem.rhs[row] = hi; + } else if (hi == inf) { + user_problem.row_sense[row] = 'G'; + user_problem.rhs[row] = lo; + } else { + user_problem.row_sense[row] = 'E'; + user_problem.rhs[row] = lo; + user_problem.range_rows.push_back(row); + user_problem.range_value.push_back(hi - lo); + } + } + user_problem.num_range_rows = user_problem.range_rows.size(); + + // Rebuild the constraint matrix and column data in place, dropping the + // slack/artificial columns (each contributes exactly one nonzero). var_types + // is filtered to the kept (structural) columns in the same pass. + csc_matrix_t& user_A = user_problem.A; + user_A.m = m; + user_A.n = new_n; + user_A.nz_max = new_nnz; + user_A.col_start.assign(new_n + 1, 0); + user_A.i.assign(new_nnz, 0); + user_A.x.assign(new_nnz, 0); + + i_t nz = 0; + i_t new_j = 0; + for (i_t j = 0; j < n; ++j) { + if (is_slack[j]) { continue; } + user_A.col_start[new_j] = nz; + for (i_t p = simplex_A.col_start[j]; p < simplex_A.col_start[j + 1]; ++p) { + user_A.i[nz] = simplex_A.i[p]; + user_A.x[nz] = simplex_A.x[p]; + ++nz; + } + user_problem.objective.push_back(simplex_problem.objective[j]); + user_problem.lower.push_back(simplex_problem.lower[j]); + user_problem.upper.push_back(simplex_problem.upper[j]); + if (!var_types.empty()) { user_problem.var_types.push_back(var_types[j]); } + ++new_j; + } + user_A.col_start[new_n] = nz; + assert(new_j == new_n); + assert(nz == new_nnz); + + user_problem.obj_scale = simplex_problem.obj_scale; + user_problem.obj_constant = simplex_problem.obj_constant; + user_problem.objective_is_integral = simplex_problem.objective_is_integral; + user_problem.objective_step = simplex_problem.objective_step; + + // LP/MIP only: no quadratic objective and no second-order cones. + user_problem.cone_var_start = 0; + user_problem.second_order_cone_dims.clear(); +} + template void convert_user_problem(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, @@ -1892,6 +2024,13 @@ template void convert_user_problem( std::vector& new_slacks, dualize_info_t& dualize_info); +template void convert_simplex_problem( + const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const std::vector& new_slacks, + user_problem_t& user_problem); + template void convert_user_lp_with_guess( const user_problem_t& user_problem, const lp_solution_t& initial_solution, diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 5dd1df0304..3607063327 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -58,13 +58,20 @@ struct lp_problem_t { f_t max_abs_obj_coeff = 0; f_t min_abs_obj_coeff = 0; - void write_mps(const std::string& path) const + // Dump the problem to an MPS file. When `var_types` is provided (length up to num_cols; any + // trailing columns such as slacks are treated as continuous), integer columns are wrapped in + // 'MARKER' 'INTORG'/'INTEND' pairs so the file round-trips as a MIP. With an empty `var_types` + // (the default) the output is the pure continuous LP as before. + void write_mps(const std::string& path, const std::vector& var_types = {}) const { std::ofstream mps_file(path); if (!mps_file.is_open()) { printf("Failed to open file %s\n", path.c_str()); return; } + auto is_integer_col = [&](i_t j) { + return j < static_cast(var_types.size()) && var_types[j] != variable_type_t::CONTINUOUS; + }; mps_file << std::setprecision(std::numeric_limits::max_digits10); mps_file << "NAME " << "cuopt_lp_problem_t" << "\n"; mps_file << "ROWS\n"; @@ -73,7 +80,17 @@ struct lp_problem_t { mps_file << " E R" << i << "\n"; } mps_file << "COLUMNS\n"; + bool in_integer_block = false; + i_t marker_id = 0; for (i_t j = 0; j < num_cols; j++) { + const bool integer_col = is_integer_col(j); + if (integer_col && !in_integer_block) { + mps_file << " MARKER" << marker_id++ << " 'MARKER' 'INTORG'\n"; + in_integer_block = true; + } else if (!integer_col && in_integer_block) { + mps_file << " MARKER" << marker_id++ << " 'MARKER' 'INTEND'\n"; + in_integer_block = false; + } const i_t col_start = A.col_start[j]; const i_t col_end = A.col_start[j + 1]; mps_file << " " << "C" << j << " OBJ " << objective[j] << "\n"; @@ -85,6 +102,10 @@ struct lp_problem_t { mps_file << " " << col_name << " " << row_name << " " << x << "\n"; } } + if (in_integer_block) { + mps_file << " MARKER" << marker_id++ << " 'MARKER' 'INTEND'\n"; + in_integer_block = false; + } mps_file << "RHS\n"; for (i_t i = 0; i < num_rows; i++) { mps_file << " RHS1 R" << i << " " << rhs[i] << "\n"; @@ -106,6 +127,10 @@ struct lp_problem_t { } if (ub != std::numeric_limits::infinity()) { mps_file << " UP BOUND1 " << col_name << " " << ub << "\n"; + } else if (is_integer_col(j)) { + // An integer column inside an INTORG/INTEND block with no explicit upper bound defaults + // to [0, 1] in most MPS readers (HiGHS, SCIP). Emit PL to keep it unbounded above. + mps_file << " PL BOUND1 " << col_name << "\n"; } } } @@ -205,6 +230,18 @@ void convert_user_problem(const user_problem_t& user_problem, std::vector& new_slacks, dualize_info_t& dualize_info); +// Inverse of convert_user_problem (LP/MIP only): reconstruct a range-form +// user_problem from a simplex standard-form lp_problem, dropping the +// slack/artificial columns recorded in new_slacks. var_types spans the full +// simplex problem (structural + slack columns) and is filtered to the +// structural columns. +template +void convert_simplex_problem(const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const std::vector& new_slacks, + user_problem_t& user_problem); + template void convert_user_problem_with_guess(const user_problem_t& user_problem, const std::vector& guess, diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index d5834f323f..eed239bcbd 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -26,6 +26,43 @@ struct benchmark_info_t; namespace cuopt::mathematical_optimization::simplex { +struct submip_settings_t { + // Enable or disable (recursive) RINS + int enable_rins = -1; + + // Base for calculating the target fix rate for the RINS neighbourhood. Actual target value is + // determined automatically according to the success and infeasible rate. + double base_target_fixrate = 0.6; + + // Minimum fix rate for accepting the RINS neighbourhood. + double min_fixrate = 0.25; + + // Hard cap for the minimum fix rate for solving a sub-MIP. + double min_fixrate_cap = 0.1; + + // MIP gap for the sub-MIP (unless the MIP gap from the B&B is lower) + double target_mip_gap = 0.01; + + // The base node limit for the sub-MIP + int node_limit_base = 200; + + // The current level in the recursion + int level = 0; + + // Maximum recursion level + int max_level = 10; + + // Presolve sub-MIP with Papilo before solving it + bool presolve = true; + + // Limit the number of simplex iterations spent in the submip. Set as a factor of the total + // number of simplex iteration from the parent B&B. + double iteration_limit_ratio = 0.8; + + // Run CPU FJ over the sub-MIP + bool enable_cpufj = true; +}; + template struct simplex_solver_settings_t { public: @@ -34,6 +71,7 @@ struct simplex_solver_settings_t { node_limit(std::numeric_limits::max()), time_limit(std::numeric_limits::infinity()), work_limit(std::numeric_limits::infinity()), + bnb_iteration_limit(std::numeric_limits::max()), absolute_mip_gap_tol(0.0), relative_mip_gap_tol(1e-3), integer_tol(1e-5), @@ -103,7 +141,7 @@ struct simplex_solver_settings_t { bnb_max_steal_attempts(-1), reliability_branching(-1), inside_mip(0), - sub_mip(0), + inside_submip(0), solution_callback(nullptr), heuristic_preemption_callback(nullptr), dual_simplex_objective_callback(nullptr), @@ -119,15 +157,16 @@ struct simplex_solver_settings_t { i_t node_limit; f_t time_limit; f_t work_limit; - f_t absolute_mip_gap_tol; // Tolerance on mip gap to declare optimal - f_t relative_mip_gap_tol; // Tolerance on mip gap to declare optimal - f_t integer_tol; // Tolerance on integralitiy violation - f_t primal_tol; // Absolute primal infeasibility tolerance - f_t dual_tol; // Absolute dual infeasibility tolerance - f_t pivot_tol; // Simplex pivot tolerance - f_t tight_tol; // A tight tolerance used to check for infeasibility - f_t fixed_tol; // If l <= x <= u with u - l < fixed_tol a variable is consider fixed - f_t zero_tol; // Values below this tolerance are considered numerically zero + int64_t bnb_iteration_limit; // Limit of the total number of simplex iterations in B&B + f_t absolute_mip_gap_tol; // Tolerance on mip gap to declare optimal + f_t relative_mip_gap_tol; // Tolerance on mip gap to declare optimal + f_t integer_tol; // Tolerance on integralitiy violation + f_t primal_tol; // Absolute primal infeasibility tolerance + f_t dual_tol; // Absolute dual infeasibility tolerance + f_t pivot_tol; // Simplex pivot tolerance + f_t tight_tol; // A tight tolerance used to check for infeasibility + f_t fixed_tol; // If l <= x <= u with u - l < fixed_tol a variable is consider fixed + f_t zero_tol; // Values below this tolerance are considered numerically zero f_t barrier_relative_feasibility_tol; // Relative feasibility tolerance for barrier method f_t barrier_relative_optimality_tol; // Relative optimality tolerance for barrier method f_t @@ -216,10 +255,11 @@ struct simplex_solver_settings_t { i_t reliability_branching; i_t inside_mip; // 0 if outside MIP, 1 if inside MIP at root node, 2 if inside MIP at leaf node - i_t sub_mip; // 0 if in regular MIP solve, 1 if in sub-MIP solve + i_t inside_submip; // 0 if in regular MIP solve, 1 if in sub-MIP solve + + submip_settings_t submip_settings; std::function&, f_t)> solution_callback; - std::function&, f_t)> node_processed_callback; std::function heuristic_preemption_callback; std::function&, std::vector&, f_t)> set_simplex_solution_callback; std::function dual_simplex_objective_callback; // Called with current dual obj @@ -227,7 +267,7 @@ struct simplex_solver_settings_t { std::atomic* concurrent_halt; // if nullptr ignored, if !nullptr, 0 if solver should // continue, 1 if solver should halt // Optional non-owning pointer to run-level benchmark stats. - cuopt::mathematical_optimization::benchmark_info_t* benchmark_info_ptr = nullptr; + benchmark_info_t* benchmark_info_ptr = nullptr; }; } // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 9d5ef320f2..7705465512 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -26,7 +26,6 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/diversity/diversity_manager.cu ${CMAKE_CURRENT_SOURCE_DIR}/diversity/multi_armed_bandit.cu ${CMAKE_CURRENT_SOURCE_DIR}/diversity/population.cu - ${CMAKE_CURRENT_SOURCE_DIR}/diversity/lns/rins.cu ${CMAKE_CURRENT_SOURCE_DIR}/relaxed_lp/relaxed_lp.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/local_search.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/bounds_repair.cu diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index 57ba659384..00c6d39443 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -55,7 +55,6 @@ diversity_manager_t::diversity_manager_t(mip_solver_context_tn_constraints, context.problem_ptr->handle_ptr->get_stream()), ls(context, lp_optimal_solution), - rins(context, *this), timer(diversity_config.default_time_limit), bound_prop_recombiner(context, context.problem_ptr->n_variables, @@ -193,7 +192,7 @@ void diversity_manager_t::add_user_given_solutions( if (has_papilo) { if ((i_t)init_sol_assignment.size() != papilo_orig_n) { CUOPT_LOG_ERROR( - "Error cannot add the provided initial solution! Initial solution %zu has %zu vars, " + "add the provided initial solution! Initial solution %zu has %zu vars, " "expected %d; skipping", sol_idx, init_sol_assignment.size(), @@ -654,8 +653,6 @@ solution_t diversity_manager_t::run_solver() return sol; } - if (omp_get_num_threads() > CUOPT_MIP_RINS_REQUIRED_THREAD_COUNT) { rins.enable(); } - generate_solution(timer.remaining_time(), false); if (timer.check_time_limit()) { population.add_external_solutions_to_population(); diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cuh b/cpp/src/mip_heuristics/diversity/diversity_manager.cuh index 51cb7e6c44..7655a58eb8 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cuh +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cuh @@ -21,7 +21,6 @@ #include #include -#include #include #include #include @@ -102,8 +101,6 @@ class diversity_manager_t { // atomic for signalling pdlp to stop std::atomic global_concurrent_halt{0}; - rins_t rins; - bool run_only_ls_recombiner{false}; bool run_only_bp_recombiner{false}; bool run_only_fp_recombiner{false}; diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu deleted file mode 100644 index 7232e517b4..0000000000 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ /dev/null @@ -1,343 +0,0 @@ -/* - * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. - * SPDX-License-Identifier: Apache-2.0 - */ - -#include - -#include -#include -#include -#include - -#include -#include -#include - -namespace cuopt::mathematical_optimization::mip { -template -rins_t::rins_t(mip_solver_context_t& context_, - diversity_manager_t& dm_, - rins_settings_t settings_) - : context(context_), problem_ptr(context.problem_ptr), dm(dm_), settings(settings_) -{ - fixrate = context.settings.heuristic_params.rins_fix_rate; - time_limit = context.settings.heuristic_params.rins_time_limit; -} - -template -void rins_t::new_best_incumbent_callback(const std::vector& solution) -{ - node_count_at_last_improvement = node_count.load(); -} - -template -void rins_t::node_callback(const std::vector& solution, f_t objective) -{ - if (!enabled) return; - node_count++; - - if (node_count - node_count_at_last_improvement < settings.nodes_after_later_improvement) return; - if (node_count - node_count_at_last_rins > settings.node_freq) { - // opportunistic early test w/ atomic to avoid having to take the lock - if (!launch_new_task.exchange(false)) return; - - bool population_ready = false; - { - std::lock_guard pop_lock(dm.population.write_mutex); - population_ready = dm.population.current_size() > 0 && dm.population.is_feasible(); - } - - if (population_ready) { - lp_optimal_solution = solution; - - CUOPT_LOG_DEBUG("Launching RINS task"); -#pragma omp task default(none) priority(CUOPT_MEDIUM_TASK_PRIORITY) - run_rins(); - } else { - launch_new_task = true; - } - } -} - -template -void rins_t::enable() -{ - seed = cuopt::seed_generator::get_seed(); - problem_ptr->handle_ptr->sync_stream(); - problem_copy = std::make_unique>(*problem_ptr, &rins_handle); - enabled = true; -} - -template -void rins_t::run_rins() -{ - raft::common::nvtx::range fun_scope("Running RINS"); - scope_guard guard([this]() { this->launch_new_task = true; }); - - RAFT_CUDA_TRY(cudaSetDevice(context.handle_ptr->get_device())); - cuopt_assert(lp_optimal_solution.size() == problem_copy->n_variables, "Assignment size mismatch"); - cuopt_assert(problem_copy->handle_ptr == &rins_handle, "Handle mismatch"); - // Do not make assertions based on problem_ptr. The original problem may have been modified within - // the FP loop relaxing integers cuopt_assert(problem_copy->n_variables == - // problem_ptr->n_variables, "Problem size mismatch"); cuopt_assert(problem_copy->n_constraints == - // problem_ptr->n_constraints, "Problem size mismatch"); cuopt_assert(problem_copy->n_integer_vars - // == problem_ptr->n_integer_vars, - // "Problem size mismatch"); - // cuopt_assert(problem_copy->n_binary_vars == problem_ptr->n_binary_vars, "Problem size - // mismatch"); - - solution_t best_sol(*problem_copy); - rins_handle.sync_stream(); - // copy the best from the population into a solution_t in the RINS stream - { - std::lock_guard lock(dm.population.write_mutex); - if (!dm.population.is_feasible()) return; - cuopt_assert(dm.population.current_size() > 0, "No solutions in population"); - auto& best_feasible_ref = dm.population.best_feasible(); - cuopt_assert(best_feasible_ref.assignment.size() == best_sol.assignment.size(), - "Assignment size mismatch"); - cuopt_assert(best_feasible_ref.get_feasible(), "Best feasible is not feasible"); - expand_device_copy(best_sol.assignment, best_feasible_ref.assignment, rins_handle.get_stream()); - best_sol.handle_ptr = &rins_handle; - best_sol.problem_ptr = problem_copy.get(); - best_sol.compute_feasibility(); - } - cuopt_assert(best_sol.handle_ptr == &rins_handle, "Handle mismatch"); - - cuopt_assert(best_sol.get_feasible(), "Best solution is not feasible"); - if (!best_sol.get_feasible()) { return; } - - i_t sol_size_before_rins = best_sol.assignment.size(); - auto lp_opt_device = cuopt::device_copy(this->lp_optimal_solution, rins_handle.get_stream()); - cuopt_assert(lp_opt_device.size() == problem_copy->n_variables, "Assignment size mismatch"); - cuopt_assert(best_sol.assignment.size() == problem_copy->n_variables, "Assignment size mismatch"); - - rmm::device_uvector vars_to_fix(problem_copy->n_integer_vars, rins_handle.get_stream()); - auto end = - thrust::copy_if(rins_handle.get_thrust_policy(), - problem_copy->integer_indices.begin(), - problem_copy->integer_indices.end(), - vars_to_fix.begin(), - [lpopt = lp_opt_device.data(), - pb = problem_copy->view(), - incumbent = best_sol.assignment.data()] __device__(i_t var_idx) -> bool { - return pb.integer_equal(lpopt[var_idx], incumbent[var_idx]); - }); - vars_to_fix.resize(end - vars_to_fix.begin(), rins_handle.get_stream()); - f_t fractional_ratio = (f_t)(vars_to_fix.size()) / (f_t)problem_copy->n_integer_vars; - - // abort if the fractional ratio is too low - if (fractional_ratio < settings.min_fractional_ratio) { - CUOPT_LOG_TRACE("RINS fractional ratio too low, aborting"); - return; - } - - thrust::default_random_engine g(seed + node_count); - - // shuffle fixing order - thrust::shuffle(rins_handle.get_thrust_policy(), vars_to_fix.begin(), vars_to_fix.end(), g); - - // fix n first according to fractional ratio - f_t rins_ratio = fixrate; - i_t n_to_fix = std::max((int)(vars_to_fix.size() * rins_ratio), 0); - vars_to_fix.resize(n_to_fix, rins_handle.get_stream()); - thrust::sort(rins_handle.get_thrust_policy(), vars_to_fix.begin(), vars_to_fix.end()); - - cuopt_assert(thrust::all_of(rins_handle.get_thrust_policy(), - vars_to_fix.begin(), - vars_to_fix.end(), - [pb = problem_copy->view()] __device__(i_t var_idx) -> bool { - return pb.is_integer_var(var_idx); - }), - "All variables to fix must be integer variables"); - - if (n_to_fix == 0) { - CUOPT_LOG_DEBUG("RINS no variables to fix"); - return; - } - - total_calls++; - node_count_at_last_rins = node_count.load(); - time_limit = std::min(time_limit, static_cast(dm.timer.remaining_time())); - CUOPT_LOG_DEBUG("Running RINS on solution with objective %g, fixing %d/%d", - best_sol.get_user_objective(), - vars_to_fix.size(), - problem_copy->n_integer_vars); - CUOPT_LOG_DEBUG("RINS fixrate %g time limit %g", fixrate, time_limit); - CUOPT_LOG_DEBUG("RINS fractional ratio %g%%", fractional_ratio * 100); - - f_t prev_obj = best_sol.get_user_objective(); - - auto [fixed_problem, fixed_assignment, variable_map] = best_sol.fix_variables(vars_to_fix); - CUOPT_LOG_DEBUG( - "new var count %d var_count %d", fixed_problem.n_variables, problem_copy->n_integer_vars); - - // should probably just do an spmv to get the objective instead. ugly mess of copies - solution_t best_sol_fixed_space(fixed_problem); - cuopt_assert(best_sol_fixed_space.handle_ptr == &rins_handle, "Handle mismatch"); - best_sol_fixed_space.copy_new_assignment( - cuopt::host_copy(fixed_assignment, rins_handle.get_stream())); - best_sol_fixed_space.compute_feasibility(); - CUOPT_LOG_DEBUG("RINS best sol fixed space objective %g", - best_sol_fixed_space.get_user_objective()); - - if (settings.objective_cut) { - f_t objective_cut = - best_sol_fixed_space.get_objective() - - std::max(std::abs(0.001 * best_sol_fixed_space.get_objective()), OBJECTIVE_EPSILON); - fixed_problem.add_cutting_plane_at_objective(objective_cut); - } - - fixed_problem.presolve_data.reset_additional_vars(fixed_problem, &rins_handle); - fixed_problem.presolve_data.initialize_var_mapping(fixed_problem, &rins_handle); - trivial_presolve(fixed_problem); - fixed_problem.check_problem_representation(true); - - std::vector> rins_solution_queue; - - mip_solver_context_t fj_context(&rins_handle, &fixed_problem, context.settings); - fj_t fj(fj_context); - solution_t fj_solution(fixed_problem); - fj_solution.copy_new_assignment(cuopt::host_copy(fixed_assignment, rins_handle.get_stream())); - std::vector default_weights(fixed_problem.n_constraints, 1.); - - std::unique_ptr> fj_cpu = - fj.create_cpu_climber(fj_solution, - default_weights, - default_weights, - 0., - context.preempt_heuristic_solver_, - fj_settings_t{}, - true); - fj_cpu->log_prefix = "[RINS] "; - - CUOPT_LOG_DEBUG("Launching CPUFJ (RINS) task"); -#pragma omp task shared(fj_cpu) firstprivate(time_limit) \ - priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) - cpufj_solve(fj_cpu.get(), time_limit); - - f_t lower_bound = context.branch_and_bound_ptr ? context.branch_and_bound_ptr->get_lower_bound() - : -std::numeric_limits::infinity(); - f_t current_mip_gap = compute_rel_mip_gap(prev_obj, lower_bound); - - // run sub-mip - namespace simplex = cuopt::mathematical_optimization::simplex; - simplex::user_problem_t branch_and_bound_problem(&rins_handle); - simplex::simplex_solver_settings_t branch_and_bound_settings; - simplex::mip_solution_t branch_and_bound_solution(1); - mip::mip_status_t branch_and_bound_status = mip::mip_status_t::UNSET; - fixed_problem.get_host_user_problem(branch_and_bound_problem); - branch_and_bound_solution.resize(branch_and_bound_problem.num_cols); - // Fill in the settings for branch and bound - branch_and_bound_settings.time_limit = time_limit; - // branch_and_bound_settings.node_limit = 5000 + node_count / 100; // try harder as time goes - // on - branch_and_bound_settings.print_presolve_stats = false; - branch_and_bound_settings.absolute_mip_gap_tol = context.settings.tolerances.absolute_mip_gap; - branch_and_bound_settings.relative_mip_gap_tol = - std::min(current_mip_gap, (f_t)settings.target_mip_gap); - branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; - branch_and_bound_settings.num_threads = 1; - branch_and_bound_settings.reliability_branching = 0; - branch_and_bound_settings.max_cut_passes = 0; - branch_and_bound_settings.clique_cuts = 0; - branch_and_bound_settings.zero_half_cuts = 0; - branch_and_bound_settings.sub_mip = 1; - branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; - branch_and_bound_settings.log.log = false; - branch_and_bound_settings.log.log_prefix = "[RINS] "; - branch_and_bound_settings.solution_callback = [&rins_solution_queue](std::vector& solution, - f_t objective) { - rins_solution_queue.push_back(solution); - }; - mip::probing_implied_bound_t empty_probing(branch_and_bound_problem.num_cols); - mip::branch_and_bound_t branch_and_bound( - branch_and_bound_problem, branch_and_bound_settings, tic(), empty_probing); - branch_and_bound.set_initial_guess(cuopt::host_copy(fixed_assignment, rins_handle.get_stream())); - branch_and_bound_status = branch_and_bound.solve(branch_and_bound_solution); - - if (!std::isnan(branch_and_bound_solution.objective)) { - CUOPT_LOG_DEBUG("RINS submip solution found. Objective %.16e. Status %d", - branch_and_bound_solution.objective, - int(branch_and_bound_status)); - // RINS submip may have just proved the initial guess is the optimal, therefore the queue might - // be empty in that case - } - if (branch_and_bound_status == mip::mip_status_t::OPTIMAL) { - CUOPT_LOG_DEBUG("RINS submip optimal"); - // do goldilocks update - fixrate = std::max(fixrate - f_t(0.05), static_cast(settings.min_fixrate)); - time_limit = std::max(time_limit - f_t(2), static_cast(settings.min_time_limit)); - } else if (branch_and_bound_status == mip::mip_status_t::TIME_LIMIT) { - CUOPT_LOG_DEBUG("RINS submip time limit"); - // do goldilocks update - fixrate = std::min(fixrate + f_t(0.05), static_cast(settings.max_fixrate)); - time_limit = std::min(time_limit + f_t(2), - static_cast(context.settings.heuristic_params.rins_max_time_limit)); - } else if (branch_and_bound_status == mip::mip_status_t::INFEASIBLE) { - CUOPT_LOG_DEBUG("RINS submip infeasible"); - // do goldilocks update, decreasing fixrate - fixrate = std::max(fixrate - f_t(0.05), static_cast(settings.min_fixrate)); - } else { - CUOPT_LOG_DEBUG("RINS solution not found"); - // do goldilocks update - fixrate = std::min(fixrate + f_t(0.05), static_cast(settings.max_fixrate)); - time_limit = std::min(time_limit + f_t(2), - static_cast(context.settings.heuristic_params.rins_max_time_limit)); - } - -#pragma omp taskwait // Wait for the CPU FJ (RINS) to finish - CUOPT_LOG_DEBUG("CPUFJ (RINS) task was stopped"); - - CUOPT_LOG_DEBUG("RINS FJ ran for %d iterations", fj_cpu->iterations); - if (fj_cpu->feasible_found) { - CUOPT_LOG_DEBUG("RINS FJ solution found. Objective %.16e", fj_cpu->h_best_objective); - rins_solution_queue.push_back(fj_cpu->h_best_assignment); - } - // Thread will be automatically terminated and joined by destructor - - bool improvement_found = false; - for (auto& fixed_sol : rins_solution_queue) { - cuopt_assert(fixed_assignment.size() == fixed_sol.size(), "Assignment size mismatch"); - rmm::device_uvector post_processed_solution(fixed_sol.size(), rins_handle.get_stream()); - raft::copy( - post_processed_solution.data(), fixed_sol.data(), fixed_sol.size(), rins_handle.get_stream()); - fixed_problem.post_process_assignment(post_processed_solution, false); - cuopt_assert(post_processed_solution.size() == fixed_assignment.size(), - "Assignment size mismatch"); - rins_handle.sync_stream(); - - rmm::device_uvector unfixed_assignment(post_processed_solution.size(), - rins_handle.get_stream()); - raft::copy(unfixed_assignment.data(), - post_processed_solution.data(), - post_processed_solution.size(), - rins_handle.get_stream()); - best_sol.unfix_variables(unfixed_assignment, variable_map); - best_sol.compute_feasibility(); - - if (best_sol.get_feasible()) { - cuopt_assert(best_sol.test_number_all_integer(), "All must be integers after RINS"); - if (best_sol.get_user_objective() < prev_obj) { improvement_found = true; } - cuopt_assert(best_sol.assignment.size() == sol_size_before_rins, "Assignment size mismatch"); - cuopt_assert(best_sol.assignment.size() == problem_copy->n_variables, - "Assignment size mismatch"); - dm.population.add_external_solution( - best_sol.get_host_assignment(), best_sol.get_objective(), solution_origin_t::RINS); - } - } - - if (improvement_found) total_success++; - CUOPT_LOG_DEBUG("RINS calls/successes %d/%d", total_calls, total_success); -} - -#if MIP_INSTANTIATE_FLOAT -template class rins_t; -#endif - -#if MIP_INSTANTIATE_DOUBLE -template class rins_t; -#endif - -} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cuh b/cpp/src/mip_heuristics/diversity/lns/rins.cuh deleted file mode 100644 index 6a0b239288..0000000000 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cuh +++ /dev/null @@ -1,75 +0,0 @@ -/* - * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. - * SPDX-License-Identifier: Apache-2.0 - */ - -#pragma once - -#include -#include - -#include - -#include - -namespace cuopt::mathematical_optimization::mip { - -// forward declare -template -class diversity_manager_t; - -struct rins_settings_t { - int node_freq = 100; - int nodes_after_later_improvement = 200; - double min_fixrate = 0.3; - double max_fixrate = 0.8; - double min_fractional_ratio = 0.3; - double min_time_limit = 3.; - double target_mip_gap = 0.03; - bool objective_cut = true; -}; - -template -class rins_t; - -template -class rins_t { - public: - rins_t(mip_solver_context_t& context, - diversity_manager_t& dm, - rins_settings_t settings = rins_settings_t()); - - void node_callback(const std::vector& solution, f_t objective); - void new_best_incumbent_callback(const std::vector& solution); - void enable(); - - void run_rins(); - - mip_solver_context_t& context; - problem_t* problem_ptr; - diversity_manager_t& dm; - rins_settings_t settings; - - // need a separate handle for RINS to operate on a separate stream and prevent graph capture - // issues - std::unique_ptr> problem_copy; - raft::handle_t rins_handle; - - std::vector lp_optimal_solution; - - f_t fixrate{0.5}; - i_t total_calls{0}; - i_t total_success{0}; - f_t time_limit{10.}; - i_t seed; - - omp_atomic_t enabled{false}; - omp_atomic_t lower_bound{0.}; - - omp_atomic_t node_count{0}; - omp_atomic_t node_count_at_last_rins{0}; - omp_atomic_t node_count_at_last_improvement{0}; - omp_atomic_t launch_new_task{true}; -}; - -} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/diversity/population.cu b/cpp/src/mip_heuristics/diversity/population.cu index b9412a5250..6fff26cbd3 100644 --- a/cpp/src/mip_heuristics/diversity/population.cu +++ b/cpp/src/mip_heuristics/diversity/population.cu @@ -304,7 +304,7 @@ void population_t::run_solution_callbacks(solution_t& sol) } CUOPT_LOG_DEBUG("Population: Found new best solution %g", sol.get_user_objective()); if (problem_ptr->branch_and_bound_callback != nullptr) { - problem_ptr->branch_and_bound_callback(sol.get_host_assignment()); + problem_ptr->branch_and_bound_callback(sol.get_host_assignment(), HEURISTICS); } for (auto callback : user_callbacks) { if (callback->get_type() == internals::base_solution_callback_type::GET_SOLUTION) { diff --git a/cpp/src/mip_heuristics/diversity/population.cuh b/cpp/src/mip_heuristics/diversity/population.cuh index 81521b95e9..593b1ddf1e 100644 --- a/cpp/src/mip_heuristics/diversity/population.cuh +++ b/cpp/src/mip_heuristics/diversity/population.cuh @@ -25,14 +25,13 @@ namespace cuopt::mathematical_optimization::mip { template class diversity_manager_t; -enum class solution_origin_t { BRANCH_AND_BOUND, CPUFJ, RINS, EXTERNAL }; +enum class solution_origin_t { BRANCH_AND_BOUND, CPUFJ, EXTERNAL }; constexpr const char* solution_origin_to_string(solution_origin_t origin) { switch (origin) { case solution_origin_t::BRANCH_AND_BOUND: return "B&B"; case solution_origin_t::CPUFJ: return "CPUFJ"; - case solution_origin_t::RINS: return "RINS"; case solution_origin_t::EXTERNAL: return "injected"; default: return "unknown"; } diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index a5b108ab86..4ed85f970e 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -112,7 +112,8 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.zero_half_cuts = 0; - branch_and_bound_settings.sub_mip = 1; + branch_and_bound_settings.inside_submip = 1; + branch_and_bound_settings.submip_settings.enable_rins = false; branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; branch_and_bound_settings.solution_callback = [this](std::vector& solution, f_t objective) { diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh b/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh deleted file mode 100644 index 6aae5a4379..0000000000 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu_fj_thread.cuh +++ /dev/null @@ -1,56 +0,0 @@ -/* 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::mathematical_optimization::mip { - -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 simplex::lp_problem_t& problem, - const std::vector& variable_types, - const std::vector& seed_assignment, - const 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::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index ad4f584e60..5534cf10cd 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -10,10 +10,10 @@ #include #include -#include "cpu_fj_thread.cuh" #include "feasibility_jump.cuh" #include "feasibility_jump_impl_common.cuh" #include "fj_cpu.cuh" +#include "fj_cpu_worker.cuh" #include @@ -1792,50 +1792,59 @@ std::unique_ptr> init_fj_cpu_standalone( } template -void fj_cpu_task_t::fj_cpu_deleter_t::operator()(fj_cpu_climber_t* ptr) const +void fj_cpu_worker_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( +void fj_cpu_worker_t::from_simplex_lp( const lp_problem_t& problem, - const std::vector& variable_types, + const std::vector& variable_types, const std::vector& seed_assignment, const 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); + auto new_climber = init_fj_cpu_from_host_lp( + problem, variable_types, seed_assignment, settings, preemption_flag, seed); + fj_cpu.reset(new_climber.release()); 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; + fj_cpu->improvement_callback = improvement_callback; } template -void run_fj_cpu_task(fj_cpu_task_t& task, f_t time_limit, double work_unit_limit) +void fj_cpu_worker_t::run_async(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); + if (!fj_cpu) 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); } template -void stop_fj_cpu_task(fj_cpu_task_t& task) +void fj_cpu_worker_t::run_sync(f_t time_limit, double work_unit_limit) { - if (task.fj_cpu) { - auto& fj_cpu = *task.fj_cpu; - fj_cpu.preemption_flag = true; - fj_cpu.halted = true; - } + if (!fj_cpu) return; + cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); + fj_cpu.reset(); +} + +template +void fj_cpu_worker_t::stop() +{ + if (!fj_cpu) return; + + fj_cpu->preemption_flag = true; + fj_cpu->halted = true; +#pragma omp taskwait depend(in : *fj_cpu) + fj_cpu.reset(); } #if MIP_INSTANTIATE_FLOAT template class fj_t; -template struct fj_cpu_task_t; +template struct fj_cpu_worker_t; template void cpufj_solve(fj_cpu_climber_t* fj_cpu, float in_time_limit, double work_unit_limit); @@ -1844,18 +1853,6 @@ template std::unique_ptr> init_fj_cpu_standalone( solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); -template std::unique_ptr> make_fj_cpu_task_from_host_lp( - const lp_problem_t& problem, - const std::vector& variable_types, - const std::vector& seed_assignment, - const 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, @@ -1867,7 +1864,7 @@ template void finalize_fj_cpu_host_initialization( #if MIP_INSTANTIATE_DOUBLE template class fj_t; -template struct fj_cpu_task_t; +template struct fj_cpu_worker_t; template void cpufj_solve(fj_cpu_climber_t* fj_cpu, double in_time_limit, double work_unit_limit); @@ -1876,18 +1873,6 @@ template std::unique_ptr> init_fj_cpu_standalone( solution_t& solution, std::atomic& preemption_flag, fj_settings_t settings); -template std::unique_ptr> make_fj_cpu_task_from_host_lp( - const lp_problem_t& problem, - const std::vector& variable_types, - const std::vector& seed_assignment, - const 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, diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index dd2512a77e..718c89615d 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -13,8 +13,8 @@ #include #include -#include #include +#include #include #include diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh new file mode 100644 index 0000000000..38a3443c06 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh @@ -0,0 +1,58 @@ +/* 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::mathematical_optimization::mip { + +template +struct fj_cpu_climber_t; + +template +struct fj_cpu_worker_t { + // Custom deleter to avoid pulling the entire fj_cpu_climber_t class here. + 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; + std::function&, double)> improvement_callback; + + // `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. + void from_simplex_lp(const simplex::lp_problem_t& problem, + const std::vector& variable_types, + const std::vector& seed_assignment, + const simplex::simplex_solver_settings_t& settings, + std::string log_prefix, + int64_t seed = -1); + + // 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()); + + // 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(), + double work_unit_limit = std::numeric_limits::infinity()); + + void stop(); +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/mip_constants.hpp b/cpp/src/mip_heuristics/mip_constants.hpp index 1d1e604d12..f3fb68343a 100644 --- a/cpp/src/mip_heuristics/mip_constants.hpp +++ b/cpp/src/mip_heuristics/mip_constants.hpp @@ -18,7 +18,6 @@ #define CUOPT_MIP_FJ_REQUIRED_THREAD_COUNT 8 #define CUOPT_MIP_EARLY_GPUFJ_REQUIRED_THREAD_COUNT 3 #define CUOPT_MIP_EARLY_CPUFJ_REQUIRED_THREAD_COUNT 2 -#define CUOPT_MIP_RINS_REQUIRED_THREAD_COUNT 4 #define CUOPT_MIP_BATCH_PDLP_REQUIRED_THREAD_COUNT 3 #define CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT 3 diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 64a25ff556..0b382cb074 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -208,6 +208,157 @@ papilo::Problem build_papilo_problem(const optimization_problem_t return problem; } +template +papilo::Problem build_papilo_problem(const simplex::user_problem_t& problem) +{ + raft::common::nvtx::range fun_scope("Build papilo problem"); + // Build a papilo problem from a (host-side) dual-simplex user_problem_t. Unlike the + // optimization_problem_t overload, all data already lives on the host and the constraint + // matrix is stored column-major (CSC), so there are no device copies and no COO step: the + // CSC matrix is converted once to CSR and handed straight to papilo's SparseStorage. + papilo::ProblemBuilder builder; + + const i_t num_cols = problem.num_cols; + const i_t num_rows = problem.num_rows; + const i_t nnz = problem.A.nnz(); + + builder.reserve(nnz, num_rows, num_cols); + + const std::vector& obj_coeffs = problem.objective; + const std::vector& var_lb = problem.lower; + const std::vector& var_ub = problem.upper; + const std::vector& var_types = problem.var_types; + const std::vector& row_sense = problem.row_sense; + const std::vector& rhs = problem.rhs; + + // Range rows carry an extra width and are listed separately; mark them so the row-bound + // derivation below matches convert_user_problem in dual_simplex/presolve.cpp. papilo + // represents ranged rows natively as two-sided lhs <= a^T x <= rhs, so (unlike the simplex + // path) we do not add slack columns. + std::vector range_of_row(num_rows, 0); + std::vector is_range_row(num_rows, false); + for (i_t k = 0; k < problem.num_range_rows; ++k) { + const i_t row = problem.range_rows[k]; + is_range_row[row] = true; + range_of_row[row] = problem.range_value[k]; + } + + // Derive two-sided row bounds [lhs, rhs] from the row sense. + std::vector h_constr_lb(num_rows); + std::vector h_constr_ub(num_rows); + for (i_t i = 0; i < num_rows; ++i) { + const f_t b = rhs[i]; + if (is_range_row[i]) { + const f_t r = std::abs(range_of_row[i]); + if (row_sense[i] == 'L') { + h_constr_lb[i] = b - r; + h_constr_ub[i] = b; + } else if (row_sense[i] == 'G') { + h_constr_lb[i] = b; + h_constr_ub[i] = b + r; + } else { // 'E' with a range becomes a two-sided row + h_constr_lb[i] = range_of_row[i] > 0 ? b : b - r; + h_constr_ub[i] = range_of_row[i] > 0 ? b + r : b; + } + } else if (row_sense[i] == 'L') { + h_constr_lb[i] = -std::numeric_limits::infinity(); + h_constr_ub[i] = b; + } else if (row_sense[i] == 'G') { + h_constr_lb[i] = b; + h_constr_ub[i] = std::numeric_limits::infinity(); + } else { // 'E' + h_constr_lb[i] = b; + h_constr_ub[i] = b; + } + } + + builder.setNumCols(num_cols); + builder.setNumRows(num_rows); + + // user_problem_t stores the objective already in minimization sense (obj_scale carries the + // original min/max direction for reporting only), so no sign flip is needed here. + builder.setObjAll(obj_coeffs); + builder.setObjOffset(problem.obj_constant); + + if (!var_lb.empty() && !var_ub.empty()) { + builder.setColLbAll(var_lb); + builder.setColUbAll(var_ub); + if (static_cast(problem.col_names.size()) == num_cols) { + builder.setColNameAll(problem.col_names); + } + } + + for (i_t j = 0; j < num_cols; ++j) { + builder.setColIntegral(j, var_types[j] != simplex::variable_type_t::CONTINUOUS); + } + + // Row bounds + infinity flags, set on the builder so build() materializes the constraint + // matrix directly. build() also sets RowFlag::kEquation where a finite lhs == rhs. + if (num_rows > 0) { + builder.setRowLhsAll(h_constr_lb); + builder.setRowRhsAll(h_constr_ub); + } + // Per-row inf flags (mirrors the optimization_problem_t overload). The zeroed lhs/rhs and the + // flags are handed to setConstraintMatrix below. + std::vector h_row_flags(num_rows); + for (i_t i = 0; i < num_rows; ++i) { + const bool lhs_inf = h_constr_lb[i] == -std::numeric_limits::infinity(); + const bool rhs_inf = h_constr_ub[i] == std::numeric_limits::infinity(); + if (lhs_inf) { + h_row_flags[i].set(papilo::RowFlag::kLhsInf); + h_constr_lb[i] = 0; + } + if (rhs_inf) { + h_row_flags[i].set(papilo::RowFlag::kRhsInf); + h_constr_ub[i] = 0; + } + } + + for (i_t j = 0; j < num_cols; ++j) { + builder.setColLbInf(j, var_lb[j] == -std::numeric_limits::infinity()); + builder.setColUbInf(j, var_ub[j] == std::numeric_limits::infinity()); + if (var_lb[j] == -std::numeric_limits::infinity()) { builder.setColLb(j, 0); } + if (var_ub[j] == std::numeric_limits::infinity()) { builder.setColUb(j, 0); } + } + + // Assemble COO entries (row, col, value) from the CSC storage and hand the matrix to papilo via + // SparseStorage with the MIP fill-in headroom, exactly like the optimization_problem_t overload. + // The default ProblemBuilder path (addColEntries) omits that headroom, which leaves papilo's + // in-place presolve in a state where DualInfer can assert on a row it reduced. + std::vector> h_entries; + h_entries.reserve(nnz); + const std::vector& col_start = problem.A.col_start; + const std::vector& row_index = problem.A.i; + const std::vector& values = problem.A.x; + for (i_t j = 0; j < num_cols; ++j) { + for (i_t p = col_start[j]; p < col_start[j + 1]; ++p) { + h_entries.push_back(std::make_tuple(row_index[p], j, values[p])); + } + } + + auto papilo_problem = builder.build(); + if (!h_entries.empty()) { + // CSC iteration is column-major, so entries are not row-sorted; let papilo sort them. + constexpr bool sorted_entries = false; + // MIP reductions like clique merging and substitution require more fillin. + const double spare_ratio = 10.0; + const int min_inter_row_space = 30; + auto csr_storage = papilo::SparseStorage( + h_entries, num_rows, num_cols, sorted_entries, spare_ratio, min_inter_row_space); + papilo_problem.setConstraintMatrix(csr_storage, h_constr_lb, h_constr_ub, h_row_flags); + + papilo::ConstraintMatrix& matrix = papilo_problem.getConstraintMatrix(); + for (i_t i = 0; i < papilo_problem.getNRows(); ++i) { + papilo::RowFlags rowFlag = matrix.getRowFlags()[i]; + if (!rowFlag.test(papilo::RowFlag::kRhsInf) && !rowFlag.test(papilo::RowFlag::kLhsInf) && + matrix.getLeftHandSides()[i] == matrix.getRightHandSides()[i]) + matrix.getRowFlags()[i].set(papilo::RowFlag::kEquation); + } + } + + return papilo_problem; +} + struct PSLPContext { Presolver* presolver = nullptr; Settings* settings = nullptr; @@ -482,6 +633,111 @@ optimization_problem_t build_optimization_problem( return op_problem; } +// Read a reduced (presolved) papilo problem back into a host-side dual-simplex user_problem_t, +// overwriting `problem` in place. +template +void build_user_problem(papilo::Problem const& papilo_problem, + simplex::user_problem_t& problem) +{ + raft::common::nvtx::range fun_scope("Build user problem"); + + const i_t reduced_rows = papilo_problem.getNRows(); + const i_t reduced_cols = papilo_problem.getNCols(); + auto const& constraint_matrix = papilo_problem.getConstraintMatrix(); + + // Objective (already minimization sense). + auto const& obj = papilo_problem.getObjective(); + problem.objective.assign(obj.coefficients.begin(), obj.coefficients.end()); + problem.obj_constant = obj.offset; + + // Column bounds and integrality. + auto const& col_lower = papilo_problem.getLowerBounds(); + auto const& col_upper = papilo_problem.getUpperBounds(); + auto const& col_flags = papilo_problem.getColFlags(); + problem.lower.resize(reduced_cols); + problem.upper.resize(reduced_cols); + problem.var_types.resize(reduced_cols); + for (i_t j = 0; j < reduced_cols; ++j) { + problem.lower[j] = col_flags[j].test(papilo::ColFlag::kLbInf) + ? -std::numeric_limits::infinity() + : col_lower[j]; + problem.upper[j] = col_flags[j].test(papilo::ColFlag::kUbInf) + ? std::numeric_limits::infinity() + : col_upper[j]; + problem.var_types[j] = col_flags[j].test(papilo::ColFlag::kIntegral) + ? simplex::variable_type_t::INTEGER + : simplex::variable_type_t::CONTINUOUS; + } + + // Row sense / rhs / ranges -- inverse of the derivation in build_papilo_problem_mip. + auto const& lhs = constraint_matrix.getLeftHandSides(); + auto const& rhs_v = constraint_matrix.getRightHandSides(); + auto const& row_flags = constraint_matrix.getRowFlags(); + problem.row_sense.assign(reduced_rows, 'E'); + problem.rhs.assign(reduced_rows, f_t{0}); + problem.range_rows.clear(); + problem.range_value.clear(); + for (i_t r = 0; r < reduced_rows; ++r) { + const bool lhs_inf = row_flags[r].test(papilo::RowFlag::kLhsInf); + const bool rhs_inf = row_flags[r].test(papilo::RowFlag::kRhsInf); + const bool eq = row_flags[r].test(papilo::RowFlag::kEquation); + if (eq || (!lhs_inf && !rhs_inf && lhs[r] == rhs_v[r])) { + problem.row_sense[r] = 'E'; + problem.rhs[r] = rhs_v[r]; + } else if (lhs_inf && !rhs_inf) { + problem.row_sense[r] = 'L'; + problem.rhs[r] = rhs_v[r]; + } else if (!lhs_inf && rhs_inf) { + problem.row_sense[r] = 'G'; + problem.rhs[r] = lhs[r]; + } else if (!lhs_inf && !rhs_inf) { + // lhs <= a^T x <= rhs : a range row. Match the convention used by get_host_user_problem and + // convert_simplex_problem (row_sense 'E' anchored at the lower bound, width in range_value). + // Storing it as 'L' left range rows out of convert_user_problem's equality_rows, which broke + // add_artifical_variables' `range_rows subset of equality_rows` invariant (j != num_cols). + problem.row_sense[r] = 'E'; + problem.rhs[r] = lhs[r]; + problem.range_rows.push_back(r); + problem.range_value.push_back(rhs_v[r] - lhs[r]); + } else { + // Free row (both sides infinite). + problem.row_sense[r] = 'L'; + problem.rhs[r] = std::numeric_limits::infinity(); + } + } + problem.num_range_rows = static_cast(problem.range_rows.size()); + + // Constraint matrix: read papilo's column-major (CSC) transpose straight into A, packing out + // the spare gaps SparseStorage leaves between columns. + const i_t reduced_nnz = constraint_matrix.getNnz(); + problem.A.resize(reduced_rows, reduced_cols, reduced_nnz); + auto const& csc = constraint_matrix.getMatrixTranspose(); + const auto* col_ranges = csc.getRowRanges(); // per-column [start, end) + const int* row_indices = csc.getColumns(); // transpose columns == original rows + const f_t* values = csc.getValues(); + i_t pos = 0; + for (i_t j = 0; j < reduced_cols; ++j) { + problem.A.col_start[j] = pos; + for (i_t p = col_ranges[j].start; p < col_ranges[j].end; ++p) { + problem.A.i[pos] = row_indices[p]; + problem.A.x[pos] = values[p]; + ++pos; + } + } + problem.A.col_start[reduced_cols] = pos; + cuopt_assert(pos == reduced_nnz, "papilo CSC nonzero count mismatch"); + + problem.num_cols = reduced_cols; + problem.num_rows = reduced_rows; +} + +template +void papilo_round_trip(simplex::user_problem_t& problem) +{ + papilo::Problem papilo_problem = build_papilo_problem(problem); + build_user_problem(papilo_problem, problem); +} + void check_presolve_status(const papilo::PresolveStatus& status) { switch (status) { @@ -538,7 +794,7 @@ void check_postsolve_status(const papilo::PostsolveStatus& status) switch (status) { case papilo::PostsolveStatus::kOk: CUOPT_LOG_DEBUG("Post-solve status: succeeded"); break; case papilo::PostsolveStatus::kFailed: - CUOPT_LOG_INFO( + CUOPT_LOG_DEBUG( "Post-solve status: Post solved solution violates constraints. This is most likely due to " "different tolerances."); break; @@ -774,6 +1030,91 @@ third_party_presolve_result_t third_party_presolve_t::apply( original_to_reduced_map_}; } +template +third_party_presolve_status_t third_party_presolve_t::apply( + simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings, + f_t time_limit, + i_t num_threads) +{ + const bool dual_postsolve = false; + presolver_ = Papilo; + // build_papilo_problem_mip keeps the objective in minimization sense (user_problem_t carries + // the direction in obj_scale), so the read-back must not flip signs either. + maximize_ = false; + + // Capture original dimensions before the problem is overwritten in place. + const i_t orig_cols = problem.num_cols; + const i_t orig_rows = problem.num_rows; + const i_t orig_nnz = problem.A.nnz(); + + papilo::Problem papilo_problem = build_papilo_problem(problem); + + settings.log.printf("Sub-MIP presolve input: %d constraints, %d variables, %d nonzeros", + papilo_problem.getNRows(), + papilo_problem.getNCols(), + papilo_problem.getConstraintMatrix().getNnz()); + + papilo::Presolve papilo_presolver; + set_presolve_methods(papilo_presolver, problem_category_t::MIP, dual_postsolve); + set_presolve_options(papilo_presolver, + problem_category_t::MIP, + settings.primal_tol, + settings.dual_tol, + time_limit, + dual_postsolve, + num_threads); + set_presolve_parameters(papilo_presolver, problem_category_t::MIP, orig_rows, orig_cols); + + // Disable papilo logs + papilo_presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet); + + auto result = papilo_presolver.apply(papilo_problem); + auto status = convert_papilo_presolve_status_to_third_party_presolve_status(result.status); + + // Infeasible / unbounded: leave `problem` untouched; the caller branches on the status. + if (result.status == papilo::PresolveStatus::kInfeasible || + result.status == papilo::PresolveStatus::kUnbndOrInfeas || + result.status == papilo::PresolveStatus::kUnbounded) { + return status; + } + + papilo_post_solve_storage_.reset(new papilo::PostsolveStorage(result.postsolve)); + + const i_t reduced_rows = papilo_problem.getNRows(); + const i_t reduced_cols = papilo_problem.getNCols(); + const i_t reduced_nnz = papilo_problem.getConstraintMatrix().getNnz(); + settings.log.printf("Sub-MIP presolve removed: %d constraints, %d variables, %d nonzeros", + orig_rows - reduced_rows, + orig_cols - reduced_cols, + orig_nnz - reduced_nnz); + + // Presolve fully solved the problem. + if (reduced_rows == 0 && reduced_cols == 0) { status = third_party_presolve_status_t::OPTIMAL; } + + // Rebuild `problem` in place from the reduced papilo problem. + build_user_problem(papilo_problem, problem); + + // Presolve changes the dimensions, so the original row/column names no longer line up with the + // reduced problem. They are not needed for the sub-MIP solve, so clear them and let downstream + // size checks skip them. + problem.col_names.clear(); + problem.row_names.clear(); + + // Column maps for postsolve (reduced -> original and its inverse). + auto const& col_map = result.postsolve.origcol_mapping; + reduced_to_original_map_.assign(col_map.begin(), col_map.end()); + original_to_reduced_map_.assign(orig_cols, -1); + for (size_t i = 0; i < reduced_to_original_map_.size(); ++i) { + auto original_idx = reduced_to_original_map_[i]; + if (original_idx >= 0 && original_idx < original_to_reduced_map_.size()) { + original_to_reduced_map_[original_idx] = i; + } + } + + return status; +} + template void third_party_presolve_t::undo(rmm::device_uvector& primal_solution, rmm::device_uvector& dual_solution, @@ -1145,11 +1486,13 @@ void papilo_postsolve_deleter::operator()(papilo::PostsolveStorage* pt #if MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT template struct papilo_postsolve_deleter; template class third_party_presolve_t; +template void papilo_round_trip(simplex::user_problem_t&); #endif #if MIP_INSTANTIATE_DOUBLE template struct papilo_postsolve_deleter; template class third_party_presolve_t; +template void papilo_round_trip(simplex::user_problem_t&); #endif } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index a5fa089d71..5cae77a480 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -12,6 +12,8 @@ #include #include +#include +#include #include @@ -58,15 +60,21 @@ class third_party_presolve_t { third_party_presolve_t(third_party_presolve_t&&) = delete; third_party_presolve_t& operator=(third_party_presolve_t&&) = delete; - third_party_presolve_result_t apply( - optimization_problem_t const& op_problem, - problem_category_t category, - cuopt::mathematical_optimization::presolver_t presolver, - bool dual_postsolve, - f_t absolute_tolerance, - f_t relative_tolerance, - double time_limit, - i_t num_cpu_threads = 0); + third_party_presolve_result_t apply(optimization_problem_t const& op_problem, + problem_category_t category, + presolver_t presolver, + bool dual_postsolve, + f_t absolute_tolerance, + f_t relative_tolerance, + double time_limit, + i_t num_cpu_threads = 0); + + // Apply the presolve on an simplex::user_problem in-place. Used in sub MIP and (in the future) + // restarts. + third_party_presolve_status_t apply(simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings, + f_t time_limit, + i_t num_threads); void undo(rmm::device_uvector& primal_solution, rmm::device_uvector& dual_solution, @@ -105,9 +113,8 @@ class third_party_presolve_t { rmm::device_uvector& reduced_costs, rmm::cuda_stream_view stream_view); - bool maximize_ = false; - cuopt::mathematical_optimization::presolver_t presolver_ = - cuopt::mathematical_optimization::presolver_t::PSLP; + bool maximize_ = false; + presolver_t presolver_ = PSLP; // PSLP settings Settings* pslp_stgs_{nullptr}; Presolver* pslp_presolver_{nullptr}; @@ -122,4 +129,8 @@ class third_party_presolve_t { std::vector original_to_reduced_map_{}; }; +// Just for testing the conversion: user_problem -> Papilo problem -> user_problem. +template +void papilo_round_trip(simplex::user_problem_t& problem); + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/problem/problem.cuh b/cpp/src/mip_heuristics/problem/problem.cuh index 5e84514cf6..543e3bba27 100644 --- a/cpp/src/mip_heuristics/problem/problem.cuh +++ b/cpp/src/mip_heuristics/problem/problem.cuh @@ -25,6 +25,7 @@ #include +#include #include #include #include @@ -254,7 +255,7 @@ class problem_t { std::shared_ptr> integer_fixed_problem = nullptr; rmm::device_uvector integer_fixed_variable_map; - std::function&)> branch_and_bound_callback; + std::function&, worker_type_t)> branch_and_bound_callback; std::function&, const std::vector&, const std::vector&, diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 5f136cbf43..81230a58e9 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -70,7 +70,6 @@ struct branch_and_bound_solution_helper_t { void solution_callback(std::vector& solution, f_t objective) { dm->population.add_external_solution(solution, objective, solution_origin_t::BRANCH_AND_BOUND); - dm->rins.new_best_incumbent_callback(solution); } void set_simplex_solution(std::vector& solution, @@ -80,11 +79,6 @@ struct branch_and_bound_solution_helper_t { dm->set_simplex_solution(solution, dual_solution, objective); } - void node_processed_callback(const std::vector& solution, f_t objective) - { - dm->rins.node_callback(solution, objective); - } - void preempt_heuristic_solver() { dm->population.preempt_heuristic_solver(); } diversity_manager_t* dm; simplex_solver_settings_t& settings_; @@ -420,12 +414,6 @@ solution_t mip_solver_t::run_solver() std::placeholders::_1, std::placeholders::_2, std::placeholders::_3); - - branch_and_bound_settings.node_processed_callback = - std::bind(&branch_and_bound_solution_helper_t::node_processed_callback, - &solution_helper, - std::placeholders::_1, - std::placeholders::_2); } // Create the branch and bound object @@ -461,7 +449,8 @@ solution_t mip_solver_t::run_solver() context.problem_ptr->branch_and_bound_callback = std::bind(&mip::branch_and_bound_t::set_solution_from_heuristics, branch_and_bound.get(), - std::placeholders::_1); + std::placeholders::_1, + std::placeholders::_2); } else if (context.settings.determinism_mode == CUOPT_MODE_DETERMINISTIC) { branch_and_bound->set_concurrent_lp_root_solve(false); // TODO once deterministic GPU heuristics are integrated diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 2e44442599..dffbcdc41a 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp @@ -14,7 +14,6 @@ #include #include #include -#include #include #include @@ -331,5 +330,4 @@ TEST(dual_simplex, dual_variable_greater_than) EXPECT_NEAR(solution.z[0], 2.0, 1e-6); EXPECT_NEAR(solution.z[1], 0.0, 1e-6); } - } // namespace cuopt::mathematical_optimization::simplex::test diff --git a/cpp/tests/linear_programming/unit_tests/presolve_test.cu b/cpp/tests/linear_programming/unit_tests/presolve_test.cu index 0c76de192a..09c5600d69 100644 --- a/cpp/tests/linear_programming/unit_tests/presolve_test.cu +++ b/cpp/tests/linear_programming/unit_tests/presolve_test.cu @@ -27,7 +27,9 @@ #include #include #include +#include #include +#include #include namespace cuopt::mathematical_optimization::test { @@ -900,6 +902,316 @@ INSTANTIATE_TEST_SUITE_P( ); // clang-format on +class papilo_problem : public ::testing::TestWithParam {}; + +// clang-format off +INSTANTIATE_TEST_SUITE_P( + papilo_presolve, + papilo_problem, + ::testing::Values( + "mip/fiball.mps", + "mip/50v-10.mps", + "mip/drayage-25-23.mps", + "mip/neos-3004026-krka.mps", + "mip/app1-1.mps", + "mip/bnatt400.mps", + "mip/decomp2.mps", + "mip/graph20-20-1rand.mps", + "mip/neos-1582420.mps", + "mip/neos-5188808-nattai.mps", + "mip/net12.mps", + "mip/n2seq36q.mps", + "mip/seymour1.mps", + "mip/neos8.mps", + "mip/CMS750_4.mps", + "mip/cbs-cta.mps", + "mip/swath3.mps", + "mip/air05.mps", + "mip/fastxgemm-n2r6s0t2.mps", + "mip/dws008-01.mps", + "mip/neos-1445765.mps", + "mip/neos-3083819-nubu.mps", + "mip/neos-5107597-kakapo.mps", + "mip/rocI-4-11.mps" + ), + [](const ::testing::TestParamInfo& info) { + std::string name = info.param; + std::replace(name.begin(), name.end(), '/', '_'); + std::replace(name.begin(), name.end(), '.', '_'); + std::replace(name.begin(), name.end(), '-', '_'); + return name; + } +); +// clang-format on + +TEST_P(papilo_problem, round_trip) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute(GetParam()); + auto mps_data_model = cuopt::mathematical_optimization::io::read_mps(path, false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps_data_model); + auto user_problem = cuopt_optimization_problem_to_user_problem(&handle_, op_problem); + simplex::user_problem_t output_problem = user_problem; + + // papilo_round_trip is a pure conversion method + // (user_problem -> papilo::Problem -> user_problem) + mip::papilo_round_trip(output_problem); + + // Sizes. + EXPECT_EQ(output_problem.num_rows, user_problem.num_rows); + EXPECT_EQ(output_problem.num_cols, user_problem.num_cols); + EXPECT_EQ(output_problem.A.m, user_problem.A.m); + EXPECT_EQ(output_problem.A.n, user_problem.A.n); + + // Objective (minimization sense), bounds, and integrality. + EXPECT_EQ(output_problem.objective, user_problem.objective); + EXPECT_DOUBLE_EQ(output_problem.obj_constant, user_problem.obj_constant); + EXPECT_EQ(output_problem.lower, user_problem.lower); + EXPECT_EQ(output_problem.upper, user_problem.upper); + EXPECT_EQ(output_problem.var_types, user_problem.var_types); + + // Row senses / rhs / ranges recover exactly (no '>=' -> '<=' fold). + EXPECT_EQ(output_problem.row_sense, user_problem.row_sense); + ASSERT_EQ(output_problem.rhs.size(), user_problem.rhs.size()); + for (int i = 0; i < user_problem.num_rows; ++i) { + EXPECT_DOUBLE_EQ(output_problem.rhs[i], user_problem.rhs[i]) << "rhs[" << i << "]"; + } + EXPECT_EQ(output_problem.num_range_rows, user_problem.num_range_rows); + EXPECT_EQ(output_problem.range_rows, user_problem.range_rows); + ASSERT_EQ(output_problem.range_value.size(), user_problem.range_value.size()); + for (std::size_t k = 0; k < user_problem.range_value.size(); ++k) { + EXPECT_DOUBLE_EQ(output_problem.range_value[k], user_problem.range_value[k]) + << "range_value[" << k << "]"; + } + + // Constraint matrix: same per-column structure and, per column, the same + // (row, value) entries. Compare order-independently since papilo's SparseStorage + // transpose need not preserve the input CSC's within-column ordering. + ASSERT_EQ(output_problem.A.col_start, user_problem.A.col_start); + for (int j = 0; j < user_problem.num_cols; ++j) { + std::vector> in_col; + std::vector> out_col; + for (int p = user_problem.A.col_start[j]; p < user_problem.A.col_start[j + 1]; ++p) { + in_col.emplace_back(user_problem.A.i[p], user_problem.A.x[p]); + } + for (int p = output_problem.A.col_start[j]; p < output_problem.A.col_start[j + 1]; ++p) { + out_col.emplace_back(output_problem.A.i[p], output_problem.A.x[p]); + } + std::sort(in_col.begin(), in_col.end()); + std::sort(out_col.begin(), out_col.end()); + ASSERT_EQ(out_col.size(), in_col.size()) << "column " << j; + for (std::size_t p = 0; p < in_col.size(); ++p) { + EXPECT_EQ(out_col[p].first, in_col[p].first) << "column " << j << " entry " << p; + EXPECT_DOUBLE_EQ(out_col[p].second, in_col[p].second) << "column " << j << " entry " << p; + } + } + + EXPECT_EQ(output_problem.col_names, user_problem.col_names); + EXPECT_EQ(output_problem.row_names, user_problem.row_names); + EXPECT_DOUBLE_EQ(output_problem.obj_scale, user_problem.obj_scale); + EXPECT_EQ(output_problem.problem_name, user_problem.problem_name); + EXPECT_EQ(output_problem.handle_ptr, user_problem.handle_ptr); + EXPECT_EQ(output_problem.Q_values, user_problem.Q_values); + EXPECT_EQ(output_problem.second_order_cone_dims, user_problem.second_order_cone_dims); + EXPECT_EQ(output_problem.cone_var_start, user_problem.cone_var_start); +} + +// Exercises the MIP presolve path: presolver_t::apply -> third_party_presolve_t::apply +// reduces a user_problem_t in place via PaPILO. ex9 is fully solved by presolve (it collapses +// to a 0x0 problem), so this also checks the OPTIMAL status and that postsolve maps the empty +// reduced solution back to a full-dimension, objective-81 assignment. +TEST(submip_presolve, ex9_fully_reduced) +{ + const raft::handle_t handle_{}; + + auto path = make_path_absolute("mip/ex9.mps"); + auto mps_data_model = cuopt::mathematical_optimization::io::read_mps(path, false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps_data_model); + + // The MIP presolve operates on the host representation. + auto user_problem = cuopt_optimization_problem_to_user_problem(&handle_, op_problem); + + const int orig_cols = user_problem.num_cols; + const auto obj_coeffs = op_problem.get_objective_coefficients_host(); + ASSERT_GT(user_problem.num_rows, 0); + ASSERT_GT(orig_cols, 0); + + simplex::simplex_solver_settings_t settings; + + mip::third_party_presolve_t presolver; + auto status = presolver.apply(user_problem, settings, 120, 8); + + // PaPILO solves ex9 entirely during presolve -> empty reduced problem. + EXPECT_EQ(status, mip::third_party_presolve_status_t::OPTIMAL); + EXPECT_EQ(user_problem.num_rows, 0); + EXPECT_EQ(user_problem.num_cols, 0); + EXPECT_EQ(user_problem.A.nnz(), 0); + + // Postsolve reconstructs the full original assignment from the (empty) reduced solution. + std::vector reduced_solution; // no reduced columns remain + std::vector full_solution; + presolver.uncrush_primal_solution(reduced_solution, full_solution); + ASSERT_EQ(static_cast(full_solution.size()), orig_cols); + + double objective = 0.0; + for (int j = 0; j < orig_cols; ++j) { + objective += obj_coeffs[j] * full_solution[j]; + } + EXPECT_NEAR(objective, 81.0, 1e-6); +} + +class simplex_problem : public ::testing::TestWithParam {}; + +// clang-format off +INSTANTIATE_TEST_SUITE_P( + dual_simplex, + simplex_problem, + ::testing::Values( + "mip/fiball.mps", + "mip/50v-10.mps", + "mip/drayage-25-23.mps", + "mip/neos-3004026-krka.mps", + "mip/app1-1.mps", + "mip/bnatt400.mps", + "mip/decomp2.mps", + "mip/graph20-20-1rand.mps", + "mip/neos-1582420.mps", + "mip/neos-5188808-nattai.mps", + "mip/net12.mps", + "mip/n2seq36q.mps", + "mip/seymour1.mps", + "mip/neos8.mps", + "mip/CMS750_4.mps", + "mip/cbs-cta.mps", + "mip/swath3.mps", + "mip/air05.mps", + "mip/fastxgemm-n2r6s0t2.mps", + "mip/dws008-01.mps", + "mip/neos-1445765.mps", + "mip/neos-3083819-nubu.mps", + "mip/neos-5107597-kakapo.mps", + "mip/rocI-4-11.mps" + ), + [](const ::testing::TestParamInfo& info) { + std::string name = info.param; + std::replace(name.begin(), name.end(), '/', '_'); + std::replace(name.begin(), name.end(), '.', '_'); + std::replace(name.begin(), name.end(), '-', '_'); + return name; + } +); + +// Round-trip a MIP through convert_user_problem (range form -> simplex standard +// form, appending one slack/artificial column per row) and +// convert_simplex_problem (the inverse: drop the slacks and recover the row +// bounds). The problem exercises every row type: '<=', '==', '>=', and a range +// row. +// +// '>=' rows are folded into '<=' rows by negating their coefficients/rhs in the +// forward pass, so the inverse recovers them as the equivalent negated '<=' row +// rather than the original '>='. The recovered problem is therefore feasibly +// equivalent but not textually identical. We assert two things: +// 1. the directly-predictable recovered fields (sizes, row_sense, rhs, range, +// objective, bounds, var_types), and +// 2. the round-trip invariant: re-running the forward conversion on the +// recovered problem reproduces the original standard-form problem exactly. +TEST_P(simplex_problem, round_trip) +{ + init_logger_t log("", true); + raft::handle_t handle{}; + + auto path = mathematical_optimization::test::make_path_absolute(GetParam()); + auto mps = cuopt::mathematical_optimization::io::read_mps(path, false); + auto op_problem = mps_data_model_to_optimization_problem(&handle, mps); + auto user_problem = cuopt_optimization_problem_to_user_problem(&handle, op_problem); + + simplex::simplex_solver_settings_t settings; + + // Forward: range-form user_problem -> simplex standard form (A x = b, l <= x <= u), + // appending one slack/artificial column per row. + simplex::lp_problem_t simplex_problem(&handle, 0, 0, 0); + std::vector new_slacks; + simplex::dualize_info_t dualize_info; + convert_user_problem(user_problem, settings, simplex_problem, new_slacks, dualize_info); + ASSERT_FALSE(dualize_info.solving_dual) + << "the default (non-barrier) path must not dualize the problem"; + + // var_types must span the full simplex problem; the appended slack/artificial + // columns are continuous. + std::vector full_var_types = user_problem.var_types; + full_var_types.resize(simplex_problem.num_cols, simplex::variable_type_t::CONTINUOUS); + + // Inverse: recover a range-form user_problem, dropping the slack columns. + simplex::user_problem_t recovered(&handle); + convert_simplex_problem(simplex_problem, full_var_types, settings, new_slacks, recovered); + + // (1) Directly-predictable recovered fields. Sizes, objective, variable bounds, + // and var_types survive untouched: the forward pass only appends slack columns + // and the inverse drops exactly those. + EXPECT_EQ(recovered.num_rows, user_problem.num_rows); + EXPECT_EQ(recovered.num_cols, user_problem.num_cols); + EXPECT_EQ(recovered.objective, user_problem.objective); + EXPECT_EQ(recovered.lower, user_problem.lower); + EXPECT_EQ(recovered.upper, user_problem.upper); + EXPECT_EQ(recovered.var_types, user_problem.var_types); + + // row_sense / rhs: '>=' rows are folded to '<=' by negating coefficients and + // rhs in the forward pass, so they come back as 'L' with a negated rhs. '<=', + // '==', and range rows ('==' with a finite range value here) are recovered + // unchanged. + std::vector expected_row_sense = user_problem.row_sense; + std::vector expected_rhs = user_problem.rhs; + for (int i = 0; i < user_problem.num_rows; ++i) { + if (user_problem.row_sense[i] == 'G') { + expected_row_sense[i] = 'L'; + expected_rhs[i] = -user_problem.rhs[i]; + } + } + EXPECT_EQ(recovered.row_sense, expected_row_sense); + for (int i = 0; i < user_problem.num_rows; ++i) { + EXPECT_NEAR(recovered.rhs[i], expected_rhs[i], 1e-9) << "rhs[" << i << "]"; + } + + // Range rows survive as the same row indices, and each range width is |r|. + EXPECT_EQ(recovered.num_range_rows, user_problem.num_range_rows); + EXPECT_EQ(recovered.range_rows, user_problem.range_rows); + ASSERT_EQ(recovered.range_value.size(), user_problem.range_value.size()); + for (std::size_t k = 0; k < recovered.range_value.size(); ++k) { + EXPECT_NEAR(recovered.range_value[k], std::abs(user_problem.range_value[k]), 1e-9) + << "range_value[" << k << "]"; + } + + // (2) Round-trip invariant: re-running the forward conversion on the recovered + // problem reproduces the original standard-form problem exactly (same slack + // columns, same matrix, same bounds). + simplex::lp_problem_t simplex_problem2(&handle, 0, 0, 0); + std::vector new_slacks2; + simplex::dualize_info_t dualize_info2; + convert_user_problem(recovered, settings, simplex_problem2, new_slacks2, dualize_info2); + + EXPECT_EQ(new_slacks2, new_slacks); + ASSERT_EQ(simplex_problem2.num_rows, simplex_problem.num_rows); + ASSERT_EQ(simplex_problem2.num_cols, simplex_problem.num_cols); + for (int j = 0; j < simplex_problem.num_cols; ++j) { + EXPECT_DOUBLE_EQ(simplex_problem2.objective[j], simplex_problem.objective[j]) + << "objective[" << j << "]"; + EXPECT_DOUBLE_EQ(simplex_problem2.lower[j], simplex_problem.lower[j]) << "lower[" << j << "]"; + EXPECT_DOUBLE_EQ(simplex_problem2.upper[j], simplex_problem.upper[j]) << "upper[" << j << "]"; + } + for (int i = 0; i < simplex_problem.num_rows; ++i) { + EXPECT_DOUBLE_EQ(simplex_problem2.rhs[i], simplex_problem.rhs[i]) << "rhs[" << i << "]"; + } + EXPECT_EQ(simplex_problem2.A.m, simplex_problem.A.m); + EXPECT_EQ(simplex_problem2.A.n, simplex_problem.A.n); + ASSERT_EQ(simplex_problem2.A.col_start, simplex_problem.A.col_start); + const int nnz = simplex_problem.A.col_start[simplex_problem.A.n]; + for (int p = 0; p < nnz; ++p) { + EXPECT_EQ(simplex_problem2.A.i[p], simplex_problem.A.i[p]) << "A.i[" << p << "]"; + EXPECT_DOUBLE_EQ(simplex_problem2.A.x[p], simplex_problem.A.x[p]) << "A.x[" << p << "]"; + } +} + } // namespace cuopt::mathematical_optimization::test CUOPT_TEST_PROGRAM_MAIN() diff --git a/cpp/tests/mip/presolve_test.cu b/cpp/tests/mip/presolve_test.cu index 9f9b65d6e0..ee29850cee 100644 --- a/cpp/tests/mip/presolve_test.cu +++ b/cpp/tests/mip/presolve_test.cu @@ -10,8 +10,11 @@ #include #include #include +#include +#include #include #include +#include #include #include #include @@ -149,5 +152,4 @@ TEST(gf2_presolve, uses_compact_constraint_indices) EXPECT_EQ(result.status, mip::third_party_presolve_status_t::REDUCED); } - } // namespace cuopt::mathematical_optimization::test diff --git a/datasets/mip/download_miplib_test_dataset.sh b/datasets/mip/download_miplib_test_dataset.sh index 2633c35b4c..177fb8008f 100755 --- a/datasets/mip/download_miplib_test_dataset.sh +++ b/datasets/mip/download_miplib_test_dataset.sh @@ -4,6 +4,7 @@ INSTANCES=( "50v-10" + "ex9" "fiball" "gen-ip054" "sct2"