From aefe31a4ac23babf9a35829880c072834c0c58cd Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Thu, 25 Jun 2026 11:46:18 +0200 Subject: [PATCH 01/16] recursive RINS based on HiGHS. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/CMakeLists.txt | 10 +- cpp/src/branch_and_bound/branch_and_bound.cpp | 462 ++++++++++++++++-- cpp/src/branch_and_bound/branch_and_bound.hpp | 57 ++- cpp/src/branch_and_bound/constants.hpp | 20 +- .../branch_and_bound/diving_heuristics.hpp | 2 + cpp/src/branch_and_bound/pseudo_costs.cpp | 10 +- cpp/src/branch_and_bound/submip.hpp | 94 ++++ cpp/src/branch_and_bound/worker.hpp | 5 + cpp/src/cuts/cuts.cpp | 4 +- .../dual_simplex/simplex_solver_settings.hpp | 33 +- .../diversity/diversity_manager.cu | 2 +- cpp/src/mip_heuristics/diversity/lns/rins.cu | 2 +- .../diversity/recombiners/sub_mip.cuh | 2 +- 13 files changed, 619 insertions(+), 84 deletions(-) create mode 100644 cpp/src/branch_and_bound/submip.hpp diff --git a/cpp/src/branch_and_bound/CMakeLists.txt b/cpp/src/branch_and_bound/CMakeLists.txt index 1e40c1bbf1..47ca6b98f5 100644 --- a/cpp/src/branch_and_bound/CMakeLists.txt +++ b/cpp/src/branch_and_bound/CMakeLists.txt @@ -4,11 +4,11 @@ # cmake-format: on set(BRANCH_AND_BOUND_SRC_FILES - ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp - ) + ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp +) set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} - ${BRANCH_AND_BOUND_SRC_FILES} PARENT_SCOPE) + ${BRANCH_AND_BOUND_SRC_FILES} PARENT_SCOPE) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a17bc6716c..e5c112f09e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -197,6 +197,31 @@ dual_status_t convert_lp_status_to_dual_status(lp_status_t status) } } +lp_status_t convert_dual_status_to_lp_status(dual_status_t status) +{ + if (status == dual_status_t::OPTIMAL) { + return lp_status_t::OPTIMAL; + } else if (status == dual_status_t::DUAL_UNBOUNDED) { + return lp_status_t::INFEASIBLE; + } else if (status == dual_status_t::ITERATION_LIMIT) { + return lp_status_t::ITERATION_LIMIT; + } else if (status == dual_status_t::TIME_LIMIT) { + return lp_status_t::TIME_LIMIT; + } else if (status == dual_status_t::WORK_LIMIT) { + return lp_status_t::WORK_LIMIT; + } else if (status == dual_status_t::NUMERICAL) { + return lp_status_t::NUMERICAL_ISSUES; + } else if (status == dual_status_t::CUTOFF) { + return lp_status_t::CUTOFF; + } else if (status == dual_status_t::CONCURRENT_LIMIT) { + return lp_status_t::CONCURRENT_LIMIT; + } else if (status == dual_status_t::UNSET) { + return lp_status_t::UNSET; + } else { + return lp_status_t::NUMERICAL_ISSUES; + } +} + template f_t sgn(f_t x) { @@ -290,6 +315,42 @@ branch_and_bound_t::branch_and_bound_t( root_lp_current_lower_bound_ = -inf; } +template +branch_and_bound_t::branch_and_bound_t( + const branch_and_bound_t& other, + const simplex_solver_settings_t& solver_settings, + const std::vector& lower, + const std::vector& upper) + : original_problem_(other.original_problem_), + settings_(solver_settings), + probing_implied_bound_(other.probing_implied_bound_), + clique_table_(other.clique_table_), + symmetry_(nullptr), // Enabling symmetry in the submip lead to crashes in some instances + original_lp_(other.original_lp_), + Arow_(other.Arow_), + new_slacks_(other.new_slacks_), + var_types_(other.var_types_), + incumbent_(other.incumbent_), + root_vstatus_(other.root_vstatus_), + root_relax_soln_(1, 1), + root_crossover_soln_(1, 1), + pc_(other.pc_), + solver_status_(mip_status_t::UNSET) +{ + assert(lower.size() == original_lp_.num_cols); + assert(upper.size() == original_lp_.num_cols); + original_lp_.lower = lower; + original_lp_.upper = upper; + + exploration_stats_.start_time = other.exploration_stats_.start_time; + upper_bound_ = other.upper_bound_; + root_objective_ = std::numeric_limits::quiet_NaN(); + root_lp_current_lower_bound_ = -inf; + external_upper_bound_ = + other.external_upper_bound_ ? other.external_upper_bound_ : &other.upper_bound_; + root_warm_start_ = true; +} + template f_t branch_and_bound_t::get_lower_bound() { @@ -538,6 +599,42 @@ void branch_and_bound_t::set_solution_from_heuristics(const std::vecto } } +template +void branch_and_bound_t::set_solution_from_submip(const std::vector& solution, + f_t fixrate) +{ + mutex_original_lp_.lock(); + if (solution.size() != original_problem_.num_cols) { + settings_.log.printf( + "Solution size mismatch %ld %d\n", solution.size(), original_problem_.num_cols); + } + std::vector crushed_solution; + 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(); + + mutex_upper_.lock(); + if (improves_incumbent(obj)) { + settings_.log.printf("Recursive subMIP found solution with obj=%g", + compute_user_objective(original_lp_, obj)); + + 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 (settings_.solution_callback != nullptr) { + std::vector original_x; + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, original_x); + settings_.solution_callback(original_x, obj); + } + + submip_stats_.save_success(fixrate); + } + mutex_upper_.unlock(); +} + template void branch_and_bound_t::queue_external_solution_deterministic( const std::vector& solution, double work_unit_ts) @@ -722,8 +819,14 @@ 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.printf("Stopping submip solve...\n"); + return; + } + if (solver_status_ == mip_status_t::NUMERICAL) { settings_.log.printf("Numerical issue encountered. Stopping the solver...\n"); + return; } if (solver_status_ == mip_status_t::TIME_LIMIT) { @@ -773,7 +876,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 @@ -1226,6 +1329,21 @@ struct deterministic_diving_policy_t } }; +template +void branch_and_bound_t::apply_objective_step(mip_node_t* node_ptr, + f_t leaf_obj) +{ + if (original_lp_.objective_step.has_step()) { + f_t step = original_lp_.objective_step.step_size; + f_t bias = original_lp_.objective_step.bias; + // Round up to next value on the lattice: k * step + bias >= leaf_obj + f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); + node_ptr->lower_bound = k * step + bias; + } else if (original_lp_.objective_is_integral) { + node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); + } +} + template template std::pair branch_and_bound_t::update_tree_impl( @@ -1283,21 +1401,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); - } + apply_objective_step(node_ptr, leaf_obj); policy.on_optimal_callback(leaf_solution.x, leaf_obj); if (num_frac == 0) { @@ -1680,6 +1784,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 (worker->worker_id == 0) { launch_submip_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, @@ -1814,6 +1920,19 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t 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 (external_upper_bound_ && settings_.inside_submip && + external_upper_bound_->load() <= lower_bound) { + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::SUBMIP_HALT; + settings_.log.debug( + "Stopping the submip since the main B&B has a better incumbent. " + "(upper=%.g, external_upper=%.g, lower=%.g", + upper_bound_.load(), + external_upper_bound_->load(), + lower_bound); + 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 && @@ -1897,11 +2016,6 @@ void branch_and_bound_t::dive_with(diving_worker_t* worker) 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); @@ -2033,6 +2147,239 @@ bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* return false; } +template +void branch_and_bound_t::launch_submip_worker(const std::vector& sol) +{ + if (settings_.submip_settings.enable_rins == 0) return; + if (submip_worker_pool_.num_idle() == 0) return; + if (!incumbent_.has_incumbent) return; + + i_t nodes_since_last_submip = + exploration_stats_.nodes_explored - exploration_stats_.nodes_at_last_submip; + if (nodes_since_last_submip < settings_.submip_settings.node_freq) return; + + submip_worker_t* submip_worker = submip_worker_pool_.pop_idle_worker(); + + submip_worker->set_active(); + +#pragma omp task priority(CUOPT_MEDIUM_TASK_PRIORITY) affinity(submip_worker) \ + firstprivate(submip_worker, sol) + rins(submip_worker, sol); + + exploration_stats_.nodes_at_last_submip = exploration_stats_.nodes_explored; +} + +template +void branch_and_bound_t::rins(submip_worker_t* rins_worker, + const std::vector& node_solution) +{ + raft::common::nvtx::range scope("BB::rins_thread"); + if (rins_worker->orbital_fixing) { rins_worker->orbital_fixing->disable(); } + logger_t log; + log.log = false; + + i_t submip_level = settings_.submip_settings.level + 1; + std::string log_prefix = std::format("[SUBMIP {}] ", submip_level); + + ++submip_stats_.total_calls; + + bool need_to_solve_submip = false; + const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + + branch_and_bound_stats_t rins_stats; + + mip_node_t node = search_tree_.root.detach_copy(); + node.packed_vstatus = compress_vstatus(root_vstatus_); + rins_worker->leaf_vstatus = root_vstatus_; + rins_worker->leaf_problem.lower = original_lp_.lower; + rins_worker->leaf_problem.upper = original_lp_.upper; + rins_worker->leaf_solution.x = node_solution; + + std::vector& lower = rins_worker->leaf_problem.lower; + std::vector& upper = rins_worker->leaf_problem.upper; + std::vector& bounds_changed = rins_worker->bounds_changed; + std::vector& current_sol = rins_worker->leaf_solution.x; + std::vector fractional; + i_t num_frac = fractional_variables(settings_, current_sol, var_types_, fractional); + + std::vector current_incumbent; + mutex_upper_.lock(); + current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + + std::vector integer_list; + for (i_t j = 0; j < var_types_.size(); ++j) { + if (var_types_[j] == variable_type_t::CONTINUOUS) { continue; } + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + integer_list.push_back(j); + } + i_t num_integers = integer_list.size(); + + f_t max_fixrate = + submip_get_max_fixrate(submip_stats_, settings_.submip_settings, rins_worker->rng); + f_t min_fixrate = std::min(settings_.submip_settings.min_fixrate, max_fixrate); + i_t max_var_fixed = max_fixrate * num_integers; + i_t min_var_fixed = min_fixrate * num_integers; + i_t num_var_fixed = 0; + + while (solver_status_ == mip_status_t::UNSET && is_running_) { + // RINS neighbourhood 1: Fix all the integer variables where the starting solution matches the + // current incumbent, considering only the fractional values in the current node + i_t prev_num_fixed = num_var_fixed; + for (i_t j : fractional) { + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + if (std::abs(node_solution[j] - current_incumbent[j]) <= settings_.integer_tol) { + f_t fixed_val = std::round(node_solution[j]); + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + } + } + + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; + break; + } + + if (prev_num_fixed == num_var_fixed) { + // RINS neighbourhood 2: Search into the entire list of integer variables where the current + // LP solution matches the current incumbent. + for (i_t j : integer_list) { + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + if (std::abs(current_sol[j] - current_incumbent[j]) <= settings_.integer_tol) { + f_t fixed_val = std::round(current_sol[j]); + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + } + } + + if (num_var_fixed >= min_var_fixed) { + settings_.log.debug( + "%sFixed %d variables (min=%d)\n", log_prefix.c_str(), num_var_fixed, min_var_fixed); + need_to_solve_submip = true; + break; + } + + if (prev_num_fixed == num_var_fixed) { + settings_.log.debug("%sFixed the maximum (%d) number of variables for RINS (min=%d)\n", + log_prefix.c_str(), + num_var_fixed, + min_var_fixed); + need_to_solve_submip = true; + break; + } + } + + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; + break; + } + + dual_status_t lp_status = solve_node_lp(&node, rins_worker, rins_stats, log); + if (lp_status != dual_status_t::OPTIMAL) { break; } + + if (lp_status == dual_status_t::OPTIMAL) { + fractional.clear(); + num_frac = fractional_variables(settings_, current_sol, var_types_, fractional); + + f_t leaf_obj = compute_objective(rins_worker->leaf_problem, current_sol); + node.lower_bound = leaf_obj; + + apply_objective_step(&node, leaf_obj); + + if (num_frac == 0) { + add_feasible_solution(leaf_obj, current_sol, -1, SUBMIP); + break; + } + + if (leaf_obj > upper_bound_.load()) { break; } + } + } + + f_t fixrate = (f_t)num_var_fixed / num_integers; + + if (need_to_solve_submip) { + bool feasible = + rins_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); + submip_worker_pool_.return_worker_to_pool(rins_worker); + 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); + + simplex_solver_settings_t submip_settings; + submip_settings.node_limit = + settings_.submip_settings.node_limit_base + exploration_stats_.nodes_explored / 20; + submip_settings.relative_mip_gap_tol = + std::min(settings_.submip_settings.target_mip_gap, rel_gap); + submip_settings.print_presolve_stats = false; + submip_settings.num_threads = 1; + submip_settings.reliability_branching = 0; + submip_settings.max_cut_passes = 0; + submip_settings.clique_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; + submip_settings.log.log_prefix = log_prefix; + submip_settings.submip_settings.enable_rins = + submip_level <= settings_.submip_settings.max_level; + submip_settings.solution_callback = [this, fixrate](std::vector& solution, f_t objective) { + this->set_solution_from_submip(solution, fixrate); + }; + + branch_and_bound_t submip_bnb(*this, submip_settings, lower, upper); + mip_solution_t fixed_solution(original_problem_.num_cols); + mip_status_t submip_status = submip_bnb.solve(fixed_solution); + + if (submip_status == mip_status_t::INFEASIBLE) { + submip_stats_.save_infeasible(fixrate); + + } else if (submip_status == mip_status_t::OPTIMAL) { + std::vector crushed_solution; + crush_primal_solution( + original_problem_, original_lp_, fixed_solution.x, new_slacks_, crushed_solution); + f_t obj = compute_objective(original_lp_, crushed_solution); + settings_.log.debug("%sFound a solution with obj=%.4g\n", + log_prefix.c_str(), + compute_user_objective(original_lp_, obj)); + + if (improves_incumbent(obj)) { + add_feasible_solution(obj, crushed_solution, -1, SUBMIP); + submip_stats_.save_success(fixrate); + } + } + } + + settings_.log.debug( + "%ssuccess=%d, infeasible=%d, calls=%d, fixrate=%.4g (%d), max_fixrate=%.4g (%d), " + "min_fixrate=%.4g (%d)\n", + log_prefix.c_str(), + submip_stats_.total_success, + submip_stats_.total_infeasible, + submip_stats_.total_calls, + fixrate, + num_var_fixed, + max_fixrate, + max_var_fixed, + min_fixrate, + min_var_fixed); + + submip_worker_pool_.return_worker_to_pool(rins_worker); +} + template lp_status_t branch_and_bound_t::solve_root_relaxation( simplex_solver_settings_t const& lp_settings, @@ -2525,33 +2872,58 @@ 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_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - } else { - settings_.log.printf("\nSolving LP root relaxation in concurrent mode\n"); - root_status = solve_root_relaxation(lp_settings, - root_relax_soln_, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - edge_norms_); + if (root_warm_start_) { + i_t node_iter = 0; + f_t lp_start_time = tic(); + + dual_status_t lp_status = dual_phase2_with_advanced_basis(2, + 0, + true, + lp_start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + node_iter, + edge_norms_); + + root_status = convert_dual_status_to_lp_status(lp_status); } + + if (root_status == lp_status_t::NUMERICAL_ISSUES || root_status == lp_status_t::UNSET) { + 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_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + + } else { + settings_.log.printf("\nSolving LP root relaxation in concurrent mode\n"); + root_status = solve_root_relaxation(lp_settings, + root_relax_soln_, + root_vstatus_, + basis_update, + basic_list, + 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); @@ -2840,7 +3212,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut pc_.resize(original_lp_.num_cols); pc_.Arow = Arow_; - { + + if (!root_warm_start_) { raft::common::nvtx::range scope_sb("BB::strong_branching"); strong_branching(original_lp_, settings_, @@ -2959,6 +3332,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut const i_t num_bfs_workers = std::max(settings_.num_threads / 2, 1); const i_t num_diving_workers = num_workers - num_bfs_workers; bfs_worker_pool_.init(num_bfs_workers, original_lp_, Arow_, var_types_, symmetry_, settings_); + submip_worker_pool_.init(1, original_lp_, Arow_, var_types_, symmetry_, settings_); if (num_diving_workers > 0) { diving_worker_pool_.init(num_diving_workers, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 18eb38de66..f1bd0b5864 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -50,17 +51,6 @@ namespace cuopt::mathematical_optimization::mip { template struct mip_symmetry_t; -enum class mip_status_t { - OPTIMAL = 0, // The optimal integer solution was found - UNBOUNDED = 1, // The problem is unbounded - INFEASIBLE = 2, // The problem is infeasible - TIME_LIMIT = 3, // The solver reached a time limit - NODE_LIMIT = 4, // The maximum number of nodes was reached (not implemented) - NUMERICAL = 5, // The solver encountered a numerical error - UNSET = 6, // The status is not set - WORK_LIMIT = 7, // The solver reached a deterministic work limit -}; - template void upper_bound_callback(f_t upper_bound); @@ -83,6 +73,11 @@ class branch_and_bound_t { std::shared_ptr> clique_table = nullptr, mip_symmetry_t* symmetry = nullptr); + branch_and_bound_t(const branch_and_bound_t& other, + const simplex::simplex_solver_settings_t& solver_settings, + const std::vector& lower, + const std::vector& upper); + // 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; } @@ -111,6 +106,9 @@ 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); + // Set a solution based on the user problem during the course of the solve + void set_solution_from_submip(const std::vector& solution, f_t fixrate); + // 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); @@ -151,6 +149,15 @@ class branch_and_bound_t { std::vector& lower_bounds, std::vector& upper_bounds); + // Whether obj should replace the stored incumbent. Must be called under mutex_upper_. + // Compares against the stored incumbent's objective, NOT against upper_bound_, because + // set_initial_upper_bound can set a tighter bound from an OG-space solution that has no + // corresponding solver-space incumbent (e.g. papilo can't crush it back). + bool improves_incumbent(f_t obj) const + { + return !incumbent_.has_incumbent || obj < incumbent_.objective; + } + // The main entry routine. Returns the solver status and populates solution with the incumbent. mip_status_t solve(simplex::mip_solution_t& solution); @@ -199,18 +206,12 @@ 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_; + // For submip + const omp_atomic_t* external_upper_bound_{nullptr}; + // Solver-space incumbent tracked directly by B&B. simplex::mip_solution_t incumbent_; - // Whether obj should replace the stored incumbent. Must be called under mutex_upper_. - // Compares against the stored incumbent's objective, NOT against upper_bound_, because - // set_initial_upper_bound can set a tighter bound from an OG-space solution that has no - // corresponding solver-space incumbent (e.g. papilo can't crush it back). - bool improves_incumbent(f_t obj) const - { - return !incumbent_.has_incumbent || obj < incumbent_.objective; - } - // Structure with the general info of the solver. branch_and_bound_stats_t exploration_stats_; @@ -233,6 +234,7 @@ class branch_and_bound_t { std::atomic root_concurrent_halt_{0}; std::atomic node_concurrent_halt_{0}; bool is_root_solution_set{false}; + bool root_warm_start_{false}; // Pseudocosts pseudo_costs_t pc_; @@ -246,6 +248,9 @@ class branch_and_bound_t { // Worker pool dedicated to diving diving_worker_pool_t diving_worker_pool_; + submip_worker_pool_t submip_worker_pool_; + submip_stats_t submip_stats_; + // Global status of the solver. omp_atomic_t solver_status_; omp_atomic_t is_running_{false}; @@ -316,6 +321,14 @@ class branch_and_bound_t { // Launch a new diving worker from a given best-first worker. bool launch_diving_worker(bfs_worker_t* bfs_worker); + // If the objective is integral or must move in steps than + // the lower bound will be different from the leaf objective. + // We use the leaf objective for RINS (on_optimal_callback) + // and if we are integer feasible (handle_integer_solution). + // We use the lower bound to decide if we should fathom the + // node or branch. + void apply_objective_step(mip_node_t* node_ptr, f_t leaf_obj); + // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); @@ -334,6 +347,10 @@ class branch_and_bound_t { // to find integer feasible solutions. void dive_with(diving_worker_t* worker); + void launch_submip_worker(const std::vector& sol); + + void rins(submip_worker_t* rins_worker, const std::vector& node_solution); + // Solve the LP relaxation of a leaf node simplex::dual_status_t solve_node_lp(mip_node_t* node_ptr, branch_and_bound_worker_t* worker, diff --git a/cpp/src/branch_and_bound/constants.hpp b/cpp/src/branch_and_bound/constants.hpp index 0194e23da4..f3e21126dd 100644 --- a/cpp/src/branch_and_bound/constants.hpp +++ b/cpp/src/branch_and_bound/constants.hpp @@ -9,7 +9,7 @@ namespace cuopt::mathematical_optimization::mip { -constexpr int num_search_strategies = 7; +constexpr int num_search_strategies = 8; // Indicate the search and variable selection algorithms used by each thread // in B&B (See [1]). @@ -26,7 +26,8 @@ enum search_strategy_t : int { GUIDED_DIVING = 3, // Guided diving (9.2.3). COEFFICIENT_DIVING = 4, // Coefficient diving (9.2.1) FARKAS_DIVING = 5, // Farkas Diving (see [2]) - VECTOR_LENGTH_DIVING = 6 // Vector Length Diving (9.2.6) + VECTOR_LENGTH_DIVING = 6, // Vector Length Diving (9.2.6) + SUBMIP = 7, // RINS/RENS (akin to a guided diving, see HiGHS) }; constexpr search_strategy_t search_strategies[] = {BEST_FIRST, @@ -35,10 +36,23 @@ constexpr search_strategy_t search_strategies[] = {BEST_FIRST, GUIDED_DIVING, COEFFICIENT_DIVING, FARKAS_DIVING, - VECTOR_LENGTH_DIVING}; + VECTOR_LENGTH_DIVING, + SUBMIP}; enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 }; enum class branch_and_bound_mode_t { PARALLEL = 0, DETERMINISTIC = 1 }; +enum class mip_status_t { + OPTIMAL = 0, // The optimal integer solution was found + UNBOUNDED = 1, // The problem is unbounded + INFEASIBLE = 2, // The problem is infeasible + TIME_LIMIT = 3, // The solver reached a time limit + NODE_LIMIT = 4, // The maximum number of nodes was reached (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 + SUBMIP_HALT = 8 // Stop submip solve since it no longer valid +}; + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/diving_heuristics.hpp b/cpp/src/branch_and_bound/diving_heuristics.hpp index f778ae38dc..2ae279f952 100644 --- a/cpp/src/branch_and_bound/diving_heuristics.hpp +++ b/cpp/src/branch_and_bound/diving_heuristics.hpp @@ -23,6 +23,7 @@ namespace cuopt::mathematical_optimization::mip { inline char feasible_solution_symbol(search_strategy_t strategy, bool log_diving_type) { if (strategy == BEST_FIRST) return 'B'; + if (strategy == SUBMIP) return 'S'; if (!log_diving_type) { return 'D'; } switch (strategy) { case COEFFICIENT_DIVING: return 'C'; @@ -47,6 +48,7 @@ bool is_search_strategy_enabled(search_strategy_t strategy, 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; + case SUBMIP: return false; } return false; diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index baea370fa4..9987cc7529 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -1046,7 +1046,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) { @@ -1059,7 +1059,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"); } @@ -1573,7 +1573,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 @@ -1593,7 +1593,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 @@ -1609,7 +1609,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/submip.hpp b/cpp/src/branch_and_bound/submip.hpp new file mode 100644 index 0000000000..a528ce5f0a --- /dev/null +++ b/cpp/src/branch_and_bound/submip.hpp @@ -0,0 +1,94 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +template +class branch_and_bound_t; + +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; } +}; + +inline double submip_get_max_fixrate(const submip_stats_t& stats, + const simplex::submip_settings_t& submip_settings, + pcgenerator_t& rng) +{ + // Adaptive fix rate based on previous successes and failures. + double low = submip_settings.base_target_fixrate; + double high = submip_settings.base_target_fixrate; + + if (stats.total_infeasible > 0) { + double infeasible_avg_fixrate = stats.average_infeasible_fixrate(); + high = 0.9 * infeasible_avg_fixrate; + low = std::min(low, high); + } + + if (stats.total_success > 0) { + double success_avg_fixrate = stats.average_success_fixrate(); + low = std::min(low, 0.9 * success_avg_fixrate); + high = std::max(high, 1.1 * success_avg_fixrate); + } + + double fixrate = high > low ? rng.uniform(low, high) : low; + return fixrate; +} + +template +class submip_worker_t : public branch_and_bound_worker_t { + public: + using Base = branch_and_bound_worker_t; + + submip_worker_t(i_t worker_id, + const simplex::lp_problem_t& original_lp, + const simplex::csr_matrix_t& Arow, + const std::vector& var_type, + const simplex::simplex_solver_settings_t& settings, + uint64_t rng_offset = 0) + : Base(worker_id, original_lp, Arow, var_type, settings, rng_offset) + { + this->search_strategy = SUBMIP; + } + + // Set this node inactive + void set_inactive() + { + if (!this->is_active.load()) { return; } + this->is_active = false; + } + + f_t get_lower_bound() { return this->lower_bound; } +}; + +template +using submip_worker_pool_t = worker_pool_t>; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index 758bbe16d8..8fdcbc59fe 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -41,6 +41,8 @@ struct branch_and_bound_stats_t { omp_atomic_t lexical_reduction_nodes = 0; omp_atomic_t lexical_reduction_fixings_applied = 0; omp_atomic_t lexical_reduction_pruned_nodes = 0; + + omp_atomic_t nodes_at_last_submip = 0; }; template @@ -116,6 +118,9 @@ class branch_and_bound_worker_t { bool set_lp_variable_bounds(mip_node_t* node_ptr, const simplex::simplex_solver_settings_t& settings) { + // In RINS/RENS, the bounds are already set in the parent method. + if (search_strategy == SUBMIP) return true; + // Reset the bound_changed markers std::fill(bounds_changed.begin(), bounds_changed.end(), false); diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 40efbe6185..fd9d325846 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -3265,7 +3265,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; @@ -3605,7 +3605,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/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 45168ced8a..1f6e80ca83 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -26,6 +26,33 @@ struct benchmark_info_t; namespace cuopt::mathematical_optimization::simplex { +struct submip_settings_t { + // Enable or disable (recursive) RINS + int enable_rins = -1; + + // Base for calculating the target fix rate for the RINS neighbourhood (adaptive according + // to the success and infeasible rate) + double base_target_fixrate = 0.6; + + // Minimum fix rate for accepting the RINS neighbourhood + double min_fixrate = 0.2; + + // MIP gap for the submip (unless the MIP gap from the B&B is lower) + double target_mip_gap = 0.01; + + // The base node limit for the submip + int node_limit_base = 200; + + // Frequency in terms of nodes when to launch a new RINS + int node_freq = 200; + + // The current level in the recursion + int level = 0; + + // Maximum recursion level + int max_level = 5; +}; + template struct simplex_solver_settings_t { public: @@ -102,7 +129,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), @@ -214,7 +241,9 @@ struct simplex_solver_settings_t { i_t reliability_branching; i_t inside_mip; // 0 if outside MIP, 1 if inside MIP at root node, 2 if inside MIP at leaf node - i_t sub_mip; // 0 if in regular MIP solve, 1 if in sub-MIP solve + i_t inside_submip; // 0 if in regular MIP solve, 1 if in sub-MIP solve + + submip_settings_t submip_settings; std::function&, f_t)> solution_callback; std::function&, f_t)> node_processed_callback; diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index 50b8ccab02..4bb35bc35d 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -627,7 +627,7 @@ solution_t diversity_manager_t::run_solver() return sol; } - if (omp_get_num_threads() > CUOPT_MIP_RINS_REQUIRED_THREAD_COUNT) { rins.enable(); } + // if (omp_get_num_threads() > CUOPT_MIP_RINS_REQUIRED_THREAD_COUNT) { rins.enable(); } generate_solution(timer.remaining_time(), false); if (timer.check_time_limit()) { diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index a670322cb0..40594b7956 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -254,7 +254,7 @@ void rins_t::run_rins() 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.sub_mip = 1; + branch_and_bound_settings.inside_submip = 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] "; diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index 0f3625d9aa..a61d23a892 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -111,7 +111,7 @@ class sub_mip_recombiner_t : public recombiner_t { 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.sub_mip = 1; + branch_and_bound_settings.inside_submip = 1; branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; branch_and_bound_settings.solution_callback = [this](std::vector& solution, f_t objective) { From ba99a217c8168bfabcd09bfa7cec902dca7abb81 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Mon, 29 Jun 2026 12:56:31 +0200 Subject: [PATCH 02/16] fix time limit violation when solving the submip Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e5c112f09e..91626aa798 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -342,7 +342,7 @@ branch_and_bound_t::branch_and_bound_t( original_lp_.lower = lower; original_lp_.upper = upper; - exploration_stats_.start_time = other.exploration_stats_.start_time; + exploration_stats_.start_time = tic(); upper_bound_ = other.upper_bound_; root_objective_ = std::numeric_limits::quiet_NaN(); root_lp_current_lower_bound_ = -inf; @@ -2318,10 +2318,17 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, 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; - submip_settings.node_limit = - settings_.submip_settings.node_limit_base + exploration_stats_.nodes_explored / 20; + submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; + + submip_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + if (submip_settings.time_limit < 0) { + submip_worker_pool_.return_worker_to_pool(rins_worker); + return; + } + submip_settings.relative_mip_gap_tol = std::min(settings_.submip_settings.target_mip_gap, rel_gap); submip_settings.print_presolve_stats = false; From 727e20e4fdf685aa8d61c5c525a33e5798a3b237 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Mon, 29 Jun 2026 13:43:14 +0200 Subject: [PATCH 03/16] clean up the code a bit. print submip symbol in the logs instead of H Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 182 ++++++++++-------- cpp/src/branch_and_bound/branch_and_bound.hpp | 3 +- .../dual_simplex/simplex_solver_settings.hpp | 11 +- 3 files changed, 116 insertions(+), 80 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 91626aa798..7ad6e3b955 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -399,7 +399,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, char symbol) { if (is_running_) { f_t lower_bound = get_lower_bound(); @@ -409,7 +409,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}", + symbol, "", // nodes explored "", // nodes unexplored user_obj, @@ -574,7 +575,7 @@ void branch_and_bound_t::set_solution_from_heuristics(const std::vecto f_t current_upper_bound = upper_bound_.load(); upper_bound_ = std::min(current_upper_bound, obj); incumbent_.set_incumbent_solution(obj, crushed_solution); - if (current_upper_bound > upper_bound_.load()) { report_heuristic(obj); } + if (current_upper_bound > upper_bound_.load()) { report_heuristic(obj, 'H'); } } else { attempt_repair = true; constexpr bool verbose = false; @@ -616,13 +617,15 @@ void branch_and_bound_t::set_solution_from_submip(const std::vector upper_bound_.load()) { report_heuristic(obj); } + if (current_upper_bound > upper_bound_.load()) { + report_heuristic(obj, feasible_solution_symbol(SUBMIP, true)); + } if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -771,7 +774,7 @@ void branch_and_bound_t::repair_heuristic_solutions() if (improves_incumbent(repaired_obj)) { upper_bound_ = std::min(upper_bound_.load(), repaired_obj); incumbent_.set_incumbent_solution(repaired_obj, repaired_solution); - report_heuristic(repaired_obj); + report_heuristic(repaired_obj, 'H'); if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -2169,6 +2172,79 @@ void branch_and_bound_t::launch_submip_worker(const std::vector& exploration_stats_.nodes_at_last_submip = exploration_stats_.nodes_explored; } +template +void branch_and_bound_t::solve_submip(submip_worker_t* rins_worker, + i_t num_var_fixed, + i_t num_integers) +{ + i_t submip_level = settings_.submip_settings.level + 1; + std::string log_prefix = std::format("[SUBMIP {}] ", submip_level); + + std::vector& lower = rins_worker->leaf_problem.lower; + std::vector& upper = rins_worker->leaf_problem.upper; + std::vector& bounds_changed = rins_worker->bounds_changed; + f_t fixrate = (f_t)num_var_fixed / num_integers; + + bool feasible = + rins_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; + submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; + + 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.print_presolve_stats = false; + submip_settings.num_threads = 1; + submip_settings.reliability_branching = 0; + submip_settings.max_cut_passes = 0; + submip_settings.clique_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; + submip_settings.log.log_prefix = log_prefix; + submip_settings.submip_settings.enable_rins = submip_level <= settings_.submip_settings.max_level; + submip_settings.solution_callback = [this, fixrate](std::vector& solution, f_t objective) { + this->set_solution_from_submip(solution, fixrate); + }; + + branch_and_bound_t submip_bnb(*this, submip_settings, lower, upper); + mip_solution_t fixed_solution(original_problem_.num_cols); + mip_status_t submip_status = submip_bnb.solve(fixed_solution); + + if (submip_status == mip_status_t::INFEASIBLE) { + submip_stats_.save_infeasible(fixrate); + + } else if (submip_status == mip_status_t::OPTIMAL) { + std::vector crushed_solution; + crush_primal_solution( + original_problem_, original_lp_, fixed_solution.x, new_slacks_, crushed_solution); + f_t obj = compute_objective(original_lp_, crushed_solution); + settings_.log.debug("%sFound a solution with obj=%.4g\n", + log_prefix.c_str(), + compute_user_objective(original_lp_, obj)); + + if (improves_incumbent(obj)) { + add_feasible_solution(obj, crushed_solution, -1, SUBMIP); + submip_stats_.save_success(fixrate); + } + } +} + template void branch_and_bound_t::rins(submip_worker_t* rins_worker, const std::vector& node_solution) @@ -2183,8 +2259,8 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, ++submip_stats_.total_calls; - bool need_to_solve_submip = false; - const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + bool has_submip = false; + const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; branch_and_bound_stats_t rins_stats; @@ -2239,6 +2315,14 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, } } + // Enough variables has been fixed + if (num_var_fixed >= min_var_fixed) { + settings_.log.debug( + "%sFixed %d variables (min=%d)\n", log_prefix.c_str(), num_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; @@ -2260,19 +2344,22 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, } } + // Enough variables were fixed if (num_var_fixed >= min_var_fixed) { settings_.log.debug( "%sFixed %d variables (min=%d)\n", log_prefix.c_str(), num_var_fixed, min_var_fixed); - need_to_solve_submip = true; + has_submip = true; break; } + // Even considering the entire integer list, we were unable to fix a single variable in this + // iteration. Try to solve the submip with what we got. if (prev_num_fixed == num_var_fixed) { settings_.log.debug("%sFixed the maximum (%d) number of variables for RINS (min=%d)\n", log_prefix.c_str(), num_var_fixed, min_var_fixed); - need_to_solve_submip = true; + has_submip = true; break; } } @@ -2295,7 +2382,8 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, apply_objective_step(&node, leaf_obj); if (num_frac == 0) { - add_feasible_solution(leaf_obj, current_sol, -1, SUBMIP); + // We found a feasible solution when fixing the variables in RINS. + add_feasible_solution(leaf_obj, current_sol, node.depth, SUBMIP); break; } @@ -2305,68 +2393,12 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, f_t fixrate = (f_t)num_var_fixed / num_integers; - if (need_to_solve_submip) { - bool feasible = - rins_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); - submip_worker_pool_.return_worker_to_pool(rins_worker); - 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; - submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; - - submip_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); - if (submip_settings.time_limit < 0) { - submip_worker_pool_.return_worker_to_pool(rins_worker); - return; - } - - submip_settings.relative_mip_gap_tol = - std::min(settings_.submip_settings.target_mip_gap, rel_gap); - submip_settings.print_presolve_stats = false; - submip_settings.num_threads = 1; - submip_settings.reliability_branching = 0; - submip_settings.max_cut_passes = 0; - submip_settings.clique_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; - submip_settings.log.log_prefix = log_prefix; - submip_settings.submip_settings.enable_rins = - submip_level <= settings_.submip_settings.max_level; - submip_settings.solution_callback = [this, fixrate](std::vector& solution, f_t objective) { - this->set_solution_from_submip(solution, fixrate); - }; - - branch_and_bound_t submip_bnb(*this, submip_settings, lower, upper); - mip_solution_t fixed_solution(original_problem_.num_cols); - mip_status_t submip_status = submip_bnb.solve(fixed_solution); - - if (submip_status == mip_status_t::INFEASIBLE) { - submip_stats_.save_infeasible(fixrate); - - } else if (submip_status == mip_status_t::OPTIMAL) { - std::vector crushed_solution; - crush_primal_solution( - original_problem_, original_lp_, fixed_solution.x, new_slacks_, crushed_solution); - f_t obj = compute_objective(original_lp_, crushed_solution); - settings_.log.debug("%sFound a solution with obj=%.4g\n", - log_prefix.c_str(), - compute_user_objective(original_lp_, obj)); - - if (improves_incumbent(obj)) { - add_feasible_solution(obj, crushed_solution, -1, SUBMIP); - submip_stats_.save_success(fixrate); - } + if (has_submip) { + // Only solve the submip if the fix rate is greater than a hard cap (0.1). Less than that + // there is no enough variable fixing in the submip to be easier to solve than the current + // problem + if (fixrate >= settings_.submip_settings.min_fixrate_cap) { + solve_submip(rins_worker, num_var_fixed, num_integers); } } @@ -4185,7 +4217,7 @@ void branch_and_bound_t::deterministic_sort_replay_events( } if (new_upper < std::numeric_limits::infinity()) { - report_heuristic(new_upper); + report_heuristic(new_upper, 'H'); 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 f1bd0b5864..2eed732b73 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -267,7 +267,7 @@ class branch_and_bound_t { std::function user_bound_callback_; void print_table_header(); - void report_heuristic(f_t obj); + void report_heuristic(f_t obj, char symbol); void report(char symbol, f_t obj, f_t lower_bound, @@ -348,6 +348,7 @@ class branch_and_bound_t { void dive_with(diving_worker_t* worker); void launch_submip_worker(const std::vector& sol); + void solve_submip(submip_worker_t* rins_worker, i_t num_var_fixed, i_t num_integers); void rins(submip_worker_t* rins_worker, const std::vector& node_solution); diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 1f6e80ca83..6848947d86 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -30,12 +30,15 @@ struct submip_settings_t { // Enable or disable (recursive) RINS int enable_rins = -1; - // Base for calculating the target fix rate for the RINS neighbourhood (adaptive according - // to the success and infeasible rate) + // Base for calculating the target fix rate for the RINS neighbourhood. Actual target value is + // determined automatically according to the success and infeasible rate. double base_target_fixrate = 0.6; - // Minimum fix rate for accepting the RINS neighbourhood - double min_fixrate = 0.2; + // Minimum fix rate for accepting the RINS neighbourhood. + double min_fixrate = 0.25; + + // Hard cap for the minimum fix rate for solving a submip. + double min_fixrate_cap = 0.1; // MIP gap for the submip (unless the MIP gap from the B&B is lower) double target_mip_gap = 0.01; From 37182d3b07313119ef272151f678282f2f697fc9 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Mon, 29 Jun 2026 15:24:44 +0200 Subject: [PATCH 04/16] fix potential race condition Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 67 ++++++++++++------- cpp/src/branch_and_bound/branch_and_bound.hpp | 8 ++- 2 files changed, 47 insertions(+), 28 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 7ad6e3b955..3febeeb16f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -317,7 +317,7 @@ branch_and_bound_t::branch_and_bound_t( template branch_and_bound_t::branch_and_bound_t( - const branch_and_bound_t& other, + branch_and_bound_t& other, const simplex_solver_settings_t& solver_settings, const std::vector& lower, const std::vector& upper) @@ -330,17 +330,23 @@ branch_and_bound_t::branch_and_bound_t( Arow_(other.Arow_), new_slacks_(other.new_slacks_), var_types_(other.var_types_), - incumbent_(other.incumbent_), + incumbent_(1), root_vstatus_(other.root_vstatus_), root_relax_soln_(1, 1), root_crossover_soln_(1, 1), pc_(other.pc_), solver_status_(mip_status_t::UNSET) { + other.mutex_original_lp_.lock(); assert(lower.size() == original_lp_.num_cols); assert(upper.size() == original_lp_.num_cols); original_lp_.lower = lower; original_lp_.upper = upper; + other.mutex_original_lp_.unlock(); + + other.mutex_upper_.lock(); + incumbent_ = other.incumbent_; + other.mutex_upper_.unlock(); exploration_stats_.start_time = tic(); upper_bound_ = other.upper_bound_; @@ -2173,20 +2179,19 @@ void branch_and_bound_t::launch_submip_worker(const std::vector& } template -void branch_and_bound_t::solve_submip(submip_worker_t* rins_worker, +void branch_and_bound_t::solve_submip(submip_worker_t* worker, i_t num_var_fixed, - i_t num_integers) + i_t num_integers, + i_t submip_level, + std::string_view log_prefix) { - i_t submip_level = settings_.submip_settings.level + 1; - std::string log_prefix = std::format("[SUBMIP {}] ", submip_level); - - 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& lower = worker->leaf_problem.lower; + std::vector& upper = worker->leaf_problem.upper; + std::vector& bounds_changed = worker->bounds_changed; f_t fixrate = (f_t)num_var_fixed / num_integers; bool feasible = - rins_worker->node_presolver.bounds_strengthening(settings_, bounds_changed, lower, upper); + 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. @@ -2215,9 +2220,12 @@ void branch_and_bound_t::solve_submip(submip_worker_t* rins_ submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.submip_settings.level = submip_level; - submip_settings.log.log = false; + submip_settings.log.log = true; submip_settings.log.log_prefix = log_prefix; - submip_settings.submip_settings.enable_rins = submip_level <= settings_.submip_settings.max_level; + + submip_settings.submip_settings.enable_rins = settings_.submip_settings.enable_rins != 0 && + submip_level > settings_.submip_settings.max_level; + submip_settings.solution_callback = [this, fixrate](std::vector& solution, f_t objective) { this->set_solution_from_submip(solution, fixrate); }; @@ -2234,9 +2242,9 @@ void branch_and_bound_t::solve_submip(submip_worker_t* rins_ crush_primal_solution( original_problem_, original_lp_, fixed_solution.x, new_slacks_, crushed_solution); f_t obj = compute_objective(original_lp_, crushed_solution); - settings_.log.debug("%sFound a solution with obj=%.4g\n", - log_prefix.c_str(), - compute_user_objective(original_lp_, obj)); + settings_.log.debug_format("{}Found a solution with obj={:.4g}\n", + log_prefix, + compute_user_objective(original_lp_, obj)); if (improves_incumbent(obj)) { add_feasible_solution(obj, crushed_solution, -1, SUBMIP); @@ -2255,7 +2263,7 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, log.log = false; i_t submip_level = settings_.submip_settings.level + 1; - std::string log_prefix = std::format("[SUBMIP {}] ", submip_level); + std::string log_prefix = std::format("[RINS {}] ", submip_level); ++submip_stats_.total_calls; @@ -2317,8 +2325,11 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // Enough variables has been fixed if (num_var_fixed >= min_var_fixed) { - settings_.log.debug( - "%sFixed %d variables (min=%d)\n", log_prefix.c_str(), num_var_fixed, min_var_fixed); + settings_.log.print_format("{}Fixed {} variables (max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } @@ -2346,8 +2357,11 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // Enough variables were fixed if (num_var_fixed >= min_var_fixed) { - settings_.log.debug( - "%sFixed %d variables (min=%d)\n", log_prefix.c_str(), num_var_fixed, min_var_fixed); + settings_.log.print_format("{}Fixed {} variables (max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } @@ -2355,10 +2369,11 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // Even considering the entire integer list, we were unable to fix a single variable in this // iteration. Try to solve the submip with what we got. if (prev_num_fixed == num_var_fixed) { - settings_.log.debug("%sFixed the maximum (%d) number of variables for RINS (min=%d)\n", - log_prefix.c_str(), - num_var_fixed, - min_var_fixed); + settings_.log.print_format("{}Fixed {} variables, the maximum for RINS. (max={}, min={})\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); has_submip = true; break; } @@ -2398,7 +2413,7 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // there is no enough variable fixing in the submip to be easier to solve than the current // problem if (fixrate >= settings_.submip_settings.min_fixrate_cap) { - solve_submip(rins_worker, num_var_fixed, num_integers); + solve_submip(rins_worker, num_var_fixed, num_integers, submip_level, log_prefix); } } diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 2eed732b73..6be3405bf1 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -73,7 +73,7 @@ class branch_and_bound_t { std::shared_ptr> clique_table = nullptr, mip_symmetry_t* symmetry = nullptr); - branch_and_bound_t(const branch_and_bound_t& other, + branch_and_bound_t(branch_and_bound_t& other, const simplex::simplex_solver_settings_t& solver_settings, const std::vector& lower, const std::vector& upper); @@ -348,7 +348,11 @@ class branch_and_bound_t { void dive_with(diving_worker_t* worker); void launch_submip_worker(const std::vector& sol); - void solve_submip(submip_worker_t* rins_worker, i_t num_var_fixed, i_t num_integers); + void solve_submip(submip_worker_t* worker, + i_t num_var_fixed, + i_t num_integers, + i_t submip_level, + std::string_view log_prefix); void rins(submip_worker_t* rins_worker, const std::vector& node_solution); From bd7777b0adbd709ce0185bc5a670171d429b6a4d Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Mon, 29 Jun 2026 17:20:09 +0200 Subject: [PATCH 05/16] allow rins to fire whenever it is idling. fix incorrect recursion check. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 10 ++-------- cpp/src/branch_and_bound/branch_and_bound.hpp | 4 +--- cpp/src/branch_and_bound/worker.hpp | 2 -- cpp/src/dual_simplex/simplex_solver_settings.hpp | 3 --- 4 files changed, 3 insertions(+), 16 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 3febeeb16f..0e34e7620c 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2160,12 +2160,8 @@ template void branch_and_bound_t::launch_submip_worker(const std::vector& sol) { if (settings_.submip_settings.enable_rins == 0) return; - if (submip_worker_pool_.num_idle() == 0) return; if (!incumbent_.has_incumbent) return; - - i_t nodes_since_last_submip = - exploration_stats_.nodes_explored - exploration_stats_.nodes_at_last_submip; - if (nodes_since_last_submip < settings_.submip_settings.node_freq) return; + if (submip_worker_pool_.num_idle() == 0) return; submip_worker_t* submip_worker = submip_worker_pool_.pop_idle_worker(); @@ -2174,8 +2170,6 @@ void branch_and_bound_t::launch_submip_worker(const std::vector& #pragma omp task priority(CUOPT_MEDIUM_TASK_PRIORITY) affinity(submip_worker) \ firstprivate(submip_worker, sol) rins(submip_worker, sol); - - exploration_stats_.nodes_at_last_submip = exploration_stats_.nodes_explored; } template @@ -2224,7 +2218,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.log.log_prefix = log_prefix; submip_settings.submip_settings.enable_rins = settings_.submip_settings.enable_rins != 0 && - submip_level > settings_.submip_settings.max_level; + submip_level <= settings_.submip_settings.max_level; submip_settings.solution_callback = [this, fixrate](std::vector& solution, f_t objective) { this->set_solution_from_submip(solution, fixrate); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 6be3405bf1..608fd40862 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -42,11 +42,9 @@ #include namespace cuopt::mathematical_optimization::mip { + template struct clique_table_t; -} - -namespace cuopt::mathematical_optimization::mip { template struct mip_symmetry_t; diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index 8fdcbc59fe..67d8d40117 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -41,8 +41,6 @@ struct branch_and_bound_stats_t { omp_atomic_t lexical_reduction_nodes = 0; omp_atomic_t lexical_reduction_fixings_applied = 0; omp_atomic_t lexical_reduction_pruned_nodes = 0; - - omp_atomic_t nodes_at_last_submip = 0; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 6848947d86..9ee07dc54b 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -46,9 +46,6 @@ struct submip_settings_t { // The base node limit for the submip int node_limit_base = 200; - // Frequency in terms of nodes when to launch a new RINS - int node_freq = 200; - // The current level in the recursion int level = 0; From c46eecd0b037030a1bd7a2d6222f49d858bd3fe0 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 30 Jun 2026 11:50:26 +0200 Subject: [PATCH 06/16] added a presolver_t class for MIP so it can be used for presolving submip (and restarted) problems. Papilo can now be applied to an user_problem_t inplace. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/CMakeLists.txt | 1 + cpp/src/branch_and_bound/presolve.cpp | 51 +++ cpp/src/branch_and_bound/presolve.hpp | 40 +++ .../presolve/third_party_presolve.cpp | 305 ++++++++++++++++++ .../presolve/third_party_presolve.hpp | 14 +- cpp/tests/mip/presolve_test.cu | 48 +++ datasets/mip/download_miplib_test_dataset.sh | 1 + 7 files changed, 457 insertions(+), 3 deletions(-) create mode 100644 cpp/src/branch_and_bound/presolve.cpp create mode 100644 cpp/src/branch_and_bound/presolve.hpp diff --git a/cpp/src/branch_and_bound/CMakeLists.txt b/cpp/src/branch_and_bound/CMakeLists.txt index 47ca6b98f5..f8ab3c2b30 100644 --- a/cpp/src/branch_and_bound/CMakeLists.txt +++ b/cpp/src/branch_and_bound/CMakeLists.txt @@ -7,6 +7,7 @@ set(BRANCH_AND_BOUND_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/diving_heuristics.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/presolve.cpp ) diff --git a/cpp/src/branch_and_bound/presolve.cpp b/cpp/src/branch_and_bound/presolve.cpp new file mode 100644 index 0000000000..a5ba8a2f30 --- /dev/null +++ b/cpp/src/branch_and_bound/presolve.cpp @@ -0,0 +1,51 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +template +third_party_presolve_status_t presolver_t::apply( + simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings) +{ + f_t presolve_time_limit = std::min(0.1 * settings.time_limit, 60.0); + return third_party_presolver_.apply(problem, settings, presolve_time_limit, 1); +} + +template +void presolver_t::uncrush(const std::vector& reduced_primal, + std::vector& full_primal) const +{ + third_party_presolver_.uncrush_primal_solution(reduced_primal, full_primal); +} + +template +const std::vector& presolver_t::reduced_to_original_map() const +{ + return third_party_presolver_.get_reduced_to_original_map(); +} + +template +const std::vector& presolver_t::original_to_reduced_map() const +{ + return third_party_presolver_.get_original_to_reduced_map(); +} + +#if MIP_INSTANTIATE_FLOAT +template class presolver_t; +#endif + +#if MIP_INSTANTIATE_DOUBLE +template class presolver_t; +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/presolve.hpp b/cpp/src/branch_and_bound/presolve.hpp new file mode 100644 index 0000000000..4ee66dc21e --- /dev/null +++ b/cpp/src/branch_and_bound/presolve.hpp @@ -0,0 +1,40 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +// Thin owner of a PaPILO presolver scoped to a single sub-MIP solve. apply() reduces the +// problem in place; the retained column maps let a reduced-space solution be mapped back to +// the original column space via uncrush(). +template +class presolver_t { + public: + // Presolve `problem` in place using PaPILO. Returns the presolve status; on + // INFEASIBLE/UNBOUNDED the problem is left untouched. + third_party_presolve_status_t apply(simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings); + + // Map a reduced-space primal solution back to the original column space. + void uncrush(const std::vector& reduced_primal, std::vector& full_primal) const; + + const std::vector& reduced_to_original_map() const; + const std::vector& original_to_reduced_map() const; + + private: + third_party_presolve_t third_party_presolver_; +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 08bd7779a0..31602a3a12 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -204,6 +204,129 @@ papilo::Problem build_papilo_problem(const optimization_problem_t return problem; } +template +papilo::Problem build_papilo_problem(const simplex::user_problem_t& problem) +{ + raft::common::nvtx::range fun_scope("Build papilo problem"); + // Build a papilo problem from a (host-side) dual-simplex user_problem_t. Unlike the + // optimization_problem_t overload, all data already lives on the host and the constraint + // matrix is stored column-major (CSC), so there are no device copies and no COO step: the + // CSC matrix is converted once to CSR and handed straight to papilo's SparseStorage. + papilo::ProblemBuilder builder; + + const i_t num_cols = problem.num_cols; + const i_t num_rows = problem.num_rows; + const i_t nnz = problem.A.nnz(); + + builder.reserve(nnz, num_rows, num_cols); + + const std::vector& obj_coeffs = problem.objective; + const std::vector& var_lb = problem.lower; + const std::vector& var_ub = problem.upper; + const std::vector& var_types = problem.var_types; + const std::vector& row_sense = problem.row_sense; + const std::vector& rhs = problem.rhs; + + // Range rows carry an extra width and are listed separately; mark them so the row-bound + // derivation below matches convert_user_problem in dual_simplex/presolve.cpp. papilo + // represents ranged rows natively as two-sided lhs <= a^T x <= rhs, so (unlike the simplex + // path) we do not add slack columns. + std::vector range_of_row(num_rows, 0); + std::vector is_range_row(num_rows, false); + for (i_t k = 0; k < problem.num_range_rows; ++k) { + const i_t row = problem.range_rows[k]; + is_range_row[row] = true; + range_of_row[row] = problem.range_value[k]; + } + + // Derive two-sided row bounds [lhs, rhs] from the row sense. + std::vector h_constr_lb(num_rows); + std::vector h_constr_ub(num_rows); + for (i_t i = 0; i < num_rows; ++i) { + const f_t b = rhs[i]; + if (is_range_row[i]) { + const f_t r = std::abs(range_of_row[i]); + if (row_sense[i] == 'L') { + h_constr_lb[i] = b - r; + h_constr_ub[i] = b; + } else if (row_sense[i] == 'G') { + h_constr_lb[i] = b; + h_constr_ub[i] = b + r; + } else { // 'E' with a range becomes a two-sided row + h_constr_lb[i] = range_of_row[i] > 0 ? b : b - r; + h_constr_ub[i] = range_of_row[i] > 0 ? b + r : b; + } + } else if (row_sense[i] == 'L') { + h_constr_lb[i] = -std::numeric_limits::infinity(); + h_constr_ub[i] = b; + } else if (row_sense[i] == 'G') { + h_constr_lb[i] = b; + h_constr_ub[i] = std::numeric_limits::infinity(); + } else { // 'E' + h_constr_lb[i] = b; + h_constr_ub[i] = b; + } + } + + builder.setNumCols(num_cols); + builder.setNumRows(num_rows); + + // user_problem_t stores the objective already in minimization sense (obj_scale carries the + // original min/max direction for reporting only), so no sign flip is needed here. + builder.setObjAll(obj_coeffs); + builder.setObjOffset(problem.obj_constant); + + if (!var_lb.empty() && !var_ub.empty()) { + builder.setColLbAll(var_lb); + builder.setColUbAll(var_ub); + if (static_cast(problem.col_names.size()) == num_cols) { + builder.setColNameAll(problem.col_names); + } + } + + for (i_t j = 0; j < num_cols; ++j) { + builder.setColIntegral(j, var_types[j] != simplex::variable_type_t::CONTINUOUS); + } + + // Row bounds + infinity flags, set on the builder so build() materializes the constraint + // matrix directly. build() also sets RowFlag::kEquation where a finite lhs == rhs. + if (num_rows > 0) { + builder.setRowLhsAll(h_constr_lb); + builder.setRowRhsAll(h_constr_ub); + } + 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(); + builder.setRowLhsInf(i, lhs_inf); + builder.setRowRhsInf(i, rhs_inf); + // Zero the placeholder value papilo stores for an infinite side. + if (lhs_inf) { builder.setRowLhs(i, 0); } + if (rhs_inf) { builder.setRowRhs(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); } + } + + // Feed the constraint matrix column-by-column straight from the CSC storage -- no COO + // triplet list and no CSC->CSR conversion. addColEntries copies into the builder's matrix + // buffer, so nothing aliases op_problem's arrays. The matrix is assumed free of explicitly + // stored zeros (matrix_buffer.addEntry asserts a nonzero value in debug builds). + 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) { + const i_t start = col_start[j]; + const i_t len = col_start[j + 1] - start; + if (len > 0) { builder.addColEntries(j, len, &row_index[start], &values[start]); } + } + + return builder.build(); +} + struct PSLPContext { Presolver* presolver = nullptr; Settings* settings = nullptr; @@ -478,6 +601,108 @@ 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. This is the inverse of build_papilo_problem_mip: the objective +// stays in minimization sense (no sign flip), two-sided papilo row bounds are mapped back to +// (row_sense, rhs, range_rows/range_value), and the constraint matrix is read straight from +// papilo's column-major (CSC) transpose into A. Metadata not produced by presolve (handle_ptr, +// obj_scale, problem_name, quadratic/cone fields) is left untouched on `problem`. +template +void build_user_problem(papilo::Problem const& papilo_problem, + simplex::user_problem_t& problem) +{ + raft::common::nvtx::range fun_scope("Build user problem"); + + const i_t reduced_rows = papilo_problem.getNRows(); + const i_t reduced_cols = papilo_problem.getNCols(); + auto const& constraint_matrix = papilo_problem.getConstraintMatrix(); + + // Objective (already minimization sense). + auto const& obj = papilo_problem.getObjective(); + problem.objective.assign(obj.coefficients.begin(), obj.coefficients.end()); + problem.obj_constant = obj.offset; + + // Column bounds and integrality. + auto const& col_lower = papilo_problem.getLowerBounds(); + auto const& col_upper = papilo_problem.getUpperBounds(); + auto const& col_flags = papilo_problem.getColFlags(); + problem.lower.resize(reduced_cols); + problem.upper.resize(reduced_cols); + problem.var_types.resize(reduced_cols); + for (i_t j = 0; j < reduced_cols; ++j) { + problem.lower[j] = col_flags[j].test(papilo::ColFlag::kLbInf) + ? -std::numeric_limits::infinity() + : col_lower[j]; + problem.upper[j] = col_flags[j].test(papilo::ColFlag::kUbInf) + ? std::numeric_limits::infinity() + : col_upper[j]; + problem.var_types[j] = col_flags[j].test(papilo::ColFlag::kIntegral) + ? simplex::variable_type_t::INTEGER + : simplex::variable_type_t::CONTINUOUS; + } + + // Row sense / rhs / ranges -- inverse of the derivation in build_papilo_problem_mip. + auto const& lhs = constraint_matrix.getLeftHandSides(); + auto const& rhs_v = constraint_matrix.getRightHandSides(); + auto const& row_flags = constraint_matrix.getRowFlags(); + problem.row_sense.assign(reduced_rows, 'E'); + problem.rhs.assign(reduced_rows, f_t{0}); + problem.range_rows.clear(); + problem.range_value.clear(); + for (i_t r = 0; r < reduced_rows; ++r) { + const bool lhs_inf = row_flags[r].test(papilo::RowFlag::kLhsInf); + const bool rhs_inf = row_flags[r].test(papilo::RowFlag::kRhsInf); + const bool eq = row_flags[r].test(papilo::RowFlag::kEquation); + if (eq || (!lhs_inf && !rhs_inf && lhs[r] == rhs_v[r])) { + problem.row_sense[r] = 'E'; + problem.rhs[r] = rhs_v[r]; + } else if (lhs_inf && !rhs_inf) { + problem.row_sense[r] = 'L'; + problem.rhs[r] = rhs_v[r]; + } else if (!lhs_inf && rhs_inf) { + problem.row_sense[r] = 'G'; + problem.rhs[r] = lhs[r]; + } else if (!lhs_inf && !rhs_inf) { + // lhs <= a^T x <= rhs : store as 'L' with a positive range width. + problem.row_sense[r] = 'L'; + problem.rhs[r] = rhs_v[r]; + problem.range_rows.push_back(r); + problem.range_value.push_back(rhs_v[r] - lhs[r]); + } else { + // Free row (both sides infinite). + problem.row_sense[r] = 'L'; + problem.rhs[r] = std::numeric_limits::infinity(); + } + } + problem.num_range_rows = static_cast(problem.range_rows.size()); + + // Constraint matrix: read papilo's column-major (CSC) transpose straight into A, packing out + // the spare gaps SparseStorage leaves between columns. + const i_t reduced_nnz = constraint_matrix.getNnz(); + problem.A.resize(reduced_rows, reduced_cols, reduced_nnz); + auto const& csc = constraint_matrix.getMatrixTranspose(); + const auto* col_ranges = csc.getRowRanges(); // per-column [start, end) + const int* row_indices = csc.getColumns(); // transpose columns == original rows + const f_t* values = csc.getValues(); + i_t pos = 0; + for (i_t j = 0; j < reduced_cols; ++j) { + problem.A.col_start[j] = pos; + for (i_t p = col_ranges[j].start; p < col_ranges[j].end; ++p) { + problem.A.i[pos] = row_indices[p]; + problem.A.x[pos] = values[p]; + ++pos; + } + } + problem.A.col_start[reduced_cols] = pos; + cuopt_assert(pos == reduced_nnz, "papilo CSC nonzero count mismatch"); + + // Names are not needed for the sub-MIP solve; clear so downstream size checks skip them. + problem.col_names.clear(); + problem.row_names.clear(); + problem.num_cols = reduced_cols; + problem.num_rows = reduced_rows; +} + void check_presolve_status(const papilo::PresolveStatus& status) { switch (status) { @@ -767,6 +992,86 @@ third_party_presolve_result_t third_party_presolve_t::apply( original_to_reduced_map_}; } +template +third_party_presolve_status_t third_party_presolve_t::apply( + simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings, + f_t time_limit, + i_t num_threads) +{ + const bool dual_postsolve = false; + presolver_ = Papilo; + // build_papilo_problem_mip keeps the objective in minimization sense (user_problem_t carries + // the direction in obj_scale), so the read-back must not flip signs either. + maximize_ = false; + + // Capture original dimensions before the problem is overwritten in place. + const i_t orig_cols = problem.num_cols; + const i_t orig_rows = problem.num_rows; + const i_t orig_nnz = problem.A.nnz(); + + papilo::Problem papilo_problem = build_papilo_problem(problem); + + settings.log.printf("Sub-MIP presolve input: %d constraints, %d variables, %d nonzeros", + papilo_problem.getNRows(), + papilo_problem.getNCols(), + papilo_problem.getConstraintMatrix().getNnz()); + + papilo::Presolve papilo_presolver; + set_presolve_methods(papilo_presolver, problem_category_t::MIP, dual_postsolve); + set_presolve_options(papilo_presolver, + problem_category_t::MIP, + settings.primal_tol, + settings.dual_tol, + time_limit, + dual_postsolve, + num_threads); + set_presolve_parameters(papilo_presolver, problem_category_t::MIP, orig_rows, orig_cols); + + // Disable papilo logs + papilo_presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet); + + auto result = papilo_presolver.apply(papilo_problem); + check_presolve_status(result.status); + 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); + + // 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, diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index 3b6b92a40d..72be05b5ec 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 @@ -68,6 +70,13 @@ class third_party_presolve_t { double time_limit, i_t num_cpu_threads = 0); + // Apply the presolve on an simplex::user_problem in-place. Used in sub MIP and (in the future) + // restarts. + third_party_presolve_status_t apply(simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings, + f_t time_limit, + i_t num_threads); + void undo(rmm::device_uvector& primal_solution, rmm::device_uvector& dual_solution, rmm::device_uvector& reduced_costs, @@ -92,9 +101,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}; diff --git a/cpp/tests/mip/presolve_test.cu b/cpp/tests/mip/presolve_test.cu index 5d82dcec0a..bb3528323d 100644 --- a/cpp/tests/mip/presolve_test.cu +++ b/cpp/tests/mip/presolve_test.cu @@ -7,11 +7,15 @@ #include "../linear_programming/utilities/pdlp_test_utilities.cuh" +#include #include #include #include +#include +#include #include #include +#include #include #include #include @@ -63,4 +67,48 @@ TEST(problem, find_implied_integers) ((int)mip::problem_t::var_flags_t::VAR_IMPLIED_INTEGER)); } +// 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::presolver_t presolver; + auto status = presolver.apply(user_problem, settings); + + // 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(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 diff --git a/datasets/mip/download_miplib_test_dataset.sh b/datasets/mip/download_miplib_test_dataset.sh index 0105729cb8..ccfacea1ef 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" From 6596b6e6d18933a1003008a030f4a01a34870137 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Fri, 19 Jun 2026 14:59:42 +0200 Subject: [PATCH 07/16] added routine for converting from standard form to range form --- cpp/src/dual_simplex/presolve.cpp | 149 ++++++++++++++++++++ cpp/src/dual_simplex/presolve.hpp | 13 ++ cpp/tests/dual_simplex/unit_tests/solve.cpp | 119 ++++++++++++++++ 3 files changed, 281 insertions(+) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 9c55185c6d..2fc27fbf98 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -676,6 +676,147 @@ i_t add_artifical_variables(lp_problem_t& problem, return 0; } +// Inverse of convert_user_problem (LP/MIP only -- no quadratic or SOC handling): +// take a problem in simplex standard form +// minimize c^T x +// subject to A x = b +// l <= x <= u +// and express it as a user_problem in range (row-bounded) form +// minimize c^T x +// subject to l_row <= A x <= u_row +// l <= x <= u +template +void convert_simplex_problem(const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const dualize_info_t& dualize_info, + const std::vector& new_slacks, + user_problem_t& user_problem) +{ + constexpr bool verbose = false; + if (verbose) { + settings.log.printf("Converting simplex problem with %d rows and %d columns and %d nonzeros\n", + simplex_problem.num_rows, + simplex_problem.num_cols, + simplex_problem.A.col_start[simplex_problem.num_cols]); + } + + // The simplex_problem is expected to be the primal in standard form. The + // dualize path of convert_user_problem hands the dual to the solver and + // stashes the primal in dualize_info.primal_problem, so a faithful inverse + // of a dualized problem must start from that primal instead. + assert(!dualize_info.solving_dual && + "convert_simplex_problem expects the primal problem; reconstruct the primal from " + "dualize_info.primal_problem before calling"); + + const i_t m = simplex_problem.num_rows; + const i_t n = simplex_problem.num_cols; + const csc_matrix_t& simplex_A = simplex_problem.A; + assert((var_types.empty() || static_cast(var_types.size()) == n) && + "var_types must span the full simplex problem (structural + slack columns)"); + + const i_t num_slacks = static_cast(new_slacks.size()); + const i_t new_n = n - num_slacks; + const i_t new_nnz = simplex_A.col_start[n] - num_slacks; + + // Reset the destination so we can populate user_problem in place. + user_problem.num_rows = m; + user_problem.num_cols = new_n; + user_problem.range_rows.clear(); + user_problem.range_value.clear(); + user_problem.objective.clear(); + user_problem.lower.clear(); + user_problem.upper.clear(); + user_problem.var_types.clear(); + user_problem.objective.reserve(new_n); + user_problem.lower.reserve(new_n); + user_problem.upper.reserve(new_n); + if (!var_types.empty()) { user_problem.var_types.reserve(new_n); } + + // Rows with no associated slack stay equalities at their rhs; the slack loop + // below overrides the rows that had one. + user_problem.row_sense.assign(m, 'E'); + user_problem.rhs = simplex_problem.rhs; + + // Recover each row's constraint from its slack. The forward equality + // a_i^T x_struct + c * s = b_i (with l_s <= s <= u_s) gives + // a_i^T x_struct in [min(b - c l_s, b - c u_s), max(b - c l_s, b - c u_s)], + // which maps back to row_sense / rhs / range. Flag the slack columns so they + // can be dropped from the matrix below. + std::vector is_slack(n, false); + for (i_t s : new_slacks) { + assert(s >= 0 && s < n && "slack index out of range"); + is_slack[s] = true; + + const i_t col_start = simplex_A.col_start[s]; + const i_t col_end = simplex_A.col_start[s + 1]; + assert(col_end - col_start == 1 && "slack/artificial column must have a single nonzero"); + const i_t row = simplex_A.i[col_start]; + const f_t c = simplex_A.x[col_start]; + const f_t b = simplex_problem.rhs[row]; + const f_t t1 = b - c * simplex_problem.lower[s]; + const f_t t2 = b - c * simplex_problem.upper[s]; + const f_t lo = std::min(t1, t2); + const f_t hi = std::max(t1, t2); + if (lo == hi) { + user_problem.row_sense[row] = 'E'; + user_problem.rhs[row] = lo; + } else if (lo == -inf) { + user_problem.row_sense[row] = 'L'; + user_problem.rhs[row] = hi; + } else if (hi == inf) { + user_problem.row_sense[row] = 'G'; + user_problem.rhs[row] = lo; + } else { + user_problem.row_sense[row] = 'E'; + user_problem.rhs[row] = lo; + user_problem.range_rows.push_back(row); + user_problem.range_value.push_back(hi - lo); + } + } + user_problem.num_range_rows = static_cast(user_problem.range_rows.size()); + + // Rebuild the constraint matrix and column data in place, dropping the + // slack/artificial columns (each contributes exactly one nonzero). var_types + // is filtered to the kept (structural) columns in the same pass. + csc_matrix_t& user_A = user_problem.A; + user_A.m = m; + user_A.n = new_n; + user_A.nz_max = new_nnz; + user_A.col_start.assign(new_n + 1, 0); + user_A.i.assign(new_nnz, 0); + user_A.x.assign(new_nnz, 0); + + i_t nz = 0; + i_t new_j = 0; + for (i_t j = 0; j < n; ++j) { + if (is_slack[j]) { continue; } + user_A.col_start[new_j] = nz; + for (i_t p = simplex_A.col_start[j]; p < simplex_A.col_start[j + 1]; ++p) { + user_A.i[nz] = simplex_A.i[p]; + user_A.x[nz] = simplex_A.x[p]; + ++nz; + } + user_problem.objective.push_back(simplex_problem.objective[j]); + user_problem.lower.push_back(simplex_problem.lower[j]); + user_problem.upper.push_back(simplex_problem.upper[j]); + if (!var_types.empty()) { user_problem.var_types.push_back(var_types[j]); } + ++new_j; + } + user_A.col_start[new_n] = nz; + assert(new_j == new_n); + assert(nz == new_nnz); + + user_problem.obj_scale = simplex_problem.obj_scale; + user_problem.obj_constant = simplex_problem.obj_constant; + user_problem.objective_is_integral = simplex_problem.objective_is_integral; + user_problem.objective_step = simplex_problem.objective_step; + + // LP/MIP only: no quadratic objective and no second-order cones. + user_problem.cone_var_start = 0; + user_problem.second_order_cone_dims.clear(); +} + template void convert_user_problem(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, @@ -1892,6 +2033,14 @@ template void convert_user_problem( std::vector& new_slacks, dualize_info_t& dualize_info); +template void convert_simplex_problem( + const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const dualize_info_t& dualize_info, + const std::vector& new_slacks, + user_problem_t& user_problem); + template void convert_user_lp_with_guess( const user_problem_t& user_problem, const lp_solution_t& initial_solution, diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 855a4e6145..10ceee4c4d 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -205,6 +205,19 @@ void convert_user_problem(const user_problem_t& user_problem, std::vector& new_slacks, dualize_info_t& dualize_info); +// Inverse of convert_user_problem (LP/MIP only): reconstruct a range-form +// user_problem from a simplex standard-form lp_problem, dropping the +// slack/artificial columns recorded in new_slacks. var_types spans the full +// simplex problem (structural + slack columns) and is filtered to the +// structural columns. +template +void convert_simplex_problem(const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const dualize_info_t& dualize_info, + const std::vector& new_slacks, + user_problem_t& user_problem); + template void convert_user_problem_with_guess(const user_problem_t& user_problem, const std::vector& guess, diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 33d5392a9b..ad74069b44 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp @@ -332,4 +332,123 @@ TEST(dual_simplex, dual_variable_greater_than) EXPECT_NEAR(solution.z[1], 0.0, 1e-6); } +// Round-trip a MIP through convert_user_problem (range form -> simplex standard +// form, appending one slack/artificial column per row) and +// convert_simplex_problem (the inverse: drop the slacks and recover the row +// bounds). The problem exercises every row type: '<=', '==', '>=', and a range +// row. +// +// '>=' rows are folded into '<=' rows by negating their coefficients/rhs in the +// forward pass, so the inverse recovers them as the equivalent negated '<=' row +// rather than the original '>='. The recovered problem is therefore feasibly +// equivalent but not textually identical. We assert two things: +// 1. the directly-predictable recovered fields (sizes, row_sense, rhs, range, +// objective, bounds, var_types), and +// 2. the round-trip invariant: re-running the forward conversion on the +// recovered problem reproduces the original standard-form problem exactly. +TEST(dual_simplex, convert_simplex_problem_mip_round_trip) +{ + cuopt::init_logger_t log("", true); + raft::handle_t handle{}; + + // minimize x0 + 2 x1 + 3 x2 + // subject to 2 x0 + x1 <= 8 (row 0, 'L') + // x0 + x2 = 4 (row 1, 'E') + // x1 + x2 >= 3 (row 2, 'G') + // 1 <= x0 + x1 <= 7 (row 3, range row stored as 'E' + range) + // x0 integer in [0, 10] + // x1 continuous in [0, inf) + // x2 integer in [0, 5] + simplex::user_problem_t original(&handle); + constexpr int m = 4; + constexpr int n = 3; + constexpr int nz = 8; + + original.num_rows = m; + original.num_cols = n; + original.objective.assign({1.0, 2.0, 3.0}); + original.A.m = m; + original.A.n = n; + original.A.nz_max = nz; + original.A.reallocate(nz); + // column 0 (x0): rows {0, 1, 3} -> {2, 1, 1} + // column 1 (x1): rows {0, 2, 3} -> {1, 1, 1} + // column 2 (x2): rows {1, 2} -> {1, 1} + original.A.col_start.assign({0, 3, 6, 8}); + original.A.i = {0, 1, 3, 0, 2, 3, 1, 2}; + original.A.x = {2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + original.rhs.assign({8.0, 4.0, 3.0, 1.0}); + original.row_sense.assign({'L', 'E', 'G', 'E'}); + original.lower.assign({0.0, 0.0, 0.0}); + original.upper.assign({10.0, dual_simplex::inf, 5.0}); + // row 3 is the range row 1 <= a^T x <= 7: stored as an 'E' range row, which + // convert_range_rows reads as [rhs, rhs + range] = [1, 7]. + original.range_rows.assign({3}); + original.range_value.assign({6.0}); + original.num_range_rows = 1; + original.obj_constant = 0.0; + original.var_types = {dual_simplex::variable_type_t::INTEGER, + dual_simplex::variable_type_t::CONTINUOUS, + dual_simplex::variable_type_t::INTEGER}; + + // Forward: range form -> simplex standard form. + simplex::simplex_solver_settings_t settings; + simplex::lp_problem_t simplex_problem( + &handle, original.num_rows, original.num_cols, original.A.col_start[original.A.n]); + std::vector new_slacks; + simplex::dualize_info_t dualize_info; + simplex::convert_user_problem(original, settings, simplex_problem, new_slacks, dualize_info); + + // Each row gets exactly one slack/artificial column, appended after the + // structural columns. + EXPECT_EQ(new_slacks.size(), static_cast(m)); + EXPECT_EQ(simplex_problem.num_cols, n + m); + + // var_types spans the full simplex problem; the appended columns are + // continuous (mirrors full_variable_types in branch_and_bound). + std::vector var_types = original.var_types; + var_types.resize(simplex_problem.num_cols, simplex::variable_type_t::CONTINUOUS); + + // Inverse: simplex standard form -> range form. + simplex::user_problem_t recovered(&handle); + simplex::convert_simplex_problem( + simplex_problem, var_types, settings, dualize_info, new_slacks, recovered); + + // (1) Directly-predictable recovered fields. + EXPECT_EQ(recovered.num_rows, m); + EXPECT_EQ(recovered.num_cols, n); + // row 0 '<=' -> 'L' rhs 8; row 1 '==' -> 'E' rhs 4; + // row 2 '>=' -> recovered as the negated 'L' (-x1 - x2 <= -3); + // row 3 range [1, 7] -> canonical 'E' range row: rhs 1 with range 6. + EXPECT_EQ(recovered.row_sense, (std::vector{'L', 'E', 'L', 'E'})); + EXPECT_EQ(recovered.rhs, (std::vector{8.0, 4.0, -3.0, 1.0})); + EXPECT_EQ(recovered.num_range_rows, 1); + EXPECT_EQ(recovered.range_rows, (std::vector{3})); + EXPECT_EQ(recovered.range_value, (std::vector{6.0})); + // Column data (objective / bounds / types) carries over unchanged. + EXPECT_EQ(recovered.objective, original.objective); + EXPECT_EQ(recovered.lower, original.lower); + EXPECT_EQ(recovered.upper, original.upper); + EXPECT_EQ(recovered.var_types, original.var_types); + + // (2) Round-trip invariant: converting the recovered problem forward again + // must reproduce the original standard-form problem exactly. + simplex::lp_problem_t simplex_again( + &handle, recovered.num_rows, recovered.num_cols, recovered.A.col_start[recovered.A.n]); + std::vector new_slacks_again; + simplex::dualize_info_t dualize_info_again; + simplex::convert_user_problem( + recovered, settings, simplex_again, new_slacks_again, dualize_info_again); + + EXPECT_EQ(simplex_again.num_rows, simplex_problem.num_rows); + EXPECT_EQ(simplex_again.num_cols, simplex_problem.num_cols); + EXPECT_EQ(simplex_again.A.col_start, simplex_problem.A.col_start); + EXPECT_EQ(simplex_again.A.i, simplex_problem.A.i); + EXPECT_EQ(simplex_again.A.x, simplex_problem.A.x); + EXPECT_EQ(simplex_again.rhs, simplex_problem.rhs); + EXPECT_EQ(simplex_again.lower, simplex_problem.lower); + EXPECT_EQ(simplex_again.upper, simplex_problem.upper); + EXPECT_EQ(new_slacks_again, new_slacks); +} + } // namespace cuopt::mathematical_optimization::simplex::test From d39e42454933efae5d94d31a5538da2ba962d83a Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 30 Jun 2026 14:19:50 +0200 Subject: [PATCH 08/16] apply presolver before solving the submip. missing warm start. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 76 ++++++++++++++----- cpp/src/branch_and_bound/presolve.cpp | 19 ++++- cpp/src/branch_and_bound/presolve.hpp | 6 +- cpp/src/dual_simplex/presolve.cpp | 10 --- cpp/src/dual_simplex/presolve.hpp | 1 - .../presolve/third_party_presolve.cpp | 1 - 6 files changed, 76 insertions(+), 37 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a5602523c9..2c86f732c1 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -44,6 +44,8 @@ #include #include +#include "presolve.hpp" + namespace cuopt::mathematical_optimization::mip { using simplex::basis_update_mpf_t; @@ -915,7 +917,8 @@ 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; @@ -2168,7 +2171,7 @@ void branch_and_bound_t::launch_submip_worker(const std::vector& submip_worker->set_active(); #pragma omp task priority(CUOPT_MEDIUM_TASK_PRIORITY) affinity(submip_worker) \ - firstprivate(submip_worker, sol) + firstprivate(submip_worker, sol) if (!settings_.inside_submip) rins(submip_worker, sol); } @@ -2211,6 +2214,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.reliability_branching = 0; submip_settings.max_cut_passes = 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; @@ -2220,30 +2224,59 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.submip_settings.enable_rins = settings_.submip_settings.enable_rins != 0 && submip_level <= settings_.submip_settings.max_level; - submip_settings.solution_callback = [this, fixrate](std::vector& solution, f_t objective) { - this->set_solution_from_submip(solution, fixrate); + user_problem_t submip_problem(original_problem_.handle_ptr); + simplex::convert_simplex_problem( + worker->leaf_problem, var_types_, settings_, new_slacks_, submip_problem); + + presolver_t presolver; + mip_status_t submip_status = presolver.apply(submip_problem, settings_); + + if (submip_status == mip_status_t::INFEASIBLE) { + submip_stats_.save_infeasible(fixrate); + return; + } + + if (submip_status == mip_status_t::OPTIMAL) { + std::vector fixed_sol(original_problem_.num_cols); + std::vector reduced_sol(submip_problem.num_cols); + presolver.uncrush(reduced_sol, fixed_sol); + set_solution_from_submip(fixed_sol, fixrate); + return; + } + + submip_settings.solution_callback = [this, fixrate, &presolver](std::vector& solution, + f_t obj) { + std::vector fixed_sol; + presolver.uncrush(solution, fixed_sol); + this->set_solution_from_submip(fixed_sol, fixrate); }; - branch_and_bound_t submip_bnb(*this, submip_settings, lower, upper); - mip_solution_t fixed_solution(original_problem_.num_cols); - mip_status_t submip_status = submip_bnb.solve(fixed_solution); + submip_settings.log.print_format("Sub-MIP: {} constraints, {} variables, {} nonzeros\n", + submip_problem.num_rows, + submip_problem.num_cols, + submip_problem.A.nnz()); + + probing_implied_bound_t empty_probing(submip_problem.num_cols); + branch_and_bound_t submip_bnb(submip_problem, submip_settings, tic(), empty_probing); + mip_solution_t submip_solution(submip_problem.num_cols); + + mutex_upper_.lock(); + std::vector current_incumbent = incumbent_.x; + mutex_upper_.unlock(); + + // submip_bnb.set_initial_guess(current_incumbent); + submip_bnb.set_initial_upper_bound(upper_bound_.load()); + submip_status = submip_bnb.solve(submip_solution); if (submip_status == mip_status_t::INFEASIBLE) { submip_stats_.save_infeasible(fixrate); + return; + } - } else if (submip_status == mip_status_t::OPTIMAL) { - std::vector crushed_solution; - crush_primal_solution( - original_problem_, original_lp_, fixed_solution.x, new_slacks_, crushed_solution); - f_t obj = compute_objective(original_lp_, crushed_solution); - settings_.log.debug_format("{}Found a solution with obj={:.4g}\n", - log_prefix, - compute_user_objective(original_lp_, obj)); - - if (improves_incumbent(obj)) { - add_feasible_solution(obj, crushed_solution, -1, SUBMIP); - submip_stats_.save_success(fixrate); - } + if (submip_solution.has_incumbent) { + std::vector fixed_sol(original_problem_.num_cols); + presolver.uncrush(submip_solution.x, fixed_sol); + set_solution_from_submip(fixed_sol, fixrate); } } @@ -2950,7 +2983,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut 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_, + root_status = solve_linear_program_with_advanced_basis(original_lp_, exploration_stats_.start_time, lp_settings, root_relax_soln_, @@ -2959,6 +2992,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nonbasic_list, root_vstatus_, edge_norms_); + root_relax_solved_by = DualSimplex; } else { settings_.log.printf("\nSolving LP root relaxation in concurrent mode\n"); diff --git a/cpp/src/branch_and_bound/presolve.cpp b/cpp/src/branch_and_bound/presolve.cpp index a5ba8a2f30..30d1f69336 100644 --- a/cpp/src/branch_and_bound/presolve.cpp +++ b/cpp/src/branch_and_bound/presolve.cpp @@ -12,13 +12,28 @@ namespace cuopt::mathematical_optimization::mip { +mip_status_t presolve_status_to_mip_status(third_party_presolve_status_t status) +{ + switch (status) { + case third_party_presolve_status_t::OPTIMAL: return mip_status_t::OPTIMAL; + case third_party_presolve_status_t::INFEASIBLE: return mip_status_t::INFEASIBLE; + case third_party_presolve_status_t::UNBOUNDED: return mip_status_t::UNBOUNDED; + case third_party_presolve_status_t::UNBNDORINFEAS: return mip_status_t::INFEASIBLE; + case third_party_presolve_status_t::REDUCED: return mip_status_t::UNSET; + case third_party_presolve_status_t::UNCHANGED: return mip_status_t::UNSET; + } + return mip_status_t::UNSET; +} + template -third_party_presolve_status_t presolver_t::apply( +mip_status_t presolver_t::apply( simplex::user_problem_t& problem, const simplex::simplex_solver_settings_t& settings) { f_t presolve_time_limit = std::min(0.1 * settings.time_limit, 60.0); - return third_party_presolver_.apply(problem, settings, presolve_time_limit, 1); + third_party_presolve_status_t status = + third_party_presolver_.apply(problem, settings, presolve_time_limit, 1); + return presolve_status_to_mip_status(status); } template diff --git a/cpp/src/branch_and_bound/presolve.hpp b/cpp/src/branch_and_bound/presolve.hpp index 4ee66dc21e..7560e3ce7e 100644 --- a/cpp/src/branch_and_bound/presolve.hpp +++ b/cpp/src/branch_and_bound/presolve.hpp @@ -14,6 +14,8 @@ #include #include +#include "constants.hpp" + namespace cuopt::mathematical_optimization::mip { // Thin owner of a PaPILO presolver scoped to a single sub-MIP solve. apply() reduces the @@ -24,8 +26,8 @@ class presolver_t { public: // Presolve `problem` in place using PaPILO. Returns the presolve status; on // INFEASIBLE/UNBOUNDED the problem is left untouched. - third_party_presolve_status_t apply(simplex::user_problem_t& problem, - const simplex::simplex_solver_settings_t& settings); + mip_status_t apply(simplex::user_problem_t& problem, + const simplex::simplex_solver_settings_t& settings); // Map a reduced-space primal solution back to the original column space. void uncrush(const std::vector& reduced_primal, std::vector& full_primal) const; diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 2fc27fbf98..55192415ab 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -689,7 +689,6 @@ template void convert_simplex_problem(const lp_problem_t& simplex_problem, const std::vector& var_types, const simplex_solver_settings_t& settings, - const dualize_info_t& dualize_info, const std::vector& new_slacks, user_problem_t& user_problem) { @@ -701,14 +700,6 @@ void convert_simplex_problem(const lp_problem_t& simplex_problem, simplex_problem.A.col_start[simplex_problem.num_cols]); } - // The simplex_problem is expected to be the primal in standard form. The - // dualize path of convert_user_problem hands the dual to the solver and - // stashes the primal in dualize_info.primal_problem, so a faithful inverse - // of a dualized problem must start from that primal instead. - assert(!dualize_info.solving_dual && - "convert_simplex_problem expects the primal problem; reconstruct the primal from " - "dualize_info.primal_problem before calling"); - const i_t m = simplex_problem.num_rows; const i_t n = simplex_problem.num_cols; const csc_matrix_t& simplex_A = simplex_problem.A; @@ -2037,7 +2028,6 @@ template void convert_simplex_problem( const lp_problem_t& simplex_problem, const std::vector& var_types, const simplex_solver_settings_t& settings, - const dualize_info_t& dualize_info, const std::vector& new_slacks, user_problem_t& user_problem); diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 10ceee4c4d..2b70273049 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -214,7 +214,6 @@ template void convert_simplex_problem(const lp_problem_t& simplex_problem, const std::vector& var_types, const simplex_solver_settings_t& settings, - const dualize_info_t& dualize_info, const std::vector& new_slacks, user_problem_t& user_problem); diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 31602a3a12..25be9af708 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -1032,7 +1032,6 @@ third_party_presolve_status_t third_party_presolve_t::apply( papilo_presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet); auto result = papilo_presolver.apply(papilo_problem); - check_presolve_status(result.status); auto status = convert_papilo_presolve_status_to_third_party_presolve_status(result.status); // Infeasible / unbounded: leave `problem` untouched; the caller branches on the status. From 3328cda2cb4f3acda4e9278f79e1a965edaa1b10 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 30 Jun 2026 16:46:01 +0200 Subject: [PATCH 09/16] pass the upper bound from the parent B&B to the children. clean up the old warm start code. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 111 ++++-------------- cpp/src/branch_and_bound/branch_and_bound.hpp | 5 + .../dual_simplex/simplex_solver_settings.hpp | 12 +- .../presolve/third_party_presolve.hpp | 28 +++-- 4 files changed, 57 insertions(+), 99 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 2c86f732c1..40ab635e30 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -317,48 +317,6 @@ branch_and_bound_t::branch_and_bound_t( root_lp_current_lower_bound_ = -inf; } -template -branch_and_bound_t::branch_and_bound_t( - branch_and_bound_t& other, - const simplex_solver_settings_t& solver_settings, - const std::vector& lower, - const std::vector& upper) - : original_problem_(other.original_problem_), - settings_(solver_settings), - probing_implied_bound_(other.probing_implied_bound_), - clique_table_(other.clique_table_), - symmetry_(nullptr), // Enabling symmetry in the submip lead to crashes in some instances - original_lp_(other.original_lp_), - Arow_(other.Arow_), - new_slacks_(other.new_slacks_), - var_types_(other.var_types_), - incumbent_(1), - root_vstatus_(other.root_vstatus_), - root_relax_soln_(1, 1), - root_crossover_soln_(1, 1), - pc_(other.pc_), - solver_status_(mip_status_t::UNSET) -{ - other.mutex_original_lp_.lock(); - assert(lower.size() == original_lp_.num_cols); - assert(upper.size() == original_lp_.num_cols); - original_lp_.lower = lower; - original_lp_.upper = upper; - other.mutex_original_lp_.unlock(); - - other.mutex_upper_.lock(); - incumbent_ = other.incumbent_; - other.mutex_upper_.unlock(); - - exploration_stats_.start_time = tic(); - upper_bound_ = other.upper_bound_; - root_objective_ = std::numeric_limits::quiet_NaN(); - root_lp_current_lower_bound_ = -inf; - external_upper_bound_ = - other.external_upper_bound_ ? other.external_upper_bound_ : &other.upper_bound_; - root_warm_start_ = true; -} - template f_t branch_and_bound_t::get_lower_bound() { @@ -2266,6 +2224,8 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke // submip_bnb.set_initial_guess(current_incumbent); submip_bnb.set_initial_upper_bound(upper_bound_.load()); + submip_bnb.set_external_upper_bound(external_upper_bound_ ? external_upper_bound_ + : &upper_bound_); submip_status = submip_bnb.solve(submip_solution); if (submip_status == mip_status_t::INFEASIBLE) { @@ -2958,52 +2918,29 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t root_relax_start_time = tic(); - if (root_warm_start_) { - i_t node_iter = 0; - f_t lp_start_time = tic(); - - dual_status_t lp_status = dual_phase2_with_advanced_basis(2, - 0, - true, - lp_start_time, - original_lp_, - lp_settings, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - root_relax_soln_, - node_iter, - edge_norms_); - - root_status = convert_dual_status_to_lp_status(lp_status); - } - - if (root_status == lp_status_t::NUMERICAL_ISSUES || root_status == lp_status_t::UNSET) { - 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_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - root_relax_solved_by = DualSimplex; + 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_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + root_relax_solved_by = DualSimplex; - } else { - settings_.log.printf("\nSolving LP root relaxation in concurrent mode\n"); - root_status = solve_root_relaxation(lp_settings, - root_relax_soln_, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - edge_norms_); - } + } else { + settings_.log.printf("\nSolving LP root relaxation in concurrent mode\n"); + root_status = solve_root_relaxation(lp_settings, + root_relax_soln_, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + edge_norms_); } solving_root_relaxation_ = false; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 608fd40862..3f8f3147b9 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -79,6 +79,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_external_upper_bound(const omp_atomic_t* external_upper_bound) + { + external_upper_bound_ = external_upper_bound; + } + // Set the root solution found by PDLP void set_root_relaxation_solution(const std::vector& primal, const std::vector& dual, diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 943cac5ce0..7aae2338b4 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -37,13 +37,13 @@ struct submip_settings_t { // Minimum fix rate for accepting the RINS neighbourhood. double min_fixrate = 0.25; - // Hard cap for the minimum fix rate for solving a submip. + // Hard cap for the minimum fix rate for solving a sub-MIP. double min_fixrate_cap = 0.1; - // MIP gap for the submip (unless the MIP gap from the B&B is lower) + // MIP gap for the sub-MIP (unless the MIP gap from the B&B is lower) double target_mip_gap = 0.01; - // The base node limit for the submip + // The base node limit for the sub-MIP int node_limit_base = 200; // The current level in the recursion @@ -51,6 +51,12 @@ struct submip_settings_t { // Maximum recursion level int max_level = 5; + + // Presolve sub-MIP with Papilo before solving it + bool presolve = true; + + // Warm start sub-MIP's B&B with the parent basis and pseudocost + bool warm_start = true; }; template diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index 72be05b5ec..ddf00c83c3 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -60,15 +60,14 @@ 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. @@ -89,6 +88,15 @@ class third_party_presolve_t { std::vector& full_primal) const; const std::vector& get_reduced_to_original_map() const { return reduced_to_original_map_; } const std::vector& get_original_to_reduced_map() const { return original_to_reduced_map_; } + const std::vector& get_reduced_to_original_row_map() const + { + return reduced_to_original_row_map_; + } + + const std::vector& get_original_to_reduced_row_map() const + { + return original_to_reduced_row_map_; + } ~third_party_presolve_t(); @@ -115,6 +123,8 @@ class third_party_presolve_t { std::vector reduced_to_original_map_{}; std::vector original_to_reduced_map_{}; + std::vector reduced_to_original_row_map_{}; + std::vector original_to_reduced_row_map_{}; }; } // namespace cuopt::mathematical_optimization::mip From a37acd4e2867a450f00290708ef75b22c4f0f224 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 30 Jun 2026 17:24:38 +0200 Subject: [PATCH 10/16] crush the incumbent and pass to the submip B&B. includes #1104 Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 10 +- cpp/src/branch_and_bound/presolve.cpp | 6 + cpp/src/branch_and_bound/presolve.hpp | 2 + cpp/src/branch_and_bound/submip.hpp | 2 +- .../diversity/diversity_manager.cu | 43 +- .../presolve/third_party_presolve.cpp | 260 ++++++++- .../presolve/third_party_presolve.hpp | 24 +- cpp/src/mip_heuristics/solve.cu | 27 +- .../unit_tests/presolve_test.cu | 516 ++++++++++++++++++ cpp/tests/mip/incumbent_callback_test.cu | 76 ++- datasets/mip/download_miplib_test_dataset.sh | 31 ++ 11 files changed, 966 insertions(+), 31 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 681091b6b3..e87abba3cd 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2181,7 +2181,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke worker->leaf_problem, var_types_, settings_, new_slacks_, submip_problem); presolver_t presolver; - mip_status_t submip_status = presolver.apply(submip_problem, settings_); + mip_status_t submip_status = presolver.apply(submip_problem, submip_settings); if (submip_status == mip_status_t::INFEASIBLE) { submip_stats_.save_infeasible(fixrate); @@ -2216,8 +2216,10 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke std::vector current_incumbent = incumbent_.x; mutex_upper_.unlock(); - // submip_bnb.set_initial_guess(current_incumbent); - submip_bnb.set_initial_upper_bound(upper_bound_.load()); + std::vector crushed_incumbent; + presolver.crush(current_incumbent, crushed_incumbent); + + submip_bnb.set_initial_guess(crushed_incumbent); submip_bnb.set_external_upper_bound(external_upper_bound_ ? external_upper_bound_ : &upper_bound_); submip_status = submip_bnb.solve(submip_solution); @@ -2868,6 +2870,8 @@ 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.debug_format("Warm start with incumbent, obj={:.4g}", computed_obj); } } diff --git a/cpp/src/branch_and_bound/presolve.cpp b/cpp/src/branch_and_bound/presolve.cpp index 30d1f69336..7375885122 100644 --- a/cpp/src/branch_and_bound/presolve.cpp +++ b/cpp/src/branch_and_bound/presolve.cpp @@ -42,6 +42,12 @@ void presolver_t::uncrush(const std::vector& reduced_primal, { third_party_presolver_.uncrush_primal_solution(reduced_primal, full_primal); } +template +void presolver_t::crush(const std::vector& full_primal, + std::vector& reduced_primal) const +{ + third_party_presolver_.crush_primal_solution(full_primal, reduced_primal); +} template const std::vector& presolver_t::reduced_to_original_map() const diff --git a/cpp/src/branch_and_bound/presolve.hpp b/cpp/src/branch_and_bound/presolve.hpp index 7560e3ce7e..2efde2732f 100644 --- a/cpp/src/branch_and_bound/presolve.hpp +++ b/cpp/src/branch_and_bound/presolve.hpp @@ -32,6 +32,8 @@ class presolver_t { // Map a reduced-space primal solution back to the original column space. void uncrush(const std::vector& reduced_primal, std::vector& full_primal) const; + void crush(const std::vector& full_primal, std::vector& reduced_primal) const; + const std::vector& reduced_to_original_map() const; const std::vector& original_to_reduced_map() const; diff --git a/cpp/src/branch_and_bound/submip.hpp b/cpp/src/branch_and_bound/submip.hpp index a528ce5f0a..72ac7b8537 100644 --- a/cpp/src/branch_and_bound/submip.hpp +++ b/cpp/src/branch_and_bound/submip.hpp @@ -69,7 +69,7 @@ class submip_worker_t : public branch_and_bound_worker_t { submip_worker_t(i_t worker_id, const simplex::lp_problem_t& original_lp, - const simplex::csr_matrix_t& Arow, + const csr_matrix_t& Arow, const std::vector& var_type, const simplex::simplex_solver_settings_t& settings, uint64_t rng_offset = 0) diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index 4bb35bc35d..c48358fa72 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -9,6 +9,7 @@ #include "diversity_manager.cuh" #include +#include #include #include @@ -182,9 +183,33 @@ void diversity_manager_t::add_user_given_solutions( std::vector>& initial_sol_vector) { raft::common::nvtx::range fun_scope("add_user_given_solutions"); - for (const auto& init_sol : context.settings.initial_solutions) { + const bool has_papilo = problem_ptr->has_papilo_presolve_data(); + const i_t papilo_orig_n = problem_ptr->get_papilo_original_num_variables(); + for (size_t sol_idx = 0; sol_idx < context.settings.initial_solutions.size(); ++sol_idx) { + const auto& init_sol = context.settings.initial_solutions[sol_idx]; solution_t sol(*problem_ptr); rmm::device_uvector init_sol_assignment(*init_sol, sol.handle_ptr->get_stream()); + + 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, " + "expected %d; skipping", + sol_idx, + init_sol_assignment.size(), + papilo_orig_n); + continue; + } + std::vector h_original = host_copy(init_sol_assignment, sol.handle_ptr->get_stream()); + std::vector h_crushed; + problem_ptr->presolve_data.papilo_presolve_ptr->crush_primal_solution(h_original, h_crushed); + init_sol_assignment = cuopt::device_copy(h_crushed, sol.handle_ptr->get_stream()); + CUOPT_LOG_DEBUG("Crushed initial solution %d through Papilo (%d -> %d vars)", + sol_idx, + papilo_orig_n, + h_crushed.size()); + } + if (problem_ptr->pre_process_assignment(init_sol_assignment)) { relaxed_lp_settings_t lp_settings; lp_settings.time_limit = std::min(60., timer.remaining_time() / 2); @@ -202,10 +227,10 @@ void diversity_manager_t::add_user_given_solutions( sol.handle_ptr->get_stream()); bool is_feasible = sol.compute_feasibility(); cuopt_func_call(sol.test_variable_bounds(true)); - CUOPT_LOG_INFO("Adding initial solution success! feas %d objective %f excess %f", - is_feasible, - sol.get_user_objective(), - sol.get_total_excess()); + CUOPT_LOG_DEBUG("Adding initial solution success! feas %d objective %f excess %f", + is_feasible, + sol.get_user_objective(), + sol.get_total_excess()); population.run_solution_callbacks(sol); initial_sol_vector.emplace_back(std::move(sol)); } else { @@ -424,6 +449,9 @@ solution_t diversity_manager_t::run_solver() CUOPT_LOG_INFO("GPU heuristics disabled via CUOPT_DISABLE_GPU_HEURISTICS=1"); population.initialize_population(); population.allocate_solutions(); + add_user_given_solutions(initial_sol_vector); + population.add_solutions_from_vec(std::move(initial_sol_vector)); + if (check_b_b_preemption()) { return population.best_feasible(); } while (!check_b_b_preemption()) { std::this_thread::sleep_for(std::chrono::milliseconds(100)); @@ -454,8 +482,9 @@ solution_t diversity_manager_t::run_solver() "The problem must not be ii"); population.initialize_population(); population.allocate_solutions(); - if (check_b_b_preemption()) { return population.best_feasible(); } add_user_given_solutions(initial_sol_vector); + population.add_solutions_from_vec(std::move(initial_sol_vector)); + if (check_b_b_preemption()) { return population.best_feasible(); } // Run CPUFJ early to find quick initial solutions ls_cpufj_raii_guard_t ls_cpufj_raii_guard(ls); // RAII to stop cpufj threads on solve stop ls.start_cpufj_scratch_threads(population); @@ -612,8 +641,6 @@ solution_t diversity_manager_t::run_solver() ls.start_cpufj_lptopt_scratch_threads(population); } - population.add_solutions_from_vec(std::move(initial_sol_vector)); - if (check_b_b_preemption()) { return population.best_feasible(); } if (context.settings.benchmark_info_ptr != nullptr) { diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 25be9af708..dc06a0dd1c 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -44,6 +44,8 @@ #include +#include + namespace cuopt::mathematical_optimization::mip { template @@ -137,8 +139,10 @@ papilo::Problem build_papilo_problem(const optimization_problem_t } } - for (size_t i = 0; i < h_var_types.size(); ++i) { - builder.setColIntegral(i, h_var_types[i] == var_t::INTEGER); + if (category == problem_category_t::MIP) { + for (size_t i = 0; i < h_var_types.size(); ++i) { + builder.setColIntegral(i, h_var_types[i] == var_t::INTEGER); + } } if (!h_constr_lb.empty() && !h_constr_ub.empty()) { @@ -788,8 +792,6 @@ void set_presolve_methods(papilo::Presolve& presolver, presolver.addPresolveMethod(uptr(new papilo::SimpleProbing())); presolver.addPresolveMethod(uptr(new papilo::ParallelRowDetection())); presolver.addPresolveMethod(uptr(new papilo::ParallelColDetection())); - - presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing())); presolver.addPresolveMethod(uptr(new papilo::DualFix())); presolver.addPresolveMethod(uptr(new papilo::SimplifyInequalities())); presolver.addPresolveMethod(uptr(new papilo::CliqueMerging())); @@ -800,6 +802,11 @@ void set_presolve_methods(papilo::Presolve& presolver, presolver.addPresolveMethod(uptr(new papilo::Probing())); if (!dual_postsolve) { + // SingletonStuffing causes dual crushing failures on: + // tr12-30, ns1208400, gmu-35-50, dws008-01, neos-1445765, + // neos-5107597-kakapo, rocI-4-11, traininstance2, traininstance6, + // radiationm18-12-05, rococoB10-011000, b1c1s1 + presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing())); presolver.addPresolveMethod(uptr(new papilo::DualInfer())); presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution())); presolver.addPresolveMethod(uptr(new papilo::Sparsify())); @@ -1181,6 +1188,251 @@ void third_party_presolve_t::uncrush_primal_solution( full_primal = std::move(full_sol.primal); } +template +void third_party_presolve_t::crush_primal_solution( + const std::vector& original_primal, std::vector& reduced_primal) const +{ + cuopt_expects(presolver_ == cuopt::mathematical_optimization::presolver_t::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); +} + +/** + * Crush an original-space primal+dual solution into the presolved (reduced) space. + * + * This is the forward counterpart of Papilo's Postsolve::undo(). It replays + * each presolve reduction in forward order to transform variable/dual values, + * then projects onto the surviving columns/rows via origcol/origrow_mapping. + * + * Only two reductions actually transform survivor coordinates: + * kParallelCol — merges x[col1] into x[col2]; survivor rc is z[col2] if + * nonzero, else z[col1] / scale (inverse of PaPILO postsolve) + * kRowBoundChangeForcedByRow — conditionally transfers y[deleted_row] → y[kept_row] + */ +template +void third_party_presolve_t::crush_primal_dual_solution( + const std::vector& x_original, + const std::vector& y_original, + std::vector& x_reduced, + std::vector& y_reduced, + const std::vector& z_original, + std::vector& z_reduced, + const std::vector& A_values, + const std::vector& A_indices, + const std::vector& A_offsets) const +{ + cuopt_expects(presolver_ == cuopt::mathematical_optimization::presolver_t::Papilo, + error_type_t::RuntimeError, + "Crushing is only supported for PaPILO presolve"); + cuopt_assert(papilo_post_solve_storage_ != nullptr, "No postsolve storage available"); + + const auto& storage = *papilo_post_solve_storage_; + const auto& types = storage.types; + const auto& indices = storage.indices; + const auto& values = storage.values; + const auto& start = storage.start; + const auto& num = storage.num; + + cuopt_assert((int)x_original.size() == (int)storage.nColsOriginal, ""); + + const bool crush_dual = !y_original.empty(); + if (crush_dual) { cuopt_assert((int)y_original.size() == (int)storage.nRowsOriginal, ""); } + + const bool crush_rc = !z_original.empty() && crush_dual; + if (crush_rc) { cuopt_assert((int)z_original.size() == (int)storage.nColsOriginal, ""); } + + std::vector x(x_original.begin(), x_original.end()); + std::vector y(y_original.begin(), y_original.end()); + std::vector z(z_original.begin(), z_original.end()); + + // Track current coefficient values for entries modified by kCoefficientChange, + // so repeated changes to the same (row, col) are handled correctly. + std::unordered_map coeff_current; + + const i_t n_cols_original = (i_t)storage.nColsOriginal; + + auto coeff_key = [&](int row, int col) -> i_t { return (i_t)row * n_cols_original + (i_t)col; }; + + auto get_coeff = [&](int row, int col) -> f_t { + auto it = coeff_current.find(coeff_key(row, col)); + if (it != coeff_current.end()) return it->second; + for (i_t p = A_offsets[row]; p < A_offsets[row + 1]; ++p) { + if (A_indices[p] == col) return A_values[p]; + } + return 0; + }; + + for (int i = 0; i < (int)types.size(); ++i) { + int first = start[i]; + + switch (types[i]) { + case ReductionType::kParallelCol: { + // Storage layout: [orig_col1, flags1, orig_col2, flags2, -1] + // [col1lb, col1ub, col2lb, col2ub, col2scale] + int col1 = indices[first]; + int col2 = indices[first + 2]; + const f_t& scale = values[first + 4]; + x[col2] += scale * x[col1]; + if (crush_rc) { + // Inverse of Postsolve::apply_parallel_col_to_original_solution reduced-cost split. + if (num.isZero(z[col2]) && !num.isZero(z[col1])) { + cuopt_assert(!num.isZero(scale), "parallel column scale must be nonzero"); + z[col2] = z[col1] / scale; + } + } + break; + } + + case ReductionType::kRowBoundChangeForcedByRow: { + if (!crush_dual) break; + cuopt_assert(i >= 1 && types[i - 1] == ReductionType::kReasonForRowBoundChangeForcedByRow, + "kRowBoundChangeForcedByRow must be preceded by its reason record"); + + bool is_lhs = indices[first] == 1; + int row = (int)values[first]; + + int reason_first = start[i - 1]; + int deleted_row = indices[reason_first + 1]; + f_t factor = values[reason_first]; + cuopt_assert(factor != 0, "parallel row factor must be nonzero"); + + // Forward rule: if the deleted row carried dual signal that the + // reverse would have attributed to the kept row, transfer it back. + f_t candidate = y[deleted_row] / factor; + bool sign_ok = is_lhs ? num.isGT(candidate, (f_t)0) : num.isLT(candidate, (f_t)0); + + if (sign_ok) { + f_t y_old = y[row]; + y[row] = candidate; + // Maintain z = c - A^T y: propagate the y change into reduced costs + if (crush_rc) { + f_t delta_y = candidate - y_old; + for (i_t p = A_offsets[row]; p < A_offsets[row + 1]; ++p) { + f_t a = get_coeff(row, A_indices[p]); + z[A_indices[p]] -= delta_y * a; + } + } + } + break; + } + + case ReductionType::kCoefficientChange: { + if (!crush_rc) break; + int row = indices[first]; + int col = indices[first + 1]; + f_t a_new = values[first]; + f_t a_old = get_coeff(row, col); + coeff_current[coeff_key(row, col)] = a_new; + z[col] += (a_old - a_new) * y[row]; + break; + } + + case ReductionType::kSubstitutedColWithDual: { + // Singleton substitution: column j is expressed via equality row k as + // x_j = (rhs_k - Σ_{l≠j} a_kl·x_l) / a_kj + // This changes the objective for every column l in row k: + // c_red[l] = c_orig[l] - (c_j / a_kj) · a_kl + // Adjust z accordingly: Δz[l] = -(a_kl / a_kj)·z[j] - a_kl·y[k] + if (!crush_rc) break; + int row_k = indices[first]; // equality row (original space) + int row_length = (int)values[first]; + // Row coefficients start at first+3 + int row_coef_start = first + 3; + // Substituted column index is after the row coefficients + int col_j = indices[row_coef_start + row_length]; + + // Find a_kj (coefficient of col j in row k) + f_t a_kj = 0; + for (int p = 0; p < row_length; ++p) { + if (indices[row_coef_start + p] == col_j) { + a_kj = values[row_coef_start + p]; + break; + } + } + if (a_kj == 0) break; // shouldn't happen + + f_t z_j = z[col_j]; + f_t y_k = y[row_k]; + + // Adjust z for each surviving column l in the equality row (l ≠ j) + for (int p = 0; p < row_length; ++p) { + int col_l = indices[row_coef_start + p]; + if (col_l == col_j) continue; + f_t a_kl = values[row_coef_start + p]; + z[col_l] -= (a_kl / a_kj) * z_j + a_kl * y_k; + } + break; + } + + case ReductionType::kFixedCol: // Handled via projection + case ReductionType::kSubstitutedCol: // Col is dropped + case ReductionType::kFixedInfCol: // Col is dropped + case ReductionType::kVarBoundChange: // Noop + case ReductionType::kRedundantRow: // Noop + case ReductionType::kRowBoundChange: // Noop + case ReductionType::kReasonForRowBoundChangeForcedByRow: // Metadata for above + case ReductionType::kSaveRow: // Metadata + case ReductionType::kReducedBoundsCost: // Noop + case ReductionType::kColumnDualValue: // Column reduced-cost only + case ReductionType::kRowDualValue: // Handled via projection + break; + // no default: case to let the compiler yell at us if a new reduction is later introduced + } + } + + const auto& col_map = storage.origcol_mapping; + const auto& row_map = storage.origrow_mapping; + + // Cancel contributions from removed rows. The original-space z was + // computed as z = c - A^T y over ALL rows. The reduced-space stationarity + // only involves surviving rows, so we must add back the terms from removed + // rows: z[j] += y[i] * a_{i,j} for every removed row i with materially nonzero y[i]. + if (crush_rc) { + std::vector row_survives((int)storage.nRowsOriginal, false); + for (size_t k = 0; k < row_map.size(); ++k) { + row_survives[row_map[k]] = true; + } + for (int i = 0; i < (int)storage.nRowsOriginal; ++i) { + if (row_survives[i] || num.isZero(y[i])) continue; + for (i_t p = A_offsets[i]; p < A_offsets[i + 1]; ++p) { + z[A_indices[p]] += y[i] * get_coeff(i, A_indices[p]); + } + } + } + + x_reduced.resize(col_map.size()); + for (size_t k = 0; k < col_map.size(); ++k) { + x_reduced[k] = x[col_map[k]]; + } + + if (crush_dual) { + y_reduced.resize(row_map.size()); + for (size_t k = 0; k < row_map.size(); ++k) { + y_reduced[k] = y[row_map[k]]; + } + } + + if (crush_rc) { + z_reduced.resize(col_map.size()); + for (size_t k = 0; k < col_map.size(); ++k) { + z_reduced[k] = z[col_map[k]]; + } + } +} + template third_party_presolve_t::~third_party_presolve_t() { diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp index ddf00c83c3..ca8c13cbd1 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.hpp @@ -86,17 +86,21 @@ class third_party_presolve_t { void uncrush_primal_solution(const std::vector& reduced_primal, std::vector& full_primal) const; + + void crush_primal_solution(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, + std::vector& y_reduced, + const std::vector& z_original, + std::vector& z_reduced, + const std::vector& A_values, + const std::vector& A_indices, + const std::vector& A_offsets) const; const std::vector& get_reduced_to_original_map() const { return reduced_to_original_map_; } const std::vector& get_original_to_reduced_map() const { return original_to_reduced_map_; } - const std::vector& get_reduced_to_original_row_map() const - { - return reduced_to_original_row_map_; - } - - const std::vector& get_original_to_reduced_row_map() const - { - return original_to_reduced_row_map_; - } ~third_party_presolve_t(); @@ -123,8 +127,6 @@ class third_party_presolve_t { std::vector reduced_to_original_map_{}; std::vector original_to_reduced_map_{}; - std::vector reduced_to_original_row_map_{}; - std::vector original_to_reduced_row_map_{}; }; } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 80f68ee500..6e6daf7a56 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -453,7 +453,6 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p op_problem, settings.get_tolerances(), settings.determinism_mode == CUOPT_MODE_DETERMINISTIC); auto run_presolve = settings.presolver != presolver_t::None; - run_presolve = run_presolve && settings.initial_solutions.size() == 0; bool has_set_solution_callback = false; for (auto callback : settings.get_mip_callbacks()) { if (callback != nullptr && @@ -473,11 +472,23 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p // Only run if presolve is enabled (gives FJ time to find solutions) // and we're not in deterministic mode + struct early_incumbent_entry_t { + f_t user_obj; + std::vector assignment; + }; + std::vector early_incumbent_pool; + // Track best incumbent found during presolve (shared across CPU and GPU FJ). // early_best_objective is in the original problem's solver-space (always minimization), // used for fast comparison in the callback. // early_best_user_obj is the corresponding user-space objective, // passed to run_mip for correct cross-space conversion. + // We attempt to crush early-heuristics solutions into the presolved space. + // This may fail if some dual reductions cut off parts of the polytope that the + // solutions lie in. We thus may hit scenarios where an excellent incumbent is found + // but is dropped due to these dual reductions, and we lose a good solution. + // This is why we still keep the solution around in original-space + // and later extract it at the end of the solve. std::atomic early_best_objective{std::numeric_limits::infinity()}; f_t early_best_user_obj{std::numeric_limits::infinity()}; std::vector early_best_user_assignment; @@ -495,6 +506,7 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p [&early_best_objective, &early_best_user_obj, &early_best_user_assignment, + &early_incumbent_pool, &early_callback_mutex, early_fj_start, mip_callbacks = settings.get_mip_callbacks(), @@ -513,6 +525,7 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p early_best_objective.store(solver_obj); early_best_user_obj = user_obj; early_best_user_assignment = assignment; + early_incumbent_pool.push_back({user_obj, assignment}); double elapsed = std::chrono::duration(std::chrono::steady_clock::now() - early_fj_start) .count(); @@ -618,6 +631,18 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p early_cpufj.reset(); } + // Add early-heuristic incumbents (original-space) to initial_solutions. + // PaPILO crushing + validation happens downstream in add_user_given_solutions(). + if (!early_incumbent_pool.empty()) { + auto stream = op_problem.get_handle_ptr()->get_stream(); + for (const auto& inc : early_incumbent_pool) { + auto d = std::make_shared>(device_copy(inc.assignment, stream)); + settings.initial_solutions.emplace_back(std::move(d)); + } + CUOPT_LOG_DEBUG("Added %zu early-heuristic incumbents to initial solutions", + early_incumbent_pool.size()); + } + if (settings.user_problem_file != "") { CUOPT_LOG_INFO("Writing user problem to file: %s", settings.user_problem_file.c_str()); op_problem.write_to_mps(settings.user_problem_file); diff --git a/cpp/tests/linear_programming/unit_tests/presolve_test.cu b/cpp/tests/linear_programming/unit_tests/presolve_test.cu index a39c6642b2..f107148d55 100644 --- a/cpp/tests/linear_programming/unit_tests/presolve_test.cu +++ b/cpp/tests/linear_programming/unit_tests/presolve_test.cu @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -23,6 +24,7 @@ #include +#include #include #include #include @@ -100,6 +102,112 @@ static void check_variable_bounds(const std::vector& solution, } } +static void compute_transpose_matvec(const std::vector& values, + const std::vector& indices, + const std::vector& offsets, + const std::vector& y, + int n_cols, + std::vector& result) +{ + assert(!offsets.empty()); + assert(values.size() == indices.size()); + assert(n_cols >= 0); + size_t n_rows = offsets.size() - 1; + assert(y.size() == n_rows); + assert((size_t)offsets.back() == values.size()); + result.assign(n_cols, 0.0); + std::vector kahan_compensation(n_cols, 0.0); + for (size_t i = 0; i < n_rows; ++i) { + for (int j = offsets[i]; j < offsets[i + 1]; ++j) { + int col = indices[j]; + assert(col >= 0 && col < n_cols); + double delta = values[j] * y[i]; + double corrected = delta - kahan_compensation[col]; + double next = result[col] + corrected; + kahan_compensation[col] = (next - result[col]) - corrected; + result[col] = next; + } + } +} + +// Complimentary slackness checks +static void check_reduced_cost_consistency(const std::vector& reduced_cost, + const std::vector& primal, + const std::vector& var_lb, + const std::vector& var_ub, + double bound_tol, + double dual_tol) +{ + constexpr double inf = std::numeric_limits::infinity(); + assert(reduced_cost.size() == primal.size()); + assert(var_lb.size() == primal.size()); + assert(var_ub.size() == primal.size()); + for (size_t j = 0; j < primal.size(); ++j) { + const bool has_lb = var_lb[j] > -inf; + const bool has_ub = var_ub[j] < inf; + const double z_j = reduced_cost[j]; + + // If a side is missing, its multiplier cannot exist. + if (!has_lb) { + ASSERT_LE(z_j, dual_tol) << "positive reduced cost requires a finite lower bound at variable " + << j; + } + if (!has_ub) { + ASSERT_GE(z_j, -dual_tol) + << "negative reduced cost requires a finite upper bound at variable " << j; + } + + // If we are strictly away from a bound, the multiplier for that side must vanish. + if (has_lb && primal[j] - var_lb[j] > bound_tol) { + ASSERT_LE(z_j, dual_tol) + << "positive reduced cost requires an active lower bound at variable " << j; + } + if (has_ub && var_ub[j] - primal[j] > bound_tol) { + ASSERT_GE(z_j, -dual_tol) + << "negative reduced cost requires an active upper bound at variable " << j; + } + } +} + +// General-form row-bound KKT for exported solver duals. +// +// Positive y_i corresponds to the lower row bound, negative y_i to the upper +// row bound. If a side is absent, the corresponding multiplier must be zero. +static void check_dual_sign_consistency(const std::vector& Ax, + const std::vector& dual, + const std::vector& con_lb, + const std::vector& con_ub, + double bound_tol, + double dual_tol) +{ + constexpr double inf = std::numeric_limits::infinity(); + assert(Ax.size() == dual.size()); + assert(con_lb.size() == dual.size()); + assert(con_ub.size() == dual.size()); + for (size_t i = 0; i < dual.size(); ++i) { + const bool has_lb = con_lb[i] > -inf; + const bool has_ub = con_ub[i] < inf; + + if (!has_lb) { + ASSERT_LE(dual[i], dual_tol) + << "positive row dual requires a finite lower bound at constraint " << i; + } + if (!has_ub) { + ASSERT_GE(dual[i], -dual_tol) + << "negative row dual requires a finite upper bound at constraint " << i; + } + + if (has_lb && Ax[i] - con_lb[i] > bound_tol) { + ASSERT_LE(dual[i], dual_tol) + << "positive row dual requires an active lower bound at constraint " << i; + } + if (has_ub && con_ub[i] - Ax[i] > bound_tol) { + ASSERT_GE(dual[i], -dual_tol) + << "negative row dual requires an active upper bound at constraint " << i; + } + } +} + // Test PSLP presolver postsolve accuracy using afiro problem TEST(pslp_presolve, postsolve_accuracy_afiro) { @@ -383,6 +491,414 @@ TEST(pslp_presolve, postsolve_multiple_problems) EXPECT_LT(rel_error, 0.01) << "Problem " << name << " objective mismatch"; } } +struct crush_test_param { + std::string mps_path; + bool use_pdlp; +}; +class dual_crush_round_trip : public ::testing::TestWithParam {}; + +// Crush an optimal original-space (x, y, z) into reduced space and verify the +// user-space general-form KKT conditions on both the original and reduced LPs. +TEST_P(dual_crush_round_trip, kkt_check) +{ + const raft::handle_t handle_{}; + auto stream = handle_.get_stream(); + + constexpr double bound_tol = 1e-5; + constexpr double kkt_tol = 9e-2; + + const auto& param = GetParam(); + auto path = make_path_absolute(param.mps_path); + auto mps = cuopt::mathematical_optimization::io::read_mps(path, false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps); + + // Step 1: Presolve with a single presolver instance (same one used for crush later) + mip::sort_csr(op_problem); + mip::third_party_presolve_t presolver; + auto result = presolver.apply(op_problem, + problem_category_t::LP, + presolver_t::Papilo, + /*dual_postsolve=*/true, + /*abs_tol=*/1e-6, + /*rel_tol=*/1e-9, + /*time_limit=*/60.0); + ASSERT_TRUE(result.status == mip::third_party_presolve_status_t::REDUCED || + result.status == mip::third_party_presolve_status_t::UNCHANGED); + + // Step 2: Solve the reduced problem (no presolve, we already did it) + auto settings = pdlp_solver_settings_t{}; + settings.presolver = presolver_t::None; + settings.dual_postsolve = true; + settings.time_limit = 60.0; + if (param.use_pdlp) { + settings.method = cuopt::mathematical_optimization::method_t::PDLP; + settings.tolerances.absolute_dual_tolerance = 1e-7; + settings.tolerances.relative_dual_tolerance = 0; + settings.tolerances.absolute_primal_tolerance = 1e-7; + settings.tolerances.relative_primal_tolerance = 0; + } else { + settings.method = cuopt::mathematical_optimization::method_t::DualSimplex; + } + + auto reduced_solution = solve_lp(result.reduced_problem, settings); + ASSERT_EQ(reduced_solution.get_termination_status(), pdlp_termination_status_t::Optimal); + + // Step 3: Postsolve to get original-space (x, y, z) using the same presolver. + // For PDLP, derive z = c_red - A_red^T y_red instead of using get_reduced_cost(), + // since PDLP's approximate solution may have inconsistent reduced costs. + auto primal_sol = cuopt::device_copy(reduced_solution.get_primal_solution(), stream); + auto dual_sol = cuopt::device_copy(reduced_solution.get_dual_solution(), stream); + rmm::device_uvector rc_sol(0, stream); + if (param.use_pdlp) { + auto red_A_vals = result.reduced_problem.get_constraint_matrix_values_host(); + auto red_A_indices = result.reduced_problem.get_constraint_matrix_indices_host(); + auto red_A_offsets = result.reduced_problem.get_constraint_matrix_offsets_host(); + auto red_c = result.reduced_problem.get_objective_coefficients_host(); + auto h_y_red = host_copy(dual_sol, stream); + int red_n_vars = result.reduced_problem.get_n_variables(); + + std::vector ATy_red; + compute_transpose_matvec( + red_A_vals, red_A_indices, red_A_offsets, h_y_red, red_n_vars, ATy_red); + std::vector z_red(red_n_vars); + for (int j = 0; j < red_n_vars; ++j) { + z_red[j] = red_c[j] - ATy_red[j]; + } + rc_sol = cuopt::device_copy(z_red, stream); + } else { + rc_sol = cuopt::device_copy(reduced_solution.get_reduced_cost(), stream); + } + + presolver.undo(primal_sol, dual_sol, rc_sol, problem_category_t::LP, false, true, stream); + + auto x_orig = host_copy(primal_sol, stream); + auto y_orig = host_copy(dual_sol, stream); + auto rc_orig = host_copy(rc_sol, stream); + + ASSERT_EQ((int)x_orig.size(), op_problem.get_n_variables()); + ASSERT_EQ((int)y_orig.size(), op_problem.get_n_constraints()); + + // Step 4: Sanity-check the postsolve output in original space before crushing. + auto orig_A_vals = op_problem.get_constraint_matrix_values_host(); + auto orig_A_indices = op_problem.get_constraint_matrix_indices_host(); + auto orig_A_offsets = op_problem.get_constraint_matrix_offsets_host(); + auto orig_c = op_problem.get_objective_coefficients_host(); + auto orig_var_lb = op_problem.get_variable_lower_bounds_host(); + auto orig_var_ub = op_problem.get_variable_upper_bounds_host(); + auto orig_con_lb = op_problem.get_constraint_lower_bounds_host(); + auto orig_con_ub = op_problem.get_constraint_upper_bounds_host(); + { + int orig_n_vars = op_problem.get_n_variables(); + + // Stationarity: z == c - A^T y + std::vector ATy_orig; + compute_transpose_matvec( + orig_A_vals, orig_A_indices, orig_A_offsets, y_orig, orig_n_vars, ATy_orig); + for (size_t j = 0; j < rc_orig.size(); ++j) { + double z_derived = orig_c[j] - ATy_orig[j]; + EXPECT_NEAR(rc_orig[j], z_derived, kkt_tol) + << "postsolve stationarity violation at original variable " << j; + } + + // Primal feasibility + check_variable_bounds(x_orig, orig_var_lb, orig_var_ub, bound_tol); + std::vector Ax_orig; + compute_constraint_residuals(orig_A_vals, orig_A_indices, orig_A_offsets, x_orig, Ax_orig); + check_constraint_satisfaction(Ax_orig, orig_con_lb, orig_con_ub, bound_tol); + + // Variable- and row-bound KKT in original space + check_reduced_cost_consistency(rc_orig, x_orig, orig_var_lb, orig_var_ub, bound_tol, kkt_tol); + check_dual_sign_consistency(Ax_orig, y_orig, orig_con_lb, orig_con_ub, bound_tol, kkt_tol); + } + + // Step 5: Crush using the same presolver that produced the postsolve storage + std::vector x_crushed, y_crushed, rc_crushed; + presolver.crush_primal_dual_solution(x_orig, + y_orig, + x_crushed, + y_crushed, + rc_orig, + rc_crushed, + orig_A_vals, + orig_A_indices, + orig_A_offsets); + + int n_vars = result.reduced_problem.get_n_variables(); + int n_cons = result.reduced_problem.get_n_constraints(); + ASSERT_EQ((int)x_crushed.size(), n_vars); + ASSERT_EQ((int)y_crushed.size(), n_cons); + ASSERT_EQ((int)rc_crushed.size(), n_vars); + + auto A_vals = result.reduced_problem.get_constraint_matrix_values_host(); + auto A_indices = result.reduced_problem.get_constraint_matrix_indices_host(); + auto A_offsets = result.reduced_problem.get_constraint_matrix_offsets_host(); + auto c_red = result.reduced_problem.get_objective_coefficients_host(); + auto var_lb = result.reduced_problem.get_variable_lower_bounds_host(); + auto var_ub = result.reduced_problem.get_variable_upper_bounds_host(); + auto con_lb = result.reduced_problem.get_constraint_lower_bounds_host(); + auto con_ub = result.reduced_problem.get_constraint_upper_bounds_host(); + + // Primal feasibility: variable bounds + check_variable_bounds(x_crushed, var_lb, var_ub, bound_tol); + + // Primal feasibility: constraint satisfaction (l_c <= Ax <= u_c) + std::vector Ax; + compute_constraint_residuals(A_vals, A_indices, A_offsets, x_crushed, Ax); + check_constraint_satisfaction(Ax, con_lb, con_ub, bound_tol); + + // Dual feasibility: compute implied reduced costs z = c - A^T y + std::vector ATy; + compute_transpose_matvec(A_vals, A_indices, A_offsets, y_crushed, n_vars, ATy); + + // Variable-bound KKT in reduced space + check_reduced_cost_consistency(rc_crushed, x_crushed, var_lb, var_ub, bound_tol, kkt_tol); + + // Row-bound KKT in reduced space + check_dual_sign_consistency(Ax, y_crushed, con_lb, con_ub, bound_tol, kkt_tol); + + // Crushed reduced costs: consistency with derived z = c - A^T y + for (size_t j = 0; j < rc_crushed.size(); ++j) { + double z_derived = c_red[j] - ATy[j]; + ASSERT_NEAR(rc_crushed[j], z_derived, kkt_tol) + << "crushed reduced cost vs derived mismatch at variable " << j; + } + + // Cross-check: crushed primal should match the solver's primal + auto x_ref = host_copy(reduced_solution.get_primal_solution(), stream); + ASSERT_EQ(x_crushed.size(), x_ref.size()); + for (size_t i = 0; i < x_crushed.size(); ++i) { + EXPECT_NEAR(x_crushed[i], x_ref[i], kkt_tol) << "primal mismatch at reduced variable " << i; + } + + // Cross-check: crushed objective should match the solver's objective + auto obj_ref = reduced_solution.get_additional_termination_information().primal_objective; + double obj_crushed = 0.0; + for (size_t i = 0; i < x_crushed.size(); ++i) { + obj_crushed += c_red[i] * x_crushed[i]; + } + obj_crushed += result.reduced_problem.get_objective_offset(); + EXPECT_NEAR(obj_crushed, obj_ref, kkt_tol * std::max(1.0, std::abs(obj_ref))) + << "crushed objective mismatch"; +} + +// clang-format off +INSTANTIATE_TEST_SUITE_P( + papilo_presolve, + dual_crush_round_trip, + ::testing::Values( + crush_test_param{"linear_programming/graph40-40/graph40-40.mps", true}, + //crush_test_param{"linear_programming/ex10/ex10.mps", true}, + //crush_test_param{"linear_programming/datt256_lp/datt256_lp.mps", true}, + //crush_test_param{"linear_programming/woodlands09/woodlands09.mps", false}, + //crush_test_param{"linear_programming/savsched1/savsched1.mps", false}, + crush_test_param{"linear_programming/nug08-3rd/nug08-3rd.mps", true}, + crush_test_param{"linear_programming/qap15/qap15.mps", true}, + //crush_test_param{"linear_programming/scpm1/scpm1.mps", true}, + //crush_test_param{"linear_programming/neos3/neos3.mps", true}, + //crush_test_param{"linear_programming/a2864/a2864.mps", true}, + crush_test_param{"linear_programming/afiro_original.mps", false}, + crush_test_param{"mip/enlight_hard.mps", false}, + crush_test_param{"mip/enlight11.mps", false}, + crush_test_param{"mip/50v-10.mps", false}, + crush_test_param{"mip/fiball.mps", false}, + crush_test_param{"mip/gen-ip054.mps", false}, + crush_test_param{"mip/sct2.mps", false}, + crush_test_param{"mip/drayage-25-23.mps", false}, + crush_test_param{"mip/tr12-30.mps", false}, + crush_test_param{"mip/neos-3004026-krka.mps", false}, + crush_test_param{"mip/ns1208400.mps", false}, + crush_test_param{"mip/gmu-35-50.mps", true}, + crush_test_param{"mip/n2seq36q.mps", false}, + crush_test_param{"mip/seymour1.mps", false}, + //crush_test_param{"mip/thor50dday.mps", false}, + crush_test_param{"mip/neos8.mps", false}, + crush_test_param{"mip/app1-1.mps", false}, + crush_test_param{"mip/bnatt400.mps", false}, + crush_test_param{"mip/bnatt500.mps", false}, + //crush_test_param{"mip/brazil3.mps", false}, + crush_test_param{"mip/cbs-cta.mps", false}, + crush_test_param{"mip/CMS750_4.mps", false}, + crush_test_param{"mip/decomp2.mps", false}, + crush_test_param{"mip/dws008-01.mps", false}, + crush_test_param{"mip/germanrr.mps", false}, + crush_test_param{"mip/graph20-20-1rand.mps", false}, + crush_test_param{"mip/milo-v12-6-r2-40-1.mps", false}, + crush_test_param{"mip/neos-1445765.mps", false}, + crush_test_param{"mip/neos-1582420.mps", false}, + crush_test_param{"mip/neos-3083819-nubu.mps", false}, + crush_test_param{"mip/neos-5107597-kakapo.mps", false}, + crush_test_param{"mip/neos-5188808-nattai.mps", false}, + crush_test_param{"mip/net12.mps", false}, + crush_test_param{"mip/rocI-4-11.mps", false}, + crush_test_param{"mip/traininstance2.mps", false}, + crush_test_param{"mip/traininstance6.mps", false}, + crush_test_param{"mip/neos-787933.mps", false}, + crush_test_param{"mip/radiationm18-12-05.mps", false}, + crush_test_param{"mip/momentum1.mps", false}, + crush_test_param{"mip/rococoB10-011000.mps", false}, + crush_test_param{"mip/b1c1s1.mps", false}, + crush_test_param{"mip/nu25-pr12.mps", false}, + crush_test_param{"mip/air05.mps", false}, + crush_test_param{"mip/seymour.mps", false}, + crush_test_param{"mip/swath3.mps", false}, + crush_test_param{"mip/neos-950242.mps", false}, + crush_test_param{"mip/fastxgemm-n2r6s0t2.mps", false} + ), + [](const ::testing::TestParamInfo& info) { + std::string name = info.param.mps_path; + std::replace(name.begin(), name.end(), '/', '_'); + std::replace(name.begin(), name.end(), '.', '_'); + std::replace(name.begin(), name.end(), '-', '_'); + if (info.param.use_pdlp) name += "_pdlp"; + return name; + } +); +// clang-format on + +// Test that crushed solutions work as warmstarts in the full solving pipeline. +// Presolve → cold PDLP solve → postsolve → crush → warmstarted PDLP solve. +// The warmstarted solve should converge in fewer iterations than the cold start. +class crush_warmstart : public ::testing::TestWithParam {}; + +TEST_P(crush_warmstart, round_trip) +{ + const raft::handle_t handle_{}; + auto stream = handle_.get_stream(); + + auto path = make_path_absolute(GetParam()); + auto mps = cuopt::mathematical_optimization::io::read_mps(path, false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps); + + // Step 1: Presolve + mip::sort_csr(op_problem); + mip::third_party_presolve_t presolver; + auto result = presolver.apply(op_problem, + problem_category_t::LP, + presolver_t::Papilo, + /*dual_postsolve=*/true, + /*abs_tol=*/1e-6, + /*rel_tol=*/1e-9, + /*time_limit=*/60.0); + ASSERT_TRUE(result.status == mip::third_party_presolve_status_t::REDUCED || + result.status == mip::third_party_presolve_status_t::UNCHANGED); + + int n_red_vars = result.reduced_problem.get_n_variables(); + int n_red_cons = result.reduced_problem.get_n_constraints(); + + // Step 2: Cold PDLP solve of the reduced problem + auto settings = pdlp_solver_settings_t{}; + settings.presolver = presolver_t::None; + settings.dual_postsolve = true; + settings.method = cuopt::mathematical_optimization::method_t::PDLP; + settings.time_limit = 60.0; + + auto cold_solution = solve_lp(result.reduced_problem, settings); + ASSERT_EQ(cold_solution.get_termination_status(), pdlp_termination_status_t::Optimal); + + auto cold_iters = cold_solution.get_additional_termination_information().number_of_steps_taken; + double cold_obj = cold_solution.get_additional_termination_information().primal_objective; + + // Step 3: Postsolve to original space. + // Recompute z via kahan summation + auto primal_sol = cuopt::device_copy(cold_solution.get_primal_solution(), stream); + auto dual_sol = cuopt::device_copy(cold_solution.get_dual_solution(), stream); + + auto red_A_vals = result.reduced_problem.get_constraint_matrix_values_host(); + auto red_A_indices = result.reduced_problem.get_constraint_matrix_indices_host(); + auto red_A_offsets = result.reduced_problem.get_constraint_matrix_offsets_host(); + auto red_c = result.reduced_problem.get_objective_coefficients_host(); + auto h_y_red = host_copy(dual_sol, stream); + + std::vector ATy_red; + compute_transpose_matvec(red_A_vals, red_A_indices, red_A_offsets, h_y_red, n_red_vars, ATy_red); + std::vector z_red(n_red_vars); + for (int j = 0; j < n_red_vars; ++j) { + z_red[j] = red_c[j] - ATy_red[j]; + } + auto rc_sol = cuopt::device_copy(z_red, stream); + + presolver.undo(primal_sol, dual_sol, rc_sol, problem_category_t::LP, false, true, stream); + + auto x_orig = host_copy(primal_sol, stream); + auto y_orig = host_copy(dual_sol, stream); + + // Step 4: Crush back to reduced space (no rc needed for warmstart) + auto orig_A_vals = op_problem.get_constraint_matrix_values_host(); + auto orig_A_indices = op_problem.get_constraint_matrix_indices_host(); + auto orig_A_offsets = op_problem.get_constraint_matrix_offsets_host(); + + std::vector x_crushed, y_crushed, rc_unused; + presolver.crush_primal_dual_solution(x_orig, + y_orig, + x_crushed, + y_crushed, + {}, + rc_unused, + orig_A_vals, + orig_A_indices, + orig_A_offsets); + + ASSERT_EQ((int)x_crushed.size(), n_red_vars); + ASSERT_EQ((int)y_crushed.size(), n_red_cons); + + // Step 5: Warmstarted PDLP solve of the reduced problem + auto warm_settings = settings; + warm_settings.set_initial_primal_solution(x_crushed.data(), n_red_vars, stream); + warm_settings.set_initial_dual_solution(y_crushed.data(), n_red_cons, stream); + + auto warm_solution = solve_lp(result.reduced_problem, warm_settings); + ASSERT_EQ(warm_solution.get_termination_status(), pdlp_termination_status_t::Optimal); + + double warm_obj = warm_solution.get_additional_termination_information().primal_objective; + auto warm_iters = warm_solution.get_additional_termination_information().number_of_steps_taken; + + double obj_tol = 1e-3 * (1.0 + std::abs(cold_obj)); + EXPECT_NEAR(warm_obj, cold_obj, obj_tol) << "warmstarted objective should match cold solve"; + + EXPECT_LT(warm_iters, cold_iters) + << "warmstarted solve should not take more iterations than cold solve" + << " (cold=" << cold_iters << ", warm=" << warm_iters << ")"; +} + +// clang-format off +INSTANTIATE_TEST_SUITE_P( + papilo_presolve, + crush_warmstart, + ::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 } // namespace cuopt::mathematical_optimization::test diff --git a/cpp/tests/mip/incumbent_callback_test.cu b/cpp/tests/mip/incumbent_callback_test.cu index f9d6b1445b..68131ed054 100644 --- a/cpp/tests/mip/incumbent_callback_test.cu +++ b/cpp/tests/mip/incumbent_callback_test.cu @@ -34,6 +34,22 @@ namespace cuopt::mathematical_optimization::test { +class scoped_env_restore_t { + public: + scoped_env_restore_t(const char* env_name, const char* new_value) : name_(env_name) + { + if (const char* prev = std::getenv(env_name)) { prev_value_ = prev; } + ::setenv(env_name, new_value, 1); + } + ~scoped_env_restore_t() { ::setenv(name_, prev_value_.c_str(), 1); } + scoped_env_restore_t(const scoped_env_restore_t&) = delete; + scoped_env_restore_t& operator=(const scoped_env_restore_t&) = delete; + + private: + const char* name_; + std::string prev_value_; +}; + class test_set_solution_callback_t : public cuopt::internals::set_solution_callback_t { public: test_set_solution_callback_t(std::vector, double>>& solutions_, @@ -161,7 +177,7 @@ TEST(mip_solve, incumbent_get_set_callback_test) // population stays empty. The fallback in solver.cu must use the OG-space incumbent. TEST(mip_solve, early_heuristic_incumbent_fallback) { - setenv("CUOPT_DISABLE_GPU_HEURISTICS", "1", 1); + scoped_env_restore_t disable_gpu_heuristics_env("CUOPT_DISABLE_GPU_HEURISTICS", "1"); const raft::handle_t handle_{}; auto path = make_path_absolute("mip/pk1.mps"); @@ -182,8 +198,6 @@ TEST(mip_solve, early_heuristic_incumbent_fallback) auto solution = solve_mip(op_problem, settings); - unsetenv("CUOPT_DISABLE_GPU_HEURISTICS"); - EXPECT_GE(get_cb.n_calls, 1) << "Early heuristics should have emitted at least one incumbent"; auto status = solution.get_termination_status(); EXPECT_TRUE(status == mip_termination_status_t::FeasibleFound || @@ -195,4 +209,60 @@ TEST(mip_solve, early_heuristic_incumbent_fallback) if (!callback_solutions.empty()) { check_solutions(get_cb, mps_problem, settings); } } +// Verify that a user-provided MIP start in original space is correctly crushed +// through PaPILO presolve and accepted into the heuristic population. +TEST(mip_solve, initial_solution_survives_papilo_crush) +{ + scoped_env_restore_t disable_gpu_heuristics_env("CUOPT_DISABLE_GPU_HEURISTICS", "1"); + + const raft::handle_t handle_{}; + auto path = make_path_absolute("mip/pk1.mps"); + cuopt::mathematical_optimization::io::mps_data_model_t mps_problem = + cuopt::mathematical_optimization::io::read_mps(path, false); + handle_.sync_stream(); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, mps_problem); + auto stream = op_problem.get_handle_ptr()->get_stream(); + + // Step 1: solve to get a reference feasible solution. Pkl is easily solved to optimality + auto settings1 = mip_solver_settings_t{}; + settings1.time_limit = 5.; + settings1.presolver = presolver_t::Papilo; + auto result1 = solve_mip(op_problem, settings1); + auto status1 = result1.get_termination_status(); + ASSERT_TRUE(status1 == mip_termination_status_t::FeasibleFound || + status1 == mip_termination_status_t::Optimal) + << "Reference solve must find a feasible solution"; + auto reference_obj = result1.get_objective_value(); + auto reference_solution = cuopt::host_copy(result1.get_solution(), stream); + ASSERT_EQ((int)reference_solution.size(), op_problem.get_n_variables()); + + // Step 2: feed the reference solution as a MIP start with presolve ON + // and GPU heuristics disabled. B&B runs with node_limit=0 so it exits + // immediately. The only way we get a good objective is if the MIP start + // was crushed through PaPILO and accepted by add_user_given_solutions. + // Early FJ is not strong enough to find the 11 optimal in the given time frame. + auto settings2 = mip_solver_settings_t{}; + settings2.time_limit = 5.; + settings2.presolver = presolver_t::Papilo; + settings2.node_limit = 0; + settings2.add_initial_solution(reference_solution.data(), reference_solution.size(), stream); + + int user_data = 0; + std::vector, double>> callback_solutions; + test_get_solution_callback_t get_cb(callback_solutions, op_problem.get_n_variables(), &user_data); + settings2.set_mip_callback(&get_cb, &user_data); + + auto result2 = solve_mip(op_problem, settings2); + + auto status2 = result2.get_termination_status(); + EXPECT_TRUE(status2 == mip_termination_status_t::FeasibleFound || + status2 == mip_termination_status_t::Optimal) + << "Crushed MIP start should yield a feasible result, got " + << mip_solution_t::get_termination_status_string(status2); + EXPECT_TRUE(std::isfinite(result2.get_objective_value())); + EXPECT_NEAR(result2.get_objective_value(), reference_obj, 1e-4); + + if (!callback_solutions.empty()) { check_solutions(get_cb, mps_problem, settings2); } +} + } // namespace cuopt::mathematical_optimization::test diff --git a/datasets/mip/download_miplib_test_dataset.sh b/datasets/mip/download_miplib_test_dataset.sh index ccfacea1ef..177fb8008f 100755 --- a/datasets/mip/download_miplib_test_dataset.sh +++ b/datasets/mip/download_miplib_test_dataset.sh @@ -27,6 +27,37 @@ INSTANCES=( "enlight11" "supportcase22" "pk1" + "app1-1" + "bnatt400" + "bnatt500" + "brazil3" + "cbs-cta" + "CMS750_4" + "decomp2" + "dws008-01" + "germanrr" + "graph20-20-1rand" + "milo-v12-6-r2-40-1" + "neos-1445765" + "neos-1582420" + "neos-3083819-nubu" + "neos-5107597-kakapo" + "neos-5188808-nattai" + "net12" + "rocI-4-11" + "traininstance2" + "traininstance6" + "neos-787933" + "radiationm18-12-05" + "momentum1" + "rococoB10-011000" + "b1c1s1" + "nu25-pr12" + "air05" + "seymour" + "swath3" + "neos-950242" + "fastxgemm-n2r6s0t2" ) BASE_URL="https://miplib.zib.de/WebData/instances" From 579dc19a5633884e7bef8783588bc80bf8004c3d Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Tue, 30 Jun 2026 18:26:46 +0200 Subject: [PATCH 11/16] warm the submip solve with the crushed pseudocost from the parent. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 38 +++++++++++------ cpp/src/branch_and_bound/branch_and_bound.hpp | 4 ++ cpp/src/branch_and_bound/pseudo_costs.hpp | 42 +++++++++++++++++-- 3 files changed, 68 insertions(+), 16 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e87abba3cd..d0200a39a9 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -340,6 +340,16 @@ void branch_and_bound_t::set_initial_upper_bound(f_t bound) upper_bound_ = bound; } +template +void branch_and_bound_t::warm_start(const pseudo_costs_t& parent_pc, + const std::vector& reduced_to_original, + i_t max_samples) +{ + pc_.resize(original_lp_.num_cols); + pc_.warm_start_from(parent_pc, reduced_to_original, max_samples); + root_warm_start_ = true; +} + template void branch_and_bound_t::print_table_header() { @@ -784,6 +794,9 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& { if (solver_status_ == mip_status_t::SUBMIP_HALT) { settings_.log.printf("Stopping submip solve...\n"); + settings_.log.print_format("lower={}, obj={}\n", + compute_user_objective(original_lp_, lower_bound), + compute_user_objective(original_lp_, upper_bound_.load())); return; } @@ -1884,18 +1897,18 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t 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 (external_upper_bound_ && settings_.inside_submip && - external_upper_bound_->load() <= lower_bound) { - node_concurrent_halt_ = 1; - solver_status_ = mip_status_t::SUBMIP_HALT; - settings_.log.debug( - "Stopping the submip since the main B&B has a better incumbent. " - "(upper=%.g, external_upper=%.g, lower=%.g", - upper_bound_.load(), - external_upper_bound_->load(), - lower_bound); - break; - } + // if (external_upper_bound_ && settings_.inside_submip && + // external_upper_bound_->load() <= lower_bound) { + // node_concurrent_halt_ = 1; + // solver_status_ = mip_status_t::SUBMIP_HALT; + // settings_.log.print_format( + // "Stopping the submip since the main B&B has a better incumbent. " + // "(upper={:g}, external_upper={:g}, lower={:g}", + // upper_bound_.load(), + // external_upper_bound_->load(), + // lower_bound); + // 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. @@ -2222,6 +2235,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_bnb.set_initial_guess(crushed_incumbent); submip_bnb.set_external_upper_bound(external_upper_bound_ ? external_upper_bound_ : &upper_bound_); + submip_bnb.warm_start(pc_, presolver.reduced_to_original_map()); submip_status = submip_bnb.solve(submip_solution); if (submip_status == mip_status_t::INFEASIBLE) { diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 0fb10aa8fd..17f37c0ab2 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -126,6 +126,10 @@ class branch_and_bound_t { // `bound` must be in B&B's internal objective space. void set_initial_upper_bound(f_t bound); + void warm_start(const pseudo_costs_t& parent_pc, + const std::vector& reduced_to_original, + i_t max_samples = 1); + f_t get_upper_bound() const { return upper_bound_.load(); } bool has_solver_space_incumbent() const { return incumbent_.has_incumbent; } diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index eddf2cf28f..28bcea6d80 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,44 @@ class pseudo_costs_t { void resize(i_t num_variables) { - pseudo_cost_sum_down.assign(num_variables, 0); - pseudo_cost_sum_up.assign(num_variables, 0); - pseudo_cost_num_down.assign(num_variables, 0); - pseudo_cost_num_up.assign(num_variables, 0); pseudo_cost_mutex_up.resize(num_variables); pseudo_cost_mutex_down.resize(num_variables); + + if (!warm_start) { + std::fill(pseudo_cost_sum_down.begin(), pseudo_cost_sum_down.end(), 0); + std::fill(pseudo_cost_sum_up.begin(), pseudo_cost_sum_up.end(), 0); + std::fill(pseudo_cost_num_down.begin(), pseudo_cost_num_down.end(), 0); + std::fill(pseudo_cost_num_up.begin(), pseudo_cost_num_up.end(), 0); + } + + pseudo_cost_sum_down.resize(num_variables, 0); + pseudo_cost_sum_up.resize(num_variables, 0); + pseudo_cost_num_down.resize(num_variables, 0); + pseudo_cost_num_up.resize(num_variables, 0); + } + + void warm_start_from(const pseudo_costs_t& parent, const std::vector& reduced_to_original) + { + warm_start = true; + + for (i_t k = 0; k < reduced_to_original.size(); ++k) { + const i_t orig = reduced_to_original[k]; + assert(orig >= 0); + + 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 +262,7 @@ class pseudo_costs_t { std::vector pseudo_cost_mutex_down; omp_atomic_t strong_branching_lp_iter = 0; + bool warm_start = false; }; template From b43caab09f7c4214dc25cfd560423311a64f01ee Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 1 Jul 2026 12:06:42 +0200 Subject: [PATCH 12/16] set an iteration limit for the submip solve Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 40 +++++++++++++------ cpp/src/branch_and_bound/branch_and_bound.hpp | 3 +- cpp/src/branch_and_bound/constants.hpp | 19 ++++----- .../dual_simplex/simplex_solver_settings.hpp | 7 +++- 4 files changed, 43 insertions(+), 26 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index d0200a39a9..910221860c 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -342,11 +342,10 @@ void branch_and_bound_t::set_initial_upper_bound(f_t bound) template void branch_and_bound_t::warm_start(const pseudo_costs_t& parent_pc, - const std::vector& reduced_to_original, - i_t max_samples) + const std::vector& reduced_to_original) { pc_.resize(original_lp_.num_cols); - pc_.warm_start_from(parent_pc, reduced_to_original, max_samples); + pc_.warm_start_from(parent_pc, reduced_to_original); root_warm_start_ = true; } @@ -585,11 +584,11 @@ void branch_and_bound_t::set_solution_from_submip(const std::vector::set_final_solution(mip_solution_t& settings_.log.printf("Node limit reached. Stopping the solver...\n"); } + if (solver_status_ == mip_status_t::ITERATION_LIMIT) { + settings_.log.printf("Simplex iteration limit reached. Stopping the solver...\n"); + } + if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } @@ -1732,6 +1735,12 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, break; } + if (exploration_stats_.total_lp_iters > settings_.bnb_iteration_limit) { + solver_status_ = mip_status_t::ITERATION_LIMIT; + stack.push_front(node_ptr); + break; + } + dual_status_t lp_status = solve_node_lp(node_ptr, worker, exploration_stats_, settings_.log); ++exploration_stats_.nodes_since_last_log; ++exploration_stats_.nodes_explored; @@ -2168,6 +2177,8 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke simplex_solver_settings_t submip_settings; submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; + submip_settings.bnb_iteration_limit = + exploration_stats_.total_lp_iters * submip_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; } @@ -2237,6 +2248,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke : &upper_bound_); submip_bnb.warm_start(pc_, presolver.reduced_to_original_map()); submip_status = submip_bnb.solve(submip_solution); + exploration_stats_.total_lp_iters += submip_solution.simplex_iterations; if (submip_status == mip_status_t::INFEASIBLE) { submip_stats_.save_infeasible(fixrate); @@ -2932,7 +2944,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (!enable_concurrent_lp_root_solve()) { // RINS/SUBMIP path - settings_.log.printf("\nSolving LP root relaxation with dual simplex\n"); + 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, @@ -2945,7 +2958,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_solved_by = DualSimplex; } 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_, @@ -3007,11 +3021,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); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 17f37c0ab2..015742ef11 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -127,8 +127,7 @@ class branch_and_bound_t { void set_initial_upper_bound(f_t bound); void warm_start(const pseudo_costs_t& parent_pc, - const std::vector& reduced_to_original, - i_t max_samples = 1); + 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; } diff --git a/cpp/src/branch_and_bound/constants.hpp b/cpp/src/branch_and_bound/constants.hpp index f3e21126dd..0f8118ce2e 100644 --- a/cpp/src/branch_and_bound/constants.hpp +++ b/cpp/src/branch_and_bound/constants.hpp @@ -44,15 +44,16 @@ enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 }; enum class branch_and_bound_mode_t { PARALLEL = 0, DETERMINISTIC = 1 }; enum class mip_status_t { - OPTIMAL = 0, // The optimal integer solution was found - UNBOUNDED = 1, // The problem is unbounded - INFEASIBLE = 2, // The problem is infeasible - TIME_LIMIT = 3, // The solver reached a time limit - NODE_LIMIT = 4, // The maximum number of nodes was reached (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 - SUBMIP_HALT = 8 // Stop submip solve since it no longer valid + 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 // Stop submip solve since it no longer valid }; } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 224242ce04..ed480a1d75 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -55,8 +55,9 @@ struct submip_settings_t { // Presolve sub-MIP with Papilo before solving it bool presolve = true; - // Warm start sub-MIP's B&B with the parent basis and pseudocost - bool warm_start = true; + // Limit the number of simplex iterations spent in the submip. Set as a factor of the total + // number of simplex iteration from the parent B&B. + double iteration_limit_ratio = 0.8; }; template @@ -67,6 +68,7 @@ struct simplex_solver_settings_t { node_limit(std::numeric_limits::max()), time_limit(std::numeric_limits::infinity()), work_limit(std::numeric_limits::infinity()), + bnb_iteration_limit(std::numeric_limits::max()), absolute_mip_gap_tol(0.0), relative_mip_gap_tol(1e-3), integer_tol(1e-5), @@ -152,6 +154,7 @@ struct simplex_solver_settings_t { i_t node_limit; f_t time_limit; f_t work_limit; + i_t bnb_iteration_limit; // Limit of the total number of simplex iterations in B&B f_t absolute_mip_gap_tol; // Tolerance on mip gap to declare optimal f_t relative_mip_gap_tol; // Tolerance on mip gap to declare optimal f_t integer_tol; // Tolerance on integralitiy violation From 234361d246bf11a452a73f43682654fdef1cce6d Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 1 Jul 2026 13:43:35 +0200 Subject: [PATCH 13/16] extend the neighbourhood to reach the min fixrate using a coefficient/root change diving procedure when the RINS neighbourhood is not large enough Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 75 ++++++++++++++++--- 1 file changed, 64 insertions(+), 11 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 910221860c..89582350f4 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1708,10 +1708,10 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, if (worker->worker_id == 0) { f_t time_since_last_log = - exploration_stats_.last_log == 0 ? 1.0 : toc(exploration_stats_.last_log); + exploration_stats_.last_log == 0 ? 50.0 : toc(exploration_stats_.last_log); i_t nodes_since_last_log = exploration_stats_.nodes_since_last_log; - if (((nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + if (((nodes_since_last_log >= 100 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && time_since_last_log >= 1) || (time_since_last_log > 30) || now > settings_.time_limit) { report(' ', upper_bound_, lower_bound, node_ptr->depth, node_ptr->integer_infeasible); @@ -2194,7 +2194,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.submip_settings.level = submip_level; - submip_settings.log.log = true; + submip_settings.log.log = false; submip_settings.log.log_prefix = log_prefix; submip_settings.submip_settings.enable_rins = settings_.submip_settings.enable_rins != 0 && @@ -2376,15 +2376,68 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, } // Even considering the entire integer list, we were unable to fix a single variable in this - // iteration. Try to solve the submip with what we got. + // 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) { - settings_.log.print_format("{}Fixed {} variables, the maximum for RINS. (max={}, min={})\n", - log_prefix, - num_var_fixed, - max_var_fixed, - min_var_fixed); - has_submip = true; - break; + std::vector> candidates; + for (i_t j : fractional) { + if (std::abs(lower[j] - upper[j]) <= settings_.fixed_tol) { continue; } + + f_t root_change = current_sol[j] - root_relax_soln_.x[j]; + f_t obj_coeff = rins_worker->leaf_problem.objective[j]; + f_t fixed_val = 0; + + if (root_change >= 0.4) { + fixed_val = std::ceil(current_sol[j]); + } else if (root_change <= -0.4) { + fixed_val = std::floor(current_sol[j]); + } else if (obj_coeff > settings_.zero_tol) { + fixed_val = std::ceil(current_sol[j]); + } else if (obj_coeff < -settings_.zero_tol) { + fixed_val = std::floor(current_sol[j]); + } else { + fixed_val = std::round(current_sol[j]); + } + + candidates.push_back(std::make_tuple(std::abs(fixed_val - current_sol[j]), j, fixed_val)); + } + + std::sort(candidates.begin(), candidates.end(), [](auto a, auto b) { + return std::get<0>(a) < std::get<0>(b); + }); + + f_t change = 0; + for (auto [dist, j, fixed_val] : candidates) { + fixed_val = std::clamp(fixed_val, lower[j], upper[j]); + lower[j] = fixed_val; + upper[j] = fixed_val; + bounds_changed[j] = true; + ++num_var_fixed; + if (num_var_fixed >= max_var_fixed) break; + + // Limit the amount of fixing to the current LP. + change += dist; + if (change >= 0.5) { break; } + } + + if (num_var_fixed >= min_var_fixed) { + settings_.log.print_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.print_format("{} Could not fix more variables\n", + log_prefix, + num_var_fixed, + max_var_fixed, + min_var_fixed); + break; + } } } From f7cd0673edd99fa1ada266d986cdfa0981681b46 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 1 Jul 2026 16:14:47 +0200 Subject: [PATCH 14/16] properly handle the upper bound changes in the main solve. fixed iteration count for the root relaxation (it should reflect the number of simplex iterations). Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 145 +++++++++++------- cpp/src/branch_and_bound/branch_and_bound.hpp | 11 +- .../deterministic_workers.hpp | 4 +- cpp/src/branch_and_bound/pseudo_costs.cpp | 2 +- cpp/src/branch_and_bound/worker.hpp | 6 +- .../dual_simplex/simplex_solver_settings.hpp | 20 +-- 6 files changed, 110 insertions(+), 78 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 89582350f4..8399edb042 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -420,7 +420,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); @@ -584,8 +584,8 @@ void branch_and_bound_t::set_solution_from_submip(const std::vector::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 || @@ -890,7 +890,7 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& } 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 @@ -1560,10 +1560,10 @@ dual_status_t branch_and_bound_t::solve_node_lp( lp_settings.scale_columns = false; if (worker->search_strategy != search_strategy_t::BEST_FIRST) { - int64_t bnb_lp_iters = exploration_stats_.total_lp_iters; + int64_t bnb_lp_iters = exploration_stats_.total_simplex_iters; f_t factor = settings_.diving_settings.iteration_limit_factor; int64_t max_iter = factor * bnb_lp_iters; - lp_settings.iteration_limit = max_iter - stats.total_lp_iters; + lp_settings.iteration_limit = max_iter - stats.total_simplex_iters; if (lp_settings.iteration_limit <= 0) { return dual_status_t::ITERATION_LIMIT; } } @@ -1634,7 +1634,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; } } @@ -1735,7 +1735,7 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, break; } - if (exploration_stats_.total_lp_iters > settings_.bnb_iteration_limit) { + if (exploration_stats_.total_simplex_iters > settings_.bnb_iteration_limit) { solver_status_ = mip_status_t::ITERATION_LIMIT; stack.push_front(node_ptr); break; @@ -1770,7 +1770,7 @@ 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 (worker->worker_id == 0) { launch_submip_worker(worker->leaf_solution.x); } + launch_submip_worker(worker->leaf_solution.x); // The stack should only contain the children of the current parent. // If the stack size is greater than 0, @@ -1906,18 +1906,22 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t 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 (external_upper_bound_ && settings_.inside_submip && - // external_upper_bound_->load() <= lower_bound) { - // node_concurrent_halt_ = 1; - // solver_status_ = mip_status_t::SUBMIP_HALT; - // settings_.log.print_format( - // "Stopping the submip since the main B&B has a better incumbent. " - // "(upper={:g}, external_upper={:g}, lower={:g}", - // upper_bound_.load(), - // external_upper_bound_->load(), - // lower_bound); - // break; - // } + if (external_upper_bound_callback_) { + const f_t external_upper_bound = external_upper_bound_callback_(); + bool maximize = original_lp_.obj_scale < 0.0; + bool stop = maximize ? user_lower < external_upper_bound : user_lower > external_upper_bound; + if (stop) { + node_concurrent_halt_ = 1; + solver_status_ = mip_status_t::SUBMIP_HALT; + settings_.log.print_format( + "Stopping the submip since the main B&B has a better incumbent. " + "(upper={:g}, external_upper={:g}, lower={:g}", + user_obj, + external_upper_bound, + 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. @@ -2141,16 +2145,18 @@ void branch_and_bound_t::launch_submip_worker(const std::vector& if (submip_worker_pool_.num_idle() == 0) return; submip_worker_t* submip_worker = submip_worker_pool_.pop_idle_worker(); + if (!submip_worker) return; submip_worker->set_active(); #pragma omp task priority(CUOPT_MEDIUM_TASK_PRIORITY) affinity(submip_worker) \ - firstprivate(submip_worker, sol) if (!settings_.inside_submip) + firstprivate(submip_worker, sol) rins(submip_worker, sol); } template void branch_and_bound_t::solve_submip(submip_worker_t* worker, + const std::vector& current_incumbent, i_t num_var_fixed, i_t num_integers, i_t submip_level, @@ -2176,15 +2182,6 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke i_t explored = exploration_stats_.nodes_explored; simplex_solver_settings_t submip_settings; - submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; - submip_settings.bnb_iteration_limit = - exploration_stats_.total_lp_iters * submip_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.print_presolve_stats = false; submip_settings.num_threads = 1; submip_settings.reliability_branching = 0; @@ -2194,12 +2191,32 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.submip_settings.level = submip_level; - submip_settings.log.log = false; + submip_settings.log.log = true; submip_settings.log.log_prefix = log_prefix; + submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; + + submip_settings.bnb_iteration_limit = + exploration_stats_.total_simplex_iters * submip_settings.submip_settings.iteration_limit_ratio; + + submip_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); + if (submip_settings.time_limit < 0) { return; } + + submip_settings.relative_mip_gap_tol = + std::min(settings_.submip_settings.target_mip_gap, rel_gap); + submip_settings.submip_settings.enable_rins = settings_.submip_settings.enable_rins != 0 && submip_level <= settings_.submip_settings.max_level; + submip_settings.log.print_format( + "Sub-MIP solve settings: time_limit={:.2f}, node_limit={}, iter_limit={} (current_iter={}), " + "tol={:g}", + submip_settings.time_limit, + submip_settings.node_limit, + submip_settings.bnb_iteration_limit, + exploration_stats_.total_simplex_iters.load(), + submip_settings.relative_mip_gap_tol); + user_problem_t submip_problem(original_problem_.handle_ptr); simplex::convert_simplex_problem( worker->leaf_problem, var_types_, settings_, new_slacks_, submip_problem); @@ -2236,19 +2253,21 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke branch_and_bound_t submip_bnb(submip_problem, submip_settings, tic(), empty_probing); mip_solution_t submip_solution(submip_problem.num_cols); - mutex_upper_.lock(); - std::vector current_incumbent = incumbent_.x; - mutex_upper_.unlock(); - std::vector crushed_incumbent; presolver.crush(current_incumbent, crushed_incumbent); submip_bnb.set_initial_guess(crushed_incumbent); - submip_bnb.set_external_upper_bound(external_upper_bound_ ? external_upper_bound_ - : &upper_bound_); + + if (external_upper_bound_callback_) { + submip_bnb.set_external_upper_bound_callback(external_upper_bound_callback_); + } else { + submip_bnb.set_external_upper_bound_callback( + [this]() { return compute_user_objective(this->original_lp_, this->upper_bound_.load()); }); + } + submip_bnb.warm_start(pc_, presolver.reduced_to_original_map()); submip_status = submip_bnb.solve(submip_solution); - exploration_stats_.total_lp_iters += submip_solution.simplex_iterations; + exploration_stats_.total_simplex_iters += submip_solution.simplex_iterations; if (submip_status == mip_status_t::INFEASIBLE) { submip_stats_.save_infeasible(fixrate); @@ -2334,7 +2353,7 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // Enough variables has been fixed if (num_var_fixed >= min_var_fixed) { - settings_.log.print_format("{}Fixed {} variables (max={}, min={})\n", + settings_.log.debug_format("{}Fixed {} variables (max={}, min={})\n", log_prefix, num_var_fixed, max_var_fixed, @@ -2366,7 +2385,7 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // Enough variables were fixed if (num_var_fixed >= min_var_fixed) { - settings_.log.print_format("{}Fixed {} variables (max={}, min={})\n", + settings_.log.debug_format("{}Fixed {} variables (max={}, min={})\n", log_prefix, num_var_fixed, max_var_fixed, @@ -2421,7 +2440,7 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, } if (num_var_fixed >= min_var_fixed) { - settings_.log.print_format("{}Fixed {} variables (max={}, min={})\n", + settings_.log.debug_format("{}Fixed {} variables (max={}, min={})\n", log_prefix, num_var_fixed, max_var_fixed, @@ -2431,7 +2450,7 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, } if (prev_num_fixed == num_var_fixed) { - settings_.log.print_format("{} Could not fix more variables\n", + settings_.log.debug_format("{} Could not fix more variables\n", log_prefix, num_var_fixed, max_var_fixed, @@ -2475,7 +2494,8 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, // there is no enough variable fixing in the submip to be easier to solve than the current // problem if (fixrate >= settings_.submip_settings.min_fixrate_cap) { - solve_submip(rins_worker, num_var_fixed, num_integers, submip_level, log_prefix); + solve_submip( + rins_worker, current_incumbent, num_var_fixed, num_integers, submip_level, log_prefix); } } @@ -2568,6 +2588,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_; @@ -2616,12 +2642,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; @@ -2810,7 +2838,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); @@ -2836,7 +2864,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()); @@ -2950,7 +2978,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut upper_bound_ = computed_obj; mutex_upper_.unlock(); - settings_.log.debug_format("Warm start with incumbent, obj={:.4g}", computed_obj); + settings_.log.print_format("Warm starting B&B with an initial guess. Obj={:.4g}", + compute_user_objective(original_lp_, computed_obj)); } } @@ -2999,7 +3028,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // RINS/SUBMIP path 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_, + root_status = solve_linear_program_with_advanced_basis(original_lp_, exploration_stats_.start_time, lp_settings, root_relax_soln_, @@ -3008,7 +3037,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nonbasic_list, root_vstatus_, edge_norms_); - root_relax_solved_by = DualSimplex; + root_relax_solved_by = DualSimplex; + exploration_stats_.total_simplex_iters = root_relax_soln_.iterations; } else { settings_.log.printf("\n"); @@ -3023,7 +3053,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } 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; @@ -4057,7 +4086,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; @@ -4144,10 +4173,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) { diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 015742ef11..5bfe7c39f0 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -79,9 +79,9 @@ 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_external_upper_bound(const omp_atomic_t* external_upper_bound) + void set_external_upper_bound_callback(std::function callback) { - external_upper_bound_ = external_upper_bound; + external_upper_bound_callback_ = std::move(callback); } // Set the root solution found by PDLP @@ -212,8 +212,10 @@ 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_; - // For submip - const omp_atomic_t* external_upper_bound_{nullptr}; + // Callback that returns the objective of the current incumbent of the main solve in the + // user space. This is used for propagating changes to the upper bound to the + // submip solves. + std::function external_upper_bound_callback_; // Solver-space incumbent tracked directly by B&B. simplex::mip_solution_t incumbent_; @@ -355,6 +357,7 @@ class branch_and_bound_t { void launch_submip_worker(const std::vector& sol); void solve_submip(submip_worker_t* worker, + const std::vector& current_incumbent, i_t num_var_fixed, i_t num_integers, i_t submip_level, 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/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index 22a246456a..d3741d8836 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -1489,7 +1489,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, diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index da27945ac4..ac953a2553 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; diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index ed480a1d75..3bdd5a2e1f 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -154,16 +154,16 @@ struct simplex_solver_settings_t { i_t node_limit; f_t time_limit; f_t work_limit; - i_t bnb_iteration_limit; // Limit of the total number of simplex iterations in B&B - f_t absolute_mip_gap_tol; // Tolerance on mip gap to declare optimal - f_t relative_mip_gap_tol; // Tolerance on mip gap to declare optimal - f_t integer_tol; // Tolerance on integralitiy violation - f_t primal_tol; // Absolute primal infeasibility tolerance - f_t dual_tol; // Absolute dual infeasibility tolerance - f_t pivot_tol; // Simplex pivot tolerance - f_t tight_tol; // A tight tolerance used to check for infeasibility - f_t fixed_tol; // If l <= x <= u with u - l < fixed_tol a variable is consider fixed - f_t zero_tol; // Values below this tolerance are considered numerically zero + int64_t bnb_iteration_limit; // Limit of the total number of simplex iterations in B&B + f_t absolute_mip_gap_tol; // Tolerance on mip gap to declare optimal + f_t relative_mip_gap_tol; // Tolerance on mip gap to declare optimal + f_t integer_tol; // Tolerance on integralitiy violation + f_t primal_tol; // Absolute primal infeasibility tolerance + f_t dual_tol; // Absolute dual infeasibility tolerance + f_t pivot_tol; // Simplex pivot tolerance + f_t tight_tol; // A tight tolerance used to check for infeasibility + f_t fixed_tol; // If l <= x <= u with u - l < fixed_tol a variable is consider fixed + f_t zero_tol; // Values below this tolerance are considered numerically zero f_t barrier_relative_feasibility_tol; // Relative feasibility tolerance for barrier method f_t barrier_relative_optimality_tol; // Relative optimality tolerance for barrier method f_t From b30df804560d2750822ea5f0c78ff62b61a71bee Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 1 Jul 2026 17:31:29 +0200 Subject: [PATCH 15/16] fix incorrect submip iteration count. moved submipt worker to worker.hpp Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 64 +++++++------ cpp/src/branch_and_bound/branch_and_bound.hpp | 3 +- cpp/src/branch_and_bound/constants.hpp | 2 +- cpp/src/branch_and_bound/submip.hpp | 94 ------------------- cpp/src/branch_and_bound/worker.hpp | 73 ++++++++++++++ cpp/src/branch_and_bound/worker_pool.hpp | 3 + .../dual_simplex/simplex_solver_settings.hpp | 2 +- 7 files changed, 112 insertions(+), 129 deletions(-) delete mode 100644 cpp/src/branch_and_bound/submip.hpp diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 8399edb042..fdff4ffecb 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -796,12 +796,10 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& settings_.log.print_format("lower={}, obj={}\n", compute_user_objective(original_lp_, lower_bound), compute_user_objective(original_lp_, upper_bound_.load())); - return; } if (solver_status_ == mip_status_t::NUMERICAL) { settings_.log.printf("Numerical issue encountered. Stopping the solver...\n"); - return; } if (solver_status_ == mip_status_t::TIME_LIMIT) { @@ -1668,6 +1666,8 @@ 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(); } @@ -1770,7 +1770,7 @@ 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) { - launch_submip_worker(worker->leaf_solution.x); + if (can_launch_rins) { can_launch_rins = !launch_submip_worker(worker->leaf_solution.x); } // The stack should only contain the children of the current parent. // If the stack size is greater than 0, @@ -1915,7 +1915,7 @@ void branch_and_bound_t::best_first_search_with(bfs_worker_t solver_status_ = mip_status_t::SUBMIP_HALT; settings_.log.print_format( "Stopping the submip since the main B&B has a better incumbent. " - "(upper={:g}, external_upper={:g}, lower={:g}", + "upper={:g}, external_upper={:g}, lower={:g}\n", user_obj, external_upper_bound, user_lower); @@ -2138,20 +2138,22 @@ bool branch_and_bound_t::launch_diving_worker(bfs_worker_t* } template -void branch_and_bound_t::launch_submip_worker(const std::vector& sol) +bool branch_and_bound_t::launch_submip_worker(const std::vector& sol) { - if (settings_.submip_settings.enable_rins == 0) return; - if (!incumbent_.has_incumbent) return; - if (submip_worker_pool_.num_idle() == 0) return; + if (settings_.submip_settings.enable_rins == 0) return false; + if (!incumbent_.has_incumbent) return false; + if (submip_worker_pool_.num_idle() == 0) return false; submip_worker_t* submip_worker = submip_worker_pool_.pop_idle_worker(); - if (!submip_worker) return; + if (!submip_worker) return false; submip_worker->set_active(); #pragma omp task priority(CUOPT_MEDIUM_TASK_PRIORITY) affinity(submip_worker) \ - firstprivate(submip_worker, sol) + firstprivate(submip_worker, sol) if (!settings_.inside_submip) rins(submip_worker, sol); + + return true; } template @@ -2191,7 +2193,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.submip_settings.level = submip_level; - submip_settings.log.log = true; + submip_settings.log.log = false; submip_settings.log.log_prefix = log_prefix; submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; @@ -2395,8 +2397,8 @@ void branch_and_bound_t::rins(submip_worker_t* rins_worker, } // 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. + // iteration. Iterate over the fractional variables again and fixing those that closest to + // an integer solution first in order to reach the fixing threshold. if (prev_num_fixed == num_var_fixed) { std::vector> candidates; for (i_t j : fractional) { @@ -3573,8 +3575,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) */ @@ -3612,23 +3614,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. */ @@ -3738,7 +3739,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; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index 5bfe7c39f0..d413736e8f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -12,7 +12,6 @@ #include #include #include -#include #include #include @@ -355,7 +354,7 @@ class branch_and_bound_t { // to find integer feasible solutions. void dive_with(diving_worker_t* worker); - void launch_submip_worker(const std::vector& sol); + bool launch_submip_worker(const std::vector& sol); void solve_submip(submip_worker_t* worker, const std::vector& current_incumbent, i_t num_var_fixed, diff --git a/cpp/src/branch_and_bound/constants.hpp b/cpp/src/branch_and_bound/constants.hpp index 0f8118ce2e..29d6711abd 100644 --- a/cpp/src/branch_and_bound/constants.hpp +++ b/cpp/src/branch_and_bound/constants.hpp @@ -53,7 +53,7 @@ enum class mip_status_t { 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 // Stop submip solve since it no longer valid + SUBMIP_HALT = 9 // Halt the search on a suboptimal sub-MIP }; } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/submip.hpp b/cpp/src/branch_and_bound/submip.hpp deleted file mode 100644 index 72ac7b8537..0000000000 --- a/cpp/src/branch_and_bound/submip.hpp +++ /dev/null @@ -1,94 +0,0 @@ -/* clang-format off */ -/* - * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. - * SPDX-License-Identifier: Apache-2.0 - */ -/* clang-format on */ - -#pragma once - -#include -#include - -namespace cuopt::mathematical_optimization::mip { - -template -class branch_and_bound_t; - -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; } -}; - -inline double submip_get_max_fixrate(const submip_stats_t& stats, - const simplex::submip_settings_t& submip_settings, - pcgenerator_t& rng) -{ - // Adaptive fix rate based on previous successes and failures. - double low = submip_settings.base_target_fixrate; - double high = submip_settings.base_target_fixrate; - - if (stats.total_infeasible > 0) { - double infeasible_avg_fixrate = stats.average_infeasible_fixrate(); - high = 0.9 * infeasible_avg_fixrate; - low = std::min(low, high); - } - - if (stats.total_success > 0) { - double success_avg_fixrate = stats.average_success_fixrate(); - low = std::min(low, 0.9 * success_avg_fixrate); - high = std::max(high, 1.1 * success_avg_fixrate); - } - - double fixrate = high > low ? rng.uniform(low, high) : low; - return fixrate; -} - -template -class submip_worker_t : public branch_and_bound_worker_t { - public: - using Base = branch_and_bound_worker_t; - - submip_worker_t(i_t worker_id, - const simplex::lp_problem_t& original_lp, - const csr_matrix_t& Arow, - const std::vector& var_type, - const simplex::simplex_solver_settings_t& settings, - uint64_t rng_offset = 0) - : Base(worker_id, original_lp, Arow, var_type, settings, rng_offset) - { - this->search_strategy = SUBMIP; - } - - // Set this node inactive - void set_inactive() - { - if (!this->is_active.load()) { return; } - this->is_active = false; - } - - f_t get_lower_bound() { return this->lower_bound; } -}; - -template -using submip_worker_pool_t = worker_pool_t>; - -} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp index ac953a2553..4ad51683b5 100644 --- a/cpp/src/branch_and_bound/worker.hpp +++ b/cpp/src/branch_and_bound/worker.hpp @@ -259,4 +259,77 @@ 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; } +}; + +inline double submip_get_max_fixrate(const submip_stats_t& stats, + const simplex::submip_settings_t& submip_settings, + pcgenerator_t& rng) +{ + // Adaptive fix rate based on previous successes and failures. + double low = submip_settings.base_target_fixrate; + double high = submip_settings.base_target_fixrate; + + if (stats.total_infeasible > 0) { + double infeasible_avg_fixrate = stats.average_infeasible_fixrate(); + high = 0.9 * infeasible_avg_fixrate; + low = std::min(low, high); + } + + if (stats.total_success > 0) { + double success_avg_fixrate = stats.average_success_fixrate(); + low = std::min(low, 0.9 * success_avg_fixrate); + high = std::max(high, 1.1 * success_avg_fixrate); + } + + double fixrate = high > low ? rng.uniform(low, high) : low; + return fixrate; +} + +template +class submip_worker_t : public branch_and_bound_worker_t { + public: + using Base = branch_and_bound_worker_t; + + submip_worker_t(i_t worker_id, + const simplex::lp_problem_t& original_lp, + const csr_matrix_t& Arow, + const std::vector& var_type, + const simplex::simplex_solver_settings_t& settings, + uint64_t rng_offset = 0) + : Base(worker_id, original_lp, Arow, var_type, settings, rng_offset) + { + this->search_strategy = SUBMIP; + } + + // Set this node inactive + void set_inactive() + { + if (!this->is_active.load()) { return; } + this->is_active = false; + } + + f_t get_lower_bound() { return this->lower_bound; } +}; + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp index e9d55bafe3..5413c6faea 100644 --- a/cpp/src/branch_and_bound/worker_pool.hpp +++ b/cpp/src/branch_and_bound/worker_pool.hpp @@ -106,4 +106,7 @@ using bfs_worker_pool_t = worker_pool_t>; template using diving_worker_pool_t = worker_pool_t>; +template +using submip_worker_pool_t = worker_pool_t>; + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 3bdd5a2e1f..f9b05a619d 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -50,7 +50,7 @@ struct submip_settings_t { int level = 0; // Maximum recursion level - int max_level = 5; + int max_level = 10; // Presolve sub-MIP with Papilo before solving it bool presolve = true; From ce211f379fbddc7148f9a21342e08fc8d9b22736 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 1 Jul 2026 18:07:06 +0200 Subject: [PATCH 16/16] fix incorrect subspace of the incumbent Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 21 +++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index fdff4ffecb..efcc0412b7 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2193,7 +2193,7 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_settings.inside_submip = 1; submip_settings.strong_branching_simplex_iteration_limit = 50; submip_settings.submip_settings.level = submip_level; - submip_settings.log.log = false; + submip_settings.log.log = true; submip_settings.log.log_prefix = log_prefix; submip_settings.node_limit = settings_.submip_settings.node_limit_base + explored / 20; @@ -2255,9 +2255,14 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke branch_and_bound_t submip_bnb(submip_problem, submip_settings, tic(), empty_probing); mip_solution_t submip_solution(submip_problem.num_cols); - std::vector crushed_incumbent; - presolver.crush(current_incumbent, crushed_incumbent); + // Map the current incumbent to the user space, removing the variables related to the slacks/cuts + // added. + std::vector uncrushed_incumbent; + uncrush_primal_solution(original_problem_, original_lp_, current_incumbent, uncrushed_incumbent); + // Crush the incumbent to presolve space. + std::vector crushed_incumbent; + presolver.crush(uncrushed_incumbent, crushed_incumbent); submip_bnb.set_initial_guess(crushed_incumbent); if (external_upper_bound_callback_) { @@ -2271,6 +2276,9 @@ void branch_and_bound_t::solve_submip(submip_worker_t* worke submip_status = submip_bnb.solve(submip_solution); exploration_stats_.total_simplex_iters += submip_solution.simplex_iterations; + submip_settings.log.print_format( + "Sub-MIP: status={}, iterations={} \n", (int)submip_status, submip_solution.simplex_iterations); + if (submip_status == mip_status_t::INFEASIBLE) { submip_stats_.save_infeasible(fixrate); return; @@ -3457,9 +3465,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (settings_.deterministic) { run_deterministic_coordinator(Arow_); } 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_workers = settings_.num_threads; + const i_t num_bfs_workers = std::max(settings_.num_threads / 2, 1); + const i_t num_diving_workers = + settings_.inside_mip ? num_search_strategies - 2 : num_workers - num_bfs_workers; bfs_worker_pool_.init(num_bfs_workers, original_lp_, Arow_, var_types_, symmetry_, settings_); submip_worker_pool_.init(1, original_lp_, Arow_, var_types_, symmetry_, settings_);