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/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 5e720a645e..f6be07aaa9 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -75,6 +75,7 @@ #define CUOPT_MIP_ZERO_HALF_CUTS "mip_zero_half_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" +#define CUOPT_MIP_RINS "mip_rins" #define CUOPT_MIP_OBJECTIVE_STEP "mip_objective_step" #define CUOPT_MIP_CUT_CHANGE_THRESHOLD "mip_cut_change_threshold" #define CUOPT_MIP_CUT_MIN_ORTHOGONALITY "mip_cut_min_orthogonality" @@ -127,9 +128,19 @@ #define CUOPT_MIP_HYPER_DIVING_NODE_LIMIT "mip_hyper_diving_node_limit" #define CUOPT_MIP_HYPER_DIVING_ITERATION_LIMIT_FACTOR "mip_hyper_diving_iteration_limit_factor" #define CUOPT_MIP_HYPER_DIVING_BACKTRACK_LIMIT "mip_hyper_diving_backtrack_limit" -/* @brief Show per-strategy diving symbol in logs (true) instead of a generic 'D' */ +/* @brief Show per-strategy diving symbol in logs instead of a generic 'D' */ #define CUOPT_MIP_HYPER_DIVING_SHOW_TYPE "mip_hyper_diving_show_type" +/* @brief Recursive sub-MIP (RINS) hyper-parameters */ +#define CUOPT_MIP_HYPER_SUBMIP_BASE_TARGET_FIXRATE "mip_hyper_submip_base_target_fixrate" +#define CUOPT_MIP_HYPER_SUBMIP_MIN_FIXRATE "mip_hyper_submip_min_fixrate" +#define CUOPT_MIP_HYPER_SUBMIP_MIN_FIXRATE_CAP "mip_hyper_submip_min_fixrate_cap" +#define CUOPT_MIP_HYPER_SUBMIP_TARGET_MIP_GAP "mip_hyper_submip_target_mip_gap" +#define CUOPT_MIP_HYPER_SUBMIP_NODE_LIMIT_BASE "mip_hyper_submip_node_limit_base" +#define CUOPT_MIP_HYPER_SUBMIP_MAX_LEVEL "mip_hyper_submip_max_level" +#define CUOPT_MIP_HYPER_SUBMIP_ITERATION_LIMIT_RATIO "mip_hyper_submip_iteration_limit_ratio" +#define CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ "mip_hyper_submip_enable_cpufj" + /* @brief MIP determinism mode constants */ #define CUOPT_MODE_OPPORTUNISTIC 0 #define CUOPT_MODE_DETERMINISTIC 1 diff --git a/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp b/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp index ec08fe48dd..8ae1eda2a2 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/diving_hyper_params.hpp @@ -10,12 +10,7 @@ namespace cuopt::mathematical_optimization { /** - * @brief Tuning knobs for the dual-simplex diving heuristics used in MIP B&B. - * - * Used directly as simplex_solver_settings_t::diving_settings and copied into - * branch_and_bound_settings.diving_settings before solve. These are registered - * in the unified parameter framework via solver_settings_t and can be loaded - * from a config file with load_parameters_from_file(). + * @brief Tuning knobs for the diving heuristics for the MIP solver. */ template struct mip_diving_hyper_params_t { diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index 6d2dc8e694..8ddcdbbb8a 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -184,8 +185,8 @@ class mip_solver_settings_t { pdlp::pdlp_hyper_params_t hyper_params; mip_heuristics_hyper_params_t heuristic_params; - mip_diving_hyper_params_t diving_params; + mip_submip_hyper_params_t submip_params; private: std::vector mip_callbacks_; diff --git a/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp b/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp new file mode 100644 index 0000000000..505a3b61ba --- /dev/null +++ b/cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp @@ -0,0 +1,46 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +/** + * @brief Tuning knobs for the recursive sub-MIP. + */ +template +struct mip_submip_hyper_params_t { + // Enable or disable (recursive) RINS: -1 automatic, 0 disabled, 1 enabled + i_t rins = -1; + + // Base for calculating the target fix rate for the neighbourhood. Actual target value is + // determined automatically according to the success and infeasible rate. + f_t base_target_fixrate = 0.6; + + // Minimum fix rate for accepting the neighbourhood. + f_t min_fixrate = 0.25; + + // Hard cap for the minimum fix rate for solving a sub-MIP. + f_t min_fixrate_cap = 0.1; + + // MIP gap for the sub-MIP (unless the MIP gap from the B&B is lower) + f_t target_mip_gap = 0.01; + + // The base node limit for the sub-MIP + i_t node_limit_base = 200; + + // The current level in the recursion. This is an internal parameter and will set automatically. + i_t level = 0; + + // Maximum recursion level + i_t max_level = 10; + + // 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. + f_t iteration_limit_ratio = 0.8; + + // Run CPU FJ over the sub-MIP + bool enable_cpufj = true; +}; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e4ce4dfd7b..68c7d249f7 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -14,9 +14,10 @@ #include // benchmark_info_t #include -#include +#include #include #include +#include #include #include @@ -32,7 +33,6 @@ #include #include #include -#include #include @@ -191,6 +191,34 @@ dual_status_t convert_lp_status_to_dual_status(lp_status_t status) } } +inline char feasible_solution_symbol(heuristics_origin_t origin) +{ + switch (origin) { + case heuristics_origin_t::SUBMIP: return 'S'; + case heuristics_origin_t::HEURISTICS: return 'H'; + } + return 'U'; +} + +// 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 show_diving) +{ + if (strategy == search_strategy_t::BEST_FIRST) return 'B'; + if (!show_diving) return 'D'; + + switch (strategy) { + case search_strategy_t::BEST_FIRST: return 'B'; + case search_strategy_t::COEFFICIENT_DIVING: return 'C'; + case search_strategy_t::LINE_SEARCH_DIVING: return 'L'; + case search_strategy_t::PSEUDOCOST_DIVING: return 'P'; + case search_strategy_t::GUIDED_DIVING: return 'G'; + case search_strategy_t::FARKAS_DIVING: return 'F'; + case search_strategy_t::VECTOR_LENGTH_DIVING: return 'V'; + default: return 'U'; + } +} + template f_t sgn(f_t x) { @@ -313,6 +341,15 @@ void branch_and_bound_t::set_initial_upper_bound(f_t bound) upper_bound_ = bound; } +template +void branch_and_bound_t::set_initial_pseudocost( + const pseudo_costs_t& parent_pc, const std::vector& reduced_to_original) +{ + pc_.resize(original_lp_.num_cols); + pc_.set_initial_pseudocost(parent_pc, reduced_to_original); + has_initial_pseudocost_ = true; +} + template void branch_and_bound_t::print_table_header() { @@ -332,7 +369,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, heuristics_origin_t origin) { if (is_running_) { f_t lower_bound = get_lower_bound(); @@ -342,7 +379,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(origin), "", // nodes explored "", // nodes unexplored user_obj, @@ -383,7 +421,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 +511,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, + heuristics_origin_t origin) { mutex_original_lp_.lock(); if (solution.size() != original_problem_.num_cols) { @@ -484,11 +523,17 @@ 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(origin), + compute_user_objective(original_lp_, obj)); - if (improves_incumbent(obj)) { + if (!incumbent_.has_incumbent || obj < incumbent_.objective) { f_t primal_err; f_t bound_err; i_t num_fractional; @@ -503,11 +548,15 @@ void branch_and_bound_t::set_solution_from_heuristics(const std::vecto original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); mutex_original_lp_.unlock(); mutex_upper_.lock(); - if (is_feasible && improves_incumbent(obj)) { + if (is_feasible && (!incumbent_.has_incumbent || obj < incumbent_.objective)) { 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, origin); + success = true; + } + } else { attempt_repair = true; constexpr bool verbose = false; @@ -530,6 +579,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 +626,50 @@ 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_origin_t::HEURISTICS); + } +} + +// We need to do this dance of uncrush methods since we are working on the presolved space of +// the augmented system (structural + slack + cuts), while the `set_solution_from_heuristics` +// expects a solution on the user space. So we go from presolved space -> augmented space -> +// user space. +template +void branch_and_bound_t::set_solution_from_submip( + const std::vector& solution, + const third_party_presolve_t& presolver, + f_t fixrate, + f_t obj) +{ + std::vector leaf_sol; + presolver.uncrush_primal_solution(solution, leaf_sol); + std::vector user_sol; + mutex_original_lp_.lock(); + uncrush_primal_solution(original_problem_, original_lp_, leaf_sol, user_sol); + mutex_original_lp_.unlock(); + settings_.log.debug_format("SubMIP found a feasible solution with obj={:.4g}", obj); + bool success = set_solution_from_heuristics(user_sol, heuristics_origin_t::SUBMIP); + if (success) submip_stats_.save_success(fixrate); +} + template bool branch_and_bound_t::repair_solution(const std::vector& edge_norms, const std::vector& potential_solution, @@ -665,10 +760,10 @@ void branch_and_bound_t::repair_heuristic_solutions() if (is_feasible) { mutex_upper_.lock(); - if (improves_incumbent(repaired_obj)) { + if (!incumbent_.has_incumbent || repaired_obj < incumbent_.objective) { 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_origin_t::HEURISTICS); if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -716,6 +811,10 @@ template void branch_and_bound_t::set_final_solution(mip_solution_t& solution, f_t lower_bound) { + if (solver_status_ == mip_status_t::SUBMIP_HALT) { + settings_.log.debug("Stopping the sub-MIP solve...\n"); + } + if (solver_status_ == mip_status_t::NUMERICAL) { settings_.log.printf("Numerical issue encountered. Stopping the solver...\n"); } @@ -732,6 +831,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.debug("Simplex iteration limit reached. Stopping the solver...\n"); + } + if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } @@ -744,7 +847,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 +870,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,11 +900,12 @@ 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 @@ -811,21 +915,17 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, search_strategy_t thread_type) { bool send_solution = false; - - const bool log_diving_type = settings_.diving_settings.show_type; settings_.log.debug("%c found a feasible solution with obj=%.10e.\n", - feasible_solution_symbol(thread_type, log_diving_type), + feasible_solution_symbol(thread_type, settings_.diving_settings.show_type), compute_user_objective(original_lp_, leaf_objective)); mutex_upper_.lock(); - if (improves_incumbent(leaf_objective)) { + if (!incumbent_.has_incumbent || leaf_objective < incumbent_.objective) { incumbent_.set_incumbent_solution(leaf_objective, leaf_solution); upper_bound_ = std::min(upper_bound_.load(), leaf_objective); - report(feasible_solution_symbol(thread_type, log_diving_type), - leaf_objective, - get_lower_bound(), - leaf_depth, - 0); + + char symbol = feasible_solution_symbol(thread_type, settings_.diving_settings.show_type); + report(symbol, leaf_objective, get_lower_bound(), leaf_depth, 0); send_solution = true; } @@ -872,7 +972,7 @@ branch_variable_t branch_and_bound_t::variable_selection( std::vector& solution = worker->leaf_solution.x; switch (worker->search_strategy) { - case BEST_FIRST: + case search_strategy_t::BEST_FIRST: if (settings_.reliability_branching != 0) { branch_var = pc_.reliable_variable_selection(node_ptr, @@ -892,33 +992,37 @@ branch_variable_t branch_and_bound_t::variable_selection( return {branch_var, round_dir}; - case COEFFICIENT_DIVING: + case search_strategy_t::COEFFICIENT_DIVING: return coefficient_diving( original_lp_, fractional, solution, var_up_locks_, var_down_locks_, log); - case LINE_SEARCH_DIVING: + case search_strategy_t::LINE_SEARCH_DIVING: return line_search_diving(fractional, solution, root_relax_soln_.x, log); - case PSEUDOCOST_DIVING: + case search_strategy_t::PSEUDOCOST_DIVING: return pseudocost_diving(pc_, fractional, solution, root_relax_soln_.x, log); - case GUIDED_DIVING: + case search_strategy_t::GUIDED_DIVING: assert(incumbent_.has_incumbent); mutex_upper_.lock(); current_incumbent = incumbent_.x; mutex_upper_.unlock(); return guided_diving(pc_, fractional, solution, current_incumbent, log); - case FARKAS_DIVING: + case search_strategy_t::FARKAS_DIVING: return farkas_diving(worker->leaf_problem, fractional, solution, settings_.zero_tol, log); - case VECTOR_LENGTH_DIVING: + case search_strategy_t::VECTOR_LENGTH_DIVING: return vector_length_diving(worker->leaf_problem, fractional, solution, log); - default: - log.debug("Unknown variable selection method: %d\n", worker->search_strategy); - return {-1, branch_direction_t::NONE}; + case search_strategy_t::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}; } + + log.debug("Unknown variable selection method: %d\n", worker->search_strategy); + return {-1, branch_direction_t::NONE}; } // ============================================================================ @@ -947,7 +1051,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 @@ -1011,16 +1114,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 +1144,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 @@ -1181,10 +1273,10 @@ struct deterministic_diving_policy_t log); } - case VECTOR_LENGTH_DIVING: + case search_strategy_t::VECTOR_LENGTH_DIVING: return vector_length_diving(this->worker.leaf_problem, fractional, x, log); - case FARKAS_DIVING: + case search_strategy_t::FARKAS_DIVING: return farkas_diving( this->worker.leaf_problem, fractional, x, this->bnb.settings_.zero_tol, log); @@ -1220,6 +1312,26 @@ struct deterministic_diving_policy_t } }; +// 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. +template +void branch_and_bound_t::snap_to_lattice(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 +1389,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); + snap_to_lattice(node_ptr, leaf_obj); if (num_frac == 0) { policy.handle_integer_solution(node_ptr, leaf_obj, leaf_solution.x); @@ -1406,7 +1503,8 @@ dual_status_t branch_and_bound_t::solve_node_lp( mip_node_t* node_ptr, branch_and_bound_worker_t* worker, branch_and_bound_stats_t& stats, - logger_t& log) + logger_t& log, + i_t iter_limit) { raft::common::nvtx::range scope("BB::solve_node"); #ifdef DEBUG_BRANCHING @@ -1444,9 +1542,6 @@ dual_status_t branch_and_bound_t::solve_node_lp( } #endif - decompress_vstatus(node_ptr->packed_vstatus, worker->leaf_problem.num_cols, worker->leaf_vstatus); - assert(worker->leaf_vstatus.size() == worker->leaf_problem.num_cols); - simplex_solver_settings_t lp_settings = settings_; lp_settings.concurrent_halt = &node_concurrent_halt_; lp_settings.set_log(false); @@ -1465,17 +1560,10 @@ dual_status_t branch_and_bound_t::solve_node_lp( } else { lp_settings.cut_off = cutoff + settings_.dual_tol; } - lp_settings.inside_mip = 2; - 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; - 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; - if (lp_settings.iteration_limit <= 0) { return dual_status_t::ITERATION_LIMIT; } - } + lp_settings.inside_mip = 2; + lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + lp_settings.scale_columns = false; + lp_settings.iteration_limit = iter_limit; #ifdef LOG_NODE_SIMPLEX lp_settings.set_log(true); @@ -1502,7 +1590,8 @@ dual_status_t branch_and_bound_t::solve_node_lp( bool feasible = worker->set_lp_variable_bounds(node_ptr, settings_); dual_status_t lp_status = dual_status_t::DUAL_UNBOUNDED; worker->leaf_edge_norms = edge_norms_; - if (worker->recompute_bounds && worker->orbital_fixing && worker->search_strategy == BEST_FIRST) { + if (worker->recompute_bounds && worker->orbital_fixing && + worker->search_strategy == search_strategy_t::BEST_FIRST) { worker->orbital_fixing->reset(symmetry_, node_ptr); } @@ -1544,7 +1633,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 +1667,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); } @@ -1645,6 +1736,18 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, break; } + if (exploration_stats_.total_simplex_iters > + settings_.branch_and_bound_simplex_iteration_limit) { + solver_status_ = mip_status_t::ITERATION_LIMIT; + stack.push_front(node_ptr); + --exploration_stats_.nodes_being_solved; + break; + } + + decompress_vstatus( + node_ptr->packed_vstatus, worker->leaf_problem.num_cols, worker->leaf_vstatus); + assert(worker->leaf_vstatus.size() == worker->leaf_problem.num_cols); + 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 +1777,8 @@ 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_rins_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 +1908,37 @@ 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 (submip_halt_callback_) { + // Stops the solver if the callback returns "true". This happens when the lower bound + // in the sub-MIP solve is greater than the upper bound of the main solve (this can + // happen if one of the worker in the main solve found a better incumbent during the + // sub-MIP solve). The sub-MIP solve also stops if the status in the main solver changed + // (i.e., the gap in the main solve is sufficiently small, it reaches time/node/work limit, + // etc.) + bool stop = submip_halt_callback_(user_obj, user_lower); + if (stop) { + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::SUBMIP_HALT; + settings_.log.debug_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 +1988,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 explored the entire tree and no 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); @@ -1922,7 +2044,16 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker) } if (dive_stats.nodes_explored >= diving_node_limit) { break; } - dual_status_t lp_status = solve_node_lp(node_ptr, worker, dive_stats, log); + 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 - dive_stats.total_simplex_iters; + if (max_iter <= 0) { break; } + + decompress_vstatus( + node_ptr->packed_vstatus, worker->leaf_problem.num_cols, worker->leaf_vstatus); + assert(worker->leaf_vstatus.size() == worker->leaf_problem.num_cols); + + dual_status_t lp_status = solve_node_lp(node_ptr, worker, dive_stats, log, max_iter); ++dive_stats.nodes_explored; if (lp_status == dual_status_t::TIME_LIMIT) { @@ -1949,8 +2080,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 +2092,18 @@ 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 != search_strategy_t::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 +2136,617 @@ 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]; + 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; - 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++; - - 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_rins_worker(const std::vector& sol) +{ + if (settings_.submip_settings.rins == 0) return false; + if (!incumbent_.has_incumbent) return false; + if (rins_worker_pool_.num_idle() == 0) return false; + + diving_worker_t* worker = rins_worker_pool_.pop_idle_worker(); + if (!worker) return false; + + worker->set_active(); + worker->search_strategy = search_strategy_t::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; - return true; + 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.branch_and_bound_simplex_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.rins = + settings_.submip_settings.rins != 0 && submip_level <= settings_.submip_settings.max_level; + + submip_settings.log.debug_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.branch_and_bound_simplex_iteration_limit, + exploration_stats_.total_simplex_iters.load(), + submip_settings.relative_mip_gap_tol); + + // The `worker->leaf_problem` is directly converted to an `user_problem_t`, meaning that + // there is only equality rows (the range row vector is empty) and it contains + // structural + slacks + cuts constraints/variables. + user_problem_t submip_problem(original_problem_.handle_ptr); + simplex::convert_lp_to_user_problem(worker->leaf_problem, var_types_, settings_, submip_problem); + + 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_to_subproblem(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.debug_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); + f_t obj = 0.0; + + 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); + } + + obj += reduced_sol[j] * c; } + + set_solution_from_submip(reduced_sol, presolver, fixrate, obj); + return; + } + + submip_settings.heuristic_preemption_callback = nullptr; + submip_settings.dual_simplex_objective_callback = nullptr; + submip_settings.set_simplex_solution_callback = nullptr; + submip_settings.solution_callback = [this, &presolver, fixrate](const std::vector& solution, + f_t obj) { + this->set_solution_from_submip(solution, presolver, fixrate, obj); + }; + + submip_settings.log.debug_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); + + // 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(submip_problem, current_incumbent, presolved_incumbent); + submip_bnb.set_initial_guess(presolved_incumbent); + submip_bnb.set_initial_upper_bound(upper_bound_.load()); + submip_bnb.set_initial_pseudocost(pc_, presolver.get_reduced_to_original_map()); + + if (submip_halt_callback_) { + // Copy the halt callback to the deeper level. + submip_bnb.set_submip_halt_callback(submip_halt_callback_); + } else { + // This should only be called by the main solver. + submip_bnb.set_submip_halt_callback([this](f_t, f_t submip_lower_bound) { + f_t user_upper = compute_user_objective(this->original_lp_, this->upper_bound_.load()); + bool is_cutoff = original_lp_.obj_scale > 0 ? submip_lower_bound > user_upper + : user_upper > submip_lower_bound; + bool is_solver_running = this->solver_status_ == mip_status_t::UNSET && this->is_running_; + return is_cutoff || !is_solver_running; + }); + } + + fj_cpu_worker_t submip_fj_cpu_worker; + + if (settings_.submip_settings.enable_cpufj) { + // Launch a CPU FJ worker on the presolved sub-MIP with a fixed budget (in terms of work units) + // to run in parallel with the cut-and-branch algorithm with the goal of finding a quick + // feasible solution for the sub-MIP problem. The CPU FJ uses the current incumbent (crushed + // into the presolved space) as initial guess. The worker is automatically stop when we go out + // of the scope. + 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& solution, double work_units) { + submip_bnb.set_solution_from_cpu_fj(obj, solution, work_units); + }; + + f_t time_limit = submip_settings.time_limit; + f_t work_limit = 1.0; + submip_fj_cpu_worker.create_worker(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); + f_t submip_time = toc(start_time); + + submip_settings.log.debug_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; } - diving_worker_pool_.return_worker_to_pool(diving_worker); - return false; + if (submip_solution.has_incumbent) { + set_solution_from_submip(submip_solution.x, presolver, fixrate, submip_solution.objective); + } + + // 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; + } +} + +template +inline f_t submip_get_max_fixrate(const submip_stats_t& stats, + const mip_submip_hyper_params_t& submip_settings, + pcgenerator_t& rng) +{ + // Adaptive fix rate based on previous successes and failures. + f_t low = submip_settings.base_target_fixrate; + f_t high = submip_settings.base_target_fixrate; + + if (stats.total_infeasible > 0) { + f_t infeasible_avg_fixrate = stats.average_infeasible_fixrate(); + high = 0.9 * infeasible_avg_fixrate; + low = std::min(low, high); + } + + if (stats.total_success > 0) { + f_t 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); + } + + f_t fixrate = high > low ? rng.uniform(low, high) : low; + return fixrate; +} + +template +void get_unfixed_integer_variables(const std::vector& lower, + const std::vector& upper, + const std::vector& var_types, + f_t fixed_tol, + 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]) <= fixed_tol) { continue; } + integer_list.push_back(j); + } +} + +template +void fix_variable(i_t j, + std::vector& lower, + std::vector& upper, + std::vector& bounds_changed, + f_t fixed_val) +{ + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; +} + +template +void apply_rins_fixings(const simplex_solver_settings_t& settings, + const std::vector& node_solution, + const std::vector& fractional, + const std::vector& current_incumbent, + i_t max_var_fixed, + std::vector& lower, + std::vector& upper, + std::vector& bounds_changed, + i_t& 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]); + fix_variable(j, lower, upper, bounds_changed, fixed_val); + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + } + } +} + +template +void extend_variable_fixings(const simplex_solver_settings_t& settings, + const std::vector& obj_coeffs, + const std::vector& lower, + const std::vector& upper, + const std::vector& fractional, + const std::vector& current_sol, + const std::vector& root_solution, + 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_solution[j]; + f_t obj_coeff = obj_coeffs[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); + }); +} + +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(); } + + 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; + + // Note that this node does not have the vstatus (it was clear at the start of B&B exploration) + mip_node_t node = search_tree_.root.detach_copy(); + rins_worker->leaf_vstatus = root_vstatus_; + rins_worker->leaf_problem.lower = original_lp_.lower; + rins_worker->leaf_problem.upper = original_lp_.upper; + rins_worker->leaf_solution.x = node_solution; + rins_worker->recompute_bounds = false; + rins_worker->recompute_basis = true; + + std::vector& lower = rins_worker->leaf_problem.lower; + std::vector& upper = rins_worker->leaf_problem.upper; + std::vector& bounds_changed = rins_worker->bounds_changed; + std::vector& current_sol = rins_worker->leaf_solution.x; + std::vector fractional; + i_t num_frac = fractional_variables(settings_, current_sol, var_types_, fractional); + + std::vector current_incumbent; + mutex_upper_.lock(); + current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + + std::vector integer_list; + get_unfixed_integer_variables(lower, upper, var_types_, settings_.fixed_tol, integer_list); + + i_t num_integers = integer_list.size(); + + f_t max_fixrate = + submip_get_max_fixrate(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; + apply_rins_fixings(settings_, + node_solution, + fractional, + current_incumbent, + max_var_fixed, + lower, + upper, + bounds_changed, + num_var_fixed); + + // 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 the entire list of integer variables where the current + // LP solution matches the current incumbent. + apply_rins_fixings(settings_, + node_solution, + integer_list, + current_incumbent, + max_var_fixed, + lower, + upper, + bounds_changed, + num_var_fixed); + + // 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; + extend_variable_fixings(settings_, + rins_worker->leaf_problem.objective, + lower, + upper, + fractional, + current_sol, + root_relax_soln_.x, + candidates); + + f_t change = 0; + for (auto [dist, j, fixed_val] : candidates) { + fix_variable(j, lower, upper, bounds_changed, fixed_val); + ++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; + } + + // After fixing the variables, re-solve the LP relaxation. We use the optimal solution + // in the next iteration to find additional variable fixings. + // We continue to do this until enough variables were fixed or no variable is left to fix. + logger_t log; + log.log = false; + dual_status_t lp_status = solve_node_lp(&node, rins_worker, rins_stats, log); + + 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; + + snap_to_lattice(&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, search_strategy_t::SUBMIP); + break; + } + + rins_worker->recompute_basis = false; + } + + f_t fixrate = (f_t)num_var_fixed / num_integers; + + if (has_submip) { + // If not 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)) { + // We need to re-populate the vstatus of the node since it was previously cleared. + rins_worker->start_node = std::move(node); + rins_worker->start_node.packed_vstatus = simplex::compress_vstatus(rins_worker->leaf_vstatus); + + rins_worker->start_lower = lower; + rins_worker->start_upper = upper; + + 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.create_worker(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); + + rins_worker_pool_.return_worker_to_pool(rins_worker); } template @@ -2099,6 +2821,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 +2875,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 +3082,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 +3108,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 +3221,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("Setting initial MIP start. Objective={:+.6e}", + compute_user_objective(original_lp_, computed_obj)); } } @@ -2530,14 +3263,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 +3281,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 +3295,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 +3347,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 +3427,11 @@ 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; + 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 +3464,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 +3484,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 +3499,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.create_worker( + original_lp_, var_types_, root_relax_soln_.x, settings_, "[RootCut CPUFJ] "); settings_.log.debug("Root cut CPUFJ problem build time after pass %d: %.6f seconds\n", cut_pass, toc(root_cut_cpufj_build_start_time)); @@ -2817,6 +3520,18 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.benchmark_info_ptr->cut_generation_time_sec = cut_generation_time; } if (cut_info.has_cuts()) { + // If the incumbent is set before or during the cut passes, it may not have the correct + // dimensions as cuts add additional constraints/variables to `original_lp_`. + mutex_upper_.lock(); + if (incumbent_.has_incumbent && incumbent_.x.size() != original_lp_.num_cols) { + std::vector uncrushed_incumbent; + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, uncrushed_incumbent); + crush_primal_solution( + original_problem_, original_lp_, uncrushed_incumbent, new_slacks_, incumbent_.x); + } + + mutex_upper_.unlock(); + settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time); settings_.log.printf("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", @@ -2835,30 +3550,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.create_worker(original_lp_, + var_types_, + root_relax_soln_.x, + settings_, + "[RootCut CPUFJ] ", + root_cut_cpufj_seed); settings_.log.debug("Root cut CPUFJ final problem build time: %.6f seconds\n", toc(root_cut_cpufj_build_start_time)); f_t remaining_time = f_t(settings_.time_limit - toc(exploration_stats_.start_time)); // Reserve at least half of the remaining time for B&B exploration; cap absolute spend // at 1s so generous budgets don't grant CPUFJ more than the historical ceiling. f_t fj_time_limit = - settings_.deterministic ? remaining_time : std::min(remaining_time * 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 (!has_initial_pseudocost_) { raft::common::nvtx::range scope_sb("BB::strong_branching"); strong_branching(original_lp_, settings_, @@ -2975,8 +3688,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_); + rins_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 +3702,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 +3807,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 +3846,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 +3887,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 +3971,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; @@ -3603,7 +4320,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; @@ -3651,7 +4368,7 @@ void branch_and_bound_t::deterministic_process_worker_solutions( deterministic_current_horizon_); bool improved = false; - if (improves_incumbent(sol->objective)) { + if (!incumbent_.has_incumbent || sol->objective < incumbent_.objective) { upper_bound_ = std::min(upper_bound_.load(), sol->objective); incumbent_.set_incumbent_solution(sol->objective, sol->solution); current_upper = sol->objective; @@ -3690,10 +4407,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) { @@ -3815,14 +4532,14 @@ void branch_and_bound_t::deterministic_sort_replay_events( // Process heuristic solution at its correct work unit timestamp position f_t new_upper = std::numeric_limits::infinity(); - if (improves_incumbent(hsol.objective)) { + if (!incumbent_.has_incumbent || hsol.objective < incumbent_.objective) { upper_bound_ = std::min(upper_bound_.load(), hsol.objective); incumbent_.set_incumbent_solution(hsol.objective, hsol.solution); new_upper = hsol.objective; } if (new_upper < std::numeric_limits::infinity()) { - report_heuristic(new_upper); + report_heuristic(new_upper, heuristics_origin_t::HEURISTICS); if (settings_.solution_callback != nullptr) { std::vector original_x; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 12c93fcd91..32e3ef8fdf 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -22,6 +22,7 @@ #include #include #include + #include #include @@ -32,6 +33,8 @@ #include +#include + #include #include @@ -41,28 +44,42 @@ #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 + 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 + SUBMIP_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::SUBMIP_HALT: return "SUBMIP_HALT"; + } + return "UNKNOWN"; +} + +template +struct clique_table_t; + template -void upper_bound_callback(f_t upper_bound); +struct mip_symmetry_t; template struct nondeterministic_policy_t; @@ -86,6 +103,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_submip_halt_callback(std::function callback) + { + submip_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 +131,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, heuristics_origin_t origin); + + // 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 +150,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 set_initial_pseudocost(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; } @@ -199,18 +227,15 @@ 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. The main use of this callback is to stop the sub-MIP solve when + // the status of the main solve has changed (optimal, time/node/work limit, etc.) or + // the sub-MIP become suboptimal (lower bound is greater than the current incumbent) + std::function submip_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 +258,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 has_initial_pseudocost_{false}; // Pseudocosts pseudo_costs_t pc_; @@ -246,6 +272,9 @@ class branch_and_bound_t { // Worker pool dedicated to diving diving_worker_pool_t diving_worker_pool_; + diving_worker_pool_t rins_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 +291,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, heuristics_origin_t origin); void report(char symbol, f_t obj, f_t lower_bound, @@ -316,6 +345,8 @@ class branch_and_bound_t { // Launch a new diving worker from a given best-first worker. bool launch_diving_worker(bfs_worker_t* bfs_worker); + void snap_to_lattice(mip_node_t* node_ptr, f_t leaf_obj); + // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); @@ -332,13 +363,35 @@ 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 worker + bool launch_rins_worker(const std::vector& sol); + void set_solution_from_submip(const std::vector& solution, + const third_party_presolve_t& presolver, + f_t fixrate, + f_t obj); + + // Solve the RINS 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); + + // Get the simplex settings for solving the LP of a single node + simplex::simplex_solver_settings_t get_node_lp_settings(); // Solve the LP relaxation of a leaf node simplex::dual_status_t solve_node_lp(mip_node_t* node_ptr, branch_and_bound_worker_t* worker, branch_and_bound_stats_t& stats, - simplex::logger_t& log); + simplex::logger_t& log, + i_t iter_limit = std::numeric_limits::max()); // Apply symmetry-based bound reductions (orbital fixing and, when // settings_.symmetry == 2, lexical reduction) to the current node. diff --git a/cpp/src/branch_and_bound/constants.hpp b/cpp/src/branch_and_bound/constants.hpp index 0194e23da4..5629360a86 100644 --- a/cpp/src/branch_and_bound/constants.hpp +++ b/cpp/src/branch_and_bound/constants.hpp @@ -9,7 +9,10 @@ namespace cuopt::mathematical_optimization::mip { -constexpr int num_search_strategies = 7; +enum class heuristics_origin_t { + SUBMIP = 1, + HEURISTICS = 2, +}; // Indicate the search and variable selection algorithms used by each thread // in B&B (See [1]). @@ -19,26 +22,20 @@ 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 { +// [3] E. Danna, E. Rothberg, and C. L. Pape, “Exploring relaxation induced neighborhoods to +// improve MIP solutions,” Math. Program., vol. 102, no. 1, pp. 71–90, Jan. 2005, +// doi: 10.1007/s10107-004-0518-7. +enum class search_strategy_t : int { BEST_FIRST = 0, // Best-First + Plunging. PSEUDOCOST_DIVING = 1, // Pseudocost diving (9.2.5) LINE_SEARCH_DIVING = 2, // Line search diving (9.2.4) GUIDED_DIVING = 3, // Guided diving (9.2.3). COEFFICIENT_DIVING = 4, // Coefficient diving (9.2.1) 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 (see [3]) }; -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 }; - } // 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..2ac7afcc04 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(); } diff --git a/cpp/src/branch_and_bound/diving_heuristics.hpp b/cpp/src/branch_and_bound/diving_heuristics.hpp index f778ae38dc..224a7e4b0c 100644 --- a/cpp/src/branch_and_bound/diving_heuristics.hpp +++ b/cpp/src/branch_and_bound/diving_heuristics.hpp @@ -18,38 +18,21 @@ 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(search_strategy_t::PSEUDOCOST_DIVING); + if (settings.line_search_diving != 0) + heuristic_list.push_back(search_strategy_t::LINE_SEARCH_DIVING); + if (settings.guided_diving != 0) heuristic_list.push_back(search_strategy_t::GUIDED_DIVING); + if (settings.coefficient_diving != 0) + heuristic_list.push_back(search_strategy_t::COEFFICIENT_DIVING); + if (settings.farkas_diving != 0) heuristic_list.push_back(search_strategy_t::FARKAS_DIVING); + if (settings.vector_length_diving != 0) + heuristic_list.push_back(search_strategy_t::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..55fc1fff35 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,46 @@ 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 (!has_initial_pseudocost) { + 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 set_initial_pseudocost(const pseudo_costs_t& parent, + const std::vector& reduced_to_original) + { + has_initial_pseudocost = 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 +264,7 @@ class pseudo_costs_t { std::vector pseudo_cost_mutex_down; omp_atomic_t strong_branching_lp_iter = 0; + bool has_initial_pseudocost = false; }; template diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index 4fd278417e..fd5255f86f 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; @@ -75,6 +75,9 @@ class branch_and_bound_worker_t { std::unique_ptr> lexical_reduction; mip_symmetry_t* symmetry_ptr = nullptr; + bool recompute_basis = true; + bool recompute_bounds = true; + void ensure_orbital_fixing() { if (orbital_fixing == nullptr && symmetry_ptr != nullptr) { @@ -86,9 +89,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, @@ -96,7 +96,7 @@ class branch_and_bound_worker_t { const simplex::simplex_solver_settings_t& settings, uint64_t rng_offset = 0) : worker_id(worker_id), - search_strategy(BEST_FIRST), + search_strategy(search_strategy_t::BEST_FIRST), is_active(false), lower_bound(-std::numeric_limits::infinity()), leaf_problem(original_lp), @@ -151,11 +151,10 @@ 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; + this->search_strategy = search_strategy_t::BEST_FIRST; + max_diving_workers = 0; + active_diving_workers = 0; + next_heuristic = worker_id; } void set_inactive() { this->is_active = false; } @@ -172,53 +171,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 + search_strategy_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 +228,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 +245,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 f2005795b4..be45ffeecd 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"); @@ -4309,7 +4331,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; @@ -4644,7 +4666,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..9308e350fa 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -467,6 +467,16 @@ i_t convert_greater_to_less(const user_problem_t& user_problem, return 0; } +template +row_bounds_t get_range_bounds_from_sense(char row_sense, f_t rhs, f_t range_value) +{ + const f_t r = std::abs(range_value); + if (row_sense == 'L') { return {rhs - r, rhs}; } + if (row_sense == 'G') { return {rhs, rhs + r}; } + // 'E' with a range becomes a two-sided row + return range_value > 0 ? row_bounds_t{rhs, rhs + r} : row_bounds_t{rhs - r, rhs}; +} + template i_t convert_range_rows(const user_problem_t& user_problem, std::vector& row_sense, @@ -497,32 +507,19 @@ i_t convert_range_rows(const user_problem_t& user_problem, i_t p = problem.A.col_start[problem.num_cols]; i_t j = problem.num_cols; for (i_t k = 0; k < num_range_rows; k++) { - const i_t i = user_problem.range_rows[k]; - const f_t r = user_problem.range_value[k]; - const f_t b = problem.rhs[i]; - f_t h; - f_t u; + const i_t i = user_problem.range_rows[k]; + const f_t r = user_problem.range_value[k]; + const f_t b = problem.rhs[i]; + auto [lower, upper] = get_range_bounds_from_sense(row_sense[i], b, r); if (row_sense[i] == 'L') { - h = b - std::abs(r); - u = b; less_rows--; equal_rows++; } else if (row_sense[i] == 'G') { - h = b; - u = b + std::abs(r); greater_rows--; equal_rows++; - } else if (row_sense[i] == 'E') { - if (r > 0) { - h = b; - u = b + std::abs(r); - } else { - h = b - std::abs(r); - u = b; - } } - problem.lower[j] = h; - problem.upper[j] = u; + problem.lower[j] = lower; + problem.upper[j] = upper; problem.objective[j] = 0.0; problem.A.i[p] = i; problem.A.x[p] = -1.0; @@ -676,6 +673,50 @@ i_t add_artifical_variables(lp_problem_t& problem, return 0; } +template +void convert_lp_to_user_problem(const lp_problem_t& lp, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + 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", + lp.num_rows, + lp.num_cols, + lp.A.col_start[lp.num_cols]); + } + + const i_t m = lp.num_rows; + const i_t n = lp.num_cols; + + user_problem.handle_ptr = lp.handle_ptr; + user_problem.num_cols = n; + user_problem.num_rows = m; + user_problem.objective = lp.objective; + + user_problem.A = lp.A; + + user_problem.rhs = lp.rhs; + user_problem.lower = lp.lower; + user_problem.upper = lp.upper; + user_problem.row_sense.assign(m, 'E'); + user_problem.range_rows.clear(); + user_problem.range_value.clear(); + user_problem.num_range_rows = 0; + + user_problem.var_types = var_types; + + user_problem.obj_scale = lp.obj_scale; + user_problem.obj_constant = lp.obj_constant; + user_problem.objective_is_integral = lp.objective_is_integral; + user_problem.objective_step = lp.objective_step; + + user_problem.Q_indices.clear(); + user_problem.Q_offsets.clear(); + user_problem.Q_values.clear(); +} + template void convert_user_problem(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, @@ -1892,6 +1933,12 @@ template void convert_user_problem( std::vector& new_slacks, dualize_info_t& dualize_info); +template void convert_lp_to_user_problem( + const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + user_problem_t& user_problem); + template void convert_user_lp_with_guess( const user_problem_t& user_problem, const lp_solution_t& initial_solution, @@ -1940,6 +1987,14 @@ template void uncrush_solution(const presolve_info_t& std::vector& uncrushed_y, std::vector& uncrushed_z); +template row_bounds_t get_range_bounds_from_sense(char, double, double); + #endif +// Emitted unconditionally: third_party_presolve.cpp always instantiates its +// variant (its guard MIP_INSTANTIATE_FLOAT || PDLP_INSTANTIATE_FLOAT is always true because +// PDLP_INSTANTIATE_FLOAT == 1), so this float symbol must exist even though the rest of this +// dual_simplex TU is double-only. +template row_bounds_t get_range_bounds_from_sense(char, float, float); + } // namespace cuopt::mathematical_optimization::simplex diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 5dd1df0304..bd62977f54 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -22,6 +22,22 @@ namespace cuopt::mathematical_optimization::simplex { +// Two-sided activity bounds [lower, upper] of a constraint row. +template +struct row_bounds_t { + f_t lower; + f_t upper; +}; + +// Compute the two-sided bounds [lower, upper] of a range row from its `row_sense`, +// right-side vector `b`, and signed range value `r: +// 'L' -> [b - |r|, b] +// 'G' -> [b, b + |r|] +// 'E', r > 0 -> [b, b + |r|] +// 'E', r <= 0 -> [b - |r|, b] +template +row_bounds_t get_range_bounds_from_sense(char row_sense, f_t rhs, f_t range_value); + template struct lp_problem_t { lp_problem_t(raft::handle_t const* handle_ptr_, i_t m, i_t n, i_t nz) @@ -58,13 +74,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 < 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 +96,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 +118,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 +143,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 +246,17 @@ 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_lp_to_user_problem(const lp_problem_t& lp, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + 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..6a69cdfcd2 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -8,6 +8,7 @@ #pragma once #include +#include #include #include @@ -34,6 +35,7 @@ struct simplex_solver_settings_t { node_limit(std::numeric_limits::max()), time_limit(std::numeric_limits::infinity()), work_limit(std::numeric_limits::infinity()), + branch_and_bound_simplex_iteration_limit(std::numeric_limits::max()), absolute_mip_gap_tol(0.0), relative_mip_gap_tol(1e-3), integer_tol(1e-5), @@ -103,7 +105,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 +121,17 @@ 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 branch_and_bound_simplex_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 +220,12 @@ 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 + + // Settings for the recursive sub-MIP + mip_submip_hyper_params_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 +233,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/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 8bfd7cd8ba..9193112d71 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -117,6 +117,12 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_SEMICONTINUOUS_BIG_M, &mip_settings.semi_continuous_big_m, f_t(1.0), std::numeric_limits::infinity(), f_t(1e10), "big-M value for semi-continuous variables with no finite upper bound"}, // Diving heuristic hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_DIVING_ITERATION_LIMIT_FACTOR, &mip_settings.diving_params.iteration_limit_factor, f_t(0.0), f_t(1.0), f_t(0.05), "fraction of best-first iterations allowed per dive"}, + // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") + {CUOPT_MIP_HYPER_SUBMIP_BASE_TARGET_FIXRATE, &mip_settings.submip_params.base_target_fixrate, f_t(0.0), f_t(1.0), f_t(0.6), "base target fix rate for the RINS neighbourhood"}, + {CUOPT_MIP_HYPER_SUBMIP_MIN_FIXRATE, &mip_settings.submip_params.min_fixrate, f_t(0.0), f_t(1.0), f_t(0.25), "minimum fix rate for accepting the RINS neighbourhood"}, + {CUOPT_MIP_HYPER_SUBMIP_MIN_FIXRATE_CAP, &mip_settings.submip_params.min_fixrate_cap, f_t(0.0), f_t(1.0), f_t(0.1), "hard cap on the minimum fix rate for solving a sub-MIP"}, + {CUOPT_MIP_HYPER_SUBMIP_TARGET_MIP_GAP, &mip_settings.submip_params.target_mip_gap, f_t(0.0), f_t(1.0), f_t(0.01), "MIP gap target for the sub-MIP"}, + {CUOPT_MIP_HYPER_SUBMIP_ITERATION_LIMIT_RATIO, &mip_settings.submip_params.iteration_limit_ratio, f_t(0.0), f_t(1.0), f_t(0.8), "sub-MIP simplex-iteration limit as a factor of parent B&B iterations"}, }; // Int parameters @@ -142,6 +148,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, + {CUOPT_MIP_RINS, &mip_settings.submip_params.rins, -1, 1, -1}, {CUOPT_MIP_OBJECTIVE_STEP, &mip_settings.objective_step, 0, 1, 1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, {CUOPT_NUM_GPUS, &mip_settings.num_gpus, 1, 2, 1}, @@ -174,6 +181,9 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_HYPER_DIVING_MIN_NODE_DEPTH, &mip_settings.diving_params.min_node_depth, 0, std::numeric_limits::max(), 10, "minimum depth at which to start diving"}, {CUOPT_MIP_HYPER_DIVING_NODE_LIMIT, &mip_settings.diving_params.node_limit, 0, std::numeric_limits::max(), 500, "maximum nodes explored per dive"}, {CUOPT_MIP_HYPER_DIVING_BACKTRACK_LIMIT, &mip_settings.diving_params.backtrack_limit, 0, std::numeric_limits::max(), 5, "maximum backtracking allowed per dive"}, + // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") + {CUOPT_MIP_HYPER_SUBMIP_NODE_LIMIT_BASE, &mip_settings.submip_params.node_limit_base, 0, std::numeric_limits::max(), 200, "base node limit for the sub-MIP"}, + {CUOPT_MIP_HYPER_SUBMIP_MAX_LEVEL, &mip_settings.submip_params.max_level, 0, std::numeric_limits::max(), 10, "maximum sub-MIP recursion level"}, }; // Bool parameters @@ -194,6 +204,8 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_PROBING, &mip_settings.probing, true}, // Diving heuristic hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_DIVING_SHOW_TYPE, &mip_settings.diving_params.show_type, false, "log diving heuristic type when it finds a new incumbent"}, + // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") + {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"}, }; // String parameters string_parameters = { 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 d4fcbfa738..61c90944f4 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -56,7 +56,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, @@ -194,7 +193,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(), @@ -690,8 +689,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..553e5d6e93 100644 --- a/cpp/src/mip_heuristics/diversity/population.cu +++ b/cpp/src/mip_heuristics/diversity/population.cu @@ -304,7 +304,8 @@ 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_origin_t::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..474becef25 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.rins = 0; 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..57a6a89479 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::create_worker( 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..ff6022c4c2 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh @@ -0,0 +1,60 @@ +/* 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; + + ~fj_cpu_worker_t() { stop(); } + + // `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 create_worker(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 6e22b67c40..c3102d8935 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -35,6 +35,7 @@ #else #pragma GCC diagnostic pop #endif +#include #include #include #include @@ -209,6 +210,149 @@ 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 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]) { + auto [lower, upper] = simplex::get_range_bounds_from_sense(row_sense[i], b, range_of_row[i]); + h_constr_lb[i] = lower; + h_constr_ub[i] = upper; + } 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; @@ -483,6 +627,113 @@ 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.clear(); + problem.rhs.clear(); + problem.row_sense.reserve(reduced_rows); + problem.rhs.reserve(reduced_rows); + 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.push_back('E'); + problem.rhs.push_back(rhs_v[r]); + } else if (lhs_inf && !rhs_inf) { + problem.row_sense.push_back('L'); + problem.rhs.push_back(rhs_v[r]); + } else if (!lhs_inf && rhs_inf) { + problem.row_sense.push_back('G'); + problem.rhs.push_back(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.push_back('E'); + problem.rhs.push_back(lhs[r]); + problem.range_rows.push_back(r); + problem.range_value.push_back(rhs_v[r] - lhs[r]); + } else { + assert(false && "Papilo should remove all the free rows"); + } + } + + problem.num_range_rows = 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) { @@ -779,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_to_subproblem( + 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, @@ -869,7 +1205,7 @@ template void third_party_presolve_t::uncrush_primal_solution( const std::vector& reduced_primal, std::vector& full_primal) const { - if (presolver_ == cuopt::mathematical_optimization::presolver_t::PSLP) { + if (presolver_ == PSLP) { cuopt_expects(false, error_type_t::RuntimeError, "This code path should be never called, as this is meant for callbacks and they " @@ -928,6 +1264,45 @@ void third_party_presolve_t::crush_primal_solution( } } +template +void third_party_presolve_t::crush_primal_solution( + const simplex::user_problem_t& reduced_problem, + const std::vector& original_primal, + std::vector& reduced_primal) const +{ + cuopt_expects(presolver_ == Papilo, + error_type_t::RuntimeError, + "Primal crushing is only supported for PaPILO presolve"); + cuopt_assert(papilo_post_solve_storage_ != nullptr, "No postsolve storage available"); + std::vector unused_y, unused_z; + std::vector empty_vals; + std::vector empty_indices, empty_offsets; + crush_primal_dual_solution(original_primal, + {}, + reduced_primal, + unused_y, + {}, + unused_z, + empty_vals, + empty_indices, + empty_offsets); + + // Dual bound strengthening (e.g. DualFix kVarBoundChange which aren't emitted in primal mode) + // can tighten a bound past a value + // that was feasible in the original polytope. A simple clamp is often enough to crush a solution + // inside the tightened polytope without breaking feasibility (dualfix seems to emit purely dual + // based bounds reductions so primality is safe when clamping) + cuopt_assert(reduced_problem.num_cols == (i_t)reduced_to_original_map_.size(), + "reduced_problem does not match this presolver's reduction"); + cuopt_assert(reduced_primal.size() == reduced_problem.lower.size() && + reduced_primal.size() == reduced_problem.upper.size(), + "reduced problem must match crush output dimension"); + for (size_t j = 0; j < reduced_primal.size(); ++j) { + reduced_primal[j] = + std::clamp(reduced_primal[j], reduced_problem.lower[j], reduced_problem.upper[j]); + } +} + /** * Crush an original-space primal+dual solution into the presolved (reduced) space. * @@ -1167,11 +1542,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 b23337d25f..51193a6e12 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,22 @@ 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_to_subproblem( + 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, @@ -83,6 +92,10 @@ class third_party_presolve_t { const std::vector& original_primal, std::vector& reduced_primal) const; + void crush_primal_solution(const simplex::user_problem_t& reduced_problem, + const std::vector& original_primal, + std::vector& reduced_primal) const; + void crush_primal_dual_solution(const std::vector& x_original, const std::vector& y_original, std::vector& x_reduced, @@ -113,9 +126,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}; @@ -134,4 +146,8 @@ class third_party_presolve_t { f_t original_objective_scaling_factor_{1}; }; +// 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..bcc3f06fc2 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&, heuristics_origin_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..b045190e9f 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_; @@ -400,6 +394,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.symmetry = context.settings.symmetry; branch_and_bound_settings.diving_settings = context.settings.diving_params; + branch_and_bound_settings.submip_settings = context.settings.submip_params; // Set the branch and bound -> primal heuristics callback branch_and_bound_settings.solution_callback = @@ -420,12 +415,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 +450,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/linear_programming/unit_tests/presolve_test.cu b/cpp/tests/linear_programming/unit_tests/presolve_test.cu index 0c76de192a..75e0565d67 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,165 @@ 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_to_subproblem(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); +} + } // namespace cuopt::mathematical_optimization::test CUOPT_TEST_PROGRAM_MAIN() 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"