From 033f1ce1cf1fd382f73022707b82d110be2b707e Mon Sep 17 00:00:00 2001 From: anandhkb Date: Sun, 8 Mar 2026 17:54:35 -0700 Subject: [PATCH 1/9] Add MIP gap to heuristic log during dual-simplex root relaxation - Add root_lp_progress_callback to simplex_solver_settings (optional) - Invoke callback in dual_phase2 with current dual objective when inside_mip=1 - B&B: store root LP lower bound via callback, use in report_heuristic - Append 'Gap: X%' to log line when heuristic finds solution during root solve Addresses issue when solver times out during root relaxation: gap can now be computed from last dual feasible solution and heuristic incumbent. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 25 +++++++++++++++---- cpp/src/branch_and_bound/branch_and_bound.hpp | 1 + cpp/src/dual_simplex/phase2.cpp | 3 +++ .../dual_simplex/simplex_solver_settings.hpp | 3 +++ 4 files changed, 27 insertions(+), 5 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 6ce9a4f4d0..f09158e7f2 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -282,8 +282,9 @@ branch_and_bound_t::branch_and_bound_t( } #endif - upper_bound_ = inf; - root_objective_ = std::numeric_limits::quiet_NaN(); + upper_bound_ = inf; + root_lp_current_lower_bound_ = -inf; + root_objective_ = std::numeric_limits::quiet_NaN(); } template @@ -318,9 +319,19 @@ void branch_and_bound_t::report_heuristic(f_t obj) user_gap.c_str(), toc(exploration_stats_.start_time)); } else { - settings_.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f\n", - compute_user_objective(original_lp_, obj), - toc(exploration_stats_.start_time)); + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t user_lower = get_lower_bound(); + if (!std::isfinite(user_lower)) { + f_t root_lower = root_lp_current_lower_bound_.load(); + if (std::isfinite(root_lower)) { user_lower = root_lower; } + } + std::string gap_str = std::isfinite(user_lower) + ? (". Gap: " + user_mip_gap(user_obj, user_lower) + "\n") + : "\n"; + settings_.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f%s", + user_obj, + toc(exploration_stats_.start_time), + gap_str.c_str()); } } @@ -1943,6 +1954,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut log.log_prefix = settings_.log.log_prefix; solver_status_ = mip_status_t::UNSET; is_running_ = false; + root_lp_current_lower_bound_ = -inf; exploration_stats_.nodes_unexplored = 0; exploration_stats_.nodes_explored = 0; original_lp_.A.to_compressed_row(Arow_); @@ -1972,6 +1984,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut lp_settings.inside_mip = 1; lp_settings.scale_columns = false; lp_settings.concurrent_halt = get_root_concurrent_halt(); + lp_settings.root_lp_progress_callback = [this](f_t user_obj) { + root_lp_current_lower_bound_.store(user_obj); + }; 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); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index a13d5cedcf..ec73792502 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -190,6 +190,7 @@ class branch_and_bound_t { lp_solution_t root_crossover_soln_; std::vector edge_norms_; std::atomic root_crossover_solution_set_{false}; + omp_atomic_t root_lp_current_lower_bound_; bool enable_concurrent_lp_root_solve_{false}; std::atomic root_concurrent_halt_{0}; bool is_root_solution_set{false}; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 426d9a7535..8c575f884f 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -3520,6 +3520,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, primal_infeasibility_squared, sum_perturb, now); + if (phase == 2 && settings.inside_mip == 1 && settings.root_lp_progress_callback) { + settings.root_lp_progress_callback(compute_user_objective(lp, obj)); + } } if (obj >= settings.cut_off) { diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 815e229232..bdf392d813 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -110,6 +110,7 @@ struct simplex_solver_settings_t { sub_mip(0), solution_callback(nullptr), heuristic_preemption_callback(nullptr), + root_lp_progress_callback(nullptr), concurrent_halt(nullptr) { } @@ -202,6 +203,8 @@ struct simplex_solver_settings_t { std::function&, f_t)> node_processed_callback; std::function heuristic_preemption_callback; std::function&, std::vector&, f_t)> set_simplex_solution_callback; + std::function + root_lp_progress_callback; // Called with current dual obj during root LP mutable logger_t log; std::atomic* concurrent_halt; // if nullptr ignored, if !nullptr, 0 if solver should // continue, 1 if solver should halt From a7e4ad3f5925776877e09a03bd551e0a672101b7 Mon Sep 17 00:00:00 2001 From: anandhkb Date: Sun, 8 Mar 2026 18:08:45 -0700 Subject: [PATCH 2/9] Address CodeRabbit review: user-space conversion, root_objective_ reset, callback every iteration - Convert get_lower_bound() to user space before user_mip_gap in report_heuristic - Reset root_objective_ at solve start to avoid stale bound on B&B reuse - Invoke root_lp_progress_callback every iteration, not just at log points --- cpp/src/branch_and_bound/branch_and_bound.cpp | 11 +++++------ cpp/src/dual_simplex/phase2.cpp | 6 +++--- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index f09158e7f2..bd4b3e2e38 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -319,12 +319,10 @@ void branch_and_bound_t::report_heuristic(f_t obj) user_gap.c_str(), toc(exploration_stats_.start_time)); } else { - f_t user_obj = compute_user_objective(original_lp_, obj); - f_t user_lower = get_lower_bound(); - if (!std::isfinite(user_lower)) { - f_t root_lower = root_lp_current_lower_bound_.load(); - if (std::isfinite(root_lower)) { user_lower = root_lower; } - } + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t lower = get_lower_bound(); + f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) + : root_lp_current_lower_bound_.load(); std::string gap_str = std::isfinite(user_lower) ? (". Gap: " + user_mip_gap(user_obj, user_lower) + "\n") : "\n"; @@ -1955,6 +1953,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut solver_status_ = mip_status_t::UNSET; is_running_ = false; root_lp_current_lower_bound_ = -inf; + root_objective_ = std::numeric_limits::quiet_NaN(); exploration_stats_.nodes_unexplored = 0; exploration_stats_.nodes_explored = 0; original_lp_.A.to_compressed_row(Arow_); diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 8c575f884f..07fc944333 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -3520,9 +3520,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, primal_infeasibility_squared, sum_perturb, now); - if (phase == 2 && settings.inside_mip == 1 && settings.root_lp_progress_callback) { - settings.root_lp_progress_callback(compute_user_objective(lp, obj)); - } + } + if (phase == 2 && settings.inside_mip == 1 && settings.root_lp_progress_callback) { + settings.root_lp_progress_callback(compute_user_objective(lp, obj)); } if (obj >= settings.cut_off) { From 7de698551765cb2757cf525f4f8866f5d0148e69 Mon Sep 17 00:00:00 2001 From: anandhkb Date: Sun, 8 Mar 2026 19:14:42 -0700 Subject: [PATCH 3/9] Use root_lp_progress_callback bound during root cut passes (CodeRabbit comment 4) When nodes_explored==0, prefer root_lp_current_lower_bound_ over get_lower_bound() when it provides a tighter bound. Fixes stale gap for heuristics during cut-pass solve. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index bd4b3e2e38..70baa3e0a0 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -308,8 +308,16 @@ template void branch_and_bound_t::report_heuristic(f_t obj) { if (is_running_) { - f_t user_obj = compute_user_objective(original_lp_, obj); - f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t lower = get_lower_bound(); + f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; + if (exploration_stats_.nodes_explored == 0) { + f_t root_progress = root_lp_current_lower_bound_.load(); + if (std::isfinite(root_progress) && + (!std::isfinite(user_lower) || root_progress > user_lower)) { + user_lower = root_progress; + } + } std::string user_gap = user_mip_gap(user_obj, user_lower); settings_.log.printf( From 971f2d1f27d13acc0e2370d3548b11598af946aa Mon Sep 17 00:00:00 2001 From: anandhkb Date: Mon, 9 Mar 2026 14:47:14 -0700 Subject: [PATCH 4/9] CodeRabbit/chris-maes: objective-sense aware bound, reset lower_bound_ceiling_, time last in heuristic log --- cpp/src/branch_and_bound/branch_and_bound.cpp | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 70baa3e0a0..93bf2c2aa6 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -312,9 +312,11 @@ void branch_and_bound_t::report_heuristic(f_t obj) f_t lower = get_lower_bound(); f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; if (exploration_stats_.nodes_explored == 0) { - f_t root_progress = root_lp_current_lower_bound_.load(); + f_t root_progress = root_lp_current_lower_bound_.load(); + const bool is_maximization = original_lp_.obj_scale < 0.0; if (std::isfinite(root_progress) && - (!std::isfinite(user_lower) || root_progress > user_lower)) { + (!std::isfinite(user_lower) || + (is_maximization ? root_progress < user_lower : root_progress > user_lower))) { user_lower = root_progress; } } @@ -327,17 +329,16 @@ void branch_and_bound_t::report_heuristic(f_t obj) user_gap.c_str(), toc(exploration_stats_.start_time)); } else { - f_t user_obj = compute_user_objective(original_lp_, obj); - f_t lower = get_lower_bound(); - f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) - : root_lp_current_lower_bound_.load(); - std::string gap_str = std::isfinite(user_lower) - ? (". Gap: " + user_mip_gap(user_obj, user_lower) + "\n") - : "\n"; - settings_.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f%s", + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t lower = get_lower_bound(); + f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) + : root_lp_current_lower_bound_.load(); + std::string gap_str = + std::isfinite(user_lower) ? (". Gap: " + user_mip_gap(user_obj, user_lower)) : ""; + settings_.log.printf("New solution from primal heuristics. Objective %+.6e%s. Time %.2f\n", user_obj, - toc(exploration_stats_.start_time), - gap_str.c_str()); + gap_str.c_str(), + toc(exploration_stats_.start_time)); } } @@ -1960,6 +1961,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut log.log_prefix = settings_.log.log_prefix; solver_status_ = mip_status_t::UNSET; is_running_ = false; + lower_bound_ceiling_ = inf; root_lp_current_lower_bound_ = -inf; root_objective_ = std::numeric_limits::quiet_NaN(); exploration_stats_.nodes_unexplored = 0; From d7ddec79e6c8031e91b0e9b8c7e5f74aff70bf8c Mon Sep 17 00:00:00 2001 From: anandhkb Date: Mon, 9 Mar 2026 15:13:02 -0700 Subject: [PATCH 5/9] Chris Maes: dual_simplex_objective_callback, solving_root_relaxation_, unify report_heuristic --- cpp/src/branch_and_bound/branch_and_bound.cpp | 67 +++++++++++-------- cpp/src/branch_and_bound/branch_and_bound.hpp | 1 + cpp/src/dual_simplex/phase2.cpp | 9 +-- .../dual_simplex/simplex_solver_settings.hpp | 5 +- 4 files changed, 48 insertions(+), 34 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 93bf2c2aa6..a1c5d775ef 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -290,6 +290,10 @@ branch_and_bound_t::branch_and_bound_t( template f_t branch_and_bound_t::get_lower_bound() { + if (solving_root_relaxation_.load()) { + f_t root = root_lp_current_lower_bound_.load(); + return std::isfinite(root) ? root : -inf; + } f_t lower_bound = lower_bound_ceiling_.load(); f_t heap_lower_bound = node_queue_.get_lower_bound(); lower_bound = std::min(heap_lower_bound, lower_bound); @@ -307,21 +311,21 @@ f_t branch_and_bound_t::get_lower_bound() template void branch_and_bound_t::report_heuristic(f_t obj) { - if (is_running_) { - f_t user_obj = compute_user_objective(original_lp_, obj); - f_t lower = get_lower_bound(); - f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; - if (exploration_stats_.nodes_explored == 0) { - f_t root_progress = root_lp_current_lower_bound_.load(); - const bool is_maximization = original_lp_.obj_scale < 0.0; - if (std::isfinite(root_progress) && - (!std::isfinite(user_lower) || - (is_maximization ? root_progress < user_lower : root_progress > user_lower))) { - user_lower = root_progress; - } + const f_t user_obj = compute_user_objective(original_lp_, obj); + f_t lower = get_lower_bound(); + f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; + if (!solving_root_relaxation_.load() && exploration_stats_.nodes_explored == 0) { + f_t root_progress = root_lp_current_lower_bound_.load(); + const bool is_maximization = original_lp_.obj_scale < 0.0; + if (std::isfinite(root_progress) && + (!std::isfinite(user_lower) || + (is_maximization ? root_progress < user_lower : root_progress > user_lower))) { + user_lower = root_progress; } - std::string user_gap = user_mip_gap(user_obj, user_lower); + } + if (is_running_) { + std::string user_gap = user_mip_gap(user_obj, user_lower); settings_.log.printf( "H %+13.6e %+10.6e %s %9.2f\n", user_obj, @@ -329,10 +333,6 @@ void branch_and_bound_t::report_heuristic(f_t obj) user_gap.c_str(), toc(exploration_stats_.start_time)); } else { - f_t user_obj = compute_user_objective(original_lp_, obj); - f_t lower = get_lower_bound(); - f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) - : root_lp_current_lower_bound_.load(); std::string gap_str = std::isfinite(user_lower) ? (". Gap: " + user_mip_gap(user_obj, user_lower)) : ""; settings_.log.printf("New solution from primal heuristics. Objective %+.6e%s. Time %.2f\n", @@ -1961,6 +1961,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut log.log_prefix = settings_.log.log_prefix; solver_status_ = mip_status_t::UNSET; is_running_ = false; + solving_root_relaxation_ = false; lower_bound_ceiling_ = inf; root_lp_current_lower_bound_ = -inf; root_objective_ = std::numeric_limits::quiet_NaN(); @@ -1988,12 +1989,13 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); - i_t original_rows = original_lp_.num_rows; - simplex_solver_settings_t lp_settings = settings_; - lp_settings.inside_mip = 1; - lp_settings.scale_columns = false; - lp_settings.concurrent_halt = get_root_concurrent_halt(); - lp_settings.root_lp_progress_callback = [this](f_t user_obj) { + solving_root_relaxation_ = true; + i_t original_rows = original_lp_.num_rows; + simplex_solver_settings_t lp_settings = settings_; + lp_settings.inside_mip = 1; + lp_settings.scale_columns = false; + lp_settings.concurrent_halt = get_root_concurrent_halt(); + lp_settings.dual_simplex_objective_callback = [this](f_t user_obj) { root_lp_current_lower_bound_.store(user_obj); }; std::vector basic_list(original_lp_.num_rows); @@ -2027,6 +2029,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time); if (root_status == lp_status_t::INFEASIBLE) { + solving_root_relaxation_ = false; settings_.log.printf("MIP Infeasible\n"); // FIXME: rarely dual simplex detects infeasible whereas it is feasible. // to add a small safety net, check if there is a primal solution already. @@ -2037,6 +2040,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::INFEASIBLE; } if (root_status == lp_status_t::UNBOUNDED) { + solving_root_relaxation_ = false; settings_.log.printf("MIP Unbounded\n"); if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); @@ -2044,19 +2048,22 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::UNBOUNDED; } if (root_status == lp_status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; + solving_root_relaxation_ = false; + solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, -inf); return solver_status_; } if (root_status == lp_status_t::WORK_LIMIT) { - solver_status_ = mip_status_t::WORK_LIMIT; + solving_root_relaxation_ = false; + solver_status_ = mip_status_t::WORK_LIMIT; set_final_solution(solution, -inf); return solver_status_; } if (root_status == lp_status_t::NUMERICAL_ISSUES) { - solver_status_ = mip_status_t::NUMERICAL; + solving_root_relaxation_ = false; + solver_status_ = mip_status_t::NUMERICAL; set_final_solution(solution, -inf); return solver_status_; } @@ -2087,6 +2094,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut cut_info_t cut_info; if (num_fractional == 0) { + solving_root_relaxation_ = false; set_solution_at_root(solution, cut_info); return mip_status_t::OPTIMAL; } @@ -2117,6 +2125,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t cut_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { if (num_fractional == 0) { + solving_root_relaxation_ = false; set_solution_at_root(solution, cut_info); return mip_status_t::OPTIMAL; } else { @@ -2278,7 +2287,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); } if (cut_status == dual::status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; + solving_root_relaxation_ = false; + solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, root_objective_); return solver_status_; } @@ -2301,6 +2311,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.total_lp_iters += root_relax_soln_.iterations; root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); } else { + solving_root_relaxation_ = false; settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); return mip_status_t::NUMERICAL; } @@ -2343,6 +2354,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); f_t abs_gap = upper_bound_.load() - root_objective_; if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { + solving_root_relaxation_ = false; set_solution_at_root(solution, cut_info); set_final_solution(solution, root_objective_); return mip_status_t::OPTIMAL; @@ -2362,6 +2374,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } } + solving_root_relaxation_ = false; print_cut_info(settings_, cut_info); if (cut_info.has_cuts()) { diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index ec73792502..8b083acd6e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -191,6 +191,7 @@ class branch_and_bound_t { std::vector edge_norms_; std::atomic root_crossover_solution_set_{false}; omp_atomic_t root_lp_current_lower_bound_; + omp_atomic_t solving_root_relaxation_{false}; bool enable_concurrent_lp_root_solve_{false}; std::atomic root_concurrent_halt_{0}; bool is_root_solution_set{false}; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 07fc944333..111500fbee 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -3510,19 +3510,20 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, if ((iter - start_iter) < settings.first_iteration_log || (iter % settings.iteration_log_frequency) == 0) { + const f_t user_obj = compute_user_objective(lp, obj); if (phase == 1 && iter == 1) { settings.log.printf(" Iter Objective Num Inf. Sum Inf. Perturb Time\n"); } settings.log.printf("%5d %+.16e %7d %.8e %.2e %.2f\n", iter, - compute_user_objective(lp, obj), + user_obj, infeasibility_indices.size(), primal_infeasibility_squared, sum_perturb, now); - } - if (phase == 2 && settings.inside_mip == 1 && settings.root_lp_progress_callback) { - settings.root_lp_progress_callback(compute_user_objective(lp, obj)); + if (phase == 2 && settings.inside_mip == 1 && settings.dual_simplex_objective_callback) { + settings.dual_simplex_objective_callback(user_obj); + } } if (obj >= settings.cut_off) { diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index bdf392d813..3f04f8b0a6 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -110,7 +110,7 @@ struct simplex_solver_settings_t { sub_mip(0), solution_callback(nullptr), heuristic_preemption_callback(nullptr), - root_lp_progress_callback(nullptr), + dual_simplex_objective_callback(nullptr), concurrent_halt(nullptr) { } @@ -203,8 +203,7 @@ struct simplex_solver_settings_t { std::function&, f_t)> node_processed_callback; std::function heuristic_preemption_callback; std::function&, std::vector&, f_t)> set_simplex_solution_callback; - std::function - root_lp_progress_callback; // Called with current dual obj during root LP + std::function dual_simplex_objective_callback; // Called with current dual obj mutable logger_t log; std::atomic* concurrent_halt; // if nullptr ignored, if !nullptr, 0 if solver should // continue, 1 if solver should halt From bc0927097c9e8abe81aad28f971a02df906d263b Mon Sep 17 00:00:00 2001 From: anandhkb Date: Tue, 10 Mar 2026 21:10:13 -0700 Subject: [PATCH 6/9] Chris Maes: RAII guard for solving_root_relaxation_, simplify report_heuristic - Add solving_root_relaxation_guard_t for single-point reset - Move solving_root_relaxation_ = true to just before root LP solve - Remove early return in get_lower_bound(); handle in report_heuristic - Refactor report_heuristic: shared user_obj/user_lower, simplified else branch --- cpp/src/branch_and_bound/branch_and_bound.cpp | 903 +++++++++--------- 1 file changed, 455 insertions(+), 448 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a1c5d775ef..08ca97ba1a 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -235,6 +235,15 @@ inline char feasible_solution_symbol(search_strategy_t strategy) } #endif +// RAII guard: sets solving_root_relaxation_ to false on destruction. +struct solving_root_relaxation_guard_t { + omp_atomic_t* flag_; + explicit solving_root_relaxation_guard_t(omp_atomic_t& flag) : flag_(&flag) {} + ~solving_root_relaxation_guard_t() { flag_->store(false); } + solving_root_relaxation_guard_t(const solving_root_relaxation_guard_t&) = delete; + solving_root_relaxation_guard_t& operator=(const solving_root_relaxation_guard_t&) = delete; +}; + } // namespace template @@ -290,10 +299,6 @@ branch_and_bound_t::branch_and_bound_t( template f_t branch_and_bound_t::get_lower_bound() { - if (solving_root_relaxation_.load()) { - f_t root = root_lp_current_lower_bound_.load(); - return std::isfinite(root) ? root : -inf; - } f_t lower_bound = lower_bound_ceiling_.load(); f_t heap_lower_bound = node_queue_.get_lower_bound(); lower_bound = std::min(heap_lower_bound, lower_bound); @@ -312,16 +317,12 @@ template void branch_and_bound_t::report_heuristic(f_t obj) { const f_t user_obj = compute_user_objective(original_lp_, obj); - f_t lower = get_lower_bound(); - f_t user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; - if (!solving_root_relaxation_.load() && exploration_stats_.nodes_explored == 0) { - f_t root_progress = root_lp_current_lower_bound_.load(); - const bool is_maximization = original_lp_.obj_scale < 0.0; - if (std::isfinite(root_progress) && - (!std::isfinite(user_lower) || - (is_maximization ? root_progress < user_lower : root_progress > user_lower))) { - user_lower = root_progress; - } + f_t user_lower; + if (solving_root_relaxation_.load()) { + user_lower = root_lp_current_lower_bound_.load(); + } else { + const f_t lower = get_lower_bound(); + user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; } if (is_running_) { @@ -333,12 +334,19 @@ void branch_and_bound_t::report_heuristic(f_t obj) user_gap.c_str(), toc(exploration_stats_.start_time)); } else { - std::string gap_str = - std::isfinite(user_lower) ? (". Gap: " + user_mip_gap(user_obj, user_lower)) : ""; - settings_.log.printf("New solution from primal heuristics. Objective %+.6e%s. Time %.2f\n", - user_obj, - gap_str.c_str(), - toc(exploration_stats_.start_time)); + if (solving_root_relaxation_.load()) { + std::string gap_str = + std::isfinite(user_lower) ? user_mip_gap(user_obj, user_lower) : " - "; + settings_.log.printf( + "New solution from primal heuristics. Objective %+.6e. Gap %s. Time %.2f\n", + user_obj, + gap_str.c_str(), + toc(exploration_stats_.start_time)); + } else { + settings_.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f\n", + user_obj, + toc(exploration_stats_.start_time)); + } } } @@ -1989,7 +1997,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); - solving_root_relaxation_ = true; i_t original_rows = original_lp_.num_rows; simplex_solver_settings_t lp_settings = settings_; lp_settings.inside_mip = 1; @@ -2003,491 +2010,491 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basis_update_mpf_t basis_update(original_lp_.num_rows, settings_.refactor_frequency); lp_status_t root_status; - 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_); - } - exploration_stats_.total_lp_iters = root_relax_soln_.iterations; - exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time); - - if (root_status == lp_status_t::INFEASIBLE) { - solving_root_relaxation_ = false; - settings_.log.printf("MIP Infeasible\n"); - // FIXME: rarely dual simplex detects infeasible whereas it is feasible. - // to add a small safety net, check if there is a primal solution already. - // Uncomment this if the issue with cost266-UUE is resolved - // if (settings.heuristic_preemption_callback != nullptr) { - // settings.heuristic_preemption_callback(); - // } - return mip_status_t::INFEASIBLE; - } - if (root_status == lp_status_t::UNBOUNDED) { - solving_root_relaxation_ = false; - settings_.log.printf("MIP Unbounded\n"); - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); + { + solving_root_relaxation_ = true; + solving_root_relaxation_guard_t root_guard(solving_root_relaxation_); + + 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_); + } + exploration_stats_.total_lp_iters = root_relax_soln_.iterations; + exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time); + + if (root_status == lp_status_t::INFEASIBLE) { + settings_.log.printf("MIP Infeasible\n"); + // FIXME: rarely dual simplex detects infeasible whereas it is feasible. + // to add a small safety net, check if there is a primal solution already. + // Uncomment this if the issue with cost266-UUE is resolved + // if (settings.heuristic_preemption_callback != nullptr) { + // settings.heuristic_preemption_callback(); + // } + return mip_status_t::INFEASIBLE; + } + if (root_status == lp_status_t::UNBOUNDED) { + settings_.log.printf("MIP Unbounded\n"); + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + return mip_status_t::UNBOUNDED; + } + if (root_status == lp_status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, -inf); + return solver_status_; } - return mip_status_t::UNBOUNDED; - } - if (root_status == lp_status_t::TIME_LIMIT) { - solving_root_relaxation_ = false; - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, -inf); - return solver_status_; - } - if (root_status == lp_status_t::WORK_LIMIT) { - solving_root_relaxation_ = false; - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, -inf); - return solver_status_; - } + if (root_status == lp_status_t::WORK_LIMIT) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, -inf); + return solver_status_; + } - if (root_status == lp_status_t::NUMERICAL_ISSUES) { - solving_root_relaxation_ = false; - solver_status_ = mip_status_t::NUMERICAL; - set_final_solution(solution, -inf); - return solver_status_; - } + if (root_status == lp_status_t::NUMERICAL_ISSUES) { + solver_status_ = mip_status_t::NUMERICAL; + set_final_solution(solution, -inf); + return solver_status_; + } - assert(root_vstatus_.size() == original_lp_.num_cols); - set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); + assert(root_vstatus_.size() == original_lp_.num_cols); + set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - if (settings_.set_simplex_solution_callback != nullptr) { - std::vector original_x; - uncrush_primal_solution(original_problem_, original_lp_, root_relax_soln_.x, original_x); - std::vector original_dual; - std::vector original_z; - uncrush_dual_solution(original_problem_, - original_lp_, - root_relax_soln_.y, - root_relax_soln_.z, - original_dual, - original_z); - settings_.set_simplex_solution_callback( - original_x, original_dual, compute_user_objective(original_lp_, root_objective_)); - } + if (settings_.set_simplex_solution_callback != nullptr) { + std::vector original_x; + uncrush_primal_solution(original_problem_, original_lp_, root_relax_soln_.x, original_x); + std::vector original_dual; + std::vector original_z; + uncrush_dual_solution(original_problem_, + original_lp_, + root_relax_soln_.y, + root_relax_soln_.z, + original_dual, + original_z); + settings_.set_simplex_solution_callback( + original_x, original_dual, compute_user_objective(original_lp_, root_objective_)); + } - std::vector fractional; - i_t num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + std::vector fractional; + i_t num_fractional = + fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - cut_info_t cut_info; + cut_info_t cut_info; - if (num_fractional == 0) { - solving_root_relaxation_ = false; - set_solution_at_root(solution, cut_info); - return mip_status_t::OPTIMAL; - } + if (num_fractional == 0) { + set_solution_at_root(solution, cut_info); + return mip_status_t::OPTIMAL; + } - is_running_ = true; - lower_bound_ceiling_ = inf; + is_running_ = true; + lower_bound_ceiling_ = inf; - if (num_fractional != 0 && settings_.max_cut_passes > 0) { - settings_.log.printf( - " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node | " - "Gap " - "| Time |\n"); - } + if (num_fractional != 0 && settings_.max_cut_passes > 0) { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node | " + " " + "Gap " + "| Time |\n"); + } - cut_pool_t cut_pool(original_lp_.num_cols, settings_); - cut_generation_t cut_generation( - cut_pool, original_lp_, settings_, Arow_, new_slacks_, var_types_); + cut_pool_t cut_pool(original_lp_.num_cols, settings_); + cut_generation_t cut_generation( + cut_pool, original_lp_, settings_, Arow_, new_slacks_, var_types_); - std::vector saved_solution; + std::vector saved_solution; #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - read_saved_solution_for_cut_verification(original_lp_, settings_, saved_solution); + read_saved_solution_for_cut_verification(original_lp_, settings_, saved_solution); #endif - f_t last_upper_bound = std::numeric_limits::infinity(); - f_t last_objective = root_objective_; - f_t root_relax_objective = root_objective_; + f_t last_upper_bound = std::numeric_limits::infinity(); + f_t last_objective = root_objective_; + f_t root_relax_objective = root_objective_; - i_t cut_pool_size = 0; - for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { - if (num_fractional == 0) { - solving_root_relaxation_ = false; - set_solution_at_root(solution, cut_info); - return mip_status_t::OPTIMAL; - } else { + i_t cut_pool_size = 0; + for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { + if (num_fractional == 0) { + set_solution_at_root(solution, cut_info); + return mip_status_t::OPTIMAL; + } else { #ifdef PRINT_FRACTIONAL_INFO - settings_.log.printf( - "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); - for (i_t j : fractional) { - settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", - j, - original_lp_.lower[j], - root_relax_soln_.x[j], - original_lp_.upper[j]); - } + settings_.log.printf( + "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); + for (i_t j : fractional) { + settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", + j, + original_lp_.lower[j], + root_relax_soln_.x[j], + original_lp_.upper[j]); + } #endif - // Generate cuts and add them to the cut pool - f_t cut_start_time = tic(); - cut_generation.generate_cuts(original_lp_, - settings_, - Arow_, - new_slacks_, - var_types_, - basis_update, - root_relax_soln_.x, - basic_list, - nonbasic_list); - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); - } - // Score the cuts - f_t score_start_time = tic(); - cut_pool.score_cuts(root_relax_soln_.x); - f_t score_time = toc(score_start_time); - if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } - // Get the best cuts from the cut pool - csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); - std::vector cut_rhs; - std::vector cut_types; - i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); - if (num_cuts == 0) { break; } - cut_info.record_cut_types(cut_types); + // Generate cuts and add them to the cut pool + f_t cut_start_time = tic(); + cut_generation.generate_cuts(original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + basis_update, + root_relax_soln_.x, + basic_list, + nonbasic_list); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); + } + // Score the cuts + f_t score_start_time = tic(); + cut_pool.score_cuts(root_relax_soln_.x); + f_t score_time = toc(score_start_time); + if (score_time > 1.0) { + settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); + } + // Get the best cuts from the cut pool + csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); + std::vector cut_rhs; + std::vector cut_types; + i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); + if (num_cuts == 0) { break; } + cut_info.record_cut_types(cut_types); #ifdef PRINT_CUT_POOL_TYPES - cut_pool.print_cutpool_types(); - print_cut_types("In LP ", cut_types, settings_); - printf("Cut pool size: %d\n", cut_pool.pool_size()); + cut_pool.print_cutpool_types(); + print_cut_types("In LP ", cut_types, settings_); + printf("Cut pool size: %d\n", cut_pool.pool_size()); #endif #ifdef CHECK_CUT_MATRIX - if (cuts_to_add.check_matrix() != 0) { - settings_.log.printf("Bad cuts matrix\n"); - for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { - settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); + if (cuts_to_add.check_matrix() != 0) { + settings_.log.printf("Bad cuts matrix\n"); + for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { + settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); + } + return mip_status_t::NUMERICAL; } - return mip_status_t::NUMERICAL; - } #endif - // Check against saved solution + // Check against saved solution #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); + verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); #endif - cut_pool_size = cut_pool.pool_size(); - - // Resolve the LP with the new cuts - settings_.log.debug( - "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", - num_cuts, - cuts_to_add.row_start[cuts_to_add.m], - cut_pool.pool_size(), - cuts_to_add.m + original_lp_.num_rows); - lp_settings.log.log = false; - - f_t add_cuts_start_time = tic(); - mutex_original_lp_.lock(); - i_t add_cuts_status = add_cuts(settings_, - cuts_to_add, - cut_rhs, - original_lp_, - new_slacks_, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); - mutex_original_lp_.unlock(); - f_t add_cuts_time = toc(add_cuts_start_time); - if (add_cuts_time > 1.0) { - settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); - } - if (add_cuts_status != 0) { - settings_.log.printf("Failed to add cuts\n"); - return mip_status_t::NUMERICAL; - } + cut_pool_size = cut_pool.pool_size(); - if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { - mutex_upper_.lock(); - last_upper_bound = upper_bound_.load(); - std::vector lower_bounds; - std::vector upper_bounds; - find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - mutex_upper_.unlock(); + // Resolve the LP with the new cuts + settings_.log.debug( + "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", + num_cuts, + cuts_to_add.row_start[cuts_to_add.m], + cut_pool.pool_size(), + cuts_to_add.m + original_lp_.num_rows); + lp_settings.log.log = false; + + f_t add_cuts_start_time = tic(); mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; + i_t add_cuts_status = add_cuts(settings_, + cuts_to_add, + cut_rhs, + original_lp_, + new_slacks_, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); mutex_original_lp_.unlock(); - } + f_t add_cuts_time = toc(add_cuts_start_time); + if (add_cuts_time > 1.0) { + settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); + } + if (add_cuts_status != 0) { + settings_.log.printf("Failed to add cuts\n"); + return mip_status_t::NUMERICAL; + } - // Try to do bound strengthening - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; + if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { + mutex_upper_.lock(); + last_upper_bound = upper_bound_.load(); + std::vector lower_bounds; + std::vector upper_bounds; + find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + mutex_upper_.unlock(); + mutex_original_lp_.lock(); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; + mutex_original_lp_.unlock(); + } + + // Try to do bound strengthening + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; #ifdef CHECK_MATRICES - settings_.log.printf("Before A check\n"); - original_lp_.A.check_matrix(); + settings_.log.printf("Before A check\n"); + original_lp_.A.check_matrix(); #endif - original_lp_.A.to_compressed_row(Arow_); - - f_t node_presolve_start_time = tic(); - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - std::vector new_lower = original_lp_.lower; - std::vector new_upper = original_lp_.upper; - bool feasible = - node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); - mutex_original_lp_.lock(); - original_lp_.lower = new_lower; - original_lp_.upper = new_upper; - mutex_original_lp_.unlock(); - f_t node_presolve_time = toc(node_presolve_start_time); - if (node_presolve_time > 1.0) { - settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); - } - if (!feasible) { - settings_.log.printf("Bound strengthening detected infeasibility\n"); - return mip_status_t::INFEASIBLE; - } + original_lp_.A.to_compressed_row(Arow_); + + f_t node_presolve_start_time = tic(); + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + std::vector new_lower = original_lp_.lower; + std::vector new_upper = original_lp_.upper; + bool feasible = + node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + f_t node_presolve_time = toc(node_presolve_start_time); + if (node_presolve_time > 1.0) { + settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); + } + if (!feasible) { + settings_.log.printf("Bound strengthening detected infeasibility\n"); + return mip_status_t::INFEASIBLE; + } - i_t iter = 0; - bool initialize_basis = false; - lp_settings.concurrent_halt = NULL; - f_t dual_phase2_start_time = tic(); - dual::status_t cut_status = dual_phase2_with_advanced_basis(2, - 0, - initialize_basis, - exploration_stats_.start_time, - original_lp_, - lp_settings, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - root_relax_soln_, - iter, - edge_norms_); - exploration_stats_.total_lp_iters += iter; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - f_t dual_phase2_time = toc(dual_phase2_start_time); - if (dual_phase2_time > 1.0) { - settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); - } - if (cut_status == dual::status_t::TIME_LIMIT) { - solving_root_relaxation_ = false; - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return solver_status_; - } + i_t iter = 0; + bool initialize_basis = false; + lp_settings.concurrent_halt = NULL; + f_t dual_phase2_start_time = tic(); + dual::status_t cut_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + exploration_stats_.total_lp_iters += iter; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + f_t dual_phase2_time = toc(dual_phase2_start_time); + if (dual_phase2_time > 1.0) { + settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); + } + if (cut_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return solver_status_; + } - if (cut_status != dual::status_t::OPTIMAL) { - settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); - lp_status_t scratch_status = - solve_linear_program_with_advanced_basis(original_lp_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - if (scratch_status == lp_status_t::OPTIMAL) { - // We recovered - cut_status = convert_lp_status_to_dual_status(scratch_status); - exploration_stats_.total_lp_iters += root_relax_soln_.iterations; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - } else { - solving_root_relaxation_ = false; - settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); - return mip_status_t::NUMERICAL; + if (cut_status != dual::status_t::OPTIMAL) { + settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); + lp_status_t scratch_status = + solve_linear_program_with_advanced_basis(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + if (scratch_status == lp_status_t::OPTIMAL) { + // We recovered + cut_status = convert_lp_status_to_dual_status(scratch_status); + exploration_stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else { + settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); + return mip_status_t::NUMERICAL; + } } - } - f_t remove_cuts_start_time = tic(); - mutex_original_lp_.lock(); - remove_cuts(original_lp_, - settings_, - exploration_stats_.start_time, - Arow_, - new_slacks_, - original_rows, - var_types_, - root_vstatus_, - edge_norms_, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - basis_update); - mutex_original_lp_.unlock(); - f_t remove_cuts_time = toc(remove_cuts_start_time); - if (remove_cuts_time > 1.0) { - settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); - } - fractional.clear(); - num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + f_t remove_cuts_start_time = tic(); + mutex_original_lp_.lock(); + remove_cuts(original_lp_, + settings_, + exploration_stats_.start_time, + Arow_, + new_slacks_, + original_rows, + var_types_, + root_vstatus_, + edge_norms_, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + basis_update); + mutex_original_lp_.unlock(); + f_t remove_cuts_time = toc(remove_cuts_start_time); + if (remove_cuts_time > 1.0) { + settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); + } + fractional.clear(); + num_fractional = + fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + + if (num_fractional == 0) { + upper_bound_ = root_objective_; + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + mutex_upper_.unlock(); + } + f_t obj = upper_bound_.load(); + report(' ', obj, root_objective_, 0, num_fractional); + + f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); + f_t abs_gap = upper_bound_.load() - root_objective_; + if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { + set_solution_at_root(solution, cut_info); + set_final_solution(solution, root_objective_); + return mip_status_t::OPTIMAL; + } - if (num_fractional == 0) { - upper_bound_ = root_objective_; - mutex_upper_.lock(); - incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); - mutex_upper_.unlock(); + f_t change_in_objective = root_objective_ - last_objective; + const f_t factor = settings_.cut_change_threshold; + const f_t min_objective = 1e-3; + if (change_in_objective <= + factor * std::max(min_objective, std::abs(root_relax_objective))) { + settings_.log.debug( + "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", + change_in_objective, + root_relax_objective); + break; + } + last_objective = root_objective_; } - f_t obj = upper_bound_.load(); - report(' ', obj, root_objective_, 0, num_fractional); + } - f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); - f_t abs_gap = upper_bound_.load() - root_objective_; - if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { - solving_root_relaxation_ = false; - set_solution_at_root(solution, cut_info); - set_final_solution(solution, root_objective_); - return mip_status_t::OPTIMAL; - } + print_cut_info(settings_, cut_info); - f_t change_in_objective = root_objective_ - last_objective; - const f_t factor = settings_.cut_change_threshold; - const f_t min_objective = 1e-3; - if (change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { - settings_.log.debug( - "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", - change_in_objective, - root_relax_objective); - break; - } - last_objective = root_objective_; + if (cut_info.has_cuts()) { + settings_.log.printf("Cut pool size : %d\n", cut_pool_size); + settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", + original_lp_.num_rows, + original_lp_.num_cols, + original_lp_.A.col_start[original_lp_.A.n]); } - } - solving_root_relaxation_ = false; - print_cut_info(settings_, cut_info); + set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); - if (cut_info.has_cuts()) { - settings_.log.printf("Cut pool size : %d\n", cut_pool_size); - settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", - original_lp_.num_rows, - original_lp_.num_cols, - original_lp_.A.col_start[original_lp_.A.n]); - } + pc_.resize(original_lp_.num_cols); + { + raft::common::nvtx::range scope_sb("BB::strong_branching"); + strong_branching(original_problem_, + original_lp_, + settings_, + exploration_stats_.start_time, + var_types_, + root_relax_soln_.x, + fractional, + root_objective_, + root_vstatus_, + edge_norms_, + pc_); + } - set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return solver_status_; + } - pc_.resize(original_lp_.num_cols); - { - raft::common::nvtx::range scope_sb("BB::strong_branching"); - strong_branching(original_problem_, - original_lp_, - settings_, - exploration_stats_.start_time, - var_types_, - root_relax_soln_.x, - fractional, - root_objective_, - root_vstatus_, - edge_norms_, - pc_); - } + if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) { + std::vector lower_bounds; + std::vector upper_bounds; + i_t num_fixed = find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + if (num_fixed > 0) { + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; - if (toc(exploration_stats_.start_time) > settings_.time_limit) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return solver_status_; - } - - if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) { - std::vector lower_bounds; - std::vector upper_bounds; - i_t num_fixed = find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - if (num_fixed > 0) { - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; - - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - bool feasible = node_presolve.bounds_strengthening( - settings_, bounds_changed, original_lp_.lower, original_lp_.upper); - mutex_original_lp_.unlock(); - if (!feasible) { - settings_.log.printf("Bound strengthening failed\n"); - return mip_status_t::NUMERICAL; // We had a feasible integer solution, but bound - // strengthening thinks we are infeasible. - } - // Go through and check the fractional variables and remove any that are now fixed to their - // bounds - std::vector to_remove(fractional.size(), 0); - i_t num_to_remove = 0; - for (i_t k = 0; k < fractional.size(); k++) { - const i_t j = fractional[k]; - if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) { - to_remove[k] = 1; - num_to_remove++; + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + + mutex_original_lp_.lock(); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; + bool feasible = node_presolve.bounds_strengthening( + settings_, bounds_changed, original_lp_.lower, original_lp_.upper); + mutex_original_lp_.unlock(); + if (!feasible) { + settings_.log.printf("Bound strengthening failed\n"); + return mip_status_t::NUMERICAL; // We had a feasible integer solution, but bound + // strengthening thinks we are infeasible. } - } - if (num_to_remove > 0) { - std::vector new_fractional; - new_fractional.reserve(fractional.size() - num_to_remove); + // Go through and check the fractional variables and remove any that are now fixed to their + // bounds + std::vector to_remove(fractional.size(), 0); + i_t num_to_remove = 0; for (i_t k = 0; k < fractional.size(); k++) { - if (!to_remove[k]) { new_fractional.push_back(fractional[k]); } + const i_t j = fractional[k]; + if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) { + to_remove[k] = 1; + num_to_remove++; + } + } + if (num_to_remove > 0) { + std::vector new_fractional; + new_fractional.reserve(fractional.size() - num_to_remove); + for (i_t k = 0; k < fractional.size(); k++) { + if (!to_remove[k]) { new_fractional.push_back(fractional[k]); } + } + fractional = new_fractional; + num_fractional = fractional.size(); } - fractional = new_fractional; - num_fractional = fractional.size(); } } - } - // Choose variable to branch on - i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x, log); + // Choose variable to branch on + i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x, log); - search_tree_.root = std::move(mip_node_t(root_objective_, root_vstatus_)); - search_tree_.num_nodes = 0; - search_tree_.graphviz_node(settings_.log, &search_tree_.root, "lower bound", root_objective_); - search_tree_.branch(&search_tree_.root, - branch_var, - root_relax_soln_.x[branch_var], - num_fractional, - root_vstatus_, - original_lp_, - log); - node_queue_.push(search_tree_.root.get_down_child()); - node_queue_.push(search_tree_.root.get_up_child()); + search_tree_.root = std::move(mip_node_t(root_objective_, root_vstatus_)); + search_tree_.num_nodes = 0; + search_tree_.graphviz_node(settings_.log, &search_tree_.root, "lower bound", root_objective_); + search_tree_.branch(&search_tree_.root, + branch_var, + root_relax_soln_.x[branch_var], + num_fractional, + root_vstatus_, + original_lp_, + log); + node_queue_.push(search_tree_.root.get_down_child()); + node_queue_.push(search_tree_.root.get_up_child()); - settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads); + settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads); - exploration_stats_.nodes_explored = 0; - exploration_stats_.nodes_unexplored = 2; - exploration_stats_.nodes_since_last_log = 0; - exploration_stats_.last_log = tic(); - min_node_queue_size_ = 2 * settings_.num_threads; + exploration_stats_.nodes_explored = 0; + exploration_stats_.nodes_unexplored = 2; + exploration_stats_.nodes_since_last_log = 0; + exploration_stats_.last_log = tic(); + min_node_queue_size_ = 2 * settings_.num_threads; - if (settings_.diving_settings.coefficient_diving != 0) { - calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); - } - if (settings_.deterministic) { - settings_.log.printf( - " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " - "| Gap | Work | Time |\n"); - } else { - settings_.log.printf( - " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " - "| Gap | Time |\n"); + if (settings_.diving_settings.coefficient_diving != 0) { + calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); + } + if (settings_.deterministic) { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " + "| Gap | Work | Time |\n"); + } else { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " + "| Gap | Time |\n"); + } } if (settings_.deterministic) { From f10d720ab9961cce372e6aa87ed4151f6ebe5cb6 Mon Sep 17 00:00:00 2001 From: anandhkb Date: Wed, 11 Mar 2026 10:13:38 -0700 Subject: [PATCH 7/9] Add dual-simplex callback and reset root LP bound at solve start - Set dual_simplex_objective_callback on lp_settings before root solve to populate root_lp_current_lower_bound_ during dual-simplex iterations - Reset root_lp_current_lower_bound_ to -inf at beginning of solve() with solver_status_, is_running_, etc. to avoid stale bound on reused solver --- cpp/src/branch_and_bound/branch_and_bound.cpp | 882 +++++++++--------- 1 file changed, 428 insertions(+), 454 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 08ca97ba1a..6baa14bc6a 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -235,15 +235,6 @@ inline char feasible_solution_symbol(search_strategy_t strategy) } #endif -// RAII guard: sets solving_root_relaxation_ to false on destruction. -struct solving_root_relaxation_guard_t { - omp_atomic_t* flag_; - explicit solving_root_relaxation_guard_t(omp_atomic_t& flag) : flag_(&flag) {} - ~solving_root_relaxation_guard_t() { flag_->store(false); } - solving_root_relaxation_guard_t(const solving_root_relaxation_guard_t&) = delete; - solving_root_relaxation_guard_t& operator=(const solving_root_relaxation_guard_t&) = delete; -}; - } // namespace template @@ -292,8 +283,8 @@ branch_and_bound_t::branch_and_bound_t( #endif upper_bound_ = inf; - root_lp_current_lower_bound_ = -inf; root_objective_ = std::numeric_limits::quiet_NaN(); + root_lp_current_lower_bound_ = -inf; } template @@ -316,17 +307,11 @@ f_t branch_and_bound_t::get_lower_bound() template void branch_and_bound_t::report_heuristic(f_t obj) { - const f_t user_obj = compute_user_objective(original_lp_, obj); - f_t user_lower; - if (solving_root_relaxation_.load()) { - user_lower = root_lp_current_lower_bound_.load(); - } else { - const f_t lower = get_lower_bound(); - user_lower = std::isfinite(lower) ? compute_user_objective(original_lp_, lower) : -inf; - } - if (is_running_) { + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string user_gap = user_mip_gap(user_obj, user_lower); + settings_.log.printf( "H %+13.6e %+10.6e %s %9.2f\n", user_obj, @@ -335,16 +320,17 @@ void branch_and_bound_t::report_heuristic(f_t obj) toc(exploration_stats_.start_time)); } else { if (solving_root_relaxation_.load()) { - std::string gap_str = - std::isfinite(user_lower) ? user_mip_gap(user_obj, user_lower) : " - "; + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t user_lower = root_lp_current_lower_bound_.load(); + std::string user_gap = user_mip_gap(user_obj, user_lower); settings_.log.printf( "New solution from primal heuristics. Objective %+.6e. Gap %s. Time %.2f\n", user_obj, - gap_str.c_str(), + user_gap.c_str(), toc(exploration_stats_.start_time)); } else { settings_.log.printf("New solution from primal heuristics. Objective %+.6e. Time %.2f\n", - user_obj, + compute_user_objective(original_lp_, obj), toc(exploration_stats_.start_time)); } } @@ -1969,10 +1955,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut log.log_prefix = settings_.log.log_prefix; solver_status_ = mip_status_t::UNSET; is_running_ = false; - solving_root_relaxation_ = false; - lower_bound_ceiling_ = inf; root_lp_current_lower_bound_ = -inf; - root_objective_ = std::numeric_limits::quiet_NaN(); exploration_stats_.nodes_unexplored = 0; exploration_stats_.nodes_explored = 0; original_lp_.A.to_compressed_row(Arow_); @@ -2009,492 +1992,483 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut std::vector nonbasic_list; basis_update_mpf_t basis_update(original_lp_.num_rows, settings_.refactor_frequency); lp_status_t root_status; - - { - solving_root_relaxation_ = true; - solving_root_relaxation_guard_t root_guard(solving_root_relaxation_); - - 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_); - } - exploration_stats_.total_lp_iters = root_relax_soln_.iterations; - exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time); - - if (root_status == lp_status_t::INFEASIBLE) { - settings_.log.printf("MIP Infeasible\n"); - // FIXME: rarely dual simplex detects infeasible whereas it is feasible. - // to add a small safety net, check if there is a primal solution already. - // Uncomment this if the issue with cost266-UUE is resolved - // if (settings.heuristic_preemption_callback != nullptr) { - // settings.heuristic_preemption_callback(); - // } - return mip_status_t::INFEASIBLE; - } - if (root_status == lp_status_t::UNBOUNDED) { - settings_.log.printf("MIP Unbounded\n"); - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); - } - return mip_status_t::UNBOUNDED; - } - if (root_status == lp_status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, -inf); - return solver_status_; + solving_root_relaxation_ = true; + + 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; + exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time); + + if (root_status == lp_status_t::INFEASIBLE) { + settings_.log.printf("MIP Infeasible\n"); + // FIXME: rarely dual simplex detects infeasible whereas it is feasible. + // to add a small safety net, check if there is a primal solution already. + // Uncomment this if the issue with cost266-UUE is resolved + // if (settings.heuristic_preemption_callback != nullptr) { + // settings.heuristic_preemption_callback(); + // } + return mip_status_t::INFEASIBLE; + } + if (root_status == lp_status_t::UNBOUNDED) { + settings_.log.printf("MIP Unbounded\n"); + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); } + return mip_status_t::UNBOUNDED; + } + if (root_status == lp_status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, -inf); + return solver_status_; + } - if (root_status == lp_status_t::WORK_LIMIT) { - solver_status_ = mip_status_t::WORK_LIMIT; - set_final_solution(solution, -inf); - return solver_status_; - } + if (root_status == lp_status_t::WORK_LIMIT) { + solver_status_ = mip_status_t::WORK_LIMIT; + set_final_solution(solution, -inf); + return solver_status_; + } - if (root_status == lp_status_t::NUMERICAL_ISSUES) { - solver_status_ = mip_status_t::NUMERICAL; - set_final_solution(solution, -inf); - return solver_status_; - } + if (root_status == lp_status_t::NUMERICAL_ISSUES) { + solver_status_ = mip_status_t::NUMERICAL; + set_final_solution(solution, -inf); + return solver_status_; + } - assert(root_vstatus_.size() == original_lp_.num_cols); - set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); + assert(root_vstatus_.size() == original_lp_.num_cols); + set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - if (settings_.set_simplex_solution_callback != nullptr) { - std::vector original_x; - uncrush_primal_solution(original_problem_, original_lp_, root_relax_soln_.x, original_x); - std::vector original_dual; - std::vector original_z; - uncrush_dual_solution(original_problem_, - original_lp_, - root_relax_soln_.y, - root_relax_soln_.z, - original_dual, - original_z); - settings_.set_simplex_solution_callback( - original_x, original_dual, compute_user_objective(original_lp_, root_objective_)); - } + if (settings_.set_simplex_solution_callback != nullptr) { + std::vector original_x; + uncrush_primal_solution(original_problem_, original_lp_, root_relax_soln_.x, original_x); + std::vector original_dual; + std::vector original_z; + uncrush_dual_solution(original_problem_, + original_lp_, + root_relax_soln_.y, + root_relax_soln_.z, + original_dual, + original_z); + settings_.set_simplex_solution_callback( + original_x, original_dual, compute_user_objective(original_lp_, root_objective_)); + } - std::vector fractional; - i_t num_fractional = - fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + std::vector fractional; + i_t num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - cut_info_t cut_info; + cut_info_t cut_info; - if (num_fractional == 0) { - set_solution_at_root(solution, cut_info); - return mip_status_t::OPTIMAL; - } + if (num_fractional == 0) { + set_solution_at_root(solution, cut_info); + return mip_status_t::OPTIMAL; + } - is_running_ = true; - lower_bound_ceiling_ = inf; + is_running_ = true; + lower_bound_ceiling_ = inf; - if (num_fractional != 0 && settings_.max_cut_passes > 0) { - settings_.log.printf( - " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node | " - " " - "Gap " - "| Time |\n"); - } + if (num_fractional != 0 && settings_.max_cut_passes > 0) { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node | " + "Gap " + "| Time |\n"); + } - cut_pool_t cut_pool(original_lp_.num_cols, settings_); - cut_generation_t cut_generation( - cut_pool, original_lp_, settings_, Arow_, new_slacks_, var_types_); + cut_pool_t cut_pool(original_lp_.num_cols, settings_); + cut_generation_t cut_generation( + cut_pool, original_lp_, settings_, Arow_, new_slacks_, var_types_); - std::vector saved_solution; + std::vector saved_solution; #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - read_saved_solution_for_cut_verification(original_lp_, settings_, saved_solution); + read_saved_solution_for_cut_verification(original_lp_, settings_, saved_solution); #endif - f_t last_upper_bound = std::numeric_limits::infinity(); - f_t last_objective = root_objective_; - f_t root_relax_objective = root_objective_; + f_t last_upper_bound = std::numeric_limits::infinity(); + f_t last_objective = root_objective_; + f_t root_relax_objective = root_objective_; - i_t cut_pool_size = 0; - for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { - if (num_fractional == 0) { - set_solution_at_root(solution, cut_info); - return mip_status_t::OPTIMAL; - } else { + i_t cut_pool_size = 0; + for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { + if (num_fractional == 0) { + set_solution_at_root(solution, cut_info); + return mip_status_t::OPTIMAL; + } else { #ifdef PRINT_FRACTIONAL_INFO - settings_.log.printf( - "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); - for (i_t j : fractional) { - settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", - j, - original_lp_.lower[j], - root_relax_soln_.x[j], - original_lp_.upper[j]); - } + settings_.log.printf( + "Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); + for (i_t j : fractional) { + settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", + j, + original_lp_.lower[j], + root_relax_soln_.x[j], + original_lp_.upper[j]); + } #endif - // Generate cuts and add them to the cut pool - f_t cut_start_time = tic(); - cut_generation.generate_cuts(original_lp_, - settings_, - Arow_, - new_slacks_, - var_types_, - basis_update, - root_relax_soln_.x, - basic_list, - nonbasic_list); - f_t cut_generation_time = toc(cut_start_time); - if (cut_generation_time > 1.0) { - settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); - } - // Score the cuts - f_t score_start_time = tic(); - cut_pool.score_cuts(root_relax_soln_.x); - f_t score_time = toc(score_start_time); - if (score_time > 1.0) { - settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); - } - // Get the best cuts from the cut pool - csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); - std::vector cut_rhs; - std::vector cut_types; - i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); - if (num_cuts == 0) { break; } - cut_info.record_cut_types(cut_types); + // Generate cuts and add them to the cut pool + f_t cut_start_time = tic(); + cut_generation.generate_cuts(original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + basis_update, + root_relax_soln_.x, + basic_list, + nonbasic_list); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); + } + // Score the cuts + f_t score_start_time = tic(); + cut_pool.score_cuts(root_relax_soln_.x); + f_t score_time = toc(score_start_time); + if (score_time > 1.0) { settings_.log.debug("Cut scoring time %.2f seconds\n", score_time); } + // Get the best cuts from the cut pool + csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); + std::vector cut_rhs; + std::vector cut_types; + i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); + if (num_cuts == 0) { break; } + cut_info.record_cut_types(cut_types); #ifdef PRINT_CUT_POOL_TYPES - cut_pool.print_cutpool_types(); - print_cut_types("In LP ", cut_types, settings_); - printf("Cut pool size: %d\n", cut_pool.pool_size()); + cut_pool.print_cutpool_types(); + print_cut_types("In LP ", cut_types, settings_); + printf("Cut pool size: %d\n", cut_pool.pool_size()); #endif #ifdef CHECK_CUT_MATRIX - if (cuts_to_add.check_matrix() != 0) { - settings_.log.printf("Bad cuts matrix\n"); - for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { - settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); - } - return mip_status_t::NUMERICAL; + if (cuts_to_add.check_matrix() != 0) { + settings_.log.printf("Bad cuts matrix\n"); + for (i_t i = 0; i < static_cast(cut_types.size()); ++i) { + settings_.log.printf("row %d cut type %d\n", i, cut_types[i]); } + return mip_status_t::NUMERICAL; + } #endif - // Check against saved solution + // Check against saved solution #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION - verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); + verify_cuts_against_saved_solution(cuts_to_add, cut_rhs, saved_solution); #endif - cut_pool_size = cut_pool.pool_size(); + cut_pool_size = cut_pool.pool_size(); + + // Resolve the LP with the new cuts + settings_.log.debug( + "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", + num_cuts, + cuts_to_add.row_start[cuts_to_add.m], + cut_pool.pool_size(), + cuts_to_add.m + original_lp_.num_rows); + lp_settings.log.log = false; + + f_t add_cuts_start_time = tic(); + mutex_original_lp_.lock(); + i_t add_cuts_status = add_cuts(settings_, + cuts_to_add, + cut_rhs, + original_lp_, + new_slacks_, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + mutex_original_lp_.unlock(); + f_t add_cuts_time = toc(add_cuts_start_time); + if (add_cuts_time > 1.0) { + settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); + } + if (add_cuts_status != 0) { + settings_.log.printf("Failed to add cuts\n"); + return mip_status_t::NUMERICAL; + } - // Resolve the LP with the new cuts - settings_.log.debug( - "Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", - num_cuts, - cuts_to_add.row_start[cuts_to_add.m], - cut_pool.pool_size(), - cuts_to_add.m + original_lp_.num_rows); - lp_settings.log.log = false; - - f_t add_cuts_start_time = tic(); + if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { + mutex_upper_.lock(); + last_upper_bound = upper_bound_.load(); + std::vector lower_bounds; + std::vector upper_bounds; + find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + mutex_upper_.unlock(); mutex_original_lp_.lock(); - i_t add_cuts_status = add_cuts(settings_, - cuts_to_add, - cut_rhs, - original_lp_, - new_slacks_, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; mutex_original_lp_.unlock(); - f_t add_cuts_time = toc(add_cuts_start_time); - if (add_cuts_time > 1.0) { - settings_.log.debug("Add cuts time %.2f seconds\n", add_cuts_time); - } - if (add_cuts_status != 0) { - settings_.log.printf("Failed to add cuts\n"); - return mip_status_t::NUMERICAL; - } - - if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) { - mutex_upper_.lock(); - last_upper_bound = upper_bound_.load(); - std::vector lower_bounds; - std::vector upper_bounds; - find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - mutex_upper_.unlock(); - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - mutex_original_lp_.unlock(); - } + } - // Try to do bound strengthening - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; + // Try to do bound strengthening + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; #ifdef CHECK_MATRICES - settings_.log.printf("Before A check\n"); - original_lp_.A.check_matrix(); + settings_.log.printf("Before A check\n"); + original_lp_.A.check_matrix(); #endif - original_lp_.A.to_compressed_row(Arow_); - - f_t node_presolve_start_time = tic(); - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - std::vector new_lower = original_lp_.lower; - std::vector new_upper = original_lp_.upper; - bool feasible = - node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); - mutex_original_lp_.lock(); - original_lp_.lower = new_lower; - original_lp_.upper = new_upper; - mutex_original_lp_.unlock(); - f_t node_presolve_time = toc(node_presolve_start_time); - if (node_presolve_time > 1.0) { - settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); - } - if (!feasible) { - settings_.log.printf("Bound strengthening detected infeasibility\n"); - return mip_status_t::INFEASIBLE; - } + original_lp_.A.to_compressed_row(Arow_); + + f_t node_presolve_start_time = tic(); + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + std::vector new_lower = original_lp_.lower; + std::vector new_upper = original_lp_.upper; + bool feasible = + node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + f_t node_presolve_time = toc(node_presolve_start_time); + if (node_presolve_time > 1.0) { + settings_.log.debug("Node presolve time %.2f seconds\n", node_presolve_time); + } + if (!feasible) { + settings_.log.printf("Bound strengthening detected infeasibility\n"); + return mip_status_t::INFEASIBLE; + } - i_t iter = 0; - bool initialize_basis = false; - lp_settings.concurrent_halt = NULL; - f_t dual_phase2_start_time = tic(); - dual::status_t cut_status = dual_phase2_with_advanced_basis(2, - 0, - initialize_basis, - exploration_stats_.start_time, - original_lp_, - lp_settings, - root_vstatus_, - basis_update, - basic_list, - nonbasic_list, - root_relax_soln_, - iter, - edge_norms_); - exploration_stats_.total_lp_iters += iter; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - f_t dual_phase2_time = toc(dual_phase2_start_time); - if (dual_phase2_time > 1.0) { - settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); - } - if (cut_status == dual::status_t::TIME_LIMIT) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return solver_status_; - } + i_t iter = 0; + bool initialize_basis = false; + lp_settings.concurrent_halt = NULL; + f_t dual_phase2_start_time = tic(); + dual::status_t cut_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + exploration_stats_.total_lp_iters += iter; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + f_t dual_phase2_time = toc(dual_phase2_start_time); + if (dual_phase2_time > 1.0) { + settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); + } + if (cut_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return solver_status_; + } - if (cut_status != dual::status_t::OPTIMAL) { - settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); - lp_status_t scratch_status = - solve_linear_program_with_advanced_basis(original_lp_, - exploration_stats_.start_time, - lp_settings, - root_relax_soln_, - basis_update, - basic_list, - nonbasic_list, - root_vstatus_, - edge_norms_); - if (scratch_status == lp_status_t::OPTIMAL) { - // We recovered - cut_status = convert_lp_status_to_dual_status(scratch_status); - exploration_stats_.total_lp_iters += root_relax_soln_.iterations; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); - } else { - settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); - return mip_status_t::NUMERICAL; - } + if (cut_status != dual::status_t::OPTIMAL) { + settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); + lp_status_t scratch_status = + solve_linear_program_with_advanced_basis(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + if (scratch_status == lp_status_t::OPTIMAL) { + // We recovered + cut_status = convert_lp_status_to_dual_status(scratch_status); + exploration_stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else { + settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); + return mip_status_t::NUMERICAL; } + } - f_t remove_cuts_start_time = tic(); - mutex_original_lp_.lock(); - remove_cuts(original_lp_, - settings_, - exploration_stats_.start_time, - Arow_, - new_slacks_, - original_rows, - var_types_, - root_vstatus_, - edge_norms_, - root_relax_soln_.x, - root_relax_soln_.y, - root_relax_soln_.z, - basic_list, - nonbasic_list, - basis_update); - mutex_original_lp_.unlock(); - f_t remove_cuts_time = toc(remove_cuts_start_time); - if (remove_cuts_time > 1.0) { - settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); - } - fractional.clear(); - num_fractional = - fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - - if (num_fractional == 0) { - upper_bound_ = root_objective_; - mutex_upper_.lock(); - incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); - mutex_upper_.unlock(); - } - f_t obj = upper_bound_.load(); - report(' ', obj, root_objective_, 0, num_fractional); - - f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); - f_t abs_gap = upper_bound_.load() - root_objective_; - if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { - set_solution_at_root(solution, cut_info); - set_final_solution(solution, root_objective_); - return mip_status_t::OPTIMAL; - } + f_t remove_cuts_start_time = tic(); + mutex_original_lp_.lock(); + remove_cuts(original_lp_, + settings_, + exploration_stats_.start_time, + Arow_, + new_slacks_, + original_rows, + var_types_, + root_vstatus_, + edge_norms_, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + basis_update); + mutex_original_lp_.unlock(); + f_t remove_cuts_time = toc(remove_cuts_start_time); + if (remove_cuts_time > 1.0) { + settings_.log.debug("Remove cuts time %.2f seconds\n", remove_cuts_time); + } + fractional.clear(); + num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - f_t change_in_objective = root_objective_ - last_objective; - const f_t factor = settings_.cut_change_threshold; - const f_t min_objective = 1e-3; - if (change_in_objective <= - factor * std::max(min_objective, std::abs(root_relax_objective))) { - settings_.log.debug( - "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", - change_in_objective, - root_relax_objective); - break; - } - last_objective = root_objective_; + if (num_fractional == 0) { + upper_bound_ = root_objective_; + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + mutex_upper_.unlock(); } - } + f_t obj = upper_bound_.load(); + report(' ', obj, root_objective_, 0, num_fractional); - print_cut_info(settings_, cut_info); + f_t rel_gap = user_relative_gap(original_lp_, upper_bound_.load(), root_objective_); + f_t abs_gap = upper_bound_.load() - root_objective_; + if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) { + set_solution_at_root(solution, cut_info); + set_final_solution(solution, root_objective_); + return mip_status_t::OPTIMAL; + } - if (cut_info.has_cuts()) { - settings_.log.printf("Cut pool size : %d\n", cut_pool_size); - settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", - original_lp_.num_rows, - original_lp_.num_cols, - original_lp_.A.col_start[original_lp_.A.n]); + f_t change_in_objective = root_objective_ - last_objective; + const f_t factor = settings_.cut_change_threshold; + const f_t min_objective = 1e-3; + if (change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { + settings_.log.debug( + "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", + change_in_objective, + root_relax_objective); + break; + } + last_objective = root_objective_; } + } - set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); - - pc_.resize(original_lp_.num_cols); - { - raft::common::nvtx::range scope_sb("BB::strong_branching"); - strong_branching(original_problem_, - original_lp_, - settings_, - exploration_stats_.start_time, - var_types_, - root_relax_soln_.x, - fractional, - root_objective_, - root_vstatus_, - edge_norms_, - pc_); - } + print_cut_info(settings_, cut_info); - if (toc(exploration_stats_.start_time) > settings_.time_limit) { - solver_status_ = mip_status_t::TIME_LIMIT; - set_final_solution(solution, root_objective_); - return solver_status_; - } + if (cut_info.has_cuts()) { + settings_.log.printf("Cut pool size : %d\n", cut_pool_size); + settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", + original_lp_.num_rows, + original_lp_.num_cols, + original_lp_.A.col_start[original_lp_.A.n]); + } - if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) { - std::vector lower_bounds; - std::vector upper_bounds; - i_t num_fixed = find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); - if (num_fixed > 0) { - std::vector bounds_changed(original_lp_.num_cols, true); - std::vector row_sense; + set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); - bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + pc_.resize(original_lp_.num_cols); + { + raft::common::nvtx::range scope_sb("BB::strong_branching"); + strong_branching(original_problem_, + original_lp_, + settings_, + exploration_stats_.start_time, + var_types_, + root_relax_soln_.x, + fractional, + root_objective_, + root_vstatus_, + edge_norms_, + pc_); + } - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - bool feasible = node_presolve.bounds_strengthening( - settings_, bounds_changed, original_lp_.lower, original_lp_.upper); - mutex_original_lp_.unlock(); - if (!feasible) { - settings_.log.printf("Bound strengthening failed\n"); - return mip_status_t::NUMERICAL; // We had a feasible integer solution, but bound - // strengthening thinks we are infeasible. + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return solver_status_; + } + + if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) { + std::vector lower_bounds; + std::vector upper_bounds; + i_t num_fixed = find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds); + if (num_fixed > 0) { + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; + + bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); + + mutex_original_lp_.lock(); + original_lp_.lower = lower_bounds; + original_lp_.upper = upper_bounds; + bool feasible = node_presolve.bounds_strengthening( + settings_, bounds_changed, original_lp_.lower, original_lp_.upper); + mutex_original_lp_.unlock(); + if (!feasible) { + settings_.log.printf("Bound strengthening failed\n"); + return mip_status_t::NUMERICAL; // We had a feasible integer solution, but bound + // strengthening thinks we are infeasible. + } + // Go through and check the fractional variables and remove any that are now fixed to their + // bounds + std::vector to_remove(fractional.size(), 0); + i_t num_to_remove = 0; + for (i_t k = 0; k < fractional.size(); k++) { + const i_t j = fractional[k]; + if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) { + to_remove[k] = 1; + num_to_remove++; } - // Go through and check the fractional variables and remove any that are now fixed to their - // bounds - std::vector to_remove(fractional.size(), 0); - i_t num_to_remove = 0; + } + if (num_to_remove > 0) { + std::vector new_fractional; + new_fractional.reserve(fractional.size() - num_to_remove); for (i_t k = 0; k < fractional.size(); k++) { - const i_t j = fractional[k]; - if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) { - to_remove[k] = 1; - num_to_remove++; - } - } - if (num_to_remove > 0) { - std::vector new_fractional; - new_fractional.reserve(fractional.size() - num_to_remove); - for (i_t k = 0; k < fractional.size(); k++) { - if (!to_remove[k]) { new_fractional.push_back(fractional[k]); } - } - fractional = new_fractional; - num_fractional = fractional.size(); + if (!to_remove[k]) { new_fractional.push_back(fractional[k]); } } + fractional = new_fractional; + num_fractional = fractional.size(); } } + } - // Choose variable to branch on - i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x, log); + // Choose variable to branch on + i_t branch_var = pc_.variable_selection(fractional, root_relax_soln_.x, log); - search_tree_.root = std::move(mip_node_t(root_objective_, root_vstatus_)); - search_tree_.num_nodes = 0; - search_tree_.graphviz_node(settings_.log, &search_tree_.root, "lower bound", root_objective_); - search_tree_.branch(&search_tree_.root, - branch_var, - root_relax_soln_.x[branch_var], - num_fractional, - root_vstatus_, - original_lp_, - log); - node_queue_.push(search_tree_.root.get_down_child()); - node_queue_.push(search_tree_.root.get_up_child()); + search_tree_.root = std::move(mip_node_t(root_objective_, root_vstatus_)); + search_tree_.num_nodes = 0; + search_tree_.graphviz_node(settings_.log, &search_tree_.root, "lower bound", root_objective_); + search_tree_.branch(&search_tree_.root, + branch_var, + root_relax_soln_.x[branch_var], + num_fractional, + root_vstatus_, + original_lp_, + log); + node_queue_.push(search_tree_.root.get_down_child()); + node_queue_.push(search_tree_.root.get_up_child()); - settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads); + settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads); - exploration_stats_.nodes_explored = 0; - exploration_stats_.nodes_unexplored = 2; - exploration_stats_.nodes_since_last_log = 0; - exploration_stats_.last_log = tic(); - min_node_queue_size_ = 2 * settings_.num_threads; + exploration_stats_.nodes_explored = 0; + exploration_stats_.nodes_unexplored = 2; + exploration_stats_.nodes_since_last_log = 0; + exploration_stats_.last_log = tic(); + min_node_queue_size_ = 2 * settings_.num_threads; - if (settings_.diving_settings.coefficient_diving != 0) { - calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); - } - if (settings_.deterministic) { - settings_.log.printf( - " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " - "| Gap | Work | Time |\n"); - } else { - settings_.log.printf( - " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " - "| Gap | Time |\n"); - } + if (settings_.diving_settings.coefficient_diving != 0) { + calculate_variable_locks(original_lp_, var_up_locks_, var_down_locks_); + } + if (settings_.deterministic) { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " + "| Gap | Work | Time |\n"); + } else { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node " + "| Gap | Time |\n"); } if (settings_.deterministic) { From 2be4f5b35c0424e0f4c53633b568816996317947 Mon Sep 17 00:00:00 2001 From: anandhkb Date: Wed, 11 Mar 2026 17:32:44 -0700 Subject: [PATCH 8/9] Log root LP dual objective in final status for gap verification Add 'Root LP dual objective (last)' when finite so reviewers can confirm the heuristic gap calculation used the correct dual simplex objective. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index f9d9112763..23bb5e1e6a 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -724,6 +724,12 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& obj, is_maximization ? "Upper" : "Lower", user_bound); + { + const f_t root_lp_obj = root_lp_current_lower_bound_.load(); + if (std::isfinite(root_lp_obj)) { + settings_.log.printf("Root LP dual objective (last): %.16e\n", root_lp_obj); + } + } if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) { solver_status_ = mip_status_t::OPTIMAL; From 6a0d60beb58de6034167190b7d255390a80a909e Mon Sep 17 00:00:00 2001 From: anandhkb Date: Wed, 11 Mar 2026 17:34:02 -0700 Subject: [PATCH 9/9] Style: fix formatting in solve() --- cpp/src/branch_and_bound/branch_and_bound.cpp | 11 +++++------ 1 file changed, 5 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 23bb5e1e6a..d8fab6ec52 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1989,7 +1989,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); - if (settings_.clique_cuts != 0 && clique_table_ == nullptr) { signal_extend_cliques_.store(false, std::memory_order_release); typename ::cuopt::linear_programming::mip_solver_settings_t::tolerances_t @@ -2015,11 +2014,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut }); } - i_t original_rows = original_lp_.num_rows; - simplex_solver_settings_t lp_settings = settings_; - lp_settings.inside_mip = 1; - lp_settings.scale_columns = false; - lp_settings.concurrent_halt = get_root_concurrent_halt(); + i_t original_rows = original_lp_.num_rows; + simplex_solver_settings_t lp_settings = settings_; + lp_settings.inside_mip = 1; + lp_settings.scale_columns = false; + lp_settings.concurrent_halt = get_root_concurrent_halt(); lp_settings.dual_simplex_objective_callback = [this](f_t user_obj) { root_lp_current_lower_bound_.store(user_obj); };