diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index a376ee23d6..e8a9b5dce1 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -17,6 +17,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}/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/bounds_strengthening.cpp b/cpp/src/dual_simplex/bounds_strengthening.cpp new file mode 100644 index 0000000000..f1bf52c1e3 --- /dev/null +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -0,0 +1,290 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include +#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 +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(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 bounds_strengthening_t::bounds_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; + + 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; + } + } + } + } + + lower = lower_bounds; + 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.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.j[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 = 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 = A.i[p]; + + if (!constraint_changed[i]) { continue; } + 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]; + + 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 = 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]); + + 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 = 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 bounds_strengthening_t; +#endif + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/bounds_strengthening.hpp b/cpp/src/dual_simplex/bounds_strengthening.hpp new file mode 100644 index 0000000000..e7e218b824 --- /dev/null +++ b/cpp/src/dual_simplex/bounds_strengthening.hpp @@ -0,0 +1,42 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +class bounds_strengthening_t { + public: + // For pure LP bounds strengthening, var_types should be defaulted (i.e. left empty) + bounds_strengthening_t(const lp_problem_t& problem, + const csr_matrix_t& Arow, + const std::vector& row_sense, + const std::vector& var_types); + + 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 csr_matrix_t& Arow; + const std::vector& var_types; + + std::vector lower; + std::vector upper; + + 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/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 9f207b6a6f..eecafb14e4 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include #include @@ -190,6 +191,21 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound) } } +inline const char* feasible_solution_symbol(thread_type_t type) +{ + switch (type) { + case thread_type_t::EXPLORATION: return "B"; + case thread_type_t::DIVING: return "D"; + default: return "U"; + } +} + +inline bool has_children(node_solve_info_t status) +{ + return status == node_solve_info_t::UP_CHILD_FIRST || + status == node_solve_info_t::DOWN_CHILD_FIRST; +} + } // namespace template @@ -202,9 +218,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_); @@ -286,7 +302,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); @@ -296,11 +312,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)); } } @@ -405,11 +421,11 @@ 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(), - toc(stats_.start_time)); + toc(exploration_stats_.start_time)); if (settings_.solution_callback != nullptr) { std::vector original_x; @@ -430,16 +446,16 @@ 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) { - settings_.log.printf("Integer infeasible.\n"); - mip_status = mip_status_t::INFEASIBLE; - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); + 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) { + settings_.heuristic_preemption_callback(); + } } } @@ -489,8 +509,8 @@ mip_status_t branch_and_bound_t::set_final_solution(mip_solution_t 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; - 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_) { @@ -512,16 +532,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("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", - 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; } @@ -535,8 +556,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; @@ -547,44 +567,75 @@ 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(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char thread_type) +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, + 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) { - f_t abs_fathom_tol = settings_.absolute_mip_gap_tol / 10; + 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); - 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); - simplex_solver_settings_t lp_settings = settings_; 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); + +#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); + + // 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, node_presolver.bounds_changed); + + } else { + node_ptr->update_branched_variable_bounds( + leaf_problem.lower, leaf_problem.upper, node_presolver.bounds_changed); + } - // 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); + node_presolver.bounds_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -610,24 +661,26 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& 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; + 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) { // 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); - return node_status_t::INFEASIBLE; + search_tree.update(node_ptr, node_status_t::INFEASIBLE); + return node_solve_info_t::NO_CHILDREN; } else if (lp_status == dual::status_t::CUTOFF) { // Node was cut off. Do not branch 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); - return node_status_t::FATHOMED; + search_tree.update(node_ptr, node_status_t::FATHOMED); + return node_solve_info_t::NO_CHILDREN; } else if (lp_status == dual::status_t::OPTIMAL) { // LP was feasible @@ -650,8 +703,8 @@ 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); - return node_status_t::INTEGER_FEASIBLE; + search_tree.update(node_ptr, node_status_t::INTEGER_FEASIBLE); + return node_solve_info_t::NO_CHILDREN; } else if (leaf_objective <= upper_bound + abs_fathom_tol) { // Choose fractional variable to branch on @@ -660,22 +713,28 @@ 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; + + rounding_direction_t round_dir = child_selection(node_ptr); + + if (round_dir == rounding_direction_t::UP) { + return node_solve_info_t::UP_CHILD_FIRST; + } else { + return node_solve_info_t::DOWN_CHILD_FIRST; + } } else { search_tree.graphviz_node(log, node_ptr, "fathomed", leaf_objective); - search_tree.update_tree(node_ptr, node_status_t::FATHOMED); - return node_status_t::FATHOMED; + search_tree.update(node_ptr, node_status_t::FATHOMED); + 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); - search_tree.update_tree(node_ptr, node_status_t::TIME_LIMIT); - return node_status_t::TIME_LIMIT; + return node_solve_info_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 " @@ -687,19 +746,18 @@ 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); - return node_status_t::NUMERICAL; + search_tree.update(node_ptr, node_status_t::NUMERICAL); + return node_solve_info_t::NUMERICAL; } } 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, +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; } + 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 @@ -710,71 +768,81 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* 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); - search_tree->update_tree(node, node_status_t::FATHOMED); + search_tree->update(node, node_status_t::FATHOMED); 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; } - // Set the correct bounds for the leaf problem - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; - - node_status_t node_status = - solve_node(*search_tree, node, leaf_problem, Arow, upper_bound, settings_.log, 'B'); - - if (node_status == node_status_t::TIME_LIMIT) { - status_ = mip_exploration_status_t::TIME_LIMIT; + // 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 node_presolver(leaf_problem, Arow, row_sense, var_types_); + + 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; - } 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( - search_tree, node->get_down_child(), leaf_problem, Arow, 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(), leaf_problem, Arow, 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 @@ -788,15 +856,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, + search_tree_t& search_tree, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) + 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(); @@ -815,62 +884,71 @@ 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); - search_tree.update_tree(node_ptr, node_status_t::FATHOMED); + search_tree.update(node_ptr, node_status_t::FATHOMED); + recompute_bounds = true; 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; } - 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; } - // Set the correct bounds for the leaf problem - leaf_problem.lower = original_lp_.lower; - leaf_problem.upper = original_lp_.upper; + 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); - node_status_t node_status = - solve_node(search_tree, node_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + recompute_bounds = !has_children(status); - if (node_status == node_status_t::TIME_LIMIT) { - status_ = mip_exploration_status_t::TIME_LIMIT; + if (status == node_solve_info_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, @@ -886,8 +964,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, node_presolver.bounds_changed); + mutex_dive_queue_.lock(); - dive_queue_.emplace(node->detach_copy(), leaf_problem.lower, leaf_problem.upper); + diving_queue_.emplace(node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } @@ -896,52 +978,60 @@ 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_solve_info_t::UP_CHILD_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()); + } } } } 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, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) + const csr_matrix_t& Arow) { f_t lower_bound = -inf; f_t upper_bound = inf; f_t abs_gap = inf; f_t rel_gap = inf; - while (status_ == mip_exploration_status_t::RUNNING && abs_gap > settings_.absolute_mip_gap_tol && - rel_gap > settings_.relative_mip_gap_tol && + // 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 node_presolver(leaf_problem, Arow, row_sense, var_types_); + + 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* 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_tree(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, search_tree, node_ptr, leaf_problem, Arow); + explore_subtree(task_id, start_node, search_tree, leaf_problem, node_presolver); active_subtrees_--; } @@ -953,63 +1043,79 @@ void branch_and_bound_t::best_first_thread(i_t 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_[id] = inf; + local_lower_bounds_[task_id] = inf; } } } template -void branch_and_bound_t::diving_thread(lp_problem_t& leaf_problem, - const csc_matrix_t& Arow) +void branch_and_bound_t::diving_thread(const csr_matrix_t& Arow) { logger_t log; log.log = false; - while (status_ == mip_exploration_status_t::RUNNING && + // 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 node_presolver(leaf_problem, Arow, row_sense, var_types_); + + while (solver_status_ == mip_exploration_status_t::RUNNING && (active_subtrees_ > 0 || get_heap_size() > 0)) { 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()) { if (get_upper_bound() < start_node->node.lower_bound) { continue; } + bool recompute_bounds = true; search_tree_t subtree(std::move(start_node->node)); 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(); 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_bounds = true; continue; } - if (toc(stats_.start_time) > settings_.time_limit) { return; } + if (toc(exploration_stats_.start_time) > settings_.time_limit) { return; } - // Set the correct bounds for the leaf problem - leaf_problem.lower = start_node->lp_lower; - leaf_problem.upper = start_node->lp_upper; + 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); - node_status_t node_status = - solve_node(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); + recompute_bounds = !has_children(status); - if (node_status == node_status_t::TIME_LIMIT) { + if (status == node_solve_info_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_solve_info_t::UP_CHILD_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) { @@ -1017,11 +1123,16 @@ void branch_and_bound_t::diving_thread(lp_problem_t& leaf_pr // 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_) { - mutex_dive_queue_.lock(); + if (diving_queue_.size() < min_diving_queue_size_) { 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, node_presolver.bounds_changed); + + mutex_dive_queue_.lock(); + diving_queue_.emplace(new_node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } } @@ -1034,11 +1145,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; @@ -1060,12 +1171,17 @@ 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"); // FIXME: rarely dual simplex detects infeasible whereas it is feasible. @@ -1085,7 +1201,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); } @@ -1127,7 +1243,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); @@ -1141,7 +1257,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, @@ -1150,8 +1266,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_); } @@ -1167,33 +1283,29 @@ 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); + 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); settings_.log.printf( " | 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; - 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; #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_; - #pragma omp master { auto down_child = search_tree.root.get_down_child(); @@ -1201,27 +1313,24 @@ 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(down_child, &search_tree, Arow, initial_size); #pragma omp task - exploration_ramp_up(&search_tree, up_child, leaf_problem, Arow, initial_size); + exploration_ramp_up(up_child, &search_tree, Arow, initial_size); } #pragma omp barrier #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, leaf_problem, 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(leaf_problem, Arow); - } + diving_thread(Arow); } } } diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 5b304addd5..6ad5b3d81d 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -7,10 +7,10 @@ #pragma once +#include #include #include #include -#include #include #include #include @@ -42,69 +42,29 @@ enum class mip_exploration_status_t { COMPLETED = 5, // The solver finished exploring the tree }; -template -void upper_bound_callback(f_t upper_bound); +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 +}; -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; - } +// 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, }; -// 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; +class bounds_strengthening_t; - 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 +void upper_bound_callback(f_t upper_bound); template class branch_and_bound_t { @@ -169,7 +129,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_; @@ -193,11 +153,11 @@ 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 diving_queue_; 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. @@ -211,49 +171,48 @@ 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(); // 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, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, + 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. void explore_subtree(i_t task_id, - search_tree_t& search_tree, mip_node_t* start_node, + search_tree_t& search_tree, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow); + 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. - void best_first_thread(i_t id, + void best_first_thread(i_t task_id, search_tree_t& search_tree, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow); + 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(lp_problem_t& leaf_problem, const csc_matrix_t& Arow); + 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(search_tree_t& search_tree, - mip_node_t* node_ptr, - lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char thread_type); + 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. - 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/diving_queue.hpp b/cpp/src/dual_simplex/diving_queue.hpp new file mode 100644 index 0000000000..91d0aa3309 --- /dev/null +++ b/cpp/src/dual_simplex/diving_queue.hpp @@ -0,0 +1,73 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025 NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include +#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, std::vector&& lower, std::vector&& upper) + : node(std::move(node)), lower(std::move(lower)), upper(std::move(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 INT16_MAX, 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_ = INT16_MAX; + + 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() - 1) { buffer.pop_back(); } + } + + void emplace(mip_node_t&& node, std::vector&& lower, std::vector&& 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() - 1) { 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/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; } diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 9034bfa228..ea900a0c62 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -19,28 +19,29 @@ namespace cuopt::linear_programming::dual_simplex { enum class node_status_t : int { - ACTIVE = 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 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 : int8_t { NONE = -1, DOWN = 0, UP = 1 }; + bool inactive_status(node_status_t status); 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(-1), + branch_dir(rounding_direction_t::NONE), vstatus(basis) { children[0] = nullptr; @@ -51,10 +52,10 @@ class mip_node_t { mip_node_t* parent_node, i_t node_num, i_t branch_variable, - i_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), @@ -65,38 +66,49 @@ 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 == 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; } void get_variable_bounds(std::vector& lower, 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_branched_variable_bounds(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; + while (parent_ptr != nullptr && parent_ptr->node_id != 0) { + parent_ptr->update_branched_variable_bounds(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_branched_variable_bounds(std::vector& lower, + std::vector& upper, + std::vector& bounds_changed) const + { + 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 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 + 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(); } @@ -193,9 +205,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); @@ -213,7 +224,7 @@ class mip_node_t { i_t depth; i_t node_id; i_t branch_var; - i_t branch_dir; + rounding_direction_t branch_dir; f_t branch_var_lower; f_t branch_var_upper; f_t fractional_val; @@ -260,13 +271,12 @@ 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::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, @@ -278,17 +288,35 @@ 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)); - - // up child - auto up_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, 1, fractional_val, parent_vstatus); - - graphviz_edge(log, parent_node, up_child.get(), branch_var, 1, std::ceil(fractional_val)); + auto down_child = std::make_unique>(original_lp, + parent_node, + ++id, + branch_var, + rounding_direction_t::DOWN, + fractional_val, + parent_vstatus); + + graphviz_edge(log, + parent_node, + down_child.get(), + branch_var, + rounding_direction_t::DOWN, + std::floor(fractional_val)); + + 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, + rounding_direction_t::UP, + std::ceil(fractional_val)); assert(parent_vstatus.size() == original_lp.num_cols); parent_node->add_children(std::move(down_child), @@ -309,7 +337,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, + rounding_direction_t branch_dir, const f_t bound) { if (write_graphviz) { @@ -317,7 +345,7 @@ class search_tree_t { origin_ptr->node_id, dest_ptr->node_id, branch_var, - branch_dir == 0 ? "<=" : ">=", + branch_dir == rounding_direction_t::DOWN ? "<=" : ">=", bound); } } diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 6c55f56236..5cfe252085 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -7,6 +7,7 @@ #include +#include #include #include #include @@ -18,272 +19,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, @@ -832,15 +567,19 @@ 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) { + constexpr bool run_bounds_strengthening = false; + if (run_bounds_strengthening) { + csr_matrix_t Arow(1, 1, 1); + problem.A.to_compressed_row(Arow); + 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); + + // 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); } + 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 && @@ -1601,13 +1340,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 4f08d0b75a..d9f85fdae0 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -180,13 +180,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 diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index ca3e58041a..70150282ff 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -200,10 +200,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 == 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 == 0) { + 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 {