From 09a5e7e74091e33da22889a86d953635bcf23326 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 6 Oct 2025 15:31:50 +0200 Subject: [PATCH 01/51] child nodes now reuse the bounds of the parent. fixed incorrect bounds when detaching the node from the tree. add backtracking argument to diving. --- cpp/src/dual_simplex/branch_and_bound.cpp | 88 +++++++++++++++-------- cpp/src/dual_simplex/branch_and_bound.hpp | 55 ++++++++++---- cpp/src/dual_simplex/mip_node.hpp | 31 ++++---- 3 files changed, 119 insertions(+), 55 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index cf6fd69798..027a27e95a 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -556,9 +557,10 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& mip_node_t* node_ptr, lp_problem_t& leaf_problem, const csc_matrix_t& Arow, + const std::vector& bounds_changed, + char thread_type, f_t upper_bound, - logger_t& log, - char thread_type) + logger_t& log) { f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; @@ -566,16 +568,6 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); assert(leaf_vstatus.size() == leaf_problem.num_cols); - // Set the correct bounds for the leaf problem - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; - - std::vector bounds_changed(leaf_problem.num_cols, false); - // Technically, we can get the already strengthened bounds from the node/parent instead of - // getting it from the original problem and re-strengthening. But this requires storing - // two vectors at each node and potentially cause memory issues - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - i_t node_iter = 0; f_t lp_start_time = tic(); std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; @@ -744,8 +736,16 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* status_ = mip_exploration_status_t::TIME_LIMIT; return; } - node_status_t node_status = - solve_node(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + + std::vector bounds_changed(leaf_problem.num_cols, false); + + // Set the correct bounds for the leaf problem + leaf_problem.lower = original_lp_.lower; + leaf_problem.upper = original_lp_.upper; + node->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); + + node_status_t node_status = solve_node( + *search_tree, node, leaf_problem, Arow, bounds_changed, 'B', upper_bound, settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -780,9 +780,12 @@ void branch_and_bound_t::explore_subtree(i_t id, lp_problem_t& leaf_problem, const csc_matrix_t& Arow) { + std::vector bounds_changed(leaf_problem.num_cols); std::deque*> stack; stack.push_front(start_node); + bool recompute_bounds = true; + while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { if (omp_get_thread_num() == 0) { repair_heuristic_solutions(); } @@ -842,8 +845,20 @@ void branch_and_bound_t::explore_subtree(i_t id, return; } - node_status_t node_status = - solve_node(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + // Reset the bounds for the leaf problem to the original problem + if (recompute_bounds) { + leaf_problem.lower = original_lp_.lower; + leaf_problem.upper = original_lp_.upper; + std::fill(bounds_changed.begin(), bounds_changed.end(), false); + node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); + + } else { + node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); + } + + node_status_t node_status = solve_node( + search_tree, node_ptr, leaf_problem, Arow, bounds_changed, 'B', upper_bound, settings_.log); + recompute_bounds = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -866,7 +881,7 @@ void branch_and_bound_t::explore_subtree(i_t id, // would be better if we discard the node instead. if (get_heap_size() > settings_.num_bfs_threads) { mutex_dive_queue_.lock(); - dive_queue_.push(node->detach_copy()); + dive_queue_.emplace(node->detach_copy(), leaf_problem.lower, leaf_problem.upper); mutex_dive_queue_.unlock(); } @@ -943,26 +958,31 @@ void branch_and_bound_t::best_first_thread(i_t id, template void branch_and_bound_t::diving_thread(lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) + const csc_matrix_t& Arow, + i_t backtracking) { logger_t log; log.log = false; + std::vector bounds_changed(leaf_problem.num_cols); + while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { - std::optional> start_node; + std::optional> start_node; mutex_dive_queue_.lock(); if (dive_queue_.size() > 0) { start_node = dive_queue_.pop(); } mutex_dive_queue_.unlock(); if (start_node.has_value()) { - if (get_upper_bound() < start_node->lower_bound) { continue; } + if (get_upper_bound() < start_node->node.lower_bound) { continue; } - search_tree_t subtree(std::move(start_node.value())); + search_tree_t subtree(std::move(start_node->node)); std::deque*> stack; stack.push_front(&subtree.root); + bool recompute_bounds = true; + while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); @@ -975,8 +995,19 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (toc(stats_.start_time) > settings_.time_limit) { return; } + // Reset the bounds for the leaf problem based on the starting node + if (recompute_bounds) { + leaf_problem.lower = start_node->lp_lower; + leaf_problem.upper = start_node->lp_upper; + std::fill(bounds_changed.begin(), bounds_changed.end(), false); + node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); + } else { + node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); + } + node_status_t node_status = - solve_node(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); + solve_node(subtree, node_ptr, leaf_problem, Arow, bounds_changed, 'D', upper_bound, log); + recompute_bounds = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { return; @@ -985,16 +1016,15 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr auto [first, second] = child_selection(node_ptr); stack.push_front(second); stack.push_front(first); + } - // If the diving thread is consuming the nodes faster than the - // best first search, then we split the current subtree at the - // lowest possible point and move to the queue, so it can - // be picked by another thread. - if (dive_queue_.size() < min_diving_queue_size_) { + if (stack.size() > 1) { + if (dive_queue_.size() < min_diving_queue_size_ || + (stack.front()->depth - stack.back()->depth) > backtracking) { mutex_dive_queue_.lock(); mip_node_t* new_node = stack.back(); stack.pop_back(); - dive_queue_.push(new_node->detach_copy()); + dive_queue_.emplace(new_node->detach_copy(), leaf_problem.lower, leaf_problem.upper); mutex_dive_queue_.unlock(); } } @@ -1192,7 +1222,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(leaf_problem, Arow); + diving_thread(leaf_problem, Arow, 10); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 7b80f88fa9..b97a4206ba 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -55,36 +55,64 @@ enum class mip_exploration_status_t { template void upper_bound_callback(f_t upper_bound); +template +struct diving_root_t { + mip_node_t node; + std::vector lp_lower; + std::vector lp_upper; + + diving_root_t(mip_node_t&& node, + const std::vector& lower, + const std::vector& upper) + : node(std::move(node)), lp_upper(upper), lp_lower(lower) + { + } + + friend bool operator>(const diving_root_t& a, const diving_root_t& b) + { + return a.node.lower_bound > b.node.lower_bound; + } +}; + // A min-heap for storing the starting nodes for the dives. -// This has a maximum size of 8192, such that the container +// This has a maximum size of 256, such that the container // will discard the least promising node if the queue is full. template class dive_queue_t { private: - std::vector> buffer; - static constexpr i_t max_size_ = 2048; + std::vector> buffer; + static constexpr i_t max_size_ = 256; public: dive_queue_t() { buffer.reserve(max_size_); } - void push(mip_node_t&& node) + void push(diving_root_t&& node) { buffer.push_back(std::move(node)); - std::push_heap(buffer.begin(), buffer.end(), node_compare_t()); + std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); + if (buffer.size() > max_size()) { buffer.pop_back(); } + } + + void emplace(mip_node_t&& node, + const std::vector& lower, + const std::vector& upper) + { + buffer.emplace_back(std::move(node), lower, upper); + std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); if (buffer.size() > max_size()) { buffer.pop_back(); } } - mip_node_t pop() + diving_root_t pop() { - std::pop_heap(buffer.begin(), buffer.end(), node_compare_t()); - mip_node_t node = std::move(buffer.back()); + std::pop_heap(buffer.begin(), buffer.end(), std::greater<>()); + diving_root_t node = std::move(buffer.back()); buffer.pop_back(); return node; } i_t size() const { return buffer.size(); } constexpr i_t max_size() const { return max_size_; } - const mip_node_t& top() const { return buffer.front(); } + const diving_root_t& top() const { return buffer.front(); } void clear() { buffer.clear(); } }; @@ -222,16 +250,19 @@ class branch_and_bound_t { // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(lp_problem_t& leaf_problem, const csc_matrix_t& Arow); + void diving_thread(lp_problem_t& leaf_problem, + const csc_matrix_t& Arow, + i_t backtracking); // Solve the LP relaxation of a leaf node and update the tree. node_status_t solve_node(search_tree_t& search_tree, mip_node_t* node_ptr, lp_problem_t& leaf_problem, const csc_matrix_t& Arow, + const std::vector& bounds_changed, + char thread_type, f_t upper_bound, - logger_t& log, - char thread_type); + logger_t& log); // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection( diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index a3a34eb81d..23ea214eb2 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -87,26 +87,29 @@ class mip_node_t { std::vector& upper, std::vector& bounds_changed) const { - std::fill(bounds_changed.begin(), bounds_changed.end(), false); - // Apply the bounds at the current node - assert(lower.size() > branch_var); - assert(upper.size() > branch_var); - lower[branch_var] = branch_var_lower; - upper[branch_var] = branch_var_upper; - bounds_changed[branch_var] = true; + update_variable_bound(lower, upper, bounds_changed); + mip_node_t* parent_ptr = parent; while (parent_ptr != nullptr) { if (parent_ptr->node_id == 0) { break; } - assert(parent_ptr->branch_var >= 0); - assert(lower.size() > parent_ptr->branch_var); - assert(upper.size() > parent_ptr->branch_var); - lower[parent_ptr->branch_var] = parent_ptr->branch_var_lower; - upper[parent_ptr->branch_var] = parent_ptr->branch_var_upper; - bounds_changed[parent_ptr->branch_var] = true; - parent_ptr = parent_ptr->parent; + parent_ptr->update_variable_bound(lower, upper, bounds_changed); + parent_ptr = parent_ptr->parent; } } + void update_variable_bound(std::vector& lower, + std::vector& upper, + std::vector& bounds_changed) const + { + // Apply the bounds at the current node + assert(branch_var >= 0); + assert(lower.size() > branch_var); + assert(upper.size() > branch_var); + lower[branch_var] = branch_var_lower; + upper[branch_var] = branch_var_upper; + bounds_changed[branch_var] = true; + } + mip_node_t* get_down_child() const { return children[0].get(); } mip_node_t* get_up_child() const { return children[1].get(); } From 0a7a7ebecef44ca54301f6090e5753b6ea1fce2a Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 9 Oct 2025 12:17:44 +0200 Subject: [PATCH 02/51] refactor based on the reuse-parent-branch --- cpp/src/dual_simplex/branch_and_bound.cpp | 132 ++++++++++++++-------- cpp/src/dual_simplex/branch_and_bound.hpp | 27 +++-- cpp/src/dual_simplex/mip_node.hpp | 2 +- 3 files changed, 100 insertions(+), 61 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 027a27e95a..f4122fa1bc 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -553,28 +553,20 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) } template -node_status_t branch_and_bound_t::solve_node(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, - const std::vector& bounds_changed, - char thread_type, - f_t upper_bound, - logger_t& log) +dual::status_t branch_and_bound_t::solve_node_lp(mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + lp_solution_t& leaf_solution, + const csc_matrix_t& Arow, + const std::vector& bounds_changed, + logger_t& log) { - f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; - + std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; std::vector& leaf_vstatus = node_ptr->vstatus; - lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); assert(leaf_vstatus.size() == leaf_problem.num_cols); - i_t node_iter = 0; - f_t lp_start_time = tic(); - std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; - simplex_solver_settings_t lp_settings = settings_; lp_settings.set_log(false); - lp_settings.cut_off = upper_bound + settings_.dual_tol; + lp_settings.cut_off = get_upper_bound() + settings_.dual_tol; lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); @@ -586,7 +578,9 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; if (feasible) { - lp_status = dual_phase2(2, + i_t node_iter = 0; + f_t lp_start_time = tic(); + lp_status = dual_phase2(2, 0, lp_start_time, leaf_problem, @@ -602,16 +596,32 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); lp_status = convert_lp_status_to_dual_status(second_status); } + + stats_.total_lp_solve_time += toc(lp_start_time); + stats_.total_lp_iters += node_iter; } - stats_.total_lp_solve_time += toc(lp_start_time); - stats_.total_lp_iters += node_iter; + return lp_status; +} + +template +node_status_t branch_and_bound_t::update_tree( + search_tree_t& search_tree, + mip_node_t* node_ptr, + const lp_problem_t& leaf_problem, + const lp_solution_t& leaf_solution, + dual::status_t lp_status, + char thread_type, + logger_t& log) +{ + f_t upper_bound = get_upper_bound(); + f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; if (lp_status == dual::status_t::DUAL_UNBOUNDED) { // Node was infeasible. Do not branch node_ptr->lower_bound = inf; search_tree.graphviz_node(log, node_ptr, "infeasible", 0.0); - search_tree.update_tree(node_ptr, node_status_t::INFEASIBLE); + search_tree.update(node_ptr, node_status_t::INFEASIBLE); return node_status_t::INFEASIBLE; } else if (lp_status == dual::status_t::CUTOFF) { @@ -619,7 +629,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& node_ptr->lower_bound = upper_bound; f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); search_tree.graphviz_node(log, node_ptr, "cut off", leaf_objective); - search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.update(node_ptr, node_status_t::FATHOMED); return node_status_t::FATHOMED; } else if (lp_status == dual::status_t::OPTIMAL) { @@ -637,28 +647,30 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& // Found a integer feasible solution add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, thread_type); search_tree.graphviz_node(log, node_ptr, "integer feasible", leaf_objective); - search_tree.update_tree(node_ptr, node_status_t::INTEGER_FEASIBLE); + search_tree.update(node_ptr, node_status_t::INTEGER_FEASIBLE); return node_status_t::INTEGER_FEASIBLE; } else if (leaf_objective <= upper_bound + abs_fathom_tol) { + logger_t pc_log = log; + pc_log.log = false; + // Choose fractional variable to branch on - const i_t branch_var = - pc_.variable_selection(leaf_fractional, leaf_solution.x, lp_settings.log); + const i_t branch_var = pc_.variable_selection(leaf_fractional, leaf_solution.x, pc_log); - assert(leaf_vstatus.size() == leaf_problem.num_cols); + assert(node_ptr->vstatus.size() == leaf_problem.num_cols); search_tree.branch( - node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_, log); + node_ptr, branch_var, leaf_solution.x[branch_var], node_ptr->vstatus, original_lp_, log); node_ptr->status = node_status_t::HAS_CHILDREN; return node_status_t::HAS_CHILDREN; } else { search_tree.graphviz_node(log, node_ptr, "fathomed", leaf_objective); - search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.update(node_ptr, node_status_t::FATHOMED); return node_status_t::FATHOMED; } } else if (lp_status == dual::status_t::TIME_LIMIT) { search_tree.graphviz_node(log, node_ptr, "timeout", 0.0); - search_tree.update_tree(node_ptr, node_status_t::TIME_LIMIT); + search_tree.update(node_ptr, node_status_t::TIME_LIMIT); return node_status_t::TIME_LIMIT; } else { @@ -674,7 +686,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& } search_tree.graphviz_node(log, node_ptr, "numerical", 0.0); - search_tree.update_tree(node_ptr, node_status_t::NUMERICAL); + search_tree.update(node_ptr, node_status_t::NUMERICAL); return node_status_t::NUMERICAL; } } @@ -702,7 +714,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree->graphviz_node(settings_.log, node, "cutoff", node->lower_bound); - search_tree->update_tree(node, node_status_t::FATHOMED); + search_tree->update(node, node_status_t::FATHOMED); return; } @@ -728,6 +740,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* gap_user.c_str(), now); + stats_.last_log = tic(); stats_.nodes_since_last_log = 0; } } @@ -737,15 +750,21 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* return; } - std::vector bounds_changed(leaf_problem.num_cols, false); + const i_t m = leaf_problem.num_rows; + const i_t n = leaf_problem.num_cols; + std::vector bounds_changed(n, false); + lp_solution_t leaf_solution(m, n); // Set the correct bounds for the leaf problem leaf_problem.lower = original_lp_.lower; leaf_problem.upper = original_lp_.upper; node->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - node_status_t node_status = solve_node( - *search_tree, node, leaf_problem, Arow, bounds_changed, 'B', upper_bound, settings_.log); + dual::status_t lp_status = + solve_node_lp(node, leaf_problem, leaf_solution, Arow, bounds_changed, settings_.log); + + node_status_t node_status = + update_tree(*search_tree, node, leaf_problem, leaf_solution, lp_status, 'B', settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -780,12 +799,15 @@ void branch_and_bound_t::explore_subtree(i_t id, lp_problem_t& leaf_problem, const csc_matrix_t& Arow) { - std::vector bounds_changed(leaf_problem.num_cols); + const i_t m = leaf_problem.num_rows; + const i_t n = leaf_problem.num_cols; + bool recompute = true; + std::vector bounds_changed(n); + lp_solution_t leaf_solution(m, n); + std::deque*> stack; stack.push_front(start_node); - bool recompute_bounds = true; - while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { if (omp_get_thread_num() == 0) { repair_heuristic_solutions(); } @@ -811,7 +833,8 @@ void branch_and_bound_t::explore_subtree(i_t id, if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); - search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.update(node_ptr, node_status_t::FATHOMED); + recompute = true; continue; } @@ -845,8 +868,8 @@ void branch_and_bound_t::explore_subtree(i_t id, return; } - // Reset the bounds for the leaf problem to the original problem - if (recompute_bounds) { + // Recompute the bounds + if (recompute) { leaf_problem.lower = original_lp_.lower; leaf_problem.upper = original_lp_.upper; std::fill(bounds_changed.begin(), bounds_changed.end(), false); @@ -856,9 +879,11 @@ void branch_and_bound_t::explore_subtree(i_t id, node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); } - node_status_t node_status = solve_node( - search_tree, node_ptr, leaf_problem, Arow, bounds_changed, 'B', upper_bound, settings_.log); - recompute_bounds = node_status != node_status_t::HAS_CHILDREN; + dual::status_t lp_status = + solve_node_lp(node_ptr, leaf_problem, leaf_solution, Arow, bounds_changed, settings_.log); + node_status_t node_status = update_tree( + search_tree, node_ptr, leaf_problem, leaf_solution, lp_status, 'B', settings_.log); + recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -929,7 +954,7 @@ void branch_and_bound_t::best_first_thread(i_t id, // This node was put on the heap earlier but its lower bound is now greater than the // current upper bound search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); - search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.update(node_ptr, node_status_t::FATHOMED); active_subtrees_--; continue; } @@ -964,7 +989,8 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr logger_t log; log.log = false; - std::vector bounds_changed(leaf_problem.num_cols); + const i_t m = leaf_problem.num_rows; + const i_t n = leaf_problem.num_cols; while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { @@ -977,12 +1003,13 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (start_node.has_value()) { if (get_upper_bound() < start_node->node.lower_bound) { continue; } + bool recompute = true; + std::vector bounds_changed(n); + lp_solution_t leaf_solution(m, n); search_tree_t subtree(std::move(start_node->node)); std::deque*> stack; stack.push_front(&subtree.root); - bool recompute_bounds = true; - while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); @@ -990,13 +1017,14 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr f_t rel_gap = user_relative_gap(original_lp_, upper_bound, node_ptr->lower_bound); if (node_ptr->lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { + recompute = true; continue; } if (toc(stats_.start_time) > settings_.time_limit) { return; } - // Reset the bounds for the leaf problem based on the starting node - if (recompute_bounds) { + // Recompute the bounds + if (recompute) { leaf_problem.lower = start_node->lp_lower; leaf_problem.upper = start_node->lp_upper; std::fill(bounds_changed.begin(), bounds_changed.end(), false); @@ -1005,9 +1033,12 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); } + dual::status_t lp_status = + solve_node_lp(node_ptr, leaf_problem, leaf_solution, Arow, bounds_changed, log); + node_status_t node_status = - solve_node(subtree, node_ptr, leaf_problem, Arow, bounds_changed, 'D', upper_bound, log); - recompute_bounds = node_status != node_status_t::HAS_CHILDREN; + update_tree(subtree, node_ptr, leaf_problem, leaf_solution, lp_status, 'D', log); + recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { return; @@ -1068,6 +1099,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_, stats_.start_time, lp_settings, root_relax_soln_, root_vstatus_, edge_norms_); stats_.total_lp_iters = root_relax_soln_.iterations; stats_.total_lp_solve_time = toc(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. diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index b97a4206ba..f27d9476ff 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -216,7 +216,7 @@ class branch_and_bound_t { // Set the final solution. mip_status_t set_final_solution(mip_solution_t& solution, f_t lower_bound); - // Update the incumbent solution with the new feasible solution. + // Update the incumbent solution with the new feasible solution // found during branch and bound. void add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, @@ -254,15 +254,22 @@ class branch_and_bound_t { const csc_matrix_t& Arow, i_t backtracking); - // Solve the LP relaxation of a leaf node and update the tree. - node_status_t solve_node(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, - const std::vector& bounds_changed, - char thread_type, - f_t upper_bound, - logger_t& log); + // Solve the LP relaxation of a leaf node. + dual::status_t solve_node_lp(mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + lp_solution_t& leaf_solution, + const csc_matrix_t& Arow, + const std::vector& bounds_changed, + logger_t& log); + + // Update the tree after solving the LP relaxation. + node_status_t update_tree(search_tree_t& search_tree, + mip_node_t* node_ptr, + const lp_problem_t& leaf_problem, + const lp_solution_t& leaf_solution, + dual::status_t lp_status, + char thread_type, + logger_t& log); // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection( diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 23ea214eb2..0eb9142acf 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -273,7 +273,7 @@ class search_tree_t { search_tree_t(mip_node_t&& node) : root(std::move(node)), num_nodes(0) {} - void update_tree(mip_node_t* node_ptr, node_status_t status) + void update(mip_node_t* node_ptr, node_status_t status) { mutex.lock(); std::vector*> stack; From 7997b12952b21a9b62c259db8eec5b5d4ec53188 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 9 Oct 2025 13:07:45 +0200 Subject: [PATCH 03/51] fixed value of the bounds_change for the child --- cpp/src/dual_simplex/branch_and_bound.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index f4122fa1bc..75c2aba29d 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -868,11 +868,13 @@ void branch_and_bound_t::explore_subtree(i_t id, return; } + // Reset the bound_changed markers + std::fill(bounds_changed.begin(), bounds_changed.end(), false); + // Recompute the bounds if (recompute) { leaf_problem.lower = original_lp_.lower; leaf_problem.upper = original_lp_.upper; - std::fill(bounds_changed.begin(), bounds_changed.end(), false); node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); } else { @@ -1023,11 +1025,13 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (toc(stats_.start_time) > settings_.time_limit) { return; } + // Reset the bound_changed markers + std::fill(bounds_changed.begin(), bounds_changed.end(), false); + // Recompute the bounds if (recompute) { leaf_problem.lower = start_node->lp_lower; leaf_problem.upper = start_node->lp_upper; - std::fill(bounds_changed.begin(), bounds_changed.end(), false); node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); } else { node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); From 95de283fe790de3789516cc5832be541a6502351 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 13 Oct 2025 13:07:34 +0200 Subject: [PATCH 04/51] revert refactor --- cpp/src/dual_simplex/branch_and_bound.cpp | 61 +++++++---------------- cpp/src/dual_simplex/branch_and_bound.hpp | 24 +++------ 2 files changed, 26 insertions(+), 59 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 7f70503b6f..edfd75fa87 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -17,7 +17,6 @@ #include #include -#include #include #include #include @@ -553,19 +552,24 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) } template -dual::status_t branch_and_bound_t::solve_node_lp(mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - lp_solution_t& leaf_solution, - const csc_matrix_t& Arow, - const std::vector& bounds_changed, - logger_t& log) +node_status_t branch_and_bound_t::solve_node(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + const std::vector& bounds_changed, + const csc_matrix_t& Arow, + char thread_type, + logger_t& log) { + const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + const f_t upper_bound = get_upper_bound(); + + lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); std::vector& leaf_vstatus = node_ptr->vstatus; assert(leaf_vstatus.size() == leaf_problem.num_cols); simplex_solver_settings_t lp_settings = settings_; lp_settings.set_log(false); - lp_settings.cut_off = get_upper_bound() + settings_.dual_tol; + lp_settings.cut_off = upper_bound + settings_.dual_tol; lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); @@ -602,22 +606,6 @@ dual::status_t branch_and_bound_t::solve_node_lp(mip_node_t* stats_.total_lp_iters += node_iter; } - return lp_status; -} - -template -node_status_t branch_and_bound_t::update_tree( - search_tree_t& search_tree, - mip_node_t* node_ptr, - const lp_problem_t& leaf_problem, - const lp_solution_t& leaf_solution, - dual::status_t lp_status, - char thread_type, - logger_t& log) -{ - f_t upper_bound = get_upper_bound(); - f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; - if (lp_status == dual::status_t::DUAL_UNBOUNDED) { // Node was infeasible. Do not branch node_ptr->lower_bound = inf; @@ -658,9 +646,9 @@ node_status_t branch_and_bound_t::update_tree( // Choose fractional variable to branch on const i_t branch_var = pc_.variable_selection(leaf_fractional, leaf_solution.x, pc_log); - assert(node_ptr->vstatus.size() == leaf_problem.num_cols); + assert(leaf_vstatus.size() == leaf_problem.num_cols); search_tree.branch( - node_ptr, branch_var, leaf_solution.x[branch_var], node_ptr->vstatus, original_lp_, log); + node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_, log); node_ptr->status = node_status_t::HAS_CHILDREN; return node_status_t::HAS_CHILDREN; @@ -753,21 +741,16 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* return; } - const i_t m = leaf_problem.num_rows; const i_t n = leaf_problem.num_cols; std::vector bounds_changed(n, false); - lp_solution_t leaf_solution(m, n); // Set the correct bounds for the leaf problem leaf_problem.lower = original_lp_.lower; leaf_problem.upper = original_lp_.upper; node->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - dual::status_t lp_status = - solve_node_lp(node, leaf_problem, leaf_solution, Arow, bounds_changed, settings_.log); - node_status_t node_status = - update_tree(*search_tree, node, leaf_problem, leaf_solution, lp_status, 'B', settings_.log); + solve_node(*search_tree, node, leaf_problem, bounds_changed, Arow, 'B', settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -884,10 +867,8 @@ void branch_and_bound_t::explore_subtree(i_t task_id, node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); } - dual::status_t lp_status = - solve_node_lp(node_ptr, leaf_problem, leaf_solution, Arow, bounds_changed, settings_.log); - node_status_t node_status = update_tree( - search_tree, node_ptr, leaf_problem, leaf_solution, lp_status, 'B', settings_.log); + node_status_t node_status = + solve_node(search_tree, node_ptr, leaf_problem, bounds_changed, Arow, 'B', settings_.log); recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -993,7 +974,6 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr logger_t log; log.log = false; - const i_t m = leaf_problem.num_rows; const i_t n = leaf_problem.num_cols; while (status_ == mip_exploration_status_t::RUNNING && @@ -1009,7 +989,6 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr bool recompute = true; std::vector bounds_changed(n); - lp_solution_t leaf_solution(m, n); search_tree_t subtree(std::move(start_node->node)); std::deque*> stack; stack.push_front(&subtree.root); @@ -1039,12 +1018,8 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); } - dual::status_t lp_status = - solve_node_lp(node_ptr, leaf_problem, leaf_solution, Arow, bounds_changed, log); - node_status_t node_status = - update_tree(subtree, node_ptr, leaf_problem, leaf_solution, lp_status, 'D', log); - recompute = node_status != node_status_t::HAS_CHILDREN; + solve_node(subtree, node_ptr, leaf_problem, bounds_changed, Arow, 'D', log); if (node_status == node_status_t::TIME_LIMIT) { return; diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index bea4001448..4ff7cfbde9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -252,22 +252,14 @@ class branch_and_bound_t { // a deep dive into the subtree determined by the node. void diving_thread(lp_problem_t& leaf_problem, const csc_matrix_t& Arow); - // Solve the LP relaxation of a leaf node. - dual::status_t solve_node_lp(mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - lp_solution_t& leaf_solution, - const csc_matrix_t& Arow, - const std::vector& bounds_changed, - logger_t& log); - - // Update the tree after solving the LP relaxation. - node_status_t update_tree(search_tree_t& search_tree, - mip_node_t* node_ptr, - const lp_problem_t& leaf_problem, - const lp_solution_t& leaf_solution, - dual::status_t lp_status, - char thread_type, - logger_t& log); + // Solve the LP relaxation of a leaf node and update the tree. + node_status_t solve_node(search_tree_t& search_tree, + mip_node_t* node_ptr, + lp_problem_t& leaf_problem, + const std::vector& bounds_changed, + const csc_matrix_t& Arow, + char thread_type, + logger_t& log); // Sort the children based on the Martin's criteria. std::pair*, mip_node_t*> child_selection( From f44b134ad1e1bbbbe9e24122b266254919e5f6d4 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 13 Oct 2025 14:35:48 +0200 Subject: [PATCH 05/51] fix incorrect bounds when two nodes shares the same branch variable --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- cpp/src/dual_simplex/mip_node.hpp | 17 +++++++++++------ 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index edfd75fa87..0fd03a10a6 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -648,7 +648,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& assert(leaf_vstatus.size() == leaf_problem.num_cols); search_tree.branch( - node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_, log); + node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, leaf_problem, log); node_ptr->status = node_status_t::HAS_CHILDREN; return node_status_t::HAS_CHILDREN; diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 0eb9142acf..3f2a3e4325 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -90,21 +90,27 @@ class mip_node_t { update_variable_bound(lower, upper, bounds_changed); mip_node_t* parent_ptr = parent; - while (parent_ptr != nullptr) { - if (parent_ptr->node_id == 0) { break; } + while (parent_ptr != nullptr && parent_ptr->node_id != 0) { parent_ptr->update_variable_bound(lower, upper, bounds_changed); parent_ptr = parent_ptr->parent; } } + // Here we assume that we are traversing from the deepest node to the + // root of the tree void update_variable_bound(std::vector& lower, std::vector& upper, std::vector& bounds_changed) const { - // Apply the bounds at the current node assert(branch_var >= 0); assert(lower.size() > branch_var); assert(upper.size() > branch_var); + + // If the bounds have already been updated on another node, + // skip this node as it contains a less tight bounds. + if (bounds_changed[branch_var]) { return; } + + // Apply the bounds at the current node lower[branch_var] = branch_var_lower; upper[branch_var] = branch_var_upper; bounds_changed[branch_var] = true; @@ -206,9 +212,8 @@ class mip_node_t { } // This method creates a copy of the current node - // with its parent set to `nullptr`, `node_id = 0` - // and `depth = 0` such that it is the root - // of a separated tree. + // with its parent set to `nullptr` and `depth = 0`. + // This detaches the node from the tree. mip_node_t detach_copy() const { mip_node_t copy(lower_bound, vstatus); From bfeb6608a43bdbf9d5b9f6ba7b26cecc4e2e029a Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 14 Oct 2025 18:55:36 +0200 Subject: [PATCH 06/51] constraints_changed is only set for branched variables. refactor bound_strenghtening to persist some variables between calls. --- cpp/src/dual_simplex/CMakeLists.txt | 1 + cpp/src/dual_simplex/branch_and_bound.cpp | 100 ++++---- cpp/src/dual_simplex/branch_and_bound.hpp | 19 +- cpp/src/dual_simplex/node_presolve.cpp | 293 ++++++++++++++++++++++ cpp/src/dual_simplex/node_presolve.hpp | 49 ++++ cpp/src/dual_simplex/presolve.cpp | 274 +------------------- cpp/src/dual_simplex/presolve.hpp | 9 - 7 files changed, 415 insertions(+), 330 deletions(-) create mode 100644 cpp/src/dual_simplex/node_presolve.cpp create mode 100644 cpp/src/dual_simplex/node_presolve.hpp diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index c091ebac4f..c52121a631 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -27,6 +27,7 @@ set(DUAL_SIMPLEX_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/phase1.cpp ${CMAKE_CURRENT_SOURCE_DIR}/phase2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/node_presolve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/primal.cpp ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/right_looking_lu.cpp diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 0fd03a10a6..694aa3be03 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -555,8 +555,7 @@ template node_status_t branch_and_bound_t::solve_node(search_tree_t& search_tree, mip_node_t* node_ptr, lp_problem_t& leaf_problem, - const std::vector& bounds_changed, - const csc_matrix_t& Arow, + node_presolve_t& presolve, char thread_type, logger_t& log) { @@ -575,8 +574,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& // in B&B we only have equality constraints, leave it empty for default std::vector row_sense; - bool feasible = - bound_strengthening(row_sense, lp_settings, leaf_problem, Arow, var_types_, bounds_changed); + bool feasible = presolve.bound_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -680,6 +678,29 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& } } +template +void branch_and_bound_t::set_variable_bounds(mip_node_t* node, + std::vector& lower, + std::vector& upper, + std::vector& bounds_changed, + const std::vector& root_lower, + const std::vector& root_upper, + bool recompute) +{ + // Reset the bound_changed markers + std::fill(bounds_changed.begin(), bounds_changed.end(), false); + + // Recompute the bounds + if (recompute) { + lower = root_lower; + upper = root_upper; + node->get_variable_bounds(lower, upper, bounds_changed); + + } else { + node->update_variable_bound(lower, upper, bounds_changed); + } +} + template void branch_and_bound_t::exploration_ramp_up(search_tree_t* search_tree, mip_node_t* node, @@ -741,16 +762,20 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* return; } - const i_t n = leaf_problem.num_cols; - std::vector bounds_changed(n, false); + std::vector row_sense; + node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); // Set the correct bounds for the leaf problem - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; - node->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); + set_variable_bounds(node, + leaf_problem.lower, + leaf_problem.upper, + presolve.bounds_changed, + original_lp_.lower, + original_lp_.upper, + true); node_status_t node_status = - solve_node(*search_tree, node, leaf_problem, bounds_changed, Arow, 'B', settings_.log); + solve_node(*search_tree, node, leaf_problem, presolve, 'B', settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -785,11 +810,10 @@ void branch_and_bound_t::explore_subtree(i_t task_id, lp_problem_t& leaf_problem, const csc_matrix_t& Arow) { - const i_t m = leaf_problem.num_rows; - const i_t n = leaf_problem.num_cols; bool recompute = true; - std::vector bounds_changed(n); - lp_solution_t leaf_solution(m, n); + + std::vector row_sense; + node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); std::deque*> stack; stack.push_front(start_node); @@ -854,21 +878,17 @@ void branch_and_bound_t::explore_subtree(i_t task_id, return; } - // Reset the bound_changed markers - std::fill(bounds_changed.begin(), bounds_changed.end(), false); - - // Recompute the bounds - if (recompute) { - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - - } else { - node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); - } + // Set the correct bounds for the leaf problem + set_variable_bounds(node_ptr, + leaf_problem.lower, + leaf_problem.upper, + presolve.bounds_changed, + original_lp_.lower, + original_lp_.upper, + recompute); node_status_t node_status = - solve_node(search_tree, node_ptr, leaf_problem, bounds_changed, Arow, 'B', settings_.log); + solve_node(search_tree, node_ptr, leaf_problem, presolve, 'B', settings_.log); recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -974,7 +994,8 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr logger_t log; log.log = false; - const i_t n = leaf_problem.num_cols; + std::vector row_sense; + node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { @@ -988,7 +1009,6 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (get_upper_bound() < start_node->node.lower_bound) { continue; } bool recompute = true; - std::vector bounds_changed(n); search_tree_t subtree(std::move(start_node->node)); std::deque*> stack; stack.push_front(&subtree.root); @@ -1006,20 +1026,16 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr if (toc(stats_.start_time) > settings_.time_limit) { return; } - // Reset the bound_changed markers - std::fill(bounds_changed.begin(), bounds_changed.end(), false); - - // Recompute the bounds - if (recompute) { - leaf_problem.lower = start_node->lp_lower; - leaf_problem.upper = start_node->lp_upper; - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed); - } else { - node_ptr->update_variable_bound(leaf_problem.lower, leaf_problem.upper, bounds_changed); - } + // Set the correct bounds for the leaf problem + set_variable_bounds(node_ptr, + leaf_problem.lower, + leaf_problem.upper, + presolve.bounds_changed, + start_node->lower, + start_node->upper, + recompute); - node_status_t node_status = - solve_node(subtree, node_ptr, leaf_problem, bounds_changed, Arow, 'D', log); + node_status_t node_status = solve_node(subtree, node_ptr, leaf_problem, presolve, 'D', log); if (node_status == node_status_t::TIME_LIMIT) { return; diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 4ff7cfbde9..4f9847d500 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -19,8 +19,8 @@ #include #include +#include #include -#include #include #include #include @@ -58,13 +58,13 @@ void upper_bound_callback(f_t upper_bound); template struct diving_root_t { mip_node_t node; - std::vector lp_lower; - std::vector lp_upper; + std::vector lower; + std::vector upper; diving_root_t(mip_node_t&& node, const std::vector& lower, const std::vector& upper) - : node(std::move(node)), lp_upper(upper), lp_lower(lower) + : node(std::move(node)), upper(upper), lower(lower) { } @@ -226,6 +226,14 @@ class branch_and_bound_t { // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); + void set_variable_bounds(mip_node_t* node, + std::vector& lower, + std::vector& upper, + std::vector& bounds_changed, + const std::vector& root_lower, + const std::vector& root_upper, + bool recompute); + // Ramp-up phase of the solver, where we greedily expand the tree until // there is enough unexplored nodes. This is done recursively using OpenMP tasks. void exploration_ramp_up(search_tree_t* search_tree, @@ -256,8 +264,7 @@ class branch_and_bound_t { node_status_t solve_node(search_tree_t& search_tree, mip_node_t* node_ptr, lp_problem_t& leaf_problem, - const std::vector& bounds_changed, - const csc_matrix_t& Arow, + node_presolve_t& presolve, char thread_type, logger_t& log); diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/node_presolve.cpp new file mode 100644 index 0000000000..829ef5f0e1 --- /dev/null +++ b/cpp/src/dual_simplex/node_presolve.cpp @@ -0,0 +1,293 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +static inline f_t update_lb(f_t curr_lb, f_t coeff, f_t delta_min_act, f_t delta_max_act) +{ + auto comp_bnd = (coeff < 0.) ? delta_min_act / coeff : delta_max_act / coeff; + return std::max(curr_lb, comp_bnd); +} + +template +static inline f_t update_ub(f_t curr_ub, f_t coeff, f_t delta_min_act, f_t delta_max_act) +{ + auto comp_bnd = (coeff < 0.) ? delta_max_act / coeff : delta_min_act / coeff; + return std::min(curr_ub, comp_bnd); +} + +template +static inline bool check_infeasibility(f_t min_a, f_t max_a, f_t cnst_lb, f_t cnst_ub, f_t eps) +{ + return (min_a > cnst_ub + eps) || (max_a < cnst_lb - eps); +} + +#define DEBUG_BOUND_STRENGTHENING 0 + +template +void print_bounds_stats(const std::vector& lower, + const std::vector& upper, + const simplex_solver_settings_t& settings, + const std::string msg) +{ +#if DEBUG_BOUND_STRENGTHENING + f_t lb_norm = 0.0; + f_t ub_norm = 0.0; + + i_t sz = lower.size(); + for (i_t i = 0; i < sz; ++i) { + if (std::isfinite(lower[i])) { lb_norm += abs(lower[i]); } + if (std::isfinite(upper[i])) { ub_norm += abs(upper[i]); } + } + settings.log.printf("%s :: lb norm %e, ub norm %e\n", msg.c_str(), lb_norm, ub_norm); +#endif +} + +template +node_presolve_t::node_presolve_t(const lp_problem_t& problem, + const std::vector& row_sense, + const csc_matrix_t& Arow, + const std::vector& var_types) + : bounds_changed(problem.num_cols, false), + problem(problem), + Arow(Arow), + var_types(var_types), + delta_min_activity(problem.num_rows), + delta_max_activity(problem.num_rows), + constraint_lb(problem.num_rows), + constraint_ub(problem.num_rows) +{ + const bool is_row_sense_empty = row_sense.empty(); + if (is_row_sense_empty) { + std::copy(problem.rhs.begin(), problem.rhs.end(), constraint_lb.begin()); + std::copy(problem.rhs.begin(), problem.rhs.end(), constraint_ub.begin()); + } else { + // Set the constraint bounds + for (i_t i = 0; i < problem.num_rows; ++i) { + if (row_sense[i] == 'E') { + constraint_lb[i] = problem.rhs[i]; + constraint_ub[i] = problem.rhs[i]; + } else if (row_sense[i] == 'L') { + constraint_ub[i] = problem.rhs[i]; + constraint_lb[i] = -inf; + } else { + constraint_lb[i] = problem.rhs[i]; + constraint_ub[i] = inf; + } + } + } +} + +template +bool node_presolve_t::bound_strengthening( + std::vector& lower_bounds, + std::vector& upper_bounds, + const simplex_solver_settings_t& settings) +{ + const i_t m = problem.num_rows; + const i_t n = problem.num_cols; + + std::vector constraint_changed(m, false); + std::vector variable_changed(n, false); + std::vector constraint_changed_next(m, false); + + for (i_t i = 0; i < bounds_changed.size(); ++i) { + if (bounds_changed[i]) { + const i_t row_start = problem.A.col_start[i]; + const i_t row_end = problem.A.col_start[i + 1]; + for (i_t p = row_start; p < row_end; ++p) { + const i_t j = problem.A.i[p]; + constraint_changed[j] = true; + } + } + } + + std::vector lower = lower_bounds; + std::vector upper = upper_bounds; + print_bounds_stats(lower, upper, settings, "Initial bounds"); + + i_t iter = 0; + const i_t iter_limit = 10; + while (iter < iter_limit) { + for (i_t i = 0; i < m; ++i) { + if (!constraint_changed[i]) { continue; } + const i_t row_start = Arow.col_start[i]; + const i_t row_end = Arow.col_start[i + 1]; + + f_t min_a = 0.0; + f_t max_a = 0.0; + for (i_t p = row_start; p < row_end; ++p) { + const i_t j = Arow.i[p]; + const f_t a_ij = Arow.x[p]; + + variable_changed[j] = true; + if (a_ij > 0) { + min_a += a_ij * lower[j]; + max_a += a_ij * upper[j]; + } else if (a_ij < 0) { + min_a += a_ij * upper[j]; + max_a += a_ij * lower[j]; + } + if (upper[j] == inf && a_ij > 0) { max_a = inf; } + if (lower[j] == -inf && a_ij < 0) { max_a = inf; } + + if (lower[j] == -inf && a_ij > 0) { min_a = -inf; } + if (upper[j] == inf && a_ij < 0) { min_a = -inf; } + } + + f_t cnst_lb = constraint_lb[i]; + f_t cnst_ub = constraint_ub[i]; + bool is_infeasible = + check_infeasibility(min_a, max_a, cnst_lb, cnst_ub, settings.primal_tol); + if (is_infeasible) { + settings.log.printf( + "Iter:: %d, Infeasible constraint %d, cnst_lb %e, cnst_ub %e, min_a %e, max_a %e\n", + iter, + i, + cnst_lb, + cnst_ub, + min_a, + max_a); + return false; + } + + delta_min_activity[i] = cnst_ub - min_a; + delta_max_activity[i] = cnst_lb - max_a; + } + + i_t num_bounds_changed = 0; + + for (i_t k = 0; k < n; ++k) { + if (!variable_changed[k]) { continue; } + f_t old_lb = lower[k]; + f_t old_ub = upper[k]; + + f_t new_lb = old_lb; + f_t new_ub = old_ub; + + const i_t row_start = problem.A.col_start[k]; + const i_t row_end = problem.A.col_start[k + 1]; + for (i_t p = row_start; p < row_end; ++p) { + const i_t i = problem.A.i[p]; + + if (!constraint_changed[i]) { continue; } + const f_t a_ik = problem.A.x[p]; + + f_t delta_min_act = delta_min_activity[i]; + f_t delta_max_act = delta_max_activity[i]; + + delta_min_act += (a_ik < 0) ? a_ik * old_ub : a_ik * old_lb; + delta_max_act += (a_ik > 0) ? a_ik * old_ub : a_ik * old_lb; + + new_lb = std::max(new_lb, update_lb(old_lb, a_ik, delta_min_act, delta_max_act)); + new_ub = std::min(new_ub, update_ub(old_ub, a_ik, delta_min_act, delta_max_act)); + } + + // Integer rounding + if (!var_types.empty() && + (var_types[k] == variable_type_t::INTEGER || var_types[k] == variable_type_t::BINARY)) { + new_lb = std::ceil(new_lb - settings.integer_tol); + new_ub = std::floor(new_ub + settings.integer_tol); + } + + bool lb_updated = abs(new_lb - old_lb) > 1e3 * settings.primal_tol; + bool ub_updated = abs(new_ub - old_ub) > 1e3 * settings.primal_tol; + + new_lb = std::max(new_lb, lower_bounds[k]); + new_ub = std::min(new_ub, upper_bounds[k]); + + if (new_lb > new_ub + 1e-6) { + settings.log.printf( + "Iter:: %d, Infeasible variable after update %d, %e > %e\n", iter, k, new_lb, new_ub); + return false; + } + if (new_lb != old_lb || new_ub != old_ub) { + for (i_t p = row_start; p < row_end; ++p) { + const i_t i = problem.A.i[p]; + constraint_changed_next[i] = true; + } + } + + lower[k] = std::min(new_lb, new_ub); + upper[k] = std::max(new_lb, new_ub); + + bool bounds_changed = lb_updated || ub_updated; + if (bounds_changed) { num_bounds_changed++; } + } + + if (num_bounds_changed == 0) { break; } + + std::swap(constraint_changed, constraint_changed_next); + std::fill(constraint_changed_next.begin(), constraint_changed_next.end(), false); + std::fill(variable_changed.begin(), variable_changed.end(), false); + + iter++; + } + + // settings.log.printf("Total strengthened variables %d\n", total_strengthened_variables); + +#if DEBUG_BOUND_STRENGTHENING + f_t lb_change = 0.0; + f_t ub_change = 0.0; + int num_lb_changed = 0; + int num_ub_changed = 0; + + for (i_t i = 0; i < n; ++i) { + if (lower[i] > problem.lower[i] + settings.primal_tol || + (!std::isfinite(problem.lower[i]) && std::isfinite(lower[i]))) { + num_lb_changed++; + lb_change += + std::isfinite(problem.lower[i]) + ? (lower[i] - problem.lower[i]) / (1e-6 + std::max(abs(lower[i]), abs(problem.lower[i]))) + : 1.0; + } + if (upper[i] < problem.upper[i] - settings.primal_tol || + (!std::isfinite(problem.upper[i]) && std::isfinite(upper[i]))) { + num_ub_changed++; + ub_change += + std::isfinite(problem.upper[i]) + ? (problem.upper[i] - upper[i]) / (1e-6 + std::max(abs(problem.upper[i]), abs(upper[i]))) + : 1.0; + } + } + + if (num_lb_changed > 0 || num_ub_changed > 0) { + settings.log.printf( + "lb change %e, ub change %e, num lb changed %d, num ub changed %d, iter %d\n", + 100 * lb_change / std::max(1, num_lb_changed), + 100 * ub_change / std::max(1, num_ub_changed), + num_lb_changed, + num_ub_changed, + iter); + } + print_bounds_stats(lower, upper, settings, "Final bounds"); +#endif + + lower_bounds = lower; + upper_bounds = upper; + + return true; +} + +#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE +template class node_presolve_t; +#endif + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/node_presolve.hpp b/cpp/src/dual_simplex/node_presolve.hpp new file mode 100644 index 0000000000..54a3c6c341 --- /dev/null +++ b/cpp/src/dual_simplex/node_presolve.hpp @@ -0,0 +1,49 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +class node_presolve_t { + public: + // For pure LP bounds strengthening, var_types should be defaulted (i.e. left empty) + node_presolve_t(const lp_problem_t& problem, + const std::vector& row_sense, + const csc_matrix_t& Arow, + const std::vector& var_types); + + bool bound_strengthening(std::vector& lower_bounds, + std::vector& upper_bounds, + const simplex_solver_settings_t& settings); + + std::vector bounds_changed; + + private: + const lp_problem_t& problem; + const csc_matrix_t& Arow; + const std::vector& var_types; + + std::vector delta_min_activity; + std::vector delta_max_activity; + std::vector constraint_lb; + std::vector constraint_ub; +}; +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 8d80337c74..a84d64b369 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -28,272 +28,6 @@ namespace cuopt::linear_programming::dual_simplex { -template -static inline f_t update_lb(f_t curr_lb, f_t coeff, f_t delta_min_act, f_t delta_max_act) -{ - auto comp_bnd = (coeff < 0.) ? delta_min_act / coeff : delta_max_act / coeff; - return std::max(curr_lb, comp_bnd); -} - -template -static inline f_t update_ub(f_t curr_ub, f_t coeff, f_t delta_min_act, f_t delta_max_act) -{ - auto comp_bnd = (coeff < 0.) ? delta_max_act / coeff : delta_min_act / coeff; - return std::min(curr_ub, comp_bnd); -} - -template -static inline bool check_infeasibility(f_t min_a, f_t max_a, f_t cnst_lb, f_t cnst_ub, f_t eps) -{ - return (min_a > cnst_ub + eps) || (max_a < cnst_lb - eps); -} - -#define DEBUG_BOUND_STRENGTHENING 0 - -template -void print_bounds_stats(const std::vector& lower, - const std::vector& upper, - const simplex_solver_settings_t& settings, - const std::string msg) -{ -#if DEBUG_BOUND_STRENGTHENING - f_t lb_norm = 0.0; - f_t ub_norm = 0.0; - - i_t sz = lower.size(); - for (i_t i = 0; i < sz; ++i) { - if (std::isfinite(lower[i])) { lb_norm += abs(lower[i]); } - if (std::isfinite(upper[i])) { ub_norm += abs(upper[i]); } - } - settings.log.printf("%s :: lb norm %e, ub norm %e\n", msg.c_str(), lb_norm, ub_norm); -#endif -} - -template -bool bound_strengthening(const std::vector& row_sense, - const simplex_solver_settings_t& settings, - lp_problem_t& problem, - const csc_matrix_t& Arow, - const std::vector& var_types, - const std::vector& bounds_changed) -{ - const i_t m = problem.num_rows; - const i_t n = problem.num_cols; - - std::vector delta_min_activity(m); - std::vector delta_max_activity(m); - std::vector constraint_lb(m); - std::vector constraint_ub(m); - - // FIXME:: Instead of initializing constraint_changed to true, we can only look - // at the constraints corresponding to branched variable in branch and bound - // This is because, the parent LP already checked for feasibility of the constraints - // without the branched variable bounds - std::vector constraint_changed(m, true); - std::vector variable_changed(n, false); - std::vector constraint_changed_next(m, false); - - if (false && !bounds_changed.empty()) { - std::fill(constraint_changed.begin(), constraint_changed.end(), false); - for (i_t i = 0; i < n; ++i) { - if (bounds_changed[i]) { - const i_t row_start = problem.A.col_start[i]; - const i_t row_end = problem.A.col_start[i + 1]; - for (i_t p = row_start; p < row_end; ++p) { - const i_t j = problem.A.i[p]; - constraint_changed[j] = true; - } - } - } - } - - const bool is_row_sense_empty = row_sense.empty(); - if (is_row_sense_empty) { - std::copy(problem.rhs.begin(), problem.rhs.end(), constraint_lb.begin()); - std::copy(problem.rhs.begin(), problem.rhs.end(), constraint_ub.begin()); - } else { - // Set the constraint bounds - for (i_t i = 0; i < m; ++i) { - if (row_sense[i] == 'E') { - constraint_lb[i] = problem.rhs[i]; - constraint_ub[i] = problem.rhs[i]; - } else if (row_sense[i] == 'L') { - constraint_ub[i] = problem.rhs[i]; - constraint_lb[i] = -inf; - } else { - constraint_lb[i] = problem.rhs[i]; - constraint_ub[i] = inf; - } - } - } - - std::vector lower = problem.lower; - std::vector upper = problem.upper; - print_bounds_stats(lower, upper, settings, "Initial bounds"); - - i_t iter = 0; - const i_t iter_limit = 10; - while (iter < iter_limit) { - for (i_t i = 0; i < m; ++i) { - if (!constraint_changed[i]) { continue; } - const i_t row_start = Arow.col_start[i]; - const i_t row_end = Arow.col_start[i + 1]; - - f_t min_a = 0.0; - f_t max_a = 0.0; - for (i_t p = row_start; p < row_end; ++p) { - const i_t j = Arow.i[p]; - const f_t a_ij = Arow.x[p]; - - variable_changed[j] = true; - if (a_ij > 0) { - min_a += a_ij * lower[j]; - max_a += a_ij * upper[j]; - } else if (a_ij < 0) { - min_a += a_ij * upper[j]; - max_a += a_ij * lower[j]; - } - if (upper[j] == inf && a_ij > 0) { max_a = inf; } - if (lower[j] == -inf && a_ij < 0) { max_a = inf; } - - if (lower[j] == -inf && a_ij > 0) { min_a = -inf; } - if (upper[j] == inf && a_ij < 0) { min_a = -inf; } - } - - f_t cnst_lb = constraint_lb[i]; - f_t cnst_ub = constraint_ub[i]; - bool is_infeasible = - check_infeasibility(min_a, max_a, cnst_lb, cnst_ub, settings.primal_tol); - if (is_infeasible) { - settings.log.printf( - "Iter:: %d, Infeasible constraint %d, cnst_lb %e, cnst_ub %e, min_a %e, max_a %e\n", - iter, - i, - cnst_lb, - cnst_ub, - min_a, - max_a); - return false; - } - - delta_min_activity[i] = cnst_ub - min_a; - delta_max_activity[i] = cnst_lb - max_a; - } - - i_t num_bounds_changed = 0; - - for (i_t k = 0; k < n; ++k) { - if (!variable_changed[k]) { continue; } - f_t old_lb = lower[k]; - f_t old_ub = upper[k]; - - f_t new_lb = old_lb; - f_t new_ub = old_ub; - - const i_t row_start = problem.A.col_start[k]; - const i_t row_end = problem.A.col_start[k + 1]; - for (i_t p = row_start; p < row_end; ++p) { - const i_t i = problem.A.i[p]; - - if (!constraint_changed[i]) { continue; } - const f_t a_ik = problem.A.x[p]; - - f_t delta_min_act = delta_min_activity[i]; - f_t delta_max_act = delta_max_activity[i]; - - delta_min_act += (a_ik < 0) ? a_ik * old_ub : a_ik * old_lb; - delta_max_act += (a_ik > 0) ? a_ik * old_ub : a_ik * old_lb; - - new_lb = std::max(new_lb, update_lb(old_lb, a_ik, delta_min_act, delta_max_act)); - new_ub = std::min(new_ub, update_ub(old_ub, a_ik, delta_min_act, delta_max_act)); - } - - // Integer rounding - if (!var_types.empty() && - (var_types[k] == variable_type_t::INTEGER || var_types[k] == variable_type_t::BINARY)) { - new_lb = std::ceil(new_lb - settings.integer_tol); - new_ub = std::floor(new_ub + settings.integer_tol); - } - - bool lb_updated = abs(new_lb - old_lb) > 1e3 * settings.primal_tol; - bool ub_updated = abs(new_ub - old_ub) > 1e3 * settings.primal_tol; - - new_lb = std::max(new_lb, problem.lower[k]); - new_ub = std::min(new_ub, problem.upper[k]); - - if (new_lb > new_ub + 1e-6) { - settings.log.printf( - "Iter:: %d, Infeasible variable after update %d, %e > %e\n", iter, k, new_lb, new_ub); - return false; - } - if (new_lb != old_lb || new_ub != old_ub) { - for (i_t p = row_start; p < row_end; ++p) { - const i_t i = problem.A.i[p]; - constraint_changed_next[i] = true; - } - } - - lower[k] = std::min(new_lb, new_ub); - upper[k] = std::max(new_lb, new_ub); - - bool bounds_changed = lb_updated || ub_updated; - if (bounds_changed) { num_bounds_changed++; } - } - - if (num_bounds_changed == 0) { break; } - - std::swap(constraint_changed, constraint_changed_next); - std::fill(constraint_changed_next.begin(), constraint_changed_next.end(), false); - std::fill(variable_changed.begin(), variable_changed.end(), false); - - iter++; - } - - // settings.log.printf("Total strengthened variables %d\n", total_strengthened_variables); - -#if DEBUG_BOUND_STRENGTHENING - f_t lb_change = 0.0; - f_t ub_change = 0.0; - int num_lb_changed = 0; - int num_ub_changed = 0; - - for (i_t i = 0; i < n; ++i) { - if (lower[i] > problem.lower[i] + settings.primal_tol || - (!std::isfinite(problem.lower[i]) && std::isfinite(lower[i]))) { - num_lb_changed++; - lb_change += - std::isfinite(problem.lower[i]) - ? (lower[i] - problem.lower[i]) / (1e-6 + std::max(abs(lower[i]), abs(problem.lower[i]))) - : 1.0; - } - if (upper[i] < problem.upper[i] - settings.primal_tol || - (!std::isfinite(problem.upper[i]) && std::isfinite(upper[i]))) { - num_ub_changed++; - ub_change += - std::isfinite(problem.upper[i]) - ? (problem.upper[i] - upper[i]) / (1e-6 + std::max(abs(problem.upper[i]), abs(upper[i]))) - : 1.0; - } - } - - if (num_lb_changed > 0 || num_ub_changed > 0) { - settings.log.printf( - "lb change %e, ub change %e, num lb changed %d, num ub changed %d, iter %d\n", - 100 * lb_change / std::max(1, num_lb_changed), - 100 * ub_change / std::max(1, num_ub_changed), - num_lb_changed, - num_ub_changed, - iter); - } - print_bounds_stats(lower, upper, settings, "Final bounds"); -#endif - - problem.lower = lower; - problem.upper = upper; - - return true; -} - template i_t remove_empty_cols(lp_problem_t& problem, i_t& num_empty_cols, @@ -851,6 +585,7 @@ void convert_user_problem(const user_problem_t& user_problem, problem.A.transpose(Arow); bound_strengthening(row_sense, settings, problem, Arow); } + settings.log.debug( "equality rows %d less rows %d columns %d\n", equal_rows, less_rows, problem.num_cols); if (settings.barrier && settings.dualize != 0 && @@ -1608,13 +1343,6 @@ template void uncrush_solution(const presolve_info_t& std::vector& uncrushed_y, std::vector& uncrushed_z); -template bool bound_strengthening( - const std::vector& row_sense, - const simplex_solver_settings_t& settings, - lp_problem_t& problem, - const csc_matrix_t& Arow, - const std::vector& var_types, - const std::vector& bounds_changed); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index bf0aab8997..538ca5dffe 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -190,13 +190,4 @@ void uncrush_solution(const presolve_info_t& presolve_info, std::vector& uncrushed_y, std::vector& uncrushed_z); -// For pure LP bounds strengthening, var_types should be defaulted (i.e. left empty) -template -bool bound_strengthening(const std::vector& row_sense, - const simplex_solver_settings_t& settings, - lp_problem_t& problem, - const csc_matrix_t& Arow, - const std::vector& var_types = {}, - const std::vector& bounds_changed = {}); - } // namespace cuopt::linear_programming::dual_simplex From b5ed956903b414312ffba9be179d4c88302b8bda Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 22 Oct 2025 11:27:17 +0200 Subject: [PATCH 07/51] small refactor --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 694aa3be03..58e124b7ea 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -572,8 +572,6 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); - // in B&B we only have equality constraints, leave it empty for default - std::vector row_sense; bool feasible = presolve.bound_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; From abeb2cf06d59149e4b78c3a3bea3f36d25950294 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 22 Oct 2025 13:11:38 +0200 Subject: [PATCH 08/51] fix log spacing --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 58e124b7ea..44a381b76c 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -414,7 +414,7 @@ void branch_and_bound_t::repair_heuristic_solutions() std::string user_gap = user_mip_gap(obj, lower); settings_.log.printf( - "H %+13.6e %+10.6e %s %9.2f\n", + "H %+13.6e %+10.6e %s %9.2f\n", obj, lower, user_gap.c_str(), From 3237913698ab4887231876dd24e9526de1c45a34 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 22 Oct 2025 13:15:27 +0200 Subject: [PATCH 09/51] fix log spacing for repaired solutions --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 44a381b76c..1998e39bad 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -414,7 +414,7 @@ void branch_and_bound_t::repair_heuristic_solutions() std::string user_gap = user_mip_gap(obj, lower); settings_.log.printf( - "H %+13.6e %+10.6e %s %9.2f\n", + "H %+13.6e %+10.6e %s %9.2f\n", obj, lower, user_gap.c_str(), From f9853b2ee72747a21ed388198373e437691bbf47 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 22 Oct 2025 18:16:32 +0200 Subject: [PATCH 10/51] initialize the presolver only once per thread --- cpp/src/dual_simplex/branch_and_bound.cpp | 21 ++++++++------------- cpp/src/dual_simplex/branch_and_bound.hpp | 6 +++--- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 1998e39bad..56592ab734 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -806,13 +806,9 @@ void branch_and_bound_t::explore_subtree(i_t task_id, search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) + node_presolve_t& presolve) { bool recompute = true; - - std::vector row_sense; - node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); - std::deque*> stack; stack.push_front(start_node); @@ -932,7 +928,7 @@ template void branch_and_bound_t::best_first_thread(i_t id, search_tree_t& search_tree, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) + node_presolve_t& presolve) { f_t lower_bound = -inf; f_t upper_bound = inf; @@ -964,7 +960,7 @@ void branch_and_bound_t::best_first_thread(i_t id, } // Best-first search with plunging - explore_subtree(id, search_tree, node_ptr, leaf_problem, Arow); + explore_subtree(id, search_tree, node_ptr, leaf_problem, presolve); active_subtrees_--; } @@ -987,14 +983,11 @@ void branch_and_bound_t::best_first_thread(i_t id, template void branch_and_bound_t::diving_thread(lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) + node_presolve_t& presolve) { logger_t log; log.log = false; - std::vector row_sense; - node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); - while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { std::optional> start_node; @@ -1225,6 +1218,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut { // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; + std::vector row_sense; + node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); #pragma omp master { @@ -1247,12 +1242,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut (active_subtrees_ > 0 || get_heap_size() > 0)) { for (i_t i = 0; i < settings_.num_bfs_threads; i++) { #pragma omp task - best_first_thread(i, search_tree, leaf_problem, Arow); + best_first_thread(i, search_tree, leaf_problem, presolve); } for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(leaf_problem, Arow); + diving_thread(leaf_problem, presolve); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 4f9847d500..79c8c6650b 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -247,18 +247,18 @@ class branch_and_bound_t { search_tree_t& search_tree, mip_node_t* start_node, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow); + node_presolve_t& presolve); // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. void best_first_thread(i_t id, search_tree_t& search_tree, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow); + node_presolve_t& presolve); // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(lp_problem_t& leaf_problem, const csc_matrix_t& Arow); + void diving_thread(lp_problem_t& leaf_problem, node_presolve_t& presolve); // Solve the LP relaxation of a leaf node and update the tree. node_status_t solve_node(search_tree_t& search_tree, From 6ccb6d9f1b205e675c09d0e6a81d7d9acf006b35 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 23 Oct 2025 17:05:54 +0200 Subject: [PATCH 11/51] fixed race condition --- cpp/src/dual_simplex/branch_and_bound.cpp | 81 ++++++++++++----------- cpp/src/dual_simplex/branch_and_bound.hpp | 20 +++--- cpp/src/dual_simplex/node_presolve.cpp | 32 ++++----- cpp/src/dual_simplex/node_presolve.hpp | 10 +-- 4 files changed, 74 insertions(+), 69 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 56592ab734..a1d59cc858 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -572,7 +572,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); - bool feasible = presolve.bound_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); + bool feasible = presolve.bound_strengthening(leaf_problem.lower, leaf_problem.upper); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -702,8 +702,6 @@ void branch_and_bound_t::set_variable_bounds(mip_node_t* nod template void branch_and_bound_t::exploration_ramp_up(search_tree_t* search_tree, mip_node_t* node, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, i_t initial_heap_size) { if (status_ != mip_exploration_status_t::RUNNING) { return; } @@ -713,6 +711,10 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* // to repair the heuristic solution. repair_heuristic_solutions(); + i_t tid = omp_get_thread_num(); + auto presolve = thread_data_[tid].presolve; + auto leaf_problem = thread_data_[tid].leaf_problem; + f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); @@ -760,20 +762,17 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* return; } - std::vector row_sense; - node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); - // Set the correct bounds for the leaf problem set_variable_bounds(node, - leaf_problem.lower, - leaf_problem.upper, - presolve.bounds_changed, + leaf_problem->lower, + leaf_problem->upper, + presolve->bounds_changed, original_lp_.lower, original_lp_.upper, true); node_status_t node_status = - solve_node(*search_tree, node, leaf_problem, presolve, 'B', settings_.log); + solve_node(*search_tree, node, *leaf_problem, *presolve, 'B', settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -785,11 +784,10 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* // If we haven't generated enough nodes to keep the threads busy, continue the ramp up phase if (stats_.nodes_unexplored < initial_heap_size) { #pragma omp task - exploration_ramp_up( - search_tree, node->get_down_child(), leaf_problem, Arow, initial_heap_size); + exploration_ramp_up(search_tree, node->get_down_child(), initial_heap_size); #pragma omp task - exploration_ramp_up(search_tree, node->get_up_child(), leaf_problem, Arow, initial_heap_size); + exploration_ramp_up(search_tree, node->get_up_child(), initial_heap_size); } else { // We've generated enough nodes, push further nodes onto the heap @@ -804,14 +802,16 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* template void branch_and_bound_t::explore_subtree(i_t task_id, search_tree_t& search_tree, - mip_node_t* start_node, - lp_problem_t& leaf_problem, - node_presolve_t& presolve) + mip_node_t* start_node) { bool recompute = true; std::deque*> stack; stack.push_front(start_node); + i_t tid = omp_get_thread_num(); + auto presolve = thread_data_[tid].presolve; + auto leaf_problem = thread_data_[tid].leaf_problem; + while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { if (task_id == 0) { repair_heuristic_solutions(); } @@ -874,15 +874,15 @@ void branch_and_bound_t::explore_subtree(i_t task_id, // Set the correct bounds for the leaf problem set_variable_bounds(node_ptr, - leaf_problem.lower, - leaf_problem.upper, - presolve.bounds_changed, + leaf_problem->lower, + leaf_problem->upper, + presolve->bounds_changed, original_lp_.lower, original_lp_.upper, recompute); node_status_t node_status = - solve_node(search_tree, node_ptr, leaf_problem, presolve, 'B', settings_.log); + solve_node(search_tree, node_ptr, *leaf_problem, *presolve, 'B', settings_.log); recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -906,7 +906,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, // would be better if we discard the node instead. if (get_heap_size() > settings_.num_bfs_threads) { mutex_dive_queue_.lock(); - dive_queue_.emplace(node->detach_copy(), leaf_problem.lower, leaf_problem.upper); + dive_queue_.emplace(node->detach_copy(), leaf_problem->lower, leaf_problem->upper); mutex_dive_queue_.unlock(); } @@ -925,10 +925,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, } template -void branch_and_bound_t::best_first_thread(i_t id, - search_tree_t& search_tree, - lp_problem_t& leaf_problem, - node_presolve_t& presolve) +void branch_and_bound_t::best_first_thread(i_t id, search_tree_t& search_tree) { f_t lower_bound = -inf; f_t upper_bound = inf; @@ -960,7 +957,7 @@ void branch_and_bound_t::best_first_thread(i_t id, } // Best-first search with plunging - explore_subtree(id, search_tree, node_ptr, leaf_problem, presolve); + explore_subtree(id, search_tree, node_ptr); active_subtrees_--; } @@ -982,12 +979,15 @@ void branch_and_bound_t::best_first_thread(i_t id, } template -void branch_and_bound_t::diving_thread(lp_problem_t& leaf_problem, - node_presolve_t& presolve) +void branch_and_bound_t::diving_thread() { logger_t log; log.log = false; + i_t tid = omp_get_thread_num(); + auto presolve = thread_data_[tid].presolve; + auto leaf_problem = thread_data_[tid].leaf_problem; + while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { std::optional> start_node; @@ -1019,14 +1019,15 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr // Set the correct bounds for the leaf problem set_variable_bounds(node_ptr, - leaf_problem.lower, - leaf_problem.upper, - presolve.bounds_changed, + leaf_problem->lower, + leaf_problem->upper, + presolve->bounds_changed, start_node->lower, start_node->upper, recompute); - node_status_t node_status = solve_node(subtree, node_ptr, leaf_problem, presolve, 'D', log); + node_status_t node_status = + solve_node(subtree, node_ptr, *leaf_problem, *presolve, 'D', log); if (node_status == node_status_t::TIME_LIMIT) { return; @@ -1046,7 +1047,7 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr mutex_dive_queue_.lock(); mip_node_t* new_node = stack.back(); stack.pop_back(); - dive_queue_.emplace(new_node->detach_copy(), leaf_problem.lower, leaf_problem.upper); + dive_queue_.emplace(new_node->detach_copy(), leaf_problem->lower, leaf_problem->upper); mutex_dive_queue_.unlock(); } } @@ -1213,13 +1214,17 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; + thread_data_.resize(settings_.num_threads, {nullptr, nullptr}); #pragma omp parallel num_threads(settings_.num_threads) { // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_); + node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_, settings_); + + thread_data_[omp_get_thread_num()].presolve = &presolve; + thread_data_[omp_get_thread_num()].leaf_problem = &leaf_problem; #pragma omp master { @@ -1228,10 +1233,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t initial_size = 2 * settings_.num_threads; #pragma omp task - exploration_ramp_up(&search_tree, down_child, leaf_problem, Arow, initial_size); + exploration_ramp_up(&search_tree, down_child, initial_size); #pragma omp task - exploration_ramp_up(&search_tree, up_child, leaf_problem, Arow, initial_size); + exploration_ramp_up(&search_tree, up_child, initial_size); } #pragma omp barrier @@ -1242,12 +1247,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut (active_subtrees_ > 0 || get_heap_size() > 0)) { for (i_t i = 0; i < settings_.num_bfs_threads; i++) { #pragma omp task - best_first_thread(i, search_tree, leaf_problem, presolve); + best_first_thread(i, search_tree); } for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(leaf_problem, presolve); + diving_thread(); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 79c8c6650b..f0d1b8b4b0 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -213,6 +213,13 @@ class branch_and_bound_t { // its blocks the progression of the lower bound. omp_atomic_t lower_bound_ceiling_; + // Stores the node presolver per thread. + struct thread_data_t { + node_presolve_t* presolve; + lp_problem_t* leaf_problem; + }; + std::vector thread_data_; + // Set the final solution. mip_status_t set_final_solution(mip_solution_t& solution, f_t lower_bound); @@ -238,27 +245,20 @@ class branch_and_bound_t { // there is enough unexplored nodes. This is done recursively using OpenMP tasks. void exploration_ramp_up(search_tree_t* search_tree, mip_node_t* node, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. void explore_subtree(i_t task_id, search_tree_t& search_tree, - mip_node_t* start_node, - lp_problem_t& leaf_problem, - node_presolve_t& presolve); + mip_node_t* start_node); // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. - void best_first_thread(i_t id, - search_tree_t& search_tree, - lp_problem_t& leaf_problem, - node_presolve_t& presolve); + void best_first_thread(i_t id, search_tree_t& search_tree); // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(lp_problem_t& leaf_problem, node_presolve_t& presolve); + void diving_thread(); // Solve the LP relaxation of a leaf node and update the tree. node_status_t solve_node(search_tree_t& search_tree, diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/node_presolve.cpp index 829ef5f0e1..a8b773b914 100644 --- a/cpp/src/dual_simplex/node_presolve.cpp +++ b/cpp/src/dual_simplex/node_presolve.cpp @@ -64,11 +64,13 @@ template node_presolve_t::node_presolve_t(const lp_problem_t& problem, const std::vector& row_sense, const csc_matrix_t& Arow, - const std::vector& var_types) + const std::vector& var_types, + const simplex_solver_settings_t& settings) : bounds_changed(problem.num_cols, false), - problem(problem), + A(problem.A), Arow(Arow), var_types(var_types), + settings(settings), delta_min_activity(problem.num_rows), delta_max_activity(problem.num_rows), constraint_lb(problem.num_rows), @@ -96,13 +98,11 @@ node_presolve_t::node_presolve_t(const lp_problem_t& problem } template -bool node_presolve_t::bound_strengthening( - std::vector& lower_bounds, - std::vector& upper_bounds, - const simplex_solver_settings_t& settings) +bool node_presolve_t::bound_strengthening(std::vector& lower_bounds, + std::vector& upper_bounds) { - const i_t m = problem.num_rows; - const i_t n = problem.num_cols; + const i_t m = A.m; + const i_t n = A.n; std::vector constraint_changed(m, false); std::vector variable_changed(n, false); @@ -110,10 +110,10 @@ bool node_presolve_t::bound_strengthening( for (i_t i = 0; i < bounds_changed.size(); ++i) { if (bounds_changed[i]) { - const i_t row_start = problem.A.col_start[i]; - const i_t row_end = problem.A.col_start[i + 1]; + const i_t row_start = A.col_start[i]; + const i_t row_end = A.col_start[i + 1]; for (i_t p = row_start; p < row_end; ++p) { - const i_t j = problem.A.i[p]; + const i_t j = A.i[p]; constraint_changed[j] = true; } } @@ -182,13 +182,13 @@ bool node_presolve_t::bound_strengthening( f_t new_lb = old_lb; f_t new_ub = old_ub; - const i_t row_start = problem.A.col_start[k]; - const i_t row_end = problem.A.col_start[k + 1]; + const i_t row_start = A.col_start[k]; + const i_t row_end = A.col_start[k + 1]; for (i_t p = row_start; p < row_end; ++p) { - const i_t i = problem.A.i[p]; + const i_t i = A.i[p]; if (!constraint_changed[i]) { continue; } - const f_t a_ik = problem.A.x[p]; + const f_t a_ik = A.x[p]; f_t delta_min_act = delta_min_activity[i]; f_t delta_max_act = delta_max_activity[i]; @@ -220,7 +220,7 @@ bool node_presolve_t::bound_strengthening( } if (new_lb != old_lb || new_ub != old_ub) { for (i_t p = row_start; p < row_end; ++p) { - const i_t i = problem.A.i[p]; + const i_t i = A.i[p]; constraint_changed_next[i] = true; } } diff --git a/cpp/src/dual_simplex/node_presolve.hpp b/cpp/src/dual_simplex/node_presolve.hpp index 54a3c6c341..1ebcf750d8 100644 --- a/cpp/src/dual_simplex/node_presolve.hpp +++ b/cpp/src/dual_simplex/node_presolve.hpp @@ -28,18 +28,18 @@ class node_presolve_t { node_presolve_t(const lp_problem_t& problem, const std::vector& row_sense, const csc_matrix_t& Arow, - const std::vector& var_types); + const std::vector& var_types, + const simplex_solver_settings_t& settings); - bool bound_strengthening(std::vector& lower_bounds, - std::vector& upper_bounds, - const simplex_solver_settings_t& settings); + bool bound_strengthening(std::vector& lower_bounds, std::vector& upper_bounds); std::vector bounds_changed; private: - const lp_problem_t& problem; + const csc_matrix_t& A; const csc_matrix_t& Arow; const std::vector& var_types; + const simplex_solver_settings_t& settings; std::vector delta_min_activity; std::vector delta_max_activity; From 8f55932f527ba0295060dc0246ff31500b598c79 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 23 Oct 2025 17:20:15 +0200 Subject: [PATCH 12/51] small refactor --- cpp/src/dual_simplex/branch_and_bound.cpp | 31 ++++++++++++----------- cpp/src/dual_simplex/branch_and_bound.hpp | 4 +-- cpp/src/dual_simplex/node_presolve.cpp | 16 ++++++------ cpp/src/dual_simplex/node_presolve.hpp | 12 ++++----- 4 files changed, 32 insertions(+), 31 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index a1d59cc858..1390b0091f 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -555,7 +555,7 @@ template node_status_t branch_and_bound_t::solve_node(search_tree_t& search_tree, mip_node_t* node_ptr, lp_problem_t& leaf_problem, - node_presolve_t& presolve, + node_presolver_t& presolver, char thread_type, logger_t& log) { @@ -572,7 +572,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); - bool feasible = presolve.bound_strengthening(leaf_problem.lower, leaf_problem.upper); + bool feasible = presolver.bound_strengthening(leaf_problem.lower, leaf_problem.upper); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -712,7 +712,7 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* repair_heuristic_solutions(); i_t tid = omp_get_thread_num(); - auto presolve = thread_data_[tid].presolve; + auto presolver = thread_data_[tid].presolver; auto leaf_problem = thread_data_[tid].leaf_problem; f_t lower_bound = node->lower_bound; @@ -766,13 +766,13 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* set_variable_bounds(node, leaf_problem->lower, leaf_problem->upper, - presolve->bounds_changed, + presolver->bounds_changed, original_lp_.lower, original_lp_.upper, true); node_status_t node_status = - solve_node(*search_tree, node, *leaf_problem, *presolve, 'B', settings_.log); + solve_node(*search_tree, node, *leaf_problem, *presolver, 'B', settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -809,7 +809,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, stack.push_front(start_node); i_t tid = omp_get_thread_num(); - auto presolve = thread_data_[tid].presolve; + auto presolver = thread_data_[tid].presolver; auto leaf_problem = thread_data_[tid].leaf_problem; while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { @@ -876,13 +876,13 @@ void branch_and_bound_t::explore_subtree(i_t task_id, set_variable_bounds(node_ptr, leaf_problem->lower, leaf_problem->upper, - presolve->bounds_changed, + presolver->bounds_changed, original_lp_.lower, original_lp_.upper, recompute); node_status_t node_status = - solve_node(search_tree, node_ptr, *leaf_problem, *presolve, 'B', settings_.log); + solve_node(search_tree, node_ptr, *leaf_problem, *presolver, 'B', settings_.log); recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -985,7 +985,7 @@ void branch_and_bound_t::diving_thread() log.log = false; i_t tid = omp_get_thread_num(); - auto presolve = thread_data_[tid].presolve; + auto presolver = thread_data_[tid].presolver; auto leaf_problem = thread_data_[tid].leaf_problem; while (status_ == mip_exploration_status_t::RUNNING && @@ -1021,13 +1021,13 @@ void branch_and_bound_t::diving_thread() set_variable_bounds(node_ptr, leaf_problem->lower, leaf_problem->upper, - presolve->bounds_changed, + presolver->bounds_changed, start_node->lower, start_node->upper, recompute); node_status_t node_status = - solve_node(subtree, node_ptr, *leaf_problem, *presolve, 'D', log); + solve_node(subtree, node_ptr, *leaf_problem, *presolver, 'D', log); if (node_status == node_status_t::TIME_LIMIT) { return; @@ -1214,17 +1214,18 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; - thread_data_.resize(settings_.num_threads, {nullptr, nullptr}); + thread_data_.resize(settings_.num_threads); #pragma omp parallel num_threads(settings_.num_threads) { // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolve_t presolve(leaf_problem, row_sense, Arow, var_types_, settings_); + node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_, settings_); - thread_data_[omp_get_thread_num()].presolve = &presolve; - thread_data_[omp_get_thread_num()].leaf_problem = &leaf_problem; + i_t tid = omp_get_thread_num(); + thread_data_[tid].presolver = &presolver; + thread_data_[tid].leaf_problem = &leaf_problem; #pragma omp master { diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index f0d1b8b4b0..36e7c17a8b 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -215,7 +215,7 @@ class branch_and_bound_t { // Stores the node presolver per thread. struct thread_data_t { - node_presolve_t* presolve; + node_presolver_t* presolver; lp_problem_t* leaf_problem; }; std::vector thread_data_; @@ -264,7 +264,7 @@ class branch_and_bound_t { node_status_t solve_node(search_tree_t& search_tree, mip_node_t* node_ptr, lp_problem_t& leaf_problem, - node_presolve_t& presolve, + node_presolver_t& presolve, char thread_type, logger_t& log); diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/node_presolve.cpp index a8b773b914..c4eb5ff0fd 100644 --- a/cpp/src/dual_simplex/node_presolve.cpp +++ b/cpp/src/dual_simplex/node_presolve.cpp @@ -61,11 +61,11 @@ void print_bounds_stats(const std::vector& lower, } template -node_presolve_t::node_presolve_t(const lp_problem_t& problem, - const std::vector& row_sense, - const csc_matrix_t& Arow, - const std::vector& var_types, - const simplex_solver_settings_t& settings) +node_presolver_t::node_presolver_t(const lp_problem_t& problem, + const std::vector& row_sense, + const csc_matrix_t& Arow, + const std::vector& var_types, + const simplex_solver_settings_t& settings) : bounds_changed(problem.num_cols, false), A(problem.A), Arow(Arow), @@ -98,8 +98,8 @@ node_presolve_t::node_presolve_t(const lp_problem_t& problem } template -bool node_presolve_t::bound_strengthening(std::vector& lower_bounds, - std::vector& upper_bounds) +bool node_presolver_t::bound_strengthening(std::vector& lower_bounds, + std::vector& upper_bounds) { const i_t m = A.m; const i_t n = A.n; @@ -287,7 +287,7 @@ bool node_presolve_t::bound_strengthening(std::vector& lower_boun } #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE -template class node_presolve_t; +template class node_presolver_t; #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/node_presolve.hpp b/cpp/src/dual_simplex/node_presolve.hpp index 1ebcf750d8..910088ce3d 100644 --- a/cpp/src/dual_simplex/node_presolve.hpp +++ b/cpp/src/dual_simplex/node_presolve.hpp @@ -22,14 +22,14 @@ namespace cuopt::linear_programming::dual_simplex { template -class node_presolve_t { +class node_presolver_t { public: // For pure LP bounds strengthening, var_types should be defaulted (i.e. left empty) - node_presolve_t(const lp_problem_t& problem, - const std::vector& row_sense, - const csc_matrix_t& Arow, - const std::vector& var_types, - const simplex_solver_settings_t& settings); + node_presolver_t(const lp_problem_t& problem, + const std::vector& row_sense, + const csc_matrix_t& Arow, + const std::vector& var_types, + const simplex_solver_settings_t& settings); bool bound_strengthening(std::vector& lower_bounds, std::vector& upper_bounds); From 7f337d4c14f3171b84dc822f3b004b34d51dce95 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 23 Oct 2025 17:33:59 +0200 Subject: [PATCH 13/51] fixed missing template parameters --- cpp/src/dual_simplex/branch_and_bound.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 1390b0091f..dc7ee39675 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -1214,12 +1214,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; - thread_data_.resize(settings_.num_threads); + thread_data_.resize(settings_.num_threads, {nullptr, nullptr}); #pragma omp parallel num_threads(settings_.num_threads) { // Make a copy of the original LP. We will modify its bounds at each leaf - lp_problem_t leaf_problem = original_lp_; + lp_problem_t leaf_problem = original_lp_; std::vector row_sense; node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_, settings_); From 06ca86d87f5db91788e49fd8c89aa8bd40612c43 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 23 Oct 2025 19:08:09 +0200 Subject: [PATCH 14/51] fixed crash --- cpp/src/dual_simplex/branch_and_bound.cpp | 97 +++++++++++------------ cpp/src/dual_simplex/branch_and_bound.hpp | 28 +++---- cpp/src/dual_simplex/node_presolve.cpp | 10 +-- cpp/src/dual_simplex/node_presolve.hpp | 8 +- 4 files changed, 70 insertions(+), 73 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index dc7ee39675..dd8efd562c 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -552,8 +552,8 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) } template -node_status_t branch_and_bound_t::solve_node(search_tree_t& search_tree, - mip_node_t* node_ptr, +node_status_t branch_and_bound_t::solve_node(mip_node_t* node_ptr, + search_tree_t& search_tree, lp_problem_t& leaf_problem, node_presolver_t& presolver, char thread_type, @@ -572,7 +572,8 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); - bool feasible = presolver.bound_strengthening(leaf_problem.lower, leaf_problem.upper); + bool feasible = + presolver.bound_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -700,8 +701,9 @@ void branch_and_bound_t::set_variable_bounds(mip_node_t* nod } template -void branch_and_bound_t::exploration_ramp_up(search_tree_t* search_tree, - mip_node_t* node, +void branch_and_bound_t::exploration_ramp_up(mip_node_t* node, + search_tree_t* search_tree, + const csc_matrix_t& Arow, i_t initial_heap_size) { if (status_ != mip_exploration_status_t::RUNNING) { return; } @@ -711,9 +713,10 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* // to repair the heuristic solution. repair_heuristic_solutions(); - i_t tid = omp_get_thread_num(); - auto presolver = thread_data_[tid].presolver; - auto leaf_problem = thread_data_[tid].leaf_problem; + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + std::vector row_sense; + node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_); f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); @@ -764,15 +767,15 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* // Set the correct bounds for the leaf problem set_variable_bounds(node, - leaf_problem->lower, - leaf_problem->upper, - presolver->bounds_changed, + leaf_problem.lower, + leaf_problem.upper, + presolver.bounds_changed, original_lp_.lower, original_lp_.upper, true); node_status_t node_status = - solve_node(*search_tree, node, *leaf_problem, *presolver, 'B', settings_.log); + solve_node(node, *search_tree, leaf_problem, presolver, 'B', settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -784,10 +787,10 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* // If we haven't generated enough nodes to keep the threads busy, continue the ramp up phase if (stats_.nodes_unexplored < initial_heap_size) { #pragma omp task - exploration_ramp_up(search_tree, node->get_down_child(), initial_heap_size); + exploration_ramp_up(node->get_down_child(), search_tree, Arow, initial_heap_size); #pragma omp task - exploration_ramp_up(search_tree, node->get_up_child(), initial_heap_size); + exploration_ramp_up(node->get_up_child(), search_tree, Arow, initial_heap_size); } else { // We've generated enough nodes, push further nodes onto the heap @@ -801,17 +804,15 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* template void branch_and_bound_t::explore_subtree(i_t task_id, + mip_node_t* start_node, search_tree_t& search_tree, - mip_node_t* start_node) + lp_problem_t& leaf_problem, + node_presolver_t& presolver) { bool recompute = true; std::deque*> stack; stack.push_front(start_node); - i_t tid = omp_get_thread_num(); - auto presolver = thread_data_[tid].presolver; - auto leaf_problem = thread_data_[tid].leaf_problem; - while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { if (task_id == 0) { repair_heuristic_solutions(); } @@ -874,15 +875,15 @@ void branch_and_bound_t::explore_subtree(i_t task_id, // Set the correct bounds for the leaf problem set_variable_bounds(node_ptr, - leaf_problem->lower, - leaf_problem->upper, - presolver->bounds_changed, + leaf_problem.lower, + leaf_problem.upper, + presolver.bounds_changed, original_lp_.lower, original_lp_.upper, recompute); node_status_t node_status = - solve_node(search_tree, node_ptr, *leaf_problem, *presolver, 'B', settings_.log); + solve_node(node_ptr, search_tree, leaf_problem, presolver, 'B', settings_.log); recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -906,7 +907,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, // would be better if we discard the node instead. if (get_heap_size() > settings_.num_bfs_threads) { mutex_dive_queue_.lock(); - dive_queue_.emplace(node->detach_copy(), leaf_problem->lower, leaf_problem->upper); + dive_queue_.emplace(node->detach_copy(), leaf_problem.lower, leaf_problem.upper); mutex_dive_queue_.unlock(); } @@ -925,13 +926,20 @@ void branch_and_bound_t::explore_subtree(i_t task_id, } template -void branch_and_bound_t::best_first_thread(i_t id, search_tree_t& search_tree) +void branch_and_bound_t::best_first_thread(i_t id, + search_tree_t& search_tree, + const csc_matrix_t& Arow) { f_t lower_bound = -inf; f_t upper_bound = inf; f_t abs_gap = inf; f_t rel_gap = inf; + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + std::vector row_sense; + node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_); + while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && (active_subtrees_ > 0 || get_heap_size() > 0)) { @@ -957,7 +965,7 @@ void branch_and_bound_t::best_first_thread(i_t id, search_tree_t::best_first_thread(i_t id, search_tree_t -void branch_and_bound_t::diving_thread() +void branch_and_bound_t::diving_thread(const csc_matrix_t& Arow) { logger_t log; log.log = false; - i_t tid = omp_get_thread_num(); - auto presolver = thread_data_[tid].presolver; - auto leaf_problem = thread_data_[tid].leaf_problem; + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + std::vector row_sense; + node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_); while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { @@ -1019,15 +1028,15 @@ void branch_and_bound_t::diving_thread() // Set the correct bounds for the leaf problem set_variable_bounds(node_ptr, - leaf_problem->lower, - leaf_problem->upper, - presolver->bounds_changed, + leaf_problem.lower, + leaf_problem.upper, + presolver.bounds_changed, start_node->lower, start_node->upper, recompute); node_status_t node_status = - solve_node(subtree, node_ptr, *leaf_problem, *presolver, 'D', log); + solve_node(node_ptr, subtree, leaf_problem, presolver, 'D', log); if (node_status == node_status_t::TIME_LIMIT) { return; @@ -1047,7 +1056,7 @@ void branch_and_bound_t::diving_thread() mutex_dive_queue_.lock(); mip_node_t* new_node = stack.back(); stack.pop_back(); - dive_queue_.emplace(new_node->detach_copy(), leaf_problem->lower, leaf_problem->upper); + dive_queue_.emplace(new_node->detach_copy(), leaf_problem.lower, leaf_problem.upper); mutex_dive_queue_.unlock(); } } @@ -1214,19 +1223,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; - thread_data_.resize(settings_.num_threads, {nullptr, nullptr}); #pragma omp parallel num_threads(settings_.num_threads) { - // Make a copy of the original LP. We will modify its bounds at each leaf - lp_problem_t leaf_problem = original_lp_; - std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_, settings_); - - i_t tid = omp_get_thread_num(); - thread_data_[tid].presolver = &presolver; - thread_data_[tid].leaf_problem = &leaf_problem; - #pragma omp master { auto down_child = search_tree.root.get_down_child(); @@ -1234,10 +1233,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t initial_size = 2 * settings_.num_threads; #pragma omp task - exploration_ramp_up(&search_tree, down_child, initial_size); + exploration_ramp_up(down_child, &search_tree, Arow, initial_size); #pragma omp task - exploration_ramp_up(&search_tree, up_child, initial_size); + exploration_ramp_up(up_child, &search_tree, Arow, initial_size); } #pragma omp barrier @@ -1248,12 +1247,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut (active_subtrees_ > 0 || get_heap_size() > 0)) { for (i_t i = 0; i < settings_.num_bfs_threads; i++) { #pragma omp task - best_first_thread(i, search_tree); + best_first_thread(i, search_tree, Arow); } for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(); + diving_thread(Arow); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 36e7c17a8b..d92a62dc91 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -213,13 +213,6 @@ class branch_and_bound_t { // its blocks the progression of the lower bound. omp_atomic_t lower_bound_ceiling_; - // Stores the node presolver per thread. - struct thread_data_t { - node_presolver_t* presolver; - lp_problem_t* leaf_problem; - }; - std::vector thread_data_; - // Set the final solution. mip_status_t set_final_solution(mip_solution_t& solution, f_t lower_bound); @@ -243,28 +236,33 @@ class branch_and_bound_t { // Ramp-up phase of the solver, where we greedily expand the tree until // there is enough unexplored nodes. This is done recursively using OpenMP tasks. - void exploration_ramp_up(search_tree_t* search_tree, - mip_node_t* node, + void exploration_ramp_up(mip_node_t* node, + search_tree_t* search_tree, + const csc_matrix_t& Arow, i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. void explore_subtree(i_t task_id, + mip_node_t* start_node, search_tree_t& search_tree, - mip_node_t* start_node); + lp_problem_t& leaf_problem, + node_presolver_t& presolver); // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. - void best_first_thread(i_t id, search_tree_t& search_tree); + void best_first_thread(i_t id, + search_tree_t& search_tree, + const csc_matrix_t& Arow); // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(); + void diving_thread(const csc_matrix_t& Arow); // Solve the LP relaxation of a leaf node and update the tree. - node_status_t solve_node(search_tree_t& search_tree, - mip_node_t* node_ptr, + node_status_t solve_node(mip_node_t* node_ptr, + search_tree_t& search_tree, lp_problem_t& leaf_problem, - node_presolver_t& presolve, + node_presolver_t& presolver, char thread_type, logger_t& log); diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/node_presolve.cpp index c4eb5ff0fd..f05bfab5a5 100644 --- a/cpp/src/dual_simplex/node_presolve.cpp +++ b/cpp/src/dual_simplex/node_presolve.cpp @@ -64,13 +64,11 @@ template node_presolver_t::node_presolver_t(const lp_problem_t& problem, const std::vector& row_sense, const csc_matrix_t& Arow, - const std::vector& var_types, - const simplex_solver_settings_t& settings) + const std::vector& var_types) : bounds_changed(problem.num_cols, false), A(problem.A), Arow(Arow), var_types(var_types), - settings(settings), delta_min_activity(problem.num_rows), delta_max_activity(problem.num_rows), constraint_lb(problem.num_rows), @@ -98,8 +96,10 @@ node_presolver_t::node_presolver_t(const lp_problem_t& probl } template -bool node_presolver_t::bound_strengthening(std::vector& lower_bounds, - std::vector& upper_bounds) +bool node_presolver_t::bound_strengthening( + std::vector& lower_bounds, + std::vector& upper_bounds, + const simplex_solver_settings_t& settings) { const i_t m = A.m; const i_t n = A.n; diff --git a/cpp/src/dual_simplex/node_presolve.hpp b/cpp/src/dual_simplex/node_presolve.hpp index 910088ce3d..17483d4ff4 100644 --- a/cpp/src/dual_simplex/node_presolve.hpp +++ b/cpp/src/dual_simplex/node_presolve.hpp @@ -28,10 +28,11 @@ class node_presolver_t { node_presolver_t(const lp_problem_t& problem, const std::vector& row_sense, const csc_matrix_t& Arow, - const std::vector& var_types, - const simplex_solver_settings_t& settings); + const std::vector& var_types); - bool bound_strengthening(std::vector& lower_bounds, std::vector& upper_bounds); + bool bound_strengthening(std::vector& lower_bounds, + std::vector& upper_bounds, + const simplex_solver_settings_t& settings); std::vector bounds_changed; @@ -39,7 +40,6 @@ class node_presolver_t { const csc_matrix_t& A; const csc_matrix_t& Arow; const std::vector& var_types; - const simplex_solver_settings_t& settings; std::vector delta_min_activity; std::vector delta_max_activity; From b5fc52c1c760ce8cd87a5c9b2f0e9e232da4ad2c Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 23 Oct 2025 19:24:47 +0200 Subject: [PATCH 15/51] added asserts --- cpp/src/dual_simplex/branch_and_bound.cpp | 7 +++++++ cpp/src/dual_simplex/branch_and_bound.hpp | 6 ++++-- cpp/src/dual_simplex/presolve.cpp | 12 ++---------- 3 files changed, 13 insertions(+), 12 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index dd8efd562c..88c218a1a9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -686,6 +687,12 @@ void branch_and_bound_t::set_variable_bounds(mip_node_t* nod const std::vector& root_upper, bool recompute) { + assert(bounds_changed.size() == original_lp_.lower.size()); + assert(lower.size() == original_lp_.lower.size()); + assert(upper.size() == original_lp_.upper.size()); + assert(root_lower.size() == original_lp_.lower.size()); + assert(root_upper.size() == original_lp_.upper.size()); + // Reset the bound_changed markers std::fill(bounds_changed.begin(), bounds_changed.end(), false); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index d92a62dc91..8d43300739 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -19,7 +19,6 @@ #include #include -#include #include #include #include @@ -52,6 +51,9 @@ enum class mip_exploration_status_t { COMPLETED = 5, // The solver finished exploring the tree }; +template +class node_presolver_t; + template void upper_bound_callback(f_t upper_bound); @@ -64,7 +66,7 @@ struct diving_root_t { diving_root_t(mip_node_t&& node, const std::vector& lower, const std::vector& upper) - : node(std::move(node)), upper(upper), lower(lower) + : node(std::move(node)), lower(lower), upper(upper) { } diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index a84d64b369..1682e7a9c4 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -576,15 +576,7 @@ void convert_user_problem(const user_problem_t& user_problem, convert_greater_to_less(user_problem, row_sense, problem, greater_rows, less_rows); } - // At this point the problem representation is in the form: A*x {<=, =} b - // This is the time to run bound strengthening - constexpr bool run_bound_strengthening = false; - if constexpr (run_bound_strengthening) { - settings.log.printf("Running bound strengthening\n"); - csc_matrix_t Arow(1, 1, 1); - problem.A.transpose(Arow); - bound_strengthening(row_sense, settings, problem, Arow); - } + // bounds strenghtning was moved to node_presolve.hpp settings.log.debug( "equality rows %d less rows %d columns %d\n", equal_rows, less_rows, problem.num_cols); @@ -605,7 +597,7 @@ void convert_user_problem(const user_problem_t& user_problem, } if (problem.upper[j] < inf) { num_upper_bounds++; - vars_with_upper_bounds.push_back(j); + vars_with_upper_bounds.push_back(j) } } From f7a38bcfa44e7484b6c5b5d2a2ebc9b1dd5667a3 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 23 Oct 2025 19:39:35 +0200 Subject: [PATCH 16/51] fix compilation error --- cpp/src/dual_simplex/presolve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 1682e7a9c4..f2439471af 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -597,7 +597,7 @@ void convert_user_problem(const user_problem_t& user_problem, } if (problem.upper[j] < inf) { num_upper_bounds++; - vars_with_upper_bounds.push_back(j) + vars_with_upper_bounds.push_back(j); } } From 649c8caa3054c07de8053ca9d48deec0df53b047 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Fri, 24 Oct 2025 11:00:30 +0200 Subject: [PATCH 17/51] Update cpp/src/dual_simplex/mip_node.hpp Co-authored-by: Rajesh Gandham --- cpp/src/dual_simplex/mip_node.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 3f2a3e4325..439832bec3 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -98,7 +98,7 @@ class mip_node_t { // Here we assume that we are traversing from the deepest node to the // root of the tree - void update_variable_bound(std::vector& lower, + void update_branched_variable_bounds(std::vector& lower, std::vector& upper, std::vector& bounds_changed) const { From 4afb71d8e58d8b05a415a4d40b2e30501873beeb Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 11:18:57 +0200 Subject: [PATCH 18/51] addresses reviewer's comments --- cpp/src/dual_simplex/branch_and_bound.cpp | 60 ++++++++-------- cpp/src/dual_simplex/branch_and_bound.hpp | 84 ++++------------------ cpp/src/dual_simplex/diving_queue.hpp | 88 +++++++++++++++++++++++ cpp/src/dual_simplex/node_presolve.cpp | 13 ++-- cpp/src/dual_simplex/node_presolve.hpp | 8 ++- cpp/src/dual_simplex/presolve.cpp | 2 +- cpp/src/dual_simplex/presolve.hpp | 2 + 7 files changed, 151 insertions(+), 106 deletions(-) create mode 100644 cpp/src/dual_simplex/diving_queue.hpp diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 88c218a1a9..34604b0096 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -201,6 +201,15 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound) } } +inline const char* thread_type_symbol(thread_type_t type) +{ + switch (type) { + case thread_type_t::EXPLORATION: return "B "; + case thread_type_t::DIVING: return "D "; + default: return "U"; + } +} + } // namespace template @@ -497,7 +506,7 @@ template void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, - char thread_type) + thread_type_t thread_type) { bool send_solution = false; i_t nodes_explored = stats_.nodes_explored; @@ -510,8 +519,8 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t lower_bound = get_lower_bound(); f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - thread_type, + settings_.log.printf("%s%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + thread_type_symbol(thread_type), nodes_explored, nodes_unexplored, obj, @@ -557,7 +566,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod search_tree_t& search_tree, lp_problem_t& leaf_problem, node_presolver_t& presolver, - char thread_type, + thread_type_t thread_type, logger_t& log) { const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; @@ -661,7 +670,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod return node_status_t::TIME_LIMIT; } else { - if (thread_type == 'B') { + if (thread_type == thread_type_t::EXPLORATION) { lower_bound_ceiling_.fetch_min(node_ptr->lower_bound); log.printf( "LP returned status %d on node %d. This indicates a numerical issue. The best bound is set " @@ -703,14 +712,13 @@ void branch_and_bound_t::set_variable_bounds(mip_node_t* nod node->get_variable_bounds(lower, upper, bounds_changed); } else { - node->update_variable_bound(lower, upper, bounds_changed); + node->update_branched_variable_bounds(lower, upper, bounds_changed); } } template void branch_and_bound_t::exploration_ramp_up(mip_node_t* node, search_tree_t* search_tree, - const csc_matrix_t& Arow, i_t initial_heap_size) { if (status_ != mip_exploration_status_t::RUNNING) { return; } @@ -723,7 +731,7 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_); + node_presolver_t presolver(leaf_problem, row_sense, var_types_); f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); @@ -781,8 +789,8 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod original_lp_.upper, true); - node_status_t node_status = - solve_node(node, *search_tree, leaf_problem, presolver, 'B', settings_.log); + node_status_t node_status = solve_node( + node, *search_tree, leaf_problem, presolver, thread_type_t::EXPLORATION, settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -794,10 +802,10 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod // If we haven't generated enough nodes to keep the threads busy, continue the ramp up phase if (stats_.nodes_unexplored < initial_heap_size) { #pragma omp task - exploration_ramp_up(node->get_down_child(), search_tree, Arow, initial_heap_size); + exploration_ramp_up(node->get_down_child(), search_tree, initial_heap_size); #pragma omp task - exploration_ramp_up(node->get_up_child(), search_tree, Arow, initial_heap_size); + exploration_ramp_up(node->get_up_child(), search_tree, initial_heap_size); } else { // We've generated enough nodes, push further nodes onto the heap @@ -889,8 +897,8 @@ void branch_and_bound_t::explore_subtree(i_t task_id, original_lp_.upper, recompute); - node_status_t node_status = - solve_node(node_ptr, search_tree, leaf_problem, presolver, 'B', settings_.log); + node_status_t node_status = solve_node( + node_ptr, search_tree, leaf_problem, presolver, thread_type_t::EXPLORATION, settings_.log); recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -933,9 +941,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, } template -void branch_and_bound_t::best_first_thread(i_t id, - search_tree_t& search_tree, - const csc_matrix_t& Arow) +void branch_and_bound_t::best_first_thread(i_t id, search_tree_t& search_tree) { f_t lower_bound = -inf; f_t upper_bound = inf; @@ -945,7 +951,7 @@ void branch_and_bound_t::best_first_thread(i_t id, // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_); + node_presolver_t presolver(leaf_problem, row_sense, var_types_); while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && @@ -994,7 +1000,7 @@ void branch_and_bound_t::best_first_thread(i_t id, } template -void branch_and_bound_t::diving_thread(const csc_matrix_t& Arow) +void branch_and_bound_t::diving_thread() { logger_t log; log.log = false; @@ -1002,7 +1008,7 @@ void branch_and_bound_t::diving_thread(const csc_matrix_t& A // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, Arow, var_types_); + node_presolver_t presolver(leaf_problem, row_sense, var_types_); while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { @@ -1043,7 +1049,7 @@ void branch_and_bound_t::diving_thread(const csc_matrix_t& A recompute); node_status_t node_status = - solve_node(node_ptr, subtree, leaf_problem, presolver, 'D', log); + solve_node(node_ptr, subtree, leaf_problem, presolver, thread_type_t::DIVING, log); if (node_status == node_status_t::TIME_LIMIT) { return; @@ -1219,9 +1225,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " "| Time |\n"); - csc_matrix_t Arow(1, 1, 1); - original_lp_.A.transpose(Arow); - stats_.nodes_explored = 0; stats_.nodes_unexplored = 2; stats_.nodes_since_last_log = 0; @@ -1230,6 +1233,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; + original_lp_.A.transpose(original_lp_.Arow); #pragma omp parallel num_threads(settings_.num_threads) { @@ -1240,10 +1244,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t initial_size = 2 * settings_.num_threads; #pragma omp task - exploration_ramp_up(down_child, &search_tree, Arow, initial_size); + exploration_ramp_up(down_child, &search_tree, initial_size); #pragma omp task - exploration_ramp_up(up_child, &search_tree, Arow, initial_size); + exploration_ramp_up(up_child, &search_tree, initial_size); } #pragma omp barrier @@ -1254,12 +1258,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut (active_subtrees_ > 0 || get_heap_size() > 0)) { for (i_t i = 0; i < settings_.num_bfs_threads; i++) { #pragma omp task - best_first_thread(i, search_tree, Arow); + best_first_thread(i, search_tree); } for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(Arow); + diving_thread(); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 8d43300739..0a1a556ec8 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -17,6 +17,7 @@ #pragma once +#include #include #include #include @@ -51,73 +52,21 @@ enum class mip_exploration_status_t { COMPLETED = 5, // The solver finished exploring the tree }; +// Indicate the search and variable selection algorithms used by the thread (See [1]). +// +// [1] T. Achterberg, “Constraint Integer Programming,” PhD, Technischen Universität Berlin, +// Berlin, 2007. doi: 10.14279/depositonce-1634. +enum class thread_type_t { + EXPLORATION = 0, // Best-First + Plunging. Pseudocost branching + Martin's criteria. + DIVING = 1, +}; + template class node_presolver_t; template void upper_bound_callback(f_t upper_bound); -template -struct diving_root_t { - mip_node_t node; - std::vector lower; - std::vector upper; - - diving_root_t(mip_node_t&& node, - const std::vector& lower, - const std::vector& upper) - : node(std::move(node)), lower(lower), upper(upper) - { - } - - friend bool operator>(const diving_root_t& a, const diving_root_t& b) - { - return a.node.lower_bound > b.node.lower_bound; - } -}; - -// A min-heap for storing the starting nodes for the dives. -// This has a maximum size of 256, such that the container -// will discard the least promising node if the queue is full. -template -class dive_queue_t { - private: - std::vector> buffer; - static constexpr i_t max_size_ = 256; - - public: - dive_queue_t() { buffer.reserve(max_size_); } - - void push(diving_root_t&& node) - { - buffer.push_back(std::move(node)); - std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); - if (buffer.size() > max_size()) { buffer.pop_back(); } - } - - void emplace(mip_node_t&& node, - const std::vector& lower, - const std::vector& upper) - { - buffer.emplace_back(std::move(node), lower, upper); - std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); - if (buffer.size() > max_size()) { buffer.pop_back(); } - } - - diving_root_t pop() - { - std::pop_heap(buffer.begin(), buffer.end(), std::greater<>()); - diving_root_t node = std::move(buffer.back()); - buffer.pop_back(); - return node; - } - - i_t size() const { return buffer.size(); } - constexpr i_t max_size() const { return max_size_; } - const diving_root_t& top() const { return buffer.front(); } - void clear() { buffer.clear(); } -}; - template class branch_and_bound_t { public: @@ -205,7 +154,7 @@ class branch_and_bound_t { // Queue for storing the promising node for performing dives. omp_mutex_t mutex_dive_queue_; - dive_queue_t dive_queue_; + diving_queue_t dive_queue_; i_t min_diving_queue_size_; // Global status of the solver. @@ -223,7 +172,7 @@ class branch_and_bound_t { void add_feasible_solution(f_t leaf_objective, const std::vector& leaf_solution, i_t leaf_depth, - char thread_type); + thread_type_t thread_type); // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); @@ -240,7 +189,6 @@ class branch_and_bound_t { // there is enough unexplored nodes. This is done recursively using OpenMP tasks. void exploration_ramp_up(mip_node_t* node, search_tree_t* search_tree, - const csc_matrix_t& Arow, i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. @@ -252,20 +200,18 @@ class branch_and_bound_t { // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. - void best_first_thread(i_t id, - search_tree_t& search_tree, - const csc_matrix_t& Arow); + void best_first_thread(i_t id, search_tree_t& search_tree); // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(const csc_matrix_t& Arow); + void diving_thread(); // Solve the LP relaxation of a leaf node and update the tree. node_status_t solve_node(mip_node_t* node_ptr, search_tree_t& search_tree, lp_problem_t& leaf_problem, node_presolver_t& presolver, - char thread_type, + thread_type_t thread_type, logger_t& log); // Sort the children based on the Martin's criteria. diff --git a/cpp/src/dual_simplex/diving_queue.hpp b/cpp/src/dual_simplex/diving_queue.hpp new file mode 100644 index 0000000000..d86e5c79ca --- /dev/null +++ b/cpp/src/dual_simplex/diving_queue.hpp @@ -0,0 +1,88 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#include +#include + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +struct diving_root_t { + mip_node_t node; + std::vector lower; + std::vector upper; + + diving_root_t(mip_node_t&& node, + const std::vector& lower, + const std::vector& upper) + : node(std::move(node)), lower(lower), upper(upper) + { + } + + friend bool operator>(const diving_root_t& a, const diving_root_t& b) + { + return a.node.lower_bound > b.node.lower_bound; + } +}; + +// A min-heap for storing the starting nodes for the dives. +// This has a maximum size of 256, such that the container +// will discard the least promising node if the queue is full. +template +class diving_queue_t { + private: + std::vector> buffer; + static constexpr i_t max_size_ = 256; + + public: + diving_queue_t() { buffer.reserve(max_size_); } + + void push(diving_root_t&& node) + { + buffer.push_back(std::move(node)); + std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); + if (buffer.size() > max_size()) { buffer.pop_back(); } + } + + void emplace(mip_node_t&& node, + const std::vector& lower, + const std::vector& upper) + { + buffer.emplace_back(std::move(node), lower, upper); + std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); + if (buffer.size() > max_size()) { buffer.pop_back(); } + } + + diving_root_t pop() + { + std::pop_heap(buffer.begin(), buffer.end(), std::greater<>()); + diving_root_t node = std::move(buffer.back()); + buffer.pop_back(); + return node; + } + + i_t size() const { return buffer.size(); } + constexpr i_t max_size() const { return max_size_; } + const diving_root_t& top() const { return buffer.front(); } + void clear() { buffer.clear(); } +}; + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/node_presolve.cpp index f05bfab5a5..a81cbad406 100644 --- a/cpp/src/dual_simplex/node_presolve.cpp +++ b/cpp/src/dual_simplex/node_presolve.cpp @@ -63,11 +63,10 @@ void print_bounds_stats(const std::vector& lower, template node_presolver_t::node_presolver_t(const lp_problem_t& problem, const std::vector& row_sense, - const csc_matrix_t& Arow, const std::vector& var_types) : bounds_changed(problem.num_cols, false), A(problem.A), - Arow(Arow), + Arow(problem.Arow), var_types(var_types), delta_min_activity(problem.num_rows), delta_max_activity(problem.num_rows), @@ -104,9 +103,9 @@ bool node_presolver_t::bound_strengthening( const i_t m = A.m; const i_t n = A.n; - std::vector constraint_changed(m, false); - std::vector variable_changed(n, false); - std::vector constraint_changed_next(m, false); + constraint_changed.assign(m, false); + variable_changed.assign(n, false); + constraint_changed_next.assign(m, false); for (i_t i = 0; i < bounds_changed.size(); ++i) { if (bounds_changed[i]) { @@ -119,8 +118,8 @@ bool node_presolver_t::bound_strengthening( } } - std::vector lower = lower_bounds; - std::vector upper = upper_bounds; + lower = lower_bounds; + upper = upper_bounds; print_bounds_stats(lower, upper, settings, "Initial bounds"); i_t iter = 0; diff --git a/cpp/src/dual_simplex/node_presolve.hpp b/cpp/src/dual_simplex/node_presolve.hpp index 17483d4ff4..1bb7fb899f 100644 --- a/cpp/src/dual_simplex/node_presolve.hpp +++ b/cpp/src/dual_simplex/node_presolve.hpp @@ -27,7 +27,6 @@ class node_presolver_t { // For pure LP bounds strengthening, var_types should be defaulted (i.e. left empty) node_presolver_t(const lp_problem_t& problem, const std::vector& row_sense, - const csc_matrix_t& Arow, const std::vector& var_types); bool bound_strengthening(std::vector& lower_bounds, @@ -41,6 +40,13 @@ class node_presolver_t { const csc_matrix_t& Arow; const std::vector& var_types; + std::vector constraint_changed; + std::vector variable_changed; + std::vector constraint_changed_next; + + std::vector lower; + std::vector upper; + std::vector delta_min_activity; std::vector delta_max_activity; std::vector constraint_lb; diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index f2439471af..29e63d8676 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -576,7 +576,7 @@ void convert_user_problem(const user_problem_t& user_problem, convert_greater_to_less(user_problem, row_sense, problem, greater_rows, less_rows); } - // bounds strenghtning was moved to node_presolve.hpp + // bounds strengthening was moved to node_presolve.hpp settings.log.debug( "equality rows %d less rows %d columns %d\n", equal_rows, less_rows, problem.num_cols); diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 538ca5dffe..b0d7b0aa23 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -33,6 +33,7 @@ struct lp_problem_t { num_cols(n), objective(n), A(m, n, nz), + Arow(1, 1, 1), rhs(m), lower(n), upper(n), @@ -44,6 +45,7 @@ struct lp_problem_t { i_t num_cols; std::vector objective; csc_matrix_t A; + csc_matrix_t Arow; std::vector rhs; std::vector lower; std::vector upper; From 1c92230eb7d927c16bcad6795a3d26fb771a43e5 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 11:29:00 +0200 Subject: [PATCH 19/51] fix style --- cpp/src/dual_simplex/mip_node.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 439832bec3..39544180fd 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -99,8 +99,8 @@ class mip_node_t { // Here we assume that we are traversing from the deepest node to the // root of the tree void update_branched_variable_bounds(std::vector& lower, - std::vector& upper, - std::vector& bounds_changed) const + std::vector& upper, + std::vector& bounds_changed) const { assert(branch_var >= 0); assert(lower.size() > branch_var); From bbd410b11ea2d503e82cf434a843f528f53989a9 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 12:04:08 +0200 Subject: [PATCH 20/51] small refactor. added missing recompute flag to diving. --- cpp/src/dual_simplex/branch_and_bound.cpp | 109 ++++++++++------------ cpp/src/dual_simplex/branch_and_bound.hpp | 11 +-- 2 files changed, 50 insertions(+), 70 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 34604b0096..61aebdacfb 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -567,6 +567,9 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod lp_problem_t& leaf_problem, node_presolver_t& presolver, thread_type_t thread_type, + bool recompute, + const std::vector& root_lower, + const std::vector& root_upper, logger_t& log) { const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; @@ -582,6 +585,20 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); + // Reset the bound_changed markers + std::fill(presolver.bounds_changed.begin(), presolver.bounds_changed.end(), false); + + // Set the correct bounds for the leaf problem + if (recompute) { + leaf_problem.lower = root_lower; + leaf_problem.upper = root_upper; + node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, presolver.bounds_changed); + + } else { + node_ptr->update_branched_variable_bounds( + leaf_problem.lower, leaf_problem.upper, presolver.bounds_changed); + } + bool feasible = presolver.bound_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); @@ -687,35 +704,6 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod } } -template -void branch_and_bound_t::set_variable_bounds(mip_node_t* node, - std::vector& lower, - std::vector& upper, - std::vector& bounds_changed, - const std::vector& root_lower, - const std::vector& root_upper, - bool recompute) -{ - assert(bounds_changed.size() == original_lp_.lower.size()); - assert(lower.size() == original_lp_.lower.size()); - assert(upper.size() == original_lp_.upper.size()); - assert(root_lower.size() == original_lp_.lower.size()); - assert(root_upper.size() == original_lp_.upper.size()); - - // Reset the bound_changed markers - std::fill(bounds_changed.begin(), bounds_changed.end(), false); - - // Recompute the bounds - if (recompute) { - lower = root_lower; - upper = root_upper; - node->get_variable_bounds(lower, upper, bounds_changed); - - } else { - node->update_branched_variable_bounds(lower, upper, bounds_changed); - } -} - template void branch_and_bound_t::exploration_ramp_up(mip_node_t* node, search_tree_t* search_tree, @@ -780,17 +768,15 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod return; } - // Set the correct bounds for the leaf problem - set_variable_bounds(node, - leaf_problem.lower, - leaf_problem.upper, - presolver.bounds_changed, - original_lp_.lower, - original_lp_.upper, - true); - - node_status_t node_status = solve_node( - node, *search_tree, leaf_problem, presolver, thread_type_t::EXPLORATION, settings_.log); + node_status_t node_status = solve_node(node, + *search_tree, + leaf_problem, + presolver, + thread_type_t::EXPLORATION, + true, + original_lp_.lower, + original_lp_.upper, + settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -888,17 +874,16 @@ void branch_and_bound_t::explore_subtree(i_t task_id, return; } - // Set the correct bounds for the leaf problem - set_variable_bounds(node_ptr, - leaf_problem.lower, - leaf_problem.upper, - presolver.bounds_changed, - original_lp_.lower, - original_lp_.upper, - recompute); - - node_status_t node_status = solve_node( - node_ptr, search_tree, leaf_problem, presolver, thread_type_t::EXPLORATION, settings_.log); + node_status_t node_status = solve_node(node_ptr, + search_tree, + leaf_problem, + presolver, + thread_type_t::EXPLORATION, + recompute, + original_lp_.lower, + original_lp_.upper, + settings_.log); + recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { @@ -1039,17 +1024,17 @@ void branch_and_bound_t::diving_thread() if (toc(stats_.start_time) > settings_.time_limit) { return; } - // Set the correct bounds for the leaf problem - set_variable_bounds(node_ptr, - leaf_problem.lower, - leaf_problem.upper, - presolver.bounds_changed, - start_node->lower, - start_node->upper, - recompute); - - node_status_t node_status = - solve_node(node_ptr, subtree, leaf_problem, presolver, thread_type_t::DIVING, log); + node_status_t node_status = solve_node(node_ptr, + subtree, + leaf_problem, + presolver, + thread_type_t::DIVING, + recompute, + start_node->lower, + start_node->upper, + log); + + recompute = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { return; diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 0a1a556ec8..fd4721b98c 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -177,14 +177,6 @@ class branch_and_bound_t { // Repairs low-quality solutions from the heuristics, if it is applicable. void repair_heuristic_solutions(); - void set_variable_bounds(mip_node_t* node, - std::vector& lower, - std::vector& upper, - std::vector& bounds_changed, - const std::vector& root_lower, - const std::vector& root_upper, - bool recompute); - // Ramp-up phase of the solver, where we greedily expand the tree until // there is enough unexplored nodes. This is done recursively using OpenMP tasks. void exploration_ramp_up(mip_node_t* node, @@ -212,6 +204,9 @@ class branch_and_bound_t { lp_problem_t& leaf_problem, node_presolver_t& presolver, thread_type_t thread_type, + bool recompute, + const std::vector& root_lower, + const std::vector& root_upper, logger_t& log); // Sort the children based on the Martin's criteria. From 8c420c4abd70f2546ff92f0ff06d37e2f07165c4 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 13:05:46 +0200 Subject: [PATCH 21/51] small fixes --- cpp/src/dual_simplex/mip_node.hpp | 2 +- cpp/src/dual_simplex/node_presolve.cpp | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 39544180fd..bdd97ad0b5 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -91,7 +91,7 @@ class mip_node_t { mip_node_t* parent_ptr = parent; while (parent_ptr != nullptr && parent_ptr->node_id != 0) { - parent_ptr->update_variable_bound(lower, upper, bounds_changed); + parent_ptr->update_branched_variable_bounds(lower, upper, bounds_changed); parent_ptr = parent_ptr->parent; } } diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/node_presolve.cpp index a81cbad406..fb59c56932 100644 --- a/cpp/src/dual_simplex/node_presolve.cpp +++ b/cpp/src/dual_simplex/node_presolve.cpp @@ -17,6 +17,9 @@ #include +#include +#include + namespace cuopt::linear_programming::dual_simplex { template @@ -206,8 +209,8 @@ bool node_presolver_t::bound_strengthening( new_ub = std::floor(new_ub + settings.integer_tol); } - bool lb_updated = abs(new_lb - old_lb) > 1e3 * settings.primal_tol; - bool ub_updated = abs(new_ub - old_ub) > 1e3 * settings.primal_tol; + bool lb_updated = std::abs(new_lb - old_lb) > 1e3 * settings.primal_tol; + bool ub_updated = std::abs(new_ub - old_ub) > 1e3 * settings.primal_tol; new_lb = std::max(new_lb, lower_bounds[k]); new_ub = std::min(new_ub, upper_bounds[k]); From 5bd70e1811103b604fdc1f7a5c5963931f7ebe7d Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 13:11:10 +0200 Subject: [PATCH 22/51] small fix --- cpp/src/dual_simplex/mip_node.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index bdd97ad0b5..cfb491caf0 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -87,7 +87,7 @@ class mip_node_t { std::vector& upper, std::vector& bounds_changed) const { - update_variable_bound(lower, upper, bounds_changed); + update_branched_variable_bounds(lower, upper, bounds_changed); mip_node_t* parent_ptr = parent; while (parent_ptr != nullptr && parent_ptr->node_id != 0) { From 56747e5f2f3d30829c1fe2231b2c5586c7e888b6 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 13:14:40 +0200 Subject: [PATCH 23/51] fix logs --- cpp/src/dual_simplex/branch_and_bound.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 61aebdacfb..d435e21718 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -204,8 +204,8 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound) inline const char* thread_type_symbol(thread_type_t type) { switch (type) { - case thread_type_t::EXPLORATION: return "B "; - case thread_type_t::DIVING: return "D "; + case thread_type_t::EXPLORATION: return "B"; + case thread_type_t::DIVING: return "D"; default: return "U"; } } @@ -519,7 +519,7 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t lower_bound = get_lower_bound(); f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%s%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", thread_type_symbol(thread_type), nodes_explored, nodes_unexplored, From c7facc6e2e77be0a97ea19ba3b4822748d77671e Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 24 Oct 2025 13:19:03 +0200 Subject: [PATCH 24/51] fixed log --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index d435e21718..06e6389388 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -519,7 +519,7 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t lower_bound = get_lower_bound(); f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf("%s%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", thread_type_symbol(thread_type), nodes_explored, nodes_unexplored, From 8b294be5dd639f5b812f48fa2a1d37fc0b765fe3 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 27 Oct 2025 17:25:22 +0100 Subject: [PATCH 25/51] fixed starting bounds for diving nodes. replaced rounding direction with an enum. --- cpp/src/dual_simplex/branch_and_bound.cpp | 25 ++++++++---- cpp/src/dual_simplex/diving_queue.hpp | 18 ++++---- cpp/src/dual_simplex/mip_node.hpp | 50 ++++++++++++++--------- 3 files changed, 54 insertions(+), 39 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 06e6389388..275d1c358d 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -620,7 +620,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod leaf_edge_norms); if (lp_status == dual::status_t::NUMERICAL) { - log.printf("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); + log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); lp_status_t second_status = solve_linear_program_advanced( leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); lp_status = convert_lp_status_to_dual_status(second_status); @@ -906,8 +906,12 @@ void branch_and_bound_t::explore_subtree(i_t task_id, // This lead to a SIGSEGV. Although, in this case, it // would be better if we discard the node instead. if (get_heap_size() > settings_.num_bfs_threads) { + std::vector lower = original_lp_.lower; + std::vector upper = original_lp_.upper; + node->get_variable_bounds(lower, upper, presolver.bounds_changed); + mutex_dive_queue_.lock(); - dive_queue_.emplace(node->detach_copy(), leaf_problem.lower, leaf_problem.upper); + dive_queue_.emplace(node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } @@ -1051,10 +1055,15 @@ void branch_and_bound_t::diving_thread() // lowest possible point and move to the queue, so it can // be picked by another thread. if (dive_queue_.size() < min_diving_queue_size_) { - mutex_dive_queue_.lock(); mip_node_t* new_node = stack.back(); stack.pop_back(); - dive_queue_.emplace(new_node->detach_copy(), leaf_problem.lower, leaf_problem.upper); + + std::vector lower = start_node->lower; + std::vector upper = start_node->upper; + new_node->get_variable_bounds(lower, upper, presolver.bounds_changed); + + mutex_dive_queue_.lock(); + dive_queue_.emplace(new_node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } } @@ -1200,11 +1209,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_, log); - settings_.log.printf( - "Exploring the B&B tree using %d best-first threads and %d diving threads (%d threads)\n", - settings_.num_bfs_threads, - settings_.num_diving_threads, - settings_.num_threads); + settings_.log.printf("Exploring the B&B tree using %d best-first threads and %d diving threads\n", + settings_.num_bfs_threads, + settings_.num_diving_threads); settings_.log.printf( " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " diff --git a/cpp/src/dual_simplex/diving_queue.hpp b/cpp/src/dual_simplex/diving_queue.hpp index d86e5c79ca..4b4610824c 100644 --- a/cpp/src/dual_simplex/diving_queue.hpp +++ b/cpp/src/dual_simplex/diving_queue.hpp @@ -30,10 +30,8 @@ struct diving_root_t { std::vector lower; std::vector upper; - diving_root_t(mip_node_t&& node, - const std::vector& lower, - const std::vector& upper) - : node(std::move(node)), lower(lower), upper(upper) + diving_root_t(mip_node_t&& node, std::vector&& lower, std::vector&& upper) + : node(std::move(node)), lower(std::move(lower)), upper(std::move(upper)) { } @@ -50,7 +48,7 @@ template class diving_queue_t { private: std::vector> buffer; - static constexpr i_t max_size_ = 256; + static constexpr i_t max_size_ = INT16_MAX; public: diving_queue_t() { buffer.reserve(max_size_); } @@ -59,16 +57,14 @@ class diving_queue_t { { buffer.push_back(std::move(node)); std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); - if (buffer.size() > max_size()) { buffer.pop_back(); } + if (buffer.size() > max_size() - 1) { buffer.pop_back(); } } - void emplace(mip_node_t&& node, - const std::vector& lower, - const std::vector& upper) + void emplace(mip_node_t&& node, std::vector&& lower, std::vector&& upper) { - buffer.emplace_back(std::move(node), lower, upper); + buffer.emplace_back(std::move(node), std::move(lower), std::move(upper)); std::push_heap(buffer.begin(), buffer.end(), std::greater<>()); - if (buffer.size() > max_size()) { buffer.pop_back(); } + if (buffer.size() > max_size() - 1) { buffer.pop_back(); } } diving_root_t pop() diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index cfb491caf0..1c51c09d12 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -38,6 +38,8 @@ enum class node_status_t : int { TIME_LIMIT = 6 // Time out during the LP relaxation }; +enum class round_dir_t { NONE = -1, DOWN = 0, UP = 1 }; + bool inactive_status(node_status_t status); template @@ -50,7 +52,7 @@ class mip_node_t { parent(nullptr), node_id(0), branch_var(-1), - branch_dir(-1), + branch_dir(round_dir_t::NONE), vstatus(basis) { children[0] = nullptr; @@ -61,7 +63,7 @@ class mip_node_t { mip_node_t* parent_node, i_t node_num, i_t branch_variable, - i_t branch_direction, + round_dir_t branch_direction, f_t branch_var_value, const std::vector& basis) : status(node_status_t::ACTIVE), @@ -75,12 +77,12 @@ class mip_node_t { vstatus(basis) { - branch_var_lower = - branch_direction == 0 ? problem.lower[branch_var] : std::ceil(branch_var_value); - branch_var_upper = - branch_direction == 0 ? std::floor(branch_var_value) : problem.upper[branch_var]; - children[0] = nullptr; - children[1] = nullptr; + branch_var_lower = branch_direction == round_dir_t::DOWN ? problem.lower[branch_var] + : std::ceil(branch_var_value); + branch_var_upper = branch_direction == round_dir_t::DOWN ? std::floor(branch_var_value) + : problem.upper[branch_var]; + children[0] = nullptr; + children[1] = nullptr; } void get_variable_bounds(std::vector& lower, @@ -105,6 +107,7 @@ class mip_node_t { assert(branch_var >= 0); assert(lower.size() > branch_var); assert(upper.size() > branch_var); + assert(bounds_changed.size() > branch_var); // If the bounds have already been updated on another node, // skip this node as it contains a less tight bounds. @@ -231,7 +234,7 @@ class mip_node_t { i_t depth; i_t node_id; i_t branch_var; - i_t branch_dir; + round_dir_t branch_dir; f_t branch_var_lower; f_t branch_var_upper; f_t fractional_val; @@ -296,17 +299,26 @@ class search_tree_t { { i_t id = num_nodes.fetch_add(2); - // down child - auto down_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, 0, fractional_val, parent_vstatus); - - graphviz_edge(log, parent_node, down_child.get(), branch_var, 0, std::floor(fractional_val)); + auto down_child = std::make_unique>(original_lp, + parent_node, + ++id, + branch_var, + round_dir_t::DOWN, + fractional_val, + parent_vstatus); + + graphviz_edge(log, + parent_node, + down_child.get(), + branch_var, + round_dir_t::DOWN, + std::floor(fractional_val)); - // up child auto up_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, 1, fractional_val, parent_vstatus); + original_lp, parent_node, ++id, branch_var, round_dir_t::UP, fractional_val, parent_vstatus); - graphviz_edge(log, parent_node, up_child.get(), branch_var, 1, std::ceil(fractional_val)); + graphviz_edge( + log, parent_node, up_child.get(), branch_var, round_dir_t::UP, std::ceil(fractional_val)); assert(parent_vstatus.size() == original_lp.num_cols); parent_node->add_children(std::move(down_child), @@ -327,7 +339,7 @@ class search_tree_t { const mip_node_t* origin_ptr, const mip_node_t* dest_ptr, const i_t branch_var, - const i_t branch_dir, + round_dir_t branch_dir, const f_t bound) { if (write_graphviz) { @@ -335,7 +347,7 @@ class search_tree_t { origin_ptr->node_id, dest_ptr->node_id, branch_var, - branch_dir == 0 ? "<=" : ">=", + branch_dir == round_dir_t::DOWN ? "<=" : ">=", bound); } } From 27f9264eba1b86246a123fb43cd80c1e7795aa00 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 28 Oct 2025 10:35:13 +0100 Subject: [PATCH 26/51] fix missing enum --- cpp/src/dual_simplex/pseudo_costs.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index 4bd9590e16..1a10dc73e2 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -210,10 +210,10 @@ void pseudo_costs_t::update_pseudo_costs(mip_node_t* node_pt { mutex.lock(); const f_t change_in_obj = leaf_objective - node_ptr->lower_bound; - const f_t frac = node_ptr->branch_dir == 0 + const f_t frac = node_ptr->branch_dir == round_dir_t::DOWN ? node_ptr->fractional_val - std::floor(node_ptr->fractional_val) : std::ceil(node_ptr->fractional_val) - node_ptr->fractional_val; - if (node_ptr->branch_dir == 0) { + if (node_ptr->branch_dir == round_dir_t::DOWN) { pseudo_cost_sum_down[node_ptr->branch_var] += change_in_obj / frac; pseudo_cost_num_down[node_ptr->branch_var]++; } else { From bf142d5c9b2a5e94a01c8e57ab9eb6ff67f3abe9 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 28 Oct 2025 12:02:35 +0100 Subject: [PATCH 27/51] fixed typo --- cpp/src/dual_simplex/diving_queue.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/diving_queue.hpp b/cpp/src/dual_simplex/diving_queue.hpp index 4b4610824c..57a8bacf19 100644 --- a/cpp/src/dual_simplex/diving_queue.hpp +++ b/cpp/src/dual_simplex/diving_queue.hpp @@ -18,6 +18,7 @@ #pragma once #include +#include #include #include @@ -42,7 +43,7 @@ struct diving_root_t { }; // A min-heap for storing the starting nodes for the dives. -// This has a maximum size of 256, such that the container +// This has a maximum size of INT16_MAX, such that the container // will discard the least promising node if the queue is full. template class diving_queue_t { From cc89861cf8e73bb6bf65f98d7d2d62524405bc13 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 31 Oct 2025 18:24:11 +0100 Subject: [PATCH 28/51] renamed files and classes --- cpp/src/dual_simplex/CMakeLists.txt | 2 +- ..._presolve.cpp => bounds_strengthening.cpp} | 22 +++---- ..._presolve.hpp => bounds_strengthening.hpp} | 17 +++--- cpp/src/dual_simplex/branch_and_bound.cpp | 57 ++++++++++--------- cpp/src/dual_simplex/branch_and_bound.hpp | 13 +++-- cpp/src/dual_simplex/presolve.hpp | 2 - 6 files changed, 61 insertions(+), 52 deletions(-) rename cpp/src/dual_simplex/{node_presolve.cpp => bounds_strengthening.cpp} (94%) rename cpp/src/dual_simplex/{node_presolve.hpp => bounds_strengthening.hpp} (73%) diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index c52121a631..fb063efbe7 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -27,7 +27,7 @@ set(DUAL_SIMPLEX_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/phase1.cpp ${CMAKE_CURRENT_SOURCE_DIR}/phase2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/presolve.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/node_presolve.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/bounds_strengthening.cpp ${CMAKE_CURRENT_SOURCE_DIR}/primal.cpp ${CMAKE_CURRENT_SOURCE_DIR}/pseudo_costs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/right_looking_lu.cpp diff --git a/cpp/src/dual_simplex/node_presolve.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp similarity index 94% rename from cpp/src/dual_simplex/node_presolve.cpp rename to cpp/src/dual_simplex/bounds_strengthening.cpp index fb59c56932..fc1adb771b 100644 --- a/cpp/src/dual_simplex/node_presolve.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -15,7 +15,7 @@ * limitations under the License. */ -#include +#include #include #include @@ -64,12 +64,14 @@ void print_bounds_stats(const std::vector& lower, } template -node_presolver_t::node_presolver_t(const lp_problem_t& problem, - const std::vector& row_sense, - const std::vector& var_types) +bounds_strengthening_t::bounds_strengthening_t( + const lp_problem_t& problem, + const csr_matrix_t& Arow, + const std::vector& row_sense, + const std::vector& var_types) : bounds_changed(problem.num_cols, false), A(problem.A), - Arow(problem.Arow), + Arow(Arow), var_types(var_types), delta_min_activity(problem.num_rows), delta_max_activity(problem.num_rows), @@ -98,7 +100,7 @@ node_presolver_t::node_presolver_t(const lp_problem_t& probl } template -bool node_presolver_t::bound_strengthening( +bool bounds_strengthening_t::bounds_strengthening( std::vector& lower_bounds, std::vector& upper_bounds, const simplex_solver_settings_t& settings) @@ -130,13 +132,13 @@ bool node_presolver_t::bound_strengthening( while (iter < iter_limit) { for (i_t i = 0; i < m; ++i) { if (!constraint_changed[i]) { continue; } - const i_t row_start = Arow.col_start[i]; - const i_t row_end = Arow.col_start[i + 1]; + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; f_t min_a = 0.0; f_t max_a = 0.0; for (i_t p = row_start; p < row_end; ++p) { - const i_t j = Arow.i[p]; + const i_t j = Arow.j[p]; const f_t a_ij = Arow.x[p]; variable_changed[j] = true; @@ -289,7 +291,7 @@ bool node_presolver_t::bound_strengthening( } #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE -template class node_presolver_t; +template class bounds_strengthening_t; #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/node_presolve.hpp b/cpp/src/dual_simplex/bounds_strengthening.hpp similarity index 73% rename from cpp/src/dual_simplex/node_presolve.hpp rename to cpp/src/dual_simplex/bounds_strengthening.hpp index 1bb7fb899f..28df35d0dc 100644 --- a/cpp/src/dual_simplex/node_presolve.hpp +++ b/cpp/src/dual_simplex/bounds_strengthening.hpp @@ -22,22 +22,23 @@ namespace cuopt::linear_programming::dual_simplex { template -class node_presolver_t { +class bounds_strengthening_t { public: // For pure LP bounds strengthening, var_types should be defaulted (i.e. left empty) - node_presolver_t(const lp_problem_t& problem, - const std::vector& row_sense, - const std::vector& var_types); + bounds_strengthening_t(const lp_problem_t& problem, + const csr_matrix_t& Arow, + const std::vector& row_sense, + const std::vector& var_types); - bool bound_strengthening(std::vector& lower_bounds, - std::vector& upper_bounds, - const simplex_solver_settings_t& settings); + bool bounds_strengthening(std::vector& lower_bounds, + std::vector& upper_bounds, + const simplex_solver_settings_t& settings); std::vector bounds_changed; private: const csc_matrix_t& A; - const csc_matrix_t& Arow; + const csr_matrix_t& Arow; const std::vector& var_types; std::vector constraint_changed; diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 275d1c358d..20a3e6a994 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -17,11 +17,11 @@ #include #include +#include #include #include #include #include -#include #include #include #include @@ -565,9 +565,9 @@ template node_status_t branch_and_bound_t::solve_node(mip_node_t* node_ptr, search_tree_t& search_tree, lp_problem_t& leaf_problem, - node_presolver_t& presolver, + bounds_strengthening_t& presolver, thread_type_t thread_type, - bool recompute, + bool recompute_bounds, const std::vector& root_lower, const std::vector& root_upper, logger_t& log) @@ -589,7 +589,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod std::fill(presolver.bounds_changed.begin(), presolver.bounds_changed.end(), false); // Set the correct bounds for the leaf problem - if (recompute) { + if (recompute_bounds) { leaf_problem.lower = root_lower; leaf_problem.upper = root_upper; node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, presolver.bounds_changed); @@ -600,7 +600,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod } bool feasible = - presolver.bound_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); + presolver.bounds_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -707,6 +707,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod template void branch_and_bound_t::exploration_ramp_up(mip_node_t* node, search_tree_t* search_tree, + const csr_matrix_t& Arow, i_t initial_heap_size) { if (status_ != mip_exploration_status_t::RUNNING) { return; } @@ -719,7 +720,7 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, var_types_); + bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); @@ -788,10 +789,10 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod // If we haven't generated enough nodes to keep the threads busy, continue the ramp up phase if (stats_.nodes_unexplored < initial_heap_size) { #pragma omp task - exploration_ramp_up(node->get_down_child(), search_tree, initial_heap_size); + exploration_ramp_up(node->get_down_child(), search_tree, Arow, initial_heap_size); #pragma omp task - exploration_ramp_up(node->get_up_child(), search_tree, initial_heap_size); + exploration_ramp_up(node->get_up_child(), search_tree, Arow, initial_heap_size); } else { // We've generated enough nodes, push further nodes onto the heap @@ -808,9 +809,9 @@ void branch_and_bound_t::explore_subtree(i_t task_id, mip_node_t* start_node, search_tree_t& search_tree, lp_problem_t& leaf_problem, - node_presolver_t& presolver) + bounds_strengthening_t& presolver) { - bool recompute = true; + bool recompute_bounds = true; std::deque*> stack; stack.push_front(start_node); @@ -840,7 +841,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); search_tree.update(node_ptr, node_status_t::FATHOMED); - recompute = true; + recompute_bounds = true; continue; } @@ -879,12 +880,12 @@ void branch_and_bound_t::explore_subtree(i_t task_id, leaf_problem, presolver, thread_type_t::EXPLORATION, - recompute, + recompute_bounds, original_lp_.lower, original_lp_.upper, settings_.log); - recompute = node_status != node_status_t::HAS_CHILDREN; + recompute_bounds = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -930,7 +931,9 @@ void branch_and_bound_t::explore_subtree(i_t task_id, } template -void branch_and_bound_t::best_first_thread(i_t id, search_tree_t& search_tree) +void branch_and_bound_t::best_first_thread(i_t id, + search_tree_t& search_tree, + const csr_matrix_t& Arow) { f_t lower_bound = -inf; f_t upper_bound = inf; @@ -940,7 +943,7 @@ void branch_and_bound_t::best_first_thread(i_t id, search_tree_t leaf_problem = original_lp_; std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, var_types_); + bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && @@ -989,7 +992,7 @@ void branch_and_bound_t::best_first_thread(i_t id, search_tree_t -void branch_and_bound_t::diving_thread() +void branch_and_bound_t::diving_thread(const csr_matrix_t& Arow) { logger_t log; log.log = false; @@ -997,7 +1000,7 @@ void branch_and_bound_t::diving_thread() // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - node_presolver_t presolver(leaf_problem, row_sense, var_types_); + bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); while (status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { @@ -1010,7 +1013,7 @@ void branch_and_bound_t::diving_thread() if (start_node.has_value()) { if (get_upper_bound() < start_node->node.lower_bound) { continue; } - bool recompute = true; + bool recompute_bounds = true; search_tree_t subtree(std::move(start_node->node)); std::deque*> stack; stack.push_front(&subtree.root); @@ -1022,7 +1025,7 @@ void branch_and_bound_t::diving_thread() f_t rel_gap = user_relative_gap(original_lp_, upper_bound, node_ptr->lower_bound); if (node_ptr->lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { - recompute = true; + recompute_bounds = true; continue; } @@ -1033,12 +1036,12 @@ void branch_and_bound_t::diving_thread() leaf_problem, presolver, thread_type_t::DIVING, - recompute, + recompute_bounds, start_node->lower, start_node->upper, log); - recompute = node_status != node_status_t::HAS_CHILDREN; + recompute_bounds = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { return; @@ -1225,7 +1228,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut min_diving_queue_size_ = 4 * settings_.num_diving_threads; status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; - original_lp_.A.transpose(original_lp_.Arow); + + csr_matrix_t Arow(1, 1, 0); + original_lp_.A.to_compressed_row(Arow); #pragma omp parallel num_threads(settings_.num_threads) { @@ -1236,10 +1241,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t initial_size = 2 * settings_.num_threads; #pragma omp task - exploration_ramp_up(down_child, &search_tree, initial_size); + exploration_ramp_up(down_child, &search_tree, Arow, initial_size); #pragma omp task - exploration_ramp_up(up_child, &search_tree, initial_size); + exploration_ramp_up(up_child, &search_tree, Arow, initial_size); } #pragma omp barrier @@ -1250,12 +1255,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut (active_subtrees_ > 0 || get_heap_size() > 0)) { for (i_t i = 0; i < settings_.num_bfs_threads; i++) { #pragma omp task - best_first_thread(i, search_tree); + best_first_thread(i, search_tree, Arow); } for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(); + diving_thread(Arow); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index fd4721b98c..00da2f96df 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -62,7 +62,7 @@ enum class thread_type_t { }; template -class node_presolver_t; +class bounds_strengthening_t; template void upper_bound_callback(f_t upper_bound); @@ -181,6 +181,7 @@ class branch_and_bound_t { // there is enough unexplored nodes. This is done recursively using OpenMP tasks. void exploration_ramp_up(mip_node_t* node, search_tree_t* search_tree, + const csr_matrix_t& Arow, i_t initial_heap_size); // Explore the search tree using the best-first search with plunging strategy. @@ -188,21 +189,23 @@ class branch_and_bound_t { mip_node_t* start_node, search_tree_t& search_tree, lp_problem_t& leaf_problem, - node_presolver_t& presolver); + bounds_strengthening_t& presolver); // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. - void best_first_thread(i_t id, search_tree_t& search_tree); + void best_first_thread(i_t id, + search_tree_t& search_tree, + const csr_matrix_t& Arow); // Each diving thread pops the first node from the dive queue and then performs // a deep dive into the subtree determined by the node. - void diving_thread(); + void diving_thread(const csr_matrix_t& Arow); // Solve the LP relaxation of a leaf node and update the tree. node_status_t solve_node(mip_node_t* node_ptr, search_tree_t& search_tree, lp_problem_t& leaf_problem, - node_presolver_t& presolver, + bounds_strengthening_t& presolver, thread_type_t thread_type, bool recompute, const std::vector& root_lower, diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index b0d7b0aa23..538ca5dffe 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -33,7 +33,6 @@ struct lp_problem_t { num_cols(n), objective(n), A(m, n, nz), - Arow(1, 1, 1), rhs(m), lower(n), upper(n), @@ -45,7 +44,6 @@ struct lp_problem_t { i_t num_cols; std::vector objective; csc_matrix_t A; - csc_matrix_t Arow; std::vector rhs; std::vector lower; std::vector upper; From 874d45f31f68b5fa2b0bda159a6691b4aed3ce11 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 31 Oct 2025 18:37:19 +0100 Subject: [PATCH 29/51] fix incorrect report of infeasible solution --- cpp/src/dual_simplex/branch_and_bound.cpp | 53 +++++++++++------------ 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 20a3e6a994..b013920298 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -485,11 +485,13 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t 0 && stats_.nodes_unexplored == 0 && upper_bound == inf) { - settings_.log.printf("Integer infeasible.\n"); - mip_status = mip_status_t::INFEASIBLE; - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); + if (status_ == mip_exploration_status_t::COMPLETED) { + if (stats_.nodes_explored > 0 && stats_.nodes_unexplored == 0 && upper_bound == inf) { + settings_.log.printf("Integer infeasible.\n"); + mip_status = mip_status_t::INFEASIBLE; + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } } } @@ -717,11 +719,6 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod // to repair the heuristic solution. repair_heuristic_solutions(); - // Make a copy of the original LP. We will modify its bounds at each leaf - lp_problem_t leaf_problem = original_lp_; - std::vector row_sense; - bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); - f_t lower_bound = node->lower_bound; f_t upper_bound = get_upper_bound(); f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); @@ -769,6 +766,11 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod return; } + // Make a copy of the original LP. We will modify its bounds at each leaf + lp_problem_t leaf_problem = original_lp_; + std::vector row_sense; + bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); + node_status_t node_status = solve_node(node, *search_tree, leaf_problem, @@ -931,7 +933,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, } template -void branch_and_bound_t::best_first_thread(i_t id, +void branch_and_bound_t::best_first_thread(i_t task_id, search_tree_t& search_tree, const csr_matrix_t& Arow) { @@ -948,29 +950,29 @@ void branch_and_bound_t::best_first_thread(i_t id, while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && (active_subtrees_ > 0 || get_heap_size() > 0)) { - mip_node_t* node_ptr = nullptr; + mip_node_t* start_node = nullptr; // If there any node left in the heap, we pop the top node and explore it. mutex_heap_.lock(); if (heap_.size() > 0) { - node_ptr = heap_.top(); + start_node = heap_.top(); heap_.pop(); active_subtrees_++; } mutex_heap_.unlock(); - if (node_ptr != nullptr) { - if (get_upper_bound() < node_ptr->lower_bound) { + if (start_node != nullptr) { + if (get_upper_bound() < start_node->lower_bound) { // This node was put on the heap earlier but its lower bound is now greater than the // current upper bound - search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); - search_tree.update(node_ptr, node_status_t::FATHOMED); + search_tree.graphviz_node(settings_.log, start_node, "cutoff", start_node->lower_bound); + search_tree.update(start_node, node_status_t::FATHOMED); active_subtrees_--; continue; } // Best-first search with plunging - explore_subtree(id, node_ptr, search_tree, leaf_problem, presolver); + explore_subtree(task_id, start_node, search_tree, leaf_problem, presolver); active_subtrees_--; } @@ -986,7 +988,7 @@ void branch_and_bound_t::best_first_thread(i_t id, if (active_subtrees_ == 0) { status_ = mip_exploration_status_t::COMPLETED; } else { - local_lower_bounds_[id] = inf; + local_lower_bounds_[task_id] = inf; } } } @@ -1251,17 +1253,14 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut #pragma omp master { - if (status_ == mip_exploration_status_t::RUNNING && - (active_subtrees_ > 0 || get_heap_size() > 0)) { - for (i_t i = 0; i < settings_.num_bfs_threads; i++) { + for (i_t i = 0; i < settings_.num_bfs_threads; i++) { #pragma omp task - best_first_thread(i, search_tree, Arow); - } + best_first_thread(i, search_tree, Arow); + } - for (i_t i = 0; i < settings_.num_diving_threads; i++) { + for (i_t i = 0; i < settings_.num_diving_threads; i++) { #pragma omp task - diving_thread(Arow); - } + diving_thread(Arow); } } } From 5bfda565b5574dacac02d734d256899341311093 Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 31 Oct 2025 18:41:40 +0100 Subject: [PATCH 30/51] fixed typo --- cpp/src/dual_simplex/bounds_strengthening.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index fc1adb771b..9f92062e8f 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -15,7 +15,7 @@ * limitations under the License. */ -#include +#include #include #include From 04e3063035135bebfd3acdbbed97de219eabd6ca Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 31 Oct 2025 20:19:25 +0100 Subject: [PATCH 31/51] fixed typo --- cpp/src/dual_simplex/branch_and_bound.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index b013920298..7247116aed 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -17,7 +17,7 @@ #include #include -#include +#include #include #include #include From 1df3b4c6b4f555b862603a47603f436c1ebc8d5d Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 3 Nov 2025 11:30:29 +0100 Subject: [PATCH 32/51] rename variable --- cpp/src/dual_simplex/mip_node.hpp | 47 ++++++++++++++++----------- cpp/src/dual_simplex/pseudo_costs.cpp | 4 +-- 2 files changed, 30 insertions(+), 21 deletions(-) diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 1c51c09d12..b7d1b0a7d0 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -29,7 +29,7 @@ namespace cuopt::linear_programming::dual_simplex { enum class node_status_t : int { - ACTIVE = 0, // Node still in the tree + PENDING = 0, // Node still in the tree INTEGER_FEASIBLE = 1, // Node has an integer feasible solution INFEASIBLE = 2, // Node is infeasible FATHOMED = 3, // Node objective is greater than the upper bound @@ -38,7 +38,7 @@ enum class node_status_t : int { TIME_LIMIT = 6 // Time out during the LP relaxation }; -enum class round_dir_t { NONE = -1, DOWN = 0, UP = 1 }; +enum class rounding_direction_t { NONE = -1, DOWN = 0, UP = 1 }; bool inactive_status(node_status_t status); @@ -46,13 +46,13 @@ template class mip_node_t { public: mip_node_t(f_t root_lower_bound, const std::vector& basis) - : status(node_status_t::ACTIVE), + : status(node_status_t::PENDING), lower_bound(root_lower_bound), depth(0), parent(nullptr), node_id(0), branch_var(-1), - branch_dir(round_dir_t::NONE), + branch_dir(rounding_direction_t::NONE), vstatus(basis) { children[0] = nullptr; @@ -63,10 +63,10 @@ class mip_node_t { mip_node_t* parent_node, i_t node_num, i_t branch_variable, - round_dir_t branch_direction, + rounding_direction_t branch_direction, f_t branch_var_value, const std::vector& basis) - : status(node_status_t::ACTIVE), + : status(node_status_t::PENDING), lower_bound(parent_node->lower_bound), depth(parent_node->depth + 1), parent(parent_node), @@ -77,10 +77,10 @@ class mip_node_t { vstatus(basis) { - branch_var_lower = branch_direction == round_dir_t::DOWN ? problem.lower[branch_var] - : std::ceil(branch_var_value); - branch_var_upper = branch_direction == round_dir_t::DOWN ? std::floor(branch_var_value) - : problem.upper[branch_var]; + branch_var_lower = branch_direction == rounding_direction_t::DOWN ? problem.lower[branch_var] + : std::ceil(branch_var_value); + branch_var_upper = branch_direction == rounding_direction_t::DOWN ? std::floor(branch_var_value) + : problem.upper[branch_var]; children[0] = nullptr; children[1] = nullptr; } @@ -234,7 +234,7 @@ class mip_node_t { i_t depth; i_t node_id; i_t branch_var; - round_dir_t branch_dir; + rounding_direction_t branch_dir; f_t branch_var_lower; f_t branch_var_upper; f_t fractional_val; @@ -303,7 +303,7 @@ class search_tree_t { parent_node, ++id, branch_var, - round_dir_t::DOWN, + rounding_direction_t::DOWN, fractional_val, parent_vstatus); @@ -311,14 +311,23 @@ class search_tree_t { parent_node, down_child.get(), branch_var, - round_dir_t::DOWN, + rounding_direction_t::DOWN, std::floor(fractional_val)); - auto up_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, round_dir_t::UP, fractional_val, parent_vstatus); + auto up_child = std::make_unique>(original_lp, + parent_node, + ++id, + branch_var, + rounding_direction_t::UP, + fractional_val, + parent_vstatus); - graphviz_edge( - log, parent_node, up_child.get(), branch_var, round_dir_t::UP, std::ceil(fractional_val)); + graphviz_edge(log, + parent_node, + up_child.get(), + branch_var, + rounding_direction_t::UP, + std::ceil(fractional_val)); assert(parent_vstatus.size() == original_lp.num_cols); parent_node->add_children(std::move(down_child), @@ -339,7 +348,7 @@ class search_tree_t { const mip_node_t* origin_ptr, const mip_node_t* dest_ptr, const i_t branch_var, - round_dir_t branch_dir, + rounding_direction_t branch_dir, const f_t bound) { if (write_graphviz) { @@ -347,7 +356,7 @@ class search_tree_t { origin_ptr->node_id, dest_ptr->node_id, branch_var, - branch_dir == round_dir_t::DOWN ? "<=" : ">=", + branch_dir == rounding_direction_t::DOWN ? "<=" : ">=", bound); } } diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index 1a10dc73e2..87c11e49e5 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -210,10 +210,10 @@ void pseudo_costs_t::update_pseudo_costs(mip_node_t* node_pt { mutex.lock(); const f_t change_in_obj = leaf_objective - node_ptr->lower_bound; - const f_t frac = node_ptr->branch_dir == round_dir_t::DOWN + const f_t frac = node_ptr->branch_dir == rounding_direction_t::DOWN ? node_ptr->fractional_val - std::floor(node_ptr->fractional_val) : std::ceil(node_ptr->fractional_val) - node_ptr->fractional_val; - if (node_ptr->branch_dir == round_dir_t::DOWN) { + if (node_ptr->branch_dir == rounding_direction_t::DOWN) { pseudo_cost_sum_down[node_ptr->branch_var] += change_in_obj / frac; pseudo_cost_num_down[node_ptr->branch_var]++; } else { From 5ec4472ce5d4515118fa9f20be3c10882ea03896 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 10 Nov 2025 11:10:21 +0100 Subject: [PATCH 33/51] addresses reviewer's comments --- cpp/src/dual_simplex/branch_and_bound.cpp | 395 ++++++++++++---------- cpp/src/dual_simplex/branch_and_bound.hpp | 36 +- cpp/src/dual_simplex/mip_node.hpp | 8 +- cpp/src/dual_simplex/presolve.cpp | 11 +- 4 files changed, 251 insertions(+), 199 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 7247116aed..c060ecd0f8 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -201,7 +201,7 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound) } } -inline const char* thread_type_symbol(thread_type_t type) +inline const char* feasible_solution_symbol(thread_type_t type) { switch (type) { case thread_type_t::EXPLORATION: return "B"; @@ -210,6 +210,12 @@ inline const char* thread_type_symbol(thread_type_t type) } } +inline bool has_children(node_children_status_t status) +{ + return status == node_children_status_t::UP_CHILDREN_FIRST || + status == node_children_status_t::DOWN_CHILDREN_FIRST; +} + } // namespace template @@ -222,9 +228,9 @@ branch_and_bound_t::branch_and_bound_t( incumbent_(1), root_relax_soln_(1, 1), pc_(1), - status_(mip_exploration_status_t::UNSET) + solver_status_(mip_exploration_status_t::UNSET) { - stats_.start_time = tic(); + exploration_stats_.start_time = tic(); dualize_info_t dualize_info; convert_user_problem(original_problem_, settings_, original_lp_, new_slacks_, dualize_info); full_variable_types(original_problem_, original_lp_, var_types_); @@ -305,7 +311,7 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu mutex_upper_.unlock(); if (is_feasible) { - if (status_ == mip_exploration_status_t::RUNNING) { + if (solver_status_ == mip_exploration_status_t::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 gap = user_mip_gap(user_obj, user_lower); @@ -315,11 +321,11 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu user_obj, user_lower, gap.c_str(), - toc(stats_.start_time)); + 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(stats_.start_time)); + toc(exploration_stats_.start_time)); } } @@ -428,7 +434,7 @@ void branch_and_bound_t::repair_heuristic_solutions() obj, lower, user_gap.c_str(), - toc(stats_.start_time)); + toc(exploration_stats_.start_time)); if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -449,12 +455,12 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t::set_final_solution(mip_solution_t::set_final_solution(mip_solution_t 0 && stats_.nodes_unexplored == 0 && upper_bound == inf) { + if (solver_status_ == mip_exploration_status_t::COMPLETED) { + if (exploration_stats_.nodes_explored > 0 && exploration_stats_.nodes_unexplored == 0 && + upper_bound == inf) { settings_.log.printf("Integer infeasible.\n"); mip_status = mip_status_t::INFEASIBLE; if (settings_.heuristic_preemption_callback != nullptr) { @@ -498,8 +506,8 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t::add_feasible_solution(f_t leaf_objective, thread_type_t thread_type) { bool send_solution = false; - i_t nodes_explored = stats_.nodes_explored; - i_t nodes_unexplored = stats_.nodes_unexplored; + i_t nodes_explored = exploration_stats_.nodes_explored; + i_t nodes_unexplored = exploration_stats_.nodes_unexplored; mutex_upper_.lock(); if (leaf_objective < upper_bound_) { @@ -521,16 +529,17 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t lower_bound = get_lower_bound(); f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%s%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - thread_type_symbol(thread_type), - nodes_explored, - nodes_unexplored, - obj, - lower, - leaf_depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - user_mip_gap(obj, lower).c_str(), - toc(stats_.start_time)); + settings_.log.printf( + "%s%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + feasible_solution_symbol(thread_type), + nodes_explored, + nodes_unexplored, + obj, + lower, + leaf_depth, + nodes_explored > 0 ? exploration_stats_.total_lp_iters / nodes_explored : 0, + user_mip_gap(obj, lower).c_str(), + toc(exploration_stats_.start_time)); send_solution = true; } @@ -544,8 +553,7 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, } template -std::pair*, mip_node_t*> -branch_and_bound_t::child_selection(mip_node_t* node_ptr) +rounding_direction_t branch_and_bound_t::child_selection(mip_node_t* node_ptr) { const i_t branch_var = node_ptr->get_down_child()->branch_var; const f_t branch_var_val = node_ptr->get_down_child()->fractional_val; @@ -556,23 +564,24 @@ branch_and_bound_t::child_selection(mip_node_t* node_ptr) constexpr f_t eps = 1e-6; if (down_dist < up_dist + eps) { - return std::make_pair(node_ptr->get_down_child(), node_ptr->get_up_child()); + return rounding_direction_t::DOWN; } else { - return std::make_pair(node_ptr->get_up_child(), node_ptr->get_down_child()); + return rounding_direction_t::UP; } } template -node_status_t branch_and_bound_t::solve_node(mip_node_t* node_ptr, - search_tree_t& search_tree, - lp_problem_t& leaf_problem, - bounds_strengthening_t& presolver, - thread_type_t thread_type, - bool recompute_bounds, - const std::vector& root_lower, - const std::vector& root_upper, - logger_t& log) +node_children_status_t branch_and_bound_t::solve_node( + mip_node_t* node_ptr, + search_tree_t& search_tree, + lp_problem_t& leaf_problem, + bounds_strengthening_t& node_presolver, + thread_type_t thread_type, + bool recompute_bounds, + const std::vector& root_lower, + const std::vector& root_upper, + logger_t& log) { const f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; const f_t upper_bound = get_upper_bound(); @@ -585,24 +594,25 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod lp_settings.set_log(false); lp_settings.cut_off = upper_bound + settings_.dual_tol; lp_settings.inside_mip = 2; - lp_settings.time_limit = settings_.time_limit - toc(stats_.start_time); + lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); // Reset the bound_changed markers - std::fill(presolver.bounds_changed.begin(), presolver.bounds_changed.end(), false); + std::fill(node_presolver.bounds_changed.begin(), node_presolver.bounds_changed.end(), false); // Set the correct bounds for the leaf problem if (recompute_bounds) { leaf_problem.lower = root_lower; leaf_problem.upper = root_upper; - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, presolver.bounds_changed); + node_ptr->get_variable_bounds( + leaf_problem.lower, leaf_problem.upper, node_presolver.bounds_changed); } else { node_ptr->update_branched_variable_bounds( - leaf_problem.lower, leaf_problem.upper, presolver.bounds_changed); + leaf_problem.lower, leaf_problem.upper, node_presolver.bounds_changed); } bool feasible = - presolver.bounds_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); + node_presolver.bounds_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -622,14 +632,14 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod leaf_edge_norms); if (lp_status == dual::status_t::NUMERICAL) { - log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); + log.printf("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); lp_status_t second_status = solve_linear_program_advanced( leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); lp_status = convert_lp_status_to_dual_status(second_status); } - stats_.total_lp_solve_time += toc(lp_start_time); - stats_.total_lp_iters += node_iter; + exploration_stats_.total_lp_solve_time += toc(lp_start_time); + exploration_stats_.total_lp_iters += node_iter; } if (lp_status == dual::status_t::DUAL_UNBOUNDED) { @@ -637,7 +647,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod node_ptr->lower_bound = inf; search_tree.graphviz_node(log, node_ptr, "infeasible", 0.0); search_tree.update(node_ptr, node_status_t::INFEASIBLE); - return node_status_t::INFEASIBLE; + return node_children_status_t::NO_CHILDREN; } else if (lp_status == dual::status_t::CUTOFF) { // Node was cut off. Do not branch @@ -645,7 +655,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); search_tree.graphviz_node(log, node_ptr, "cut off", leaf_objective); search_tree.update(node_ptr, node_status_t::FATHOMED); - return node_status_t::FATHOMED; + return node_children_status_t::NO_CHILDREN; } else if (lp_status == dual::status_t::OPTIMAL) { // LP was feasible @@ -663,7 +673,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, thread_type); search_tree.graphviz_node(log, node_ptr, "integer feasible", leaf_objective); search_tree.update(node_ptr, node_status_t::INTEGER_FEASIBLE); - return node_status_t::INTEGER_FEASIBLE; + return node_children_status_t::NO_CHILDREN; } else if (leaf_objective <= upper_bound + abs_fathom_tol) { logger_t pc_log = log; @@ -676,17 +686,23 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod search_tree.branch( node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, leaf_problem, log); node_ptr->status = node_status_t::HAS_CHILDREN; - return node_status_t::HAS_CHILDREN; + + rounding_direction_t round_dir = child_selection(node_ptr); + + if (round_dir == rounding_direction_t::UP) { + return node_children_status_t::UP_CHILDREN_FIRST; + } else { + return node_children_status_t::DOWN_CHILDREN_FIRST; + } } else { search_tree.graphviz_node(log, node_ptr, "fathomed", leaf_objective); search_tree.update(node_ptr, node_status_t::FATHOMED); - return node_status_t::FATHOMED; + return node_children_status_t::NO_CHILDREN; } } else if (lp_status == dual::status_t::TIME_LIMIT) { search_tree.graphviz_node(log, node_ptr, "timeout", 0.0); - search_tree.update(node_ptr, node_status_t::TIME_LIMIT); - return node_status_t::TIME_LIMIT; + return node_children_status_t::TIME_LIMIT; } else { if (thread_type == thread_type_t::EXPLORATION) { @@ -702,7 +718,7 @@ node_status_t branch_and_bound_t::solve_node(mip_node_t* nod search_tree.graphviz_node(log, node_ptr, "numerical", 0.0); search_tree.update(node_ptr, node_status_t::NUMERICAL); - return node_status_t::NUMERICAL; + return node_children_status_t::NUMERICAL; } } @@ -712,7 +728,7 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod const csr_matrix_t& Arow, i_t initial_heap_size) { - if (status_ != mip_exploration_status_t::RUNNING) { return; } + if (solver_status_ != mip_exploration_status_t::RUNNING) { return; } // Note that we do not know which thread will execute the // `exploration_ramp_up` task, so we allow to any thread @@ -723,9 +739,9 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod f_t upper_bound = get_upper_bound(); f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); f_t abs_gap = upper_bound - lower_bound; - i_t nodes_explored = (++stats_.nodes_explored); - i_t nodes_unexplored = (--stats_.nodes_unexplored); - stats_.nodes_since_last_log++; + i_t nodes_explored = (++exploration_stats_.nodes_explored); + i_t nodes_unexplored = (--exploration_stats_.nodes_unexplored); + exploration_stats_.nodes_since_last_log++; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree->graphviz_node(settings_.log, node, "cutoff", node->lower_bound); @@ -733,63 +749,66 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod return; } - f_t now = toc(stats_.start_time); - f_t time_since_last_log = stats_.last_log == 0 ? 1.0 : toc(stats_.last_log); + f_t now = toc(exploration_stats_.start_time); + f_t time_since_last_log = + exploration_stats_.last_log == 0 ? 1.0 : toc(exploration_stats_.last_log); - if (((stats_.nodes_since_last_log >= 10 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + if (((exploration_stats_.nodes_since_last_log >= 10 || + abs_gap < 10 * settings_.absolute_mip_gap_tol) && (time_since_last_log >= 1)) || (time_since_last_log > 30) || now > settings_.time_limit) { // Check if no new node was explored until now. If this is the case, // only the last thread should report the progress - if (stats_.nodes_explored.load() == nodes_explored) { - stats_.nodes_since_last_log = 0; - stats_.last_log = tic(); + if (exploration_stats_.nodes_explored.load() == nodes_explored) { + exploration_stats_.nodes_since_last_log = 0; + exploration_stats_.last_log = tic(); f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, root_objective_); std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); + settings_.log.printf( + " %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node->depth, + nodes_explored > 0 ? exploration_stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); } } if (now > settings_.time_limit) { - status_ = mip_exploration_status_t::TIME_LIMIT; + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; } // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); - - node_status_t node_status = solve_node(node, - *search_tree, - leaf_problem, - presolver, - thread_type_t::EXPLORATION, - true, - original_lp_.lower, - original_lp_.upper, - settings_.log); - - if (node_status == node_status_t::TIME_LIMIT) { - status_ = mip_exploration_status_t::TIME_LIMIT; + bounds_strengthening_t node_presolver(leaf_problem, Arow, row_sense, var_types_); + + node_children_status_t status = solve_node(node, + *search_tree, + leaf_problem, + node_presolver, + thread_type_t::EXPLORATION, + true, + original_lp_.lower, + original_lp_.upper, + settings_.log); + + if (status == node_children_status_t::TIME_LIMIT) { + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; - } else if (node_status == node_status_t::HAS_CHILDREN) { - stats_.nodes_unexplored += 2; + } else if (has_children(status)) { + exploration_stats_.nodes_unexplored += 2; // If we haven't generated enough nodes to keep the threads busy, continue the ramp up phase - if (stats_.nodes_unexplored < initial_heap_size) { + if (exploration_stats_.nodes_unexplored < initial_heap_size) { #pragma omp task exploration_ramp_up(node->get_down_child(), search_tree, Arow, initial_heap_size); @@ -811,13 +830,13 @@ void branch_and_bound_t::explore_subtree(i_t task_id, mip_node_t* start_node, search_tree_t& search_tree, lp_problem_t& leaf_problem, - bounds_strengthening_t& presolver) + bounds_strengthening_t& node_presolver) { bool recompute_bounds = true; std::deque*> stack; stack.push_front(start_node); - while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { + while (stack.size() > 0 && solver_status_ == mip_exploration_status_t::RUNNING) { if (task_id == 0) { repair_heuristic_solutions(); } mip_node_t* node_ptr = stack.front(); @@ -836,9 +855,9 @@ void branch_and_bound_t::explore_subtree(i_t task_id, // - The lower bound of the parent is lower or equal to its children assert(task_id < local_lower_bounds_.size()); local_lower_bounds_[task_id] = lower_bound; - i_t nodes_explored = (++stats_.nodes_explored); - i_t nodes_unexplored = (--stats_.nodes_unexplored); - stats_.nodes_since_last_log++; + i_t nodes_explored = (++exploration_stats_.nodes_explored); + i_t nodes_unexplored = (--exploration_stats_.nodes_unexplored); + exploration_stats_.nodes_since_last_log++; if (lower_bound > upper_bound || rel_gap < settings_.relative_mip_gap_tol) { search_tree.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound); @@ -847,53 +866,56 @@ void branch_and_bound_t::explore_subtree(i_t task_id, continue; } - f_t now = toc(stats_.start_time); + f_t now = toc(exploration_stats_.start_time); if (task_id == 0) { - f_t time_since_last_log = stats_.last_log == 0 ? 1.0 : toc(stats_.last_log); + f_t time_since_last_log = + exploration_stats_.last_log == 0 ? 1.0 : toc(exploration_stats_.last_log); - if (((stats_.nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) && + if (((exploration_stats_.nodes_since_last_log >= 1000 || + abs_gap < 10 * settings_.absolute_mip_gap_tol) && time_since_last_log >= 1) || (time_since_last_log > 30) || now > settings_.time_limit) { f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - nodes_explored, - nodes_unexplored, - obj, - user_lower, - node_ptr->depth, - nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, - gap_user.c_str(), - now); - stats_.last_log = tic(); - stats_.nodes_since_last_log = 0; + settings_.log.printf( + " %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + nodes_explored, + nodes_unexplored, + obj, + user_lower, + node_ptr->depth, + nodes_explored > 0 ? exploration_stats_.total_lp_iters / nodes_explored : 0, + gap_user.c_str(), + now); + exploration_stats_.last_log = tic(); + exploration_stats_.nodes_since_last_log = 0; } } if (now > settings_.time_limit) { - status_ = mip_exploration_status_t::TIME_LIMIT; + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; } - node_status_t node_status = solve_node(node_ptr, - search_tree, - leaf_problem, - presolver, - thread_type_t::EXPLORATION, - recompute_bounds, - original_lp_.lower, - original_lp_.upper, - settings_.log); + node_children_status_t status = solve_node(node_ptr, + search_tree, + leaf_problem, + node_presolver, + thread_type_t::EXPLORATION, + recompute_bounds, + original_lp_.lower, + original_lp_.upper, + settings_.log); - recompute_bounds = node_status != node_status_t::HAS_CHILDREN; + recompute_bounds = !has_children(status); - if (node_status == node_status_t::TIME_LIMIT) { - status_ = mip_exploration_status_t::TIME_LIMIT; + if (status == node_children_status_t::TIME_LIMIT) { + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; - } else if (node_status == node_status_t::HAS_CHILDREN) { + } else if (has_children(status)) { // The stack should only contain the children of the current parent. // If the stack size is greater than 0, // we pop the current node from the stack and place it in the global heap, @@ -911,7 +933,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, if (get_heap_size() > settings_.num_bfs_threads) { std::vector lower = original_lp_.lower; std::vector upper = original_lp_.upper; - node->get_variable_bounds(lower, upper, presolver.bounds_changed); + node->get_variable_bounds(lower, upper, node_presolver.bounds_changed); mutex_dive_queue_.lock(); dive_queue_.emplace(node->detach_copy(), std::move(lower), std::move(upper)); @@ -923,11 +945,15 @@ void branch_and_bound_t::explore_subtree(i_t task_id, mutex_heap_.unlock(); } - stats_.nodes_unexplored += 2; + exploration_stats_.nodes_unexplored += 2; - auto [first, second] = child_selection(node_ptr); - stack.push_front(second); - stack.push_front(first); + if (status == node_children_status_t::UP_CHILDREN_FIRST) { + stack.push_front(node_ptr->get_down_child()); + stack.push_front(node_ptr->get_up_child()); + } else { + stack.push_front(node_ptr->get_up_child()); + stack.push_front(node_ptr->get_down_child()); + } } } } @@ -945,10 +971,10 @@ void branch_and_bound_t::best_first_thread(i_t task_id, // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); + bounds_strengthening_t node_presolver(leaf_problem, Arow, row_sense, var_types_); - while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && - rel_gap > settings_.relative_mip_gap_tol && + while (solver_status_ == mip_exploration_status_t::RUNNING && + abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && (active_subtrees_ > 0 || get_heap_size() > 0)) { mip_node_t* start_node = nullptr; @@ -972,7 +998,7 @@ void branch_and_bound_t::best_first_thread(i_t task_id, } // Best-first search with plunging - explore_subtree(task_id, start_node, search_tree, leaf_problem, presolver); + explore_subtree(task_id, start_node, search_tree, leaf_problem, node_presolver); active_subtrees_--; } @@ -984,9 +1010,9 @@ void branch_and_bound_t::best_first_thread(i_t task_id, // Check if it is the last thread that exited the loop and no // timeout or numerical error has happen. - if (status_ == mip_exploration_status_t::RUNNING) { + if (solver_status_ == mip_exploration_status_t::RUNNING) { if (active_subtrees_ == 0) { - status_ = mip_exploration_status_t::COMPLETED; + solver_status_ = mip_exploration_status_t::COMPLETED; } else { local_lower_bounds_[task_id] = inf; } @@ -1002,9 +1028,9 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A // Make a copy of the original LP. We will modify its bounds at each leaf lp_problem_t leaf_problem = original_lp_; std::vector row_sense; - bounds_strengthening_t presolver(leaf_problem, Arow, row_sense, var_types_); + bounds_strengthening_t node_presolver(leaf_problem, Arow, row_sense, var_types_); - while (status_ == mip_exploration_status_t::RUNNING && + while (solver_status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { std::optional> start_node; @@ -1020,7 +1046,7 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A std::deque*> stack; stack.push_front(&subtree.root); - while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { + while (stack.size() > 0 && solver_status_ == mip_exploration_status_t::RUNNING) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); f_t upper_bound = get_upper_bound(); @@ -1031,27 +1057,32 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A continue; } - if (toc(stats_.start_time) > settings_.time_limit) { return; } + if (toc(exploration_stats_.start_time) > settings_.time_limit) { return; } - node_status_t node_status = solve_node(node_ptr, - subtree, - leaf_problem, - presolver, - thread_type_t::DIVING, - recompute_bounds, - start_node->lower, - start_node->upper, - log); + node_children_status_t status = solve_node(node_ptr, + subtree, + leaf_problem, + node_presolver, + thread_type_t::DIVING, + recompute_bounds, + start_node->lower, + start_node->upper, + log); - recompute_bounds = node_status != node_status_t::HAS_CHILDREN; + recompute_bounds = has_children(status); - if (node_status == node_status_t::TIME_LIMIT) { + if (status == node_children_status_t::TIME_LIMIT) { + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; - } else if (node_status == node_status_t::HAS_CHILDREN) { - auto [first, second] = child_selection(node_ptr); - stack.push_front(second); - stack.push_front(first); + } else if (has_children(status)) { + if (status == node_children_status_t::UP_CHILDREN_FIRST) { + stack.push_front(node_ptr->get_down_child()); + stack.push_front(node_ptr->get_up_child()); + } else { + stack.push_front(node_ptr->get_up_child()); + stack.push_front(node_ptr->get_down_child()); + } } if (stack.size() > 1) { @@ -1065,7 +1096,7 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A std::vector lower = start_node->lower; std::vector upper = start_node->upper; - new_node->get_variable_bounds(lower, upper, presolver.bounds_changed); + new_node->get_variable_bounds(lower, upper, node_presolver.bounds_changed); mutex_dive_queue_.lock(); dive_queue_.emplace(new_node->detach_copy(), std::move(lower), std::move(upper)); @@ -1081,10 +1112,10 @@ template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { logger_t log; - log.log = false; - status_ = mip_exploration_status_t::UNSET; - stats_.nodes_unexplored = 0; - stats_.nodes_explored = 0; + log.log = false; + solver_status_ = mip_exploration_status_t::UNSET; + exploration_stats_.nodes_unexplored = 0; + exploration_stats_.nodes_explored = 0; if (guess_.size() != 0) { std::vector crushed_guess; @@ -1106,12 +1137,16 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); settings_.log.printf("Solving LP root relaxation\n"); - simplex_solver_settings_t lp_settings = settings_; - lp_settings.inside_mip = 1; - lp_status_t root_status = solve_linear_program_advanced( - original_lp_, stats_.start_time, lp_settings, root_relax_soln_, root_vstatus_, edge_norms_); - stats_.total_lp_iters = root_relax_soln_.iterations; - stats_.total_lp_solve_time = toc(stats_.start_time); + simplex_solver_settings_t lp_settings = settings_; + lp_settings.inside_mip = 1; + lp_status_t root_status = solve_linear_program_advanced(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + root_vstatus_, + 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"); @@ -1132,7 +1167,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (root_status == lp_status_t::TIME_LIMIT) { - status_ = mip_exploration_status_t::TIME_LIMIT; + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return set_final_solution(solution, -inf); } @@ -1174,7 +1209,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut solution.simplex_iterations = root_relax_soln_.iterations; settings_.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", compute_user_objective(original_lp_, root_objective_), - toc(stats_.start_time)); + toc(exploration_stats_.start_time)); if (settings_.solution_callback != nullptr) { settings_.solution_callback(solution.x, solution.objective); @@ -1188,7 +1223,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut pc_.resize(original_lp_.num_cols); strong_branching(original_lp_, settings_, - stats_.start_time, + exploration_stats_.start_time, var_types_, root_relax_soln_.x, fractional, @@ -1197,8 +1232,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut edge_norms_, pc_); - if (toc(stats_.start_time) > settings_.time_limit) { - status_ = mip_exploration_status_t::TIME_LIMIT; + if (toc(exploration_stats_.start_time) > settings_.time_limit) { + solver_status_ = mip_exploration_status_t::TIME_LIMIT; return set_final_solution(solution, root_objective_); } @@ -1222,14 +1257,14 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut " | Explored | Unexplored | Objective | Bound | Depth | Iter/Node | Gap " "| Time |\n"); - stats_.nodes_explored = 0; - stats_.nodes_unexplored = 2; - stats_.nodes_since_last_log = 0; - stats_.last_log = tic(); - active_subtrees_ = 0; - min_diving_queue_size_ = 4 * settings_.num_diving_threads; - status_ = mip_exploration_status_t::RUNNING; - lower_bound_ceiling_ = inf; + exploration_stats_.nodes_explored = 0; + exploration_stats_.nodes_unexplored = 2; + exploration_stats_.nodes_since_last_log = 0; + exploration_stats_.last_log = tic(); + active_subtrees_ = 0; + min_diving_queue_size_ = 4 * settings_.num_diving_threads; + solver_status_ = mip_exploration_status_t::RUNNING; + lower_bound_ceiling_ = inf; csr_matrix_t Arow(1, 1, 0); original_lp_.A.to_compressed_row(Arow); diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 00da2f96df..57ce4ed3d3 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -52,6 +52,15 @@ enum class mip_exploration_status_t { COMPLETED = 5, // The solver finished exploring the tree }; +enum class node_children_status_t { + NO_CHILDREN = 0, // The node does not produced children + UP_CHILDREN_FIRST = 1, // The up child should be explored first + DOWN_CHILDREN_FIRST = 2, // The down child should be explored first + TIME_LIMIT = 3, // The solver reached a time limit + ITERATION_LIMIT = 4, // The solver reached a iteration limit + NUMERICAL = 5 // The solver encounter a numerical error when solving the node +}; + // Indicate the search and variable selection algorithms used by the thread (See [1]). // // [1] T. Achterberg, “Constraint Integer Programming,” PhD, Technischen Universität Berlin, @@ -130,7 +139,7 @@ class branch_and_bound_t { // This should only be used by the main thread f_t last_log = 0.0; omp_atomic_t nodes_since_last_log = 0; - } stats_; + } exploration_stats_; // Mutex for repair omp_mutex_t mutex_repair_; @@ -158,7 +167,7 @@ class branch_and_bound_t { i_t min_diving_queue_size_; // Global status of the solver. - omp_atomic_t status_; + omp_atomic_t solver_status_; // In case, a best-first thread encounters a numerical issue when solving a node, // its blocks the progression of the lower bound. @@ -189,7 +198,7 @@ class branch_and_bound_t { mip_node_t* start_node, search_tree_t& search_tree, lp_problem_t& leaf_problem, - bounds_strengthening_t& presolver); + bounds_strengthening_t& node_presolver); // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. @@ -202,19 +211,18 @@ class branch_and_bound_t { void diving_thread(const csr_matrix_t& Arow); // Solve the LP relaxation of a leaf node and update the tree. - node_status_t solve_node(mip_node_t* node_ptr, - search_tree_t& search_tree, - lp_problem_t& leaf_problem, - bounds_strengthening_t& presolver, - thread_type_t thread_type, - bool recompute, - const std::vector& root_lower, - const std::vector& root_upper, - logger_t& log); + node_children_status_t solve_node(mip_node_t* node_ptr, + search_tree_t& search_tree, + lp_problem_t& leaf_problem, + bounds_strengthening_t& node_presolver, + thread_type_t thread_type, + bool recompute, + const std::vector& root_lower, + const std::vector& root_upper, + logger_t& log); // Sort the children based on the Martin's criteria. - std::pair*, mip_node_t*> child_selection( - mip_node_t* node_ptr); + rounding_direction_t child_selection(mip_node_t* node_ptr); }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index b7d1b0a7d0..3b518b7755 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -34,11 +34,10 @@ enum class node_status_t : int { INFEASIBLE = 2, // Node is infeasible FATHOMED = 3, // Node objective is greater than the upper bound HAS_CHILDREN = 4, // Node has children to explore - NUMERICAL = 5, // Encountered numerical issue when solving the LP relaxation - TIME_LIMIT = 6 // Time out during the LP relaxation + NUMERICAL = 5 // Encountered numerical issue when solving the LP relaxation }; -enum class rounding_direction_t { NONE = -1, DOWN = 0, UP = 1 }; +enum class rounding_direction_t : int8_t { NONE = -1, DOWN = 0, UP = 1 }; bool inactive_status(node_status_t status); @@ -110,7 +109,8 @@ class mip_node_t { assert(bounds_changed.size() > branch_var); // If the bounds have already been updated on another node, - // skip this node as it contains a less tight bounds. + // skip this node as it contains looser bounds, since we + // are traversing up the tree toward the root if (bounds_changed[branch_var]) { return; } // Apply the bounds at the current node diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 29e63d8676..3620667d20 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -17,6 +17,7 @@ #include +#include #include #include #include @@ -576,7 +577,15 @@ void convert_user_problem(const user_problem_t& user_problem, convert_greater_to_less(user_problem, row_sense, problem, greater_rows, less_rows); } - // bounds strengthening was moved to node_presolve.hpp + constexpr bool run_bounds_strengthening = false; + if (run_bounds_strengthening) { + csr_matrix_t Arow(1, 1, 1); + problem.A.to_compressed_row(Arow); + + // Empty var_types means that all variables are continuous + bounds_strengthening_t bounds_strenghtening(problem, Arow, row_sense, {}); + bounds_strenghtening.bounds_strengthening(problem.lower, problem.upper, settings); + } settings.log.debug( "equality rows %d less rows %d columns %d\n", equal_rows, less_rows, problem.num_cols); From 015bb00a375e4cfec02cabeba7ee64f8e053bf61 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 10 Nov 2025 11:19:16 +0100 Subject: [PATCH 34/51] fixed incorrect stats update --- cpp/src/dual_simplex/branch_and_bound.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index c060ecd0f8..b8fce1f394 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -638,8 +638,10 @@ node_children_status_t branch_and_bound_t::solve_node( lp_status = convert_lp_status_to_dual_status(second_status); } - exploration_stats_.total_lp_solve_time += toc(lp_start_time); - exploration_stats_.total_lp_iters += node_iter; + if (thread_type == thread_type_t::EXPLORATION) { + exploration_stats_.total_lp_solve_time += toc(lp_start_time); + exploration_stats_.total_lp_iters += node_iter; + } } if (lp_status == dual::status_t::DUAL_UNBOUNDED) { @@ -1069,7 +1071,7 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A start_node->upper, log); - recompute_bounds = has_children(status); + recompute_bounds = !has_children(status); if (status == node_children_status_t::TIME_LIMIT) { solver_status_ = mip_exploration_status_t::TIME_LIMIT; From ac199a32606189cc3acdd718f64143f81d3db046 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 10 Nov 2025 11:44:28 +0100 Subject: [PATCH 35/51] fix typo and styling --- cpp/src/dual_simplex/bounds_strengthening.cpp | 12 ------------ cpp/src/dual_simplex/bounds_strengthening.hpp | 12 ------------ cpp/src/dual_simplex/diving_queue.hpp | 12 ------------ cpp/src/dual_simplex/presolve.cpp | 4 ++-- 4 files changed, 2 insertions(+), 38 deletions(-) diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index 9f92062e8f..6a7d5ee62d 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -1,18 +1,6 @@ /* * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. */ #include diff --git a/cpp/src/dual_simplex/bounds_strengthening.hpp b/cpp/src/dual_simplex/bounds_strengthening.hpp index 28df35d0dc..fb237c7bb0 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.hpp +++ b/cpp/src/dual_simplex/bounds_strengthening.hpp @@ -1,18 +1,6 @@ /* * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. */ #pragma once diff --git a/cpp/src/dual_simplex/diving_queue.hpp b/cpp/src/dual_simplex/diving_queue.hpp index 57a8bacf19..91d0aa3309 100644 --- a/cpp/src/dual_simplex/diving_queue.hpp +++ b/cpp/src/dual_simplex/diving_queue.hpp @@ -1,18 +1,6 @@ /* * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. */ #pragma once diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index bac7d639be..fcc33a0e98 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -573,8 +573,8 @@ void convert_user_problem(const user_problem_t& user_problem, problem.A.to_compressed_row(Arow); // Empty var_types means that all variables are continuous - bounds_strengthening_t bounds_strenghtening(problem, Arow, row_sense, {}); - bounds_strenghtening.bounds_strengthening(problem.lower, problem.upper, settings); + bounds_strengthening_t strengthening(problem, Arow, row_sense, {}); + strengthening.bounds_strengthening(problem.lower, problem.upper, settings); } settings.log.debug( From 770adfb0fff3d6601bb948be6218961b74be9a69 Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 10 Nov 2025 15:12:20 +0100 Subject: [PATCH 36/51] small refactor --- cpp/src/dual_simplex/bounds_strengthening.cpp | 4 +++- cpp/src/dual_simplex/bounds_strengthening.hpp | 4 +++- cpp/src/dual_simplex/branch_and_bound.cpp | 23 +++++++++---------- cpp/src/dual_simplex/branch_and_bound.hpp | 4 ++-- cpp/src/dual_simplex/mip_node.hpp | 5 ++-- 5 files changed, 21 insertions(+), 19 deletions(-) diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index 6a7d5ee62d..64e745d5f5 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -1,7 +1,9 @@ +/* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ +/* clang-format on */ #include diff --git a/cpp/src/dual_simplex/bounds_strengthening.hpp b/cpp/src/dual_simplex/bounds_strengthening.hpp index fb237c7bb0..8e38cde15c 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.hpp +++ b/cpp/src/dual_simplex/bounds_strengthening.hpp @@ -1,7 +1,9 @@ +/* clang-format off */ /* - * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. * SPDX-License-Identifier: Apache-2.0 */ +/* clang-format on */ #pragma once diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index aa60379edb..6b45053b8d 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -668,11 +668,9 @@ node_children_status_t branch_and_bound_t::solve_node( return node_children_status_t::NO_CHILDREN; } else if (leaf_objective <= upper_bound + abs_fathom_tol) { - logger_t pc_log = log; - pc_log.log = false; - // Choose fractional variable to branch on - const i_t branch_var = pc_.variable_selection(leaf_fractional, leaf_solution.x, pc_log); + const i_t branch_var = + pc_.variable_selection(leaf_fractional, leaf_solution.x, lp_settings.log); assert(leaf_vstatus.size() == leaf_problem.num_cols); search_tree.branch( @@ -928,7 +926,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, node->get_variable_bounds(lower, upper, node_presolver.bounds_changed); mutex_dive_queue_.lock(); - dive_queue_.emplace(node->detach_copy(), std::move(lower), std::move(upper)); + diving_queue_.emplace(node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } @@ -1027,7 +1025,7 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A std::optional> start_node; mutex_dive_queue_.lock(); - if (dive_queue_.size() > 0) { start_node = dive_queue_.pop(); } + if (diving_queue_.size() > 0) { start_node = diving_queue_.pop(); } mutex_dive_queue_.unlock(); if (start_node.has_value()) { @@ -1082,7 +1080,7 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A // best first search, then we split the current subtree at the // lowest possible point and move to the queue, so it can // be picked by another thread. - if (dive_queue_.size() < min_diving_queue_size_) { + if (diving_queue_.size() < min_diving_queue_size_) { mip_node_t* new_node = stack.back(); stack.pop_back(); @@ -1091,7 +1089,7 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A new_node->get_variable_bounds(lower, upper, node_presolver.bounds_changed); mutex_dive_queue_.lock(); - dive_queue_.emplace(new_node->detach_copy(), std::move(lower), std::move(upper)); + diving_queue_.emplace(new_node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } } @@ -1241,7 +1239,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut original_lp_, log); - settings_.log.printf("Exploring the B&B tree using %d best-first threads and %d diving threads\n", + csr_matrix_t Arow(1, 1, 0); + original_lp_.A.to_compressed_row(Arow); + + settings_.log.printf("Exploring the B&B tree using %d threads (best-first = %d, diving = %d)\n", + settings_.num_threads, settings_.num_bfs_threads, settings_.num_diving_threads); @@ -1258,9 +1260,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut solver_status_ = mip_exploration_status_t::RUNNING; lower_bound_ceiling_ = inf; - csr_matrix_t Arow(1, 1, 0); - original_lp_.A.to_compressed_row(Arow); - #pragma omp parallel num_threads(settings_.num_threads) { #pragma omp master diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 7d50478333..3c0ffa54e9 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -153,7 +153,7 @@ class branch_and_bound_t { // Queue for storing the promising node for performing dives. omp_mutex_t mutex_dive_queue_; - diving_queue_t dive_queue_; + diving_queue_t diving_queue_; i_t min_diving_queue_size_; // Global status of the solver. @@ -192,7 +192,7 @@ class branch_and_bound_t { // Each "main" thread pops a node from the global heap and then performs a plunge // (i.e., a shallow dive) into the subtree determined by the node. - void best_first_thread(i_t id, + void best_first_thread(i_t task_id, search_tree_t& search_tree, const csr_matrix_t& Arow); diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 7b93ae37cf..ea900a0c62 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -19,7 +19,7 @@ namespace cuopt::linear_programming::dual_simplex { enum class node_status_t : int { - PENDING = 0, // Node still in the tree + PENDING = 0, // Node is still in the tree, waiting to be solved INTEGER_FEASIBLE = 1, // Node has an integer feasible solution INFEASIBLE = 2, // Node is infeasible FATHOMED = 3, // Node objective is greater than the upper bound @@ -273,11 +273,10 @@ class search_tree_t { void update(mip_node_t* node_ptr, node_status_t status) { - mutex.lock(); + std::lock_guard lock(mutex); std::vector*> stack; node_ptr->set_status(status, stack); remove_fathomed_nodes(stack); - mutex.unlock(); } void branch(mip_node_t* parent_node, From 3317eaa928a7072e8314d28818eb3fe68baf024b Mon Sep 17 00:00:00 2001 From: nicolas Date: Mon, 10 Nov 2025 15:40:50 +0100 Subject: [PATCH 37/51] fix missing initialization --- cpp/src/dual_simplex/presolve.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index fcc33a0e98..2cdb08f402 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -574,6 +574,7 @@ void convert_user_problem(const user_problem_t& user_problem, // Empty var_types means that all variables are continuous bounds_strengthening_t strengthening(problem, Arow, row_sense, {}); + std::fill(strengthening.bounds_changed.begin(), strengthening.bounds_changed.end(), true); strengthening.bounds_strengthening(problem.lower, problem.upper, settings); } From a0f8e13a2a4b3e8449242de624aebd7c38c812d9 Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 11 Nov 2025 13:48:59 +0100 Subject: [PATCH 38/51] fixed bugs after merge --- cpp/src/dual_simplex/basis_updates.cpp | 4 ++-- cpp/src/dual_simplex/branch_and_bound.cpp | 21 +++++++++++---------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 3e16411f47..2dceba2b39 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1864,7 +1864,7 @@ i_t basis_update_mpf_t::update(const std::vector& utilde, S_.x.data() + S_.col_start[S_start + 1], v_nz); - if (std::abs(mu) < 1e-13) { + if (std::abs(mu) < 1e-8) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } @@ -1928,7 +1928,7 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde S_.i.data() + S_.col_start[S_start + 1], S_.x.data() + S_.col_start[S_start + 1], S_.col_start[S_start + 2] - S_.col_start[S_start + 1]); - if (std::abs(mu) < 1e-13) { + if (std::abs(mu) < 1e-8) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index b03229ce53..9c1e3c9d34 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -455,7 +455,7 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t::set_final_solution(mip_solution_t::explore_subtree(i_t task_id, solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; } - if (stats_.nodes_explored >= settings_.node_limit) { - status_ = mip_exploration_status_t::NODE_LIMIT; + if (exploration_stats_.nodes_explored >= settings_.node_limit) { + solver_status_ = mip_exploration_status_t::NODE_LIMIT; return; } @@ -1124,11 +1125,11 @@ template mip_status_t branch_and_bound_t::solve(mip_solution_t& solution) { logger_t log; - log.log = false; - log.log_prefix = settings_.log.log_prefix; - status_ = mip_exploration_status_t::UNSET; - stats_.nodes_unexplored = 0; - stats_.nodes_explored = 0; + log.log = false; + log.log_prefix = settings_.log.log_prefix; + solver_status_ = mip_exploration_status_t::UNSET; + exploration_stats_.nodes_unexplored = 0; + exploration_stats_.nodes_explored = 0; if (guess_.size() != 0) { std::vector crushed_guess; From 7077dbc4bfe09bb1776702acce9b4980ad8602fc Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 11 Nov 2025 14:13:16 +0100 Subject: [PATCH 39/51] revert parameters to their original values --- cpp/src/dual_simplex/basis_updates.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 2dceba2b39..3e16411f47 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1864,7 +1864,7 @@ i_t basis_update_mpf_t::update(const std::vector& utilde, S_.x.data() + S_.col_start[S_start + 1], v_nz); - if (std::abs(mu) < 1e-8) { + if (std::abs(mu) < 1e-13) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } @@ -1928,7 +1928,7 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde S_.i.data() + S_.col_start[S_start + 1], S_.x.data() + S_.col_start[S_start + 1], S_.col_start[S_start + 2] - S_.col_start[S_start + 1]); - if (std::abs(mu) < 1e-8) { + if (std::abs(mu) < 1e-13) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } From d9a45cf15216763a9a28e52f0a30ff2c34fbaf0c Mon Sep 17 00:00:00 2001 From: nicolas Date: Tue, 11 Nov 2025 22:46:00 +0100 Subject: [PATCH 40/51] fixing crashes --- cpp/src/dual_simplex/bounds_strengthening.cpp | 25 +++++++++++-------- cpp/src/dual_simplex/bounds_strengthening.hpp | 4 --- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/cpp/src/dual_simplex/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp index 64e745d5f5..f1bf52c1e3 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.cpp +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -98,17 +98,20 @@ bool bounds_strengthening_t::bounds_strengthening( const i_t m = A.m; const i_t n = A.n; - constraint_changed.assign(m, false); - variable_changed.assign(n, false); - constraint_changed_next.assign(m, false); - - for (i_t i = 0; i < bounds_changed.size(); ++i) { - if (bounds_changed[i]) { - const i_t row_start = A.col_start[i]; - const i_t row_end = A.col_start[i + 1]; - for (i_t p = row_start; p < row_end; ++p) { - const i_t j = A.i[p]; - constraint_changed[j] = true; + std::vector constraint_changed(m, true); + std::vector variable_changed(n, false); + std::vector constraint_changed_next(m, false); + + if (!bounds_changed.empty()) { + std::fill(constraint_changed.begin(), constraint_changed.end(), false); + for (i_t i = 0; i < n; ++i) { + if (bounds_changed[i]) { + const i_t row_start = A.col_start[i]; + const i_t row_end = A.col_start[i + 1]; + for (i_t p = row_start; p < row_end; ++p) { + const i_t j = A.i[p]; + constraint_changed[j] = true; + } } } } diff --git a/cpp/src/dual_simplex/bounds_strengthening.hpp b/cpp/src/dual_simplex/bounds_strengthening.hpp index 8e38cde15c..e7e218b824 100644 --- a/cpp/src/dual_simplex/bounds_strengthening.hpp +++ b/cpp/src/dual_simplex/bounds_strengthening.hpp @@ -31,10 +31,6 @@ class bounds_strengthening_t { const csr_matrix_t& Arow; const std::vector& var_types; - std::vector constraint_changed; - std::vector variable_changed; - std::vector constraint_changed_next; - std::vector lower; std::vector upper; From 58c70acfe4cdc847c750b369f7937d65a32e5d9e Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 12 Nov 2025 14:45:42 +0100 Subject: [PATCH 41/51] disable RINS logs --- cpp/src/mip/diversity/lns/rins.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/src/mip/diversity/lns/rins.cu b/cpp/src/mip/diversity/lns/rins.cu index 1efc971b21..1a79bdb053 100644 --- a/cpp/src/mip/diversity/lns/rins.cu +++ b/cpp/src/mip/diversity/lns/rins.cu @@ -244,6 +244,7 @@ void rins_t::run_rins() branch_and_bound_settings.num_threads = 2; branch_and_bound_settings.num_bfs_threads = 1; branch_and_bound_settings.num_diving_threads = 1; + branch_and_bound_settings.log.log = false; branch_and_bound_settings.log.log_prefix = "[RINS] "; branch_and_bound_settings.solution_callback = [this, &rins_solution_queue]( std::vector& solution, f_t objective) { From a881b70bd4f9704ef766fd738c7629b1e0fa5ec2 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 12 Nov 2025 14:51:15 +0100 Subject: [PATCH 42/51] tighten tolerance for refactoring --- cpp/src/dual_simplex/basis_updates.cpp | 4 ++-- cpp/src/dual_simplex/basis_updates.hpp | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 3e16411f47..e583691e9a 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1864,7 +1864,7 @@ i_t basis_update_mpf_t::update(const std::vector& utilde, S_.x.data() + S_.col_start[S_start + 1], v_nz); - if (std::abs(mu) < 1e-13) { + if (std::abs(mu) < mu_tolerance_) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } @@ -1928,7 +1928,7 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde S_.i.data() + S_.col_start[S_start + 1], S_.x.data() + S_.col_start[S_start + 1], S_.col_start[S_start + 2] - S_.col_start[S_start + 1]); - if (std::abs(mu) < 1e-13) { + if (std::abs(mu) < mu_tolerance_) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index 5fc0f99eba..ca0451b919 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -424,6 +424,7 @@ class basis_update_mpf_t { void l_multiply(std::vector& inout) const; void l_transpose_multiply(std::vector& inout) const; + static constexpr f_t mu_tolerance_ = 1e-10; i_t num_updates_; // Number of rank-1 updates to L0 i_t refactor_frequency_; // Average updates before refactoring mutable csc_matrix_t L0_; // Sparse lower triangular matrix from initial factorization From 4947fe17623a5ef04573de6086d2e09ab19a0c35 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 12 Nov 2025 15:53:21 +0100 Subject: [PATCH 43/51] further tighten the tolerance --- cpp/src/dual_simplex/basis_updates.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index ca0451b919..1449145dac 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -424,7 +424,7 @@ class basis_update_mpf_t { void l_multiply(std::vector& inout) const; void l_transpose_multiply(std::vector& inout) const; - static constexpr f_t mu_tolerance_ = 1e-10; + static constexpr f_t mu_tolerance_ = 1e-8; i_t num_updates_; // Number of rank-1 updates to L0 i_t refactor_frequency_; // Average updates before refactoring mutable csc_matrix_t L0_; // Sparse lower triangular matrix from initial factorization From 5e3cd26f7e5d3be893d146cab76c0482229b6b3e Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 14 Nov 2025 12:55:23 +0100 Subject: [PATCH 44/51] changed refactor tolerances. updated logger to support different modes. --- cpp/src/dual_simplex/basis_updates.cpp | 5 +++-- cpp/src/dual_simplex/basis_updates.hpp | 1 - cpp/src/dual_simplex/branch_and_bound.cpp | 20 ++++++++++++++++++++ cpp/src/dual_simplex/logger.hpp | 4 ++-- 4 files changed, 25 insertions(+), 5 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index e583691e9a..fff33b8a58 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1864,7 +1864,7 @@ i_t basis_update_mpf_t::update(const std::vector& utilde, S_.x.data() + S_.col_start[S_start + 1], v_nz); - if (std::abs(mu) < mu_tolerance_) { + if (std::abs(mu) < 1E-8 || std::abs(mu) > 1E+8) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } @@ -1928,7 +1928,8 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde S_.i.data() + S_.col_start[S_start + 1], S_.x.data() + S_.col_start[S_start + 1], S_.col_start[S_start + 2] - S_.col_start[S_start + 1]); - if (std::abs(mu) < mu_tolerance_) { + + if (std::abs(mu) < 1E-8 || std::abs(mu) > 1E+8) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index c1c86a491f..078dfffeb3 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -424,7 +424,6 @@ class basis_update_mpf_t { void l_multiply(std::vector& inout) const; void l_transpose_multiply(std::vector& inout) const; - static constexpr f_t mu_tolerance_ = 1e-8; i_t num_updates_; // Number of rank-1 updates to L0 i_t refactor_frequency_; // Average updates before refactoring mutable csc_matrix_t L0_; // Sparse lower triangular matrix from initial factorization diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 9c1e3c9d34..110b13265b 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -599,6 +599,26 @@ node_children_status_t branch_and_bound_t::solve_node( lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); +#ifdef LOG_NODE_SIMPLEX + lp_settings.set_log(true); + std::stringstream ss; + ss << "simplex-" << std::this_thread::get_id() << ".log"; + std::string logname; + ss >> logname; + lp_settings.set_log_filename(logname); + lp_settings.log.enable_log_to_file("a+"); + lp_settings.log.log_to_console = false; + lp_settings.log.printf( + "%s node id = %d, branch var = %d, fractional val = %f, variable lower bound = %f, variable " + "upper bound = %f\n", + settings_.log.log_prefix.c_str(), + node_ptr->node_id, + node_ptr->branch_var, + node_ptr->fractional_val, + node_ptr->branch_var_lower, + node_ptr->branch_var_upper); +#endif + // Reset the bound_changed markers std::fill(node_presolver.bounds_changed.begin(), node_presolver.bounds_changed.end(), false); diff --git a/cpp/src/dual_simplex/logger.hpp b/cpp/src/dual_simplex/logger.hpp index 3b13665b1b..ac5e394f9e 100644 --- a/cpp/src/dual_simplex/logger.hpp +++ b/cpp/src/dual_simplex/logger.hpp @@ -30,10 +30,10 @@ class logger_t { { } - void enable_log_to_file() + void enable_log_to_file(std::string mode = "w") { if (log_file != nullptr) { std::fclose(log_file); } - log_file = std::fopen(log_filename.c_str(), "w"); + log_file = std::fopen(log_filename.c_str(), mode.c_str()); log_to_file = true; } From 74ef1ba6a69879fd17eb67ad8a19facc63eb6f14 Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 12 Nov 2025 14:55:39 +0100 Subject: [PATCH 45/51] tighten tolerance for refactoring the basis. disabled RINS logs. --- cpp/src/dual_simplex/basis_updates.cpp | 5 +++-- cpp/src/mip/diversity/lns/rins.cu | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 3e16411f47..fff33b8a58 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1864,7 +1864,7 @@ i_t basis_update_mpf_t::update(const std::vector& utilde, S_.x.data() + S_.col_start[S_start + 1], v_nz); - if (std::abs(mu) < 1e-13) { + if (std::abs(mu) < 1E-8 || std::abs(mu) > 1E+8) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } @@ -1928,7 +1928,8 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde S_.i.data() + S_.col_start[S_start + 1], S_.x.data() + S_.col_start[S_start + 1], S_.col_start[S_start + 2] - S_.col_start[S_start + 1]); - if (std::abs(mu) < 1e-13) { + + if (std::abs(mu) < 1E-8 || std::abs(mu) > 1E+8) { // Force a refactor. Otherwise we will get numerical issues when dividing by mu. return 1; } diff --git a/cpp/src/mip/diversity/lns/rins.cu b/cpp/src/mip/diversity/lns/rins.cu index 1efc971b21..1a79bdb053 100644 --- a/cpp/src/mip/diversity/lns/rins.cu +++ b/cpp/src/mip/diversity/lns/rins.cu @@ -244,6 +244,7 @@ void rins_t::run_rins() branch_and_bound_settings.num_threads = 2; branch_and_bound_settings.num_bfs_threads = 1; branch_and_bound_settings.num_diving_threads = 1; + branch_and_bound_settings.log.log = false; branch_and_bound_settings.log.log_prefix = "[RINS] "; branch_and_bound_settings.solution_callback = [this, &rins_solution_queue]( std::vector& solution, f_t objective) { From 1bcc57f694a5a5dbc9543ef8bb2241870b7c2afa Mon Sep 17 00:00:00 2001 From: nicolas Date: Wed, 19 Nov 2025 21:46:32 +0100 Subject: [PATCH 46/51] added missing workspace cleaning --- cpp/src/dual_simplex/basis_updates.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index 078dfffeb3..cea907074c 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -391,7 +391,11 @@ class basis_update_mpf_t { mu_values_.clear(); mu_values_.reserve(refactor_frequency_); num_updates_ = 0; + + std::fill(xi_workspace_.begin(), xi_workspace_.end(), 0); + std::fill(x_workspace_.begin(), x_workspace_.end(), 0.0); } + void grow_storage(i_t nz, i_t& S_start, i_t& S_nz); i_t index_map(i_t leaving) const; f_t u_diagonal(i_t j) const; From f687ce49cd54c1d1a99622e4cf2aef20c5f986c8 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 20 Nov 2025 10:03:55 +0100 Subject: [PATCH 47/51] small refactor --- cpp/src/dual_simplex/branch_and_bound.cpp | 90 +++++++++++------------ cpp/src/dual_simplex/branch_and_bound.hpp | 32 ++++---- cpp/src/dual_simplex/presolve.cpp | 2 + 3 files changed, 63 insertions(+), 61 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 110b13265b..eecafb14e4 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -200,10 +200,10 @@ inline const char* feasible_solution_symbol(thread_type_t type) } } -inline bool has_children(node_children_status_t status) +inline bool has_children(node_solve_info_t status) { - return status == node_children_status_t::UP_CHILDREN_FIRST || - status == node_children_status_t::DOWN_CHILDREN_FIRST; + return status == node_solve_info_t::UP_CHILD_FIRST || + status == node_solve_info_t::DOWN_CHILD_FIRST; } } // namespace @@ -575,7 +575,7 @@ rounding_direction_t branch_and_bound_t::child_selection(mip_node_t -node_children_status_t branch_and_bound_t::solve_node( +node_solve_info_t branch_and_bound_t::solve_node( mip_node_t* node_ptr, search_tree_t& search_tree, lp_problem_t& leaf_problem, @@ -672,7 +672,7 @@ node_children_status_t branch_and_bound_t::solve_node( node_ptr->lower_bound = inf; search_tree.graphviz_node(log, node_ptr, "infeasible", 0.0); search_tree.update(node_ptr, node_status_t::INFEASIBLE); - return node_children_status_t::NO_CHILDREN; + return node_solve_info_t::NO_CHILDREN; } else if (lp_status == dual::status_t::CUTOFF) { // Node was cut off. Do not branch @@ -680,7 +680,7 @@ node_children_status_t branch_and_bound_t::solve_node( f_t leaf_objective = compute_objective(leaf_problem, leaf_solution.x); search_tree.graphviz_node(log, node_ptr, "cut off", leaf_objective); search_tree.update(node_ptr, node_status_t::FATHOMED); - return node_children_status_t::NO_CHILDREN; + return node_solve_info_t::NO_CHILDREN; } else if (lp_status == dual::status_t::OPTIMAL) { // LP was feasible @@ -704,7 +704,7 @@ node_children_status_t branch_and_bound_t::solve_node( add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, thread_type); search_tree.graphviz_node(log, node_ptr, "integer feasible", leaf_objective); search_tree.update(node_ptr, node_status_t::INTEGER_FEASIBLE); - return node_children_status_t::NO_CHILDREN; + return node_solve_info_t::NO_CHILDREN; } else if (leaf_objective <= upper_bound + abs_fathom_tol) { // Choose fractional variable to branch on @@ -719,19 +719,19 @@ node_children_status_t branch_and_bound_t::solve_node( rounding_direction_t round_dir = child_selection(node_ptr); if (round_dir == rounding_direction_t::UP) { - return node_children_status_t::UP_CHILDREN_FIRST; + return node_solve_info_t::UP_CHILD_FIRST; } else { - return node_children_status_t::DOWN_CHILDREN_FIRST; + return node_solve_info_t::DOWN_CHILD_FIRST; } } else { search_tree.graphviz_node(log, node_ptr, "fathomed", leaf_objective); search_tree.update(node_ptr, node_status_t::FATHOMED); - return node_children_status_t::NO_CHILDREN; + return node_solve_info_t::NO_CHILDREN; } } else if (lp_status == dual::status_t::TIME_LIMIT) { search_tree.graphviz_node(log, node_ptr, "timeout", 0.0); - return node_children_status_t::TIME_LIMIT; + return node_solve_info_t::TIME_LIMIT; } else { if (thread_type == thread_type_t::EXPLORATION) { @@ -747,7 +747,7 @@ node_children_status_t branch_and_bound_t::solve_node( search_tree.graphviz_node(log, node_ptr, "numerical", 0.0); search_tree.update(node_ptr, node_status_t::NUMERICAL); - return node_children_status_t::NUMERICAL; + return node_solve_info_t::NUMERICAL; } } @@ -819,17 +819,17 @@ void branch_and_bound_t::exploration_ramp_up(mip_node_t* nod std::vector row_sense; bounds_strengthening_t node_presolver(leaf_problem, Arow, row_sense, var_types_); - node_children_status_t status = solve_node(node, - *search_tree, - leaf_problem, - node_presolver, - thread_type_t::EXPLORATION, - true, - original_lp_.lower, - original_lp_.upper, - settings_.log); - - if (status == node_children_status_t::TIME_LIMIT) { + node_solve_info_t status = solve_node(node, + *search_tree, + leaf_problem, + node_presolver, + thread_type_t::EXPLORATION, + true, + original_lp_.lower, + original_lp_.upper, + settings_.log); + + if (status == node_solve_info_t::TIME_LIMIT) { solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; @@ -932,19 +932,19 @@ void branch_and_bound_t::explore_subtree(i_t task_id, return; } - node_children_status_t status = solve_node(node_ptr, - search_tree, - leaf_problem, - node_presolver, - thread_type_t::EXPLORATION, - recompute_bounds, - original_lp_.lower, - original_lp_.upper, - settings_.log); + node_solve_info_t status = solve_node(node_ptr, + search_tree, + leaf_problem, + node_presolver, + thread_type_t::EXPLORATION, + recompute_bounds, + original_lp_.lower, + original_lp_.upper, + settings_.log); recompute_bounds = !has_children(status); - if (status == node_children_status_t::TIME_LIMIT) { + if (status == node_solve_info_t::TIME_LIMIT) { solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; @@ -980,7 +980,7 @@ void branch_and_bound_t::explore_subtree(i_t task_id, exploration_stats_.nodes_unexplored += 2; - if (status == node_children_status_t::UP_CHILDREN_FIRST) { + if (status == node_solve_info_t::UP_CHILD_FIRST) { stack.push_front(node_ptr->get_down_child()); stack.push_front(node_ptr->get_up_child()); } else { @@ -1092,24 +1092,24 @@ void branch_and_bound_t::diving_thread(const csr_matrix_t& A if (toc(exploration_stats_.start_time) > settings_.time_limit) { return; } - node_children_status_t status = solve_node(node_ptr, - subtree, - leaf_problem, - node_presolver, - thread_type_t::DIVING, - recompute_bounds, - start_node->lower, - start_node->upper, - log); + node_solve_info_t status = solve_node(node_ptr, + subtree, + leaf_problem, + node_presolver, + thread_type_t::DIVING, + recompute_bounds, + start_node->lower, + start_node->upper, + log); recompute_bounds = !has_children(status); - if (status == node_children_status_t::TIME_LIMIT) { + if (status == node_solve_info_t::TIME_LIMIT) { solver_status_ = mip_exploration_status_t::TIME_LIMIT; return; } else if (has_children(status)) { - if (status == node_children_status_t::UP_CHILDREN_FIRST) { + if (status == node_solve_info_t::UP_CHILD_FIRST) { stack.push_front(node_ptr->get_down_child()); stack.push_front(node_ptr->get_up_child()); } else { diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 3c0ffa54e9..6ad5b3d81d 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -42,13 +42,13 @@ enum class mip_exploration_status_t { COMPLETED = 5, // The solver finished exploring the tree }; -enum class node_children_status_t { - NO_CHILDREN = 0, // The node does not produced children - UP_CHILDREN_FIRST = 1, // The up child should be explored first - DOWN_CHILDREN_FIRST = 2, // The down child should be explored first - TIME_LIMIT = 3, // The solver reached a time limit - ITERATION_LIMIT = 4, // The solver reached a iteration limit - NUMERICAL = 5 // The solver encounter a numerical error when solving the node +enum class node_solve_info_t { + NO_CHILDREN = 0, // The node does not produced children + UP_CHILD_FIRST = 1, // The up child should be explored first + DOWN_CHILD_FIRST = 2, // The down child should be explored first + TIME_LIMIT = 3, // The solver reached a time limit + ITERATION_LIMIT = 4, // The solver reached a iteration limit + NUMERICAL = 5 // The solver encounter a numerical error when solving the node }; // Indicate the search and variable selection algorithms used by the thread (See [1]). @@ -201,15 +201,15 @@ class branch_and_bound_t { void diving_thread(const csr_matrix_t& Arow); // Solve the LP relaxation of a leaf node and update the tree. - node_children_status_t solve_node(mip_node_t* node_ptr, - search_tree_t& search_tree, - lp_problem_t& leaf_problem, - bounds_strengthening_t& node_presolver, - thread_type_t thread_type, - bool recompute, - const std::vector& root_lower, - const std::vector& root_upper, - logger_t& log); + node_solve_info_t solve_node(mip_node_t* node_ptr, + search_tree_t& search_tree, + lp_problem_t& leaf_problem, + bounds_strengthening_t& node_presolver, + thread_type_t thread_type, + bool recompute, + const std::vector& root_lower, + const std::vector& root_upper, + logger_t& log); // Sort the children based on the Martin's criteria. rounding_direction_t child_selection(mip_node_t* node_ptr); diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 10e7d6358e..5cfe252085 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -572,6 +572,8 @@ void convert_user_problem(const user_problem_t& user_problem, csr_matrix_t Arow(1, 1, 1); problem.A.to_compressed_row(Arow); + settings.log.printf("Running bound strengthening\n"); + // Empty var_types means that all variables are continuous bounds_strengthening_t strengthening(problem, Arow, row_sense, {}); std::fill(strengthening.bounds_changed.begin(), strengthening.bounds_changed.end(), true); From 1048327221206c3b2def98c7eec11dcad4b15b5c Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 20 Nov 2025 14:21:35 +0100 Subject: [PATCH 48/51] fixed crash due to selecting the incorrect last column --- cpp/src/dual_simplex/basis_updates.cpp | 6 +++++- cpp/src/dual_simplex/sparse_matrix.cpp | 6 ++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index fff33b8a58..a10c7ab42b 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1193,7 +1193,7 @@ i_t basis_update_mpf_t::scatter_into_workspace(const sparse_vector_t void basis_update_mpf_t::grow_storage(i_t nz, i_t& S_start, i_t& S_nz) { - const i_t last_S_col = num_updates_ * 2; + const i_t last_S_col = S_.n; const i_t new_last_S_col = last_S_col + 2; if (new_last_S_col >= S_.col_start.size()) { S_.col_start.resize(new_last_S_col + refactor_frequency_); @@ -1204,6 +1204,8 @@ void basis_update_mpf_t::grow_storage(i_t nz, i_t& S_start, i_t& S_nz) S_.x.resize(std::max(2 * S_nz, S_nz + nz)); } S_start = last_S_col; + assert(S_nz + nz <= S_.i.size()); + assert(S_nz + nz <= S_.x.size()); } template @@ -1915,6 +1917,8 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde i_t S_nz; grow_storage(nz + etilde.i.size(), S_start, S_nz); + if (nz + etilde.i.size() > S_.i.size()) {} + S_.append_column(nz, xi_workspace_.data() + m, x_workspace_.data()); // Gather etilde into a column of S diff --git a/cpp/src/dual_simplex/sparse_matrix.cpp b/cpp/src/dual_simplex/sparse_matrix.cpp index cdd45f720e..4e522d9daf 100644 --- a/cpp/src/dual_simplex/sparse_matrix.cpp +++ b/cpp/src/dual_simplex/sparse_matrix.cpp @@ -157,6 +157,8 @@ void csc_matrix_t::append_column(const std::vector& x) i_t nz = this->col_start[this->n]; for (i_t j = 0; j < xsz; ++j) { if (x[j] != 0.0) { + assert(nz < this->i.size()); + assert(nz < this->x.size()); this->i[nz] = j; this->x[nz] = x[j]; nz++; @@ -177,6 +179,8 @@ void csc_matrix_t::append_column(const sparse_vector_t& x) const i_t i = x.i[k]; const f_t x_val = x.x[k]; if (x_val != 0.0) { + assert(nz < this->i.size()); + assert(nz < this->x.size()); this->i[nz] = i; this->x[nz] = x_val; nz++; @@ -194,6 +198,8 @@ void csc_matrix_t::append_column(i_t x_nz, i_t* i, f_t* x) const i_t i_val = i[k]; const f_t x_val = x[i_val]; if (x_val != 0.0) { + assert(nz < this->i.size()); + assert(nz < this->x.size()); this->i[nz] = i_val; this->x[nz] = x_val; nz++; From 8b4546eb8d1a0c5720df5ee7c5a0e485c7f1d46b Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 20 Nov 2025 17:14:52 +0100 Subject: [PATCH 49/51] reverting fix --- cpp/src/dual_simplex/basis_updates.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index a10c7ab42b..2767da65de 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1193,7 +1193,8 @@ i_t basis_update_mpf_t::scatter_into_workspace(const sparse_vector_t void basis_update_mpf_t::grow_storage(i_t nz, i_t& S_start, i_t& S_nz) { - const i_t last_S_col = S_.n; + const i_t last_S_col = 2 * num_updates_; + assert(S_.n == last_S_col); const i_t new_last_S_col = last_S_col + 2; if (new_last_S_col >= S_.col_start.size()) { S_.col_start.resize(new_last_S_col + refactor_frequency_); From 3240f889645cf381b5572b8fc5e28c512ddadc52 Mon Sep 17 00:00:00 2001 From: nicolas Date: Thu, 20 Nov 2025 18:13:30 +0100 Subject: [PATCH 50/51] fixed missing vector if the matrix is rank deficient and there is no singleton --- cpp/src/dual_simplex/basis_solves.cpp | 15 ++++++++++++++- cpp/src/dual_simplex/basis_updates.cpp | 2 ++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index d0f82967d4..db24f55a25 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -566,6 +566,20 @@ i_t factorize_basis(const csc_matrix_t& A, q.resize(m); f_t fact_start = tic(); rank = right_looking_lu(A, settings, medium_tol, basic_list, q, L, U, pinv); + inverse_permutation(pinv, p); + if (rank != m) { + // Get the rank deficient columns + deficient.clear(); + deficient.resize(m - rank); + for (i_t h = rank; h < m; ++h) { + deficient[h - rank] = q[h]; + } + // Get the slacks needed + slacks_needed.resize(m - rank); + for (i_t h = rank; h < m; ++h) { + slacks_needed[h - rank] = p[h]; + } + } if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { settings.log.printf("Concurrent halt\n"); return -1; @@ -573,7 +587,6 @@ i_t factorize_basis(const csc_matrix_t& A, if (verbose) { printf("Right Lnz+Unz %d t %.3f\n", L.col_start[m] + U.col_start[m], toc(fact_start)); } - inverse_permutation(pinv, p); constexpr bool check_lu = false; if (check_lu) { csc_matrix_t C(m, m, 1); diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 2767da65de..931c9302f8 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -2100,6 +2100,8 @@ int basis_update_mpf_t::refactor_basis( #ifdef CHECK_L_FACTOR if (L0_.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); } #endif + + assert(deficient.size() > 0); return deficient.size(); } settings.log.debug("Basis repaired\n"); From a90b5addff0c25628891e99f3cf26d44a25d506f Mon Sep 17 00:00:00 2001 From: nicolas Date: Fri, 21 Nov 2025 09:12:04 +0100 Subject: [PATCH 51/51] removing debug leftover --- cpp/src/dual_simplex/basis_updates.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 200d4b219b..6b79f3c862 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1918,8 +1918,6 @@ i_t basis_update_mpf_t::update(const sparse_vector_t& utilde i_t S_nz; grow_storage(nz + etilde.i.size(), S_start, S_nz); - if (nz + etilde.i.size() > S_.i.size()) {} - S_.append_column(nz, xi_workspace_.data() + m, x_workspace_.data()); // Gather etilde into a column of S