diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index c091ebac4f..fb063efbe7 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}/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..9f92062e8f --- /dev/null +++ b/cpp/src/dual_simplex/bounds_strengthening.cpp @@ -0,0 +1,297 @@ +/* + * 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 + +#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; + + 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; + } + } + } + + 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..28df35d0dc --- /dev/null +++ b/cpp/src/dual_simplex/bounds_strengthening.hpp @@ -0,0 +1,56 @@ +/* + * 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 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 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; + 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 2ce3ee0b4e..923ffaadc9 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 @@ -200,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 @@ -279,6 +289,7 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); if (is_feasible) { upper_bound_ = obj; + incumbent_.set_incumbent_solution(obj, crushed_solution); } else { attempt_repair = true; constexpr bool verbose = false; @@ -414,7 +425,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(), @@ -449,15 +460,20 @@ mip_status_t branch_and_bound_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 (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(); + } } } - uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); + if (upper_bound != inf) { + assert(incumbent_.has_incumbent); + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); + } + solution.objective = incumbent_.objective; solution.lower_bound = lower_bound; solution.nodes_explored = stats_.nodes_explored; @@ -496,7 +518,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; @@ -509,8 +531,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, @@ -552,36 +574,54 @@ 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, - const csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char thread_type) + basis_update_mpf_t& ft, + std::vector& basic_list, + std::vector& nonbasic_list, + bounds_strengthening_t& presolver, + thread_type_t thread_type, + bool recompute_bounds_and_basis, + const std::vector& root_lower, + const std::vector& root_upper, + stats_t& stats, + 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); - // in B&B we only have equality constraints, leave it empty for default - std::vector row_sense; + // Limit the number of simplex iterations when diving. + if (thread_type == thread_type_t::DIVING) { + lp_settings.iteration_limit = 0.05 * stats_.total_lp_iters - stats.total_lp_iters; + } + + // 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_bounds_and_basis) { + 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 = - bound_strengthening(row_sense, lp_settings, leaf_problem, Arow, var_types_, bounds_changed); + presolver.bounds_strengthening(leaf_problem.lower, leaf_problem.upper, lp_settings); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; @@ -590,32 +630,44 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& f_t lp_start_time = tic(); std::vector leaf_edge_norms = edge_norms_; // = node.steepest_edge_norms; - lp_status = dual_phase2(2, - 0, - lp_start_time, - leaf_problem, - lp_settings, - leaf_vstatus, - leaf_solution, - node_iter, - leaf_edge_norms); + lp_status = dual_phase2_with_advanced_basis(2, + 0, + recompute_bounds_and_basis, + lp_start_time, + leaf_problem, + lp_settings, + leaf_vstatus, + ft, + basic_list, + nonbasic_list, + leaf_solution, + node_iter, + leaf_edge_norms); if (lp_status == dual::status_t::NUMERICAL) { - 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); + log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); + lp_status_t second_status = solve_linear_program_with_advanced_basis(leaf_problem, + lp_start_time, + lp_settings, + leaf_solution, + ft, + basic_list, + nonbasic_list, + 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; } 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) { @@ -623,7 +675,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) { @@ -641,32 +693,38 @@ 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); + + node_ptr->objective_estimate = + pc_.objective_estimate(leaf_fractional, leaf_solution.x, leaf_objective, pc_log); 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; } 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); return node_status_t::TIME_LIMIT; + } else if (lp_status == dual::status_t::ITERATION_LIMIT) { + return node_status_t::ITERATION_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 " @@ -678,16 +736,15 @@ 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; } } 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; } @@ -707,7 +764,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; } @@ -744,12 +801,29 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* 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'); + // 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_); + + const i_t m = leaf_problem.num_rows; + basis_update_mpf_t basis_update(m, settings_.refactor_frequency); + std::vector basic_list(m); + std::vector nonbasic_list; + + node_status_t node_status = solve_node(node, + *search_tree, + leaf_problem, + basis_update, + basic_list, + nonbasic_list, + presolver, + thread_type_t::EXPLORATION, + true, + original_lp_.lower, + original_lp_.upper, + stats_, + settings_.log); if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -761,11 +835,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(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 @@ -779,11 +852,15 @@ 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& presolver, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list) { + bool recompute_bounds_and_basis = true; std::deque*> stack; stack.push_front(start_node); @@ -812,7 +889,8 @@ 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_tree(node_ptr, node_status_t::FATHOMED); + search_tree.update(node_ptr, node_status_t::FATHOMED); + recompute_bounds_and_basis = true; continue; } @@ -846,12 +924,21 @@ void branch_and_bound_t::explore_subtree(i_t task_id, 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_ptr, leaf_problem, Arow, upper_bound, settings_.log, 'B'); + node_status_t node_status = solve_node(node_ptr, + search_tree, + leaf_problem, + basis_update, + basic_list, + nonbasic_list, + presolver, + thread_type_t::EXPLORATION, + recompute_bounds_and_basis, + original_lp_.lower, + original_lp_.upper, + stats_, + settings_.log); + + recompute_bounds_and_basis = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_status_t::TIME_LIMIT) { status_ = mip_exploration_status_t::TIME_LIMIT; @@ -873,8 +960,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); + diving_queue_.emplace(node->detach_copy(), std::move(lower), std::move(upper)); mutex_dive_queue_.unlock(); } @@ -893,42 +984,59 @@ 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, - 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; + // 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_); + + const i_t m = leaf_problem.num_rows; + basis_update_mpf_t basis_update(m, settings_.refactor_frequency); + std::vector basic_list(m); + std::vector nonbasic_list; + 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_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, + presolver, + basis_update, + basic_list, + nonbasic_list); + active_subtrees_--; } @@ -944,33 +1052,48 @@ 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; } } } 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; + // 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_); + + const i_t m = leaf_problem.num_rows; + basis_update_mpf_t basis_update(m, settings_.refactor_frequency); + std::vector basic_list(m); + std::vector nonbasic_list; while (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_and_basis = true; search_tree_t subtree(std::move(start_node->node)); std::deque*> stack; stack.push_front(&subtree.root); + stats_t lp_stats; + lp_stats.total_lp_iters = 0; + lp_stats.total_lp_solve_time = 0; + lp_stats.nodes_explored = 0; + lp_stats.nodes_unexplored = 0; + while (stack.size() > 0 && status_ == mip_exploration_status_t::RUNNING) { mip_node_t* node_ptr = stack.front(); stack.pop_front(); @@ -978,39 +1101,43 @@ 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_bounds_and_basis = true; continue; } if (toc(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_status_t node_status = - solve_node(subtree, node_ptr, leaf_problem, Arow, upper_bound, log, 'D'); + node_status_t node_status = solve_node(node_ptr, + subtree, + leaf_problem, + basis_update, + basic_list, + nonbasic_list, + presolver, + thread_type_t::DIVING, + recompute_bounds_and_basis, + start_node->lower, + start_node->upper, + lp_stats, + log); + + recompute_bounds_and_basis = node_status != node_status_t::HAS_CHILDREN; if (node_status == node_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 (node_status == node_status_t::ITERATION_LIMIT) { + break; - if (stack.size() > 1) { - // 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_) { - mutex_dive_queue_.lock(); + } else if (node_status == node_status_t::HAS_CHILDREN) { + if (stack.size() > 0) { mip_node_t* new_node = stack.back(); stack.pop_back(); - dive_queue_.emplace(new_node->detach_copy(), leaf_problem.lower, leaf_problem.upper); - mutex_dive_queue_.unlock(); } + + auto [first, second] = child_selection(node_ptr); + stack.push_front(second); + stack.push_front(first); } } } @@ -1052,6 +1179,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. @@ -1153,19 +1281,17 @@ 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 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 " "| 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; @@ -1174,12 +1300,10 @@ 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; + diving_queue_.set_rng_seed(settings_.random_seed); #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(); @@ -1187,27 +1311,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 23fb9eb7f8..0fbe0405b8 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -17,10 +17,10 @@ #pragma once +#include #include #include #include -#include #include #include #include @@ -52,69 +52,20 @@ enum class mip_exploration_status_t { COMPLETED = 5, // The solver finished exploring the tree }; -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; - } +// 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 { @@ -203,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 diving_queue_; i_t min_diving_queue_size_; // Global status of the solver. @@ -221,45 +172,52 @@ 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& presolver, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list); // 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, + node_status_t solve_node(mip_node_t* node_ptr, + search_tree_t& search_tree, lp_problem_t& leaf_problem, - const csc_matrix_t& Arow, - f_t upper_bound, - logger_t& log, - char thread_type); + basis_update_mpf_t& ft, + std::vector& basic_list, + std::vector& nonbasic_list, + bounds_strengthening_t& presolver, + thread_type_t thread_type, + bool recompute, + const std::vector& root_lower, + const std::vector& root_upper, + stats_t& stats, + 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/diving_queue.hpp b/cpp/src/dual_simplex/diving_queue.hpp new file mode 100644 index 0000000000..f7c4e1aa13 --- /dev/null +++ b/cpp/src/dual_simplex/diving_queue.hpp @@ -0,0 +1,100 @@ +/* + * 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 + +#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.objective_estimate > b.node.objective_estimate; + } +}; + +// 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; + PCG rng; + const double epsilon = 0.1; // Probability to grab a random node + + public: + diving_queue_t() { buffer.reserve(max_size_); } + + void set_rng_seed(uint64_t seed) { rng.set_seed(PCG::default_seed + seed); } + + 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() + { + if (rng.next() <= epsilon) { + i_t idx = rng.uniform(0, buffer.size()); + std::swap(buffer[idx], buffer.back()); + diving_root_t node = std::move(buffer.back()); + buffer.pop_back(); + std::make_heap(buffer.begin(), buffer.end(), std::greater<>()); + return node; + + } else { + 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/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index a3a34eb81d..214eefac13 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -29,28 +29,32 @@ 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 process 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 + TIME_LIMIT = 6, // Time out during the LP relaxation + ITERATION_LIMIT = 7 // Reached the iteration limit during the LP relaxation }; +enum class rounding_direction_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), + objective_estimate(inf), vstatus(basis) { children[0] = nullptr; @@ -61,10 +65,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), @@ -72,41 +76,52 @@ class mip_node_t { branch_var(branch_variable), branch_dir(branch_direction), fractional_val(branch_var_value), + objective_estimate(parent_node->objective_estimate), 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 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; + } + mip_node_t* get_down_child() const { return children[0].get(); } mip_node_t* get_up_child() const { return children[1].get(); } @@ -203,27 +218,28 @@ 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); - copy.branch_var = branch_var; - copy.branch_dir = branch_dir; - copy.branch_var_lower = branch_var_lower; - copy.branch_var_upper = branch_var_upper; - copy.fractional_val = fractional_val; - copy.node_id = node_id; + copy.branch_var = branch_var; + copy.branch_dir = branch_dir; + copy.branch_var_lower = branch_var_lower; + copy.branch_var_upper = branch_var_upper; + copy.fractional_val = fractional_val; + copy.objective_estimate = objective_estimate; + copy.node_id = node_id; return copy; } node_status_t status; f_t lower_bound; + f_t objective_estimate; 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; @@ -270,13 +286,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, @@ -288,17 +303,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), @@ -319,7 +352,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) { @@ -327,7 +360,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/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 098c3b6e2b..995a6ed581 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -2245,6 +2245,8 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, phase2::bound_info(lp, settings); if (initialize_basis) { std::vector superbasic_list; + nonbasic_list.clear(); + nonbasic_list.reserve(n - m); get_basis_from_vstatus(m, vstatus, basic_list, nonbasic_list, superbasic_list); assert(superbasic_list.size() == 0); assert(nonbasic_list.size() == n - m); diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 8d80337c74..29e63d8676 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, @@ -842,15 +576,8 @@ 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 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); if (settings.barrier && settings.dualize != 0 && @@ -1608,13 +1335,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 diff --git a/cpp/src/dual_simplex/pseudo_costs.cpp b/cpp/src/dual_simplex/pseudo_costs.cpp index 4bd9590e16..ccab371acc 100644 --- a/cpp/src/dual_simplex/pseudo_costs.cpp +++ b/cpp/src/dual_simplex/pseudo_costs.cpp @@ -208,19 +208,18 @@ template void pseudo_costs_t::update_pseudo_costs(mip_node_t* node_ptr, f_t leaf_objective) { - mutex.lock(); + std::lock_guard lock(mutex); 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 { pseudo_cost_sum_up[node_ptr->branch_var] += change_in_obj / frac; pseudo_cost_num_up[node_ptr->branch_var]++; } - mutex.unlock(); } template @@ -261,65 +260,76 @@ i_t pseudo_costs_t::variable_selection(const std::vector& fractio const std::vector& solution, logger_t& log) { - mutex.lock(); + std::lock_guard lock(mutex); - const i_t num_fractional = fractional.size(); - std::vector pseudo_cost_up(num_fractional); - std::vector pseudo_cost_down(num_fractional); - std::vector score(num_fractional); + constexpr f_t eps = 1e-6; + i_t branch_var = fractional[0]; + f_t max_score = -1; i_t num_initialized_down; i_t num_initialized_up; - f_t pseudo_cost_down_avg; - f_t pseudo_cost_up_avg; + f_t pc_down_avg; + f_t pc_up_avg; + initialized(num_initialized_down, num_initialized_up, pc_down_avg, pc_up_avg); - initialized(num_initialized_down, num_initialized_up, pseudo_cost_down_avg, pseudo_cost_up_avg); + for (auto j : fractional) { + f_t f_down = solution[j] - std::floor(solution[j]); + f_t f_up = std::ceil(solution[j]) - solution[j]; - log.printf("PC: num initialized down %d up %d avg down %e up %e\n", - num_initialized_down, - num_initialized_up, - pseudo_cost_down_avg, - pseudo_cost_up_avg); + f_t pc_down = pseudo_cost_num_down[j] != 0 ? pseudo_cost_sum_down[j] / pseudo_cost_num_down[j] + : pc_down_avg; - for (i_t k = 0; k < num_fractional; k++) { - const i_t j = fractional[k]; - if (pseudo_cost_num_down[j] != 0) { - pseudo_cost_down[k] = pseudo_cost_sum_down[j] / pseudo_cost_num_down[j]; - } else { - pseudo_cost_down[k] = pseudo_cost_down_avg; - } + f_t pc_up = + pseudo_cost_num_up[j] != 0 ? pseudo_cost_sum_up[j] / pseudo_cost_num_up[j] : pc_up_avg; - if (pseudo_cost_num_up[j] != 0) { - pseudo_cost_up[k] = pseudo_cost_sum_up[j] / pseudo_cost_num_up[j]; - } else { - pseudo_cost_up[k] = pseudo_cost_up_avg; - } - constexpr f_t eps = 1e-6; - const f_t f_down = solution[j] - std::floor(solution[j]); - const f_t f_up = std::ceil(solution[j]) - solution[j]; - score[k] = - std::max(f_down * pseudo_cost_down[k], eps) * std::max(f_up * pseudo_cost_up[k], eps); - } + f_t score = std::max(f_down * pc_down, eps) * std::max(f_up * pc_up, eps); - i_t branch_var = fractional[0]; - f_t max_score = -1; - i_t select = -1; - for (i_t k = 0; k < num_fractional; k++) { - if (score[k] > max_score) { - max_score = score[k]; - branch_var = fractional[k]; - select = k; + if (score > max_score) { + max_score = score; + branch_var = j; } } - log.printf( - "pc branching on %d. Value %e. Score %e\n", branch_var, solution[branch_var], score[select]); - - mutex.unlock(); + log.debug("Pseudocost branching: selected %d with val = %e and score = %e\n", + branch_var, + solution[branch_var], + max_score); return branch_var; } +template +f_t pseudo_costs_t::objective_estimate(const std::vector& fractional, + const std::vector& solution, + f_t lower_bound, + logger_t& log) +{ + std::lock_guard lock(mutex); + + constexpr f_t eps = 1e-6; + f_t estimate = lower_bound; + + i_t num_initialized_down; + i_t num_initialized_up; + f_t pc_down_avg; + f_t pc_up_avg; + initialized(num_initialized_down, num_initialized_up, pc_down_avg, pc_up_avg); + + for (auto j : fractional) { + f_t f_down = solution[j] - std::floor(solution[j]); + f_t f_up = std::ceil(solution[j]) - solution[j]; + + f_t pc_down = pseudo_cost_num_down[j] != 0 ? pseudo_cost_sum_down[j] / pseudo_cost_num_down[j] + : pc_down_avg; + + f_t pc_up = + pseudo_cost_num_up[j] != 0 ? pseudo_cost_sum_up[j] / pseudo_cost_num_up[j] : pc_up_avg; + estimate += std::min(std::max(pc_down * f_down, eps), std::max(pc_up * f_up, eps)); + } + + return estimate; +} + template void pseudo_costs_t::update_pseudo_costs_from_strong_branching( const std::vector& fractional, const std::vector& root_soln) diff --git a/cpp/src/dual_simplex/pseudo_costs.hpp b/cpp/src/dual_simplex/pseudo_costs.hpp index 5bd03e3fcd..4831f3783d 100644 --- a/cpp/src/dual_simplex/pseudo_costs.hpp +++ b/cpp/src/dual_simplex/pseudo_costs.hpp @@ -57,6 +57,11 @@ class pseudo_costs_t { const std::vector& solution, logger_t& log); + f_t objective_estimate(const std::vector& fractional, + const std::vector& solution, + f_t lower_bound, + logger_t& log); + void update_pseudo_costs_from_strong_branching(const std::vector& fractional, const std::vector& root_soln); std::vector pseudo_cost_sum_up; diff --git a/cpp/src/dual_simplex/solution.hpp b/cpp/src/dual_simplex/solution.hpp index e4651dd112..8783fd8e8c 100644 --- a/cpp/src/dual_simplex/solution.hpp +++ b/cpp/src/dual_simplex/solution.hpp @@ -61,7 +61,11 @@ class lp_solution_t { template class mip_solution_t { public: - mip_solution_t(i_t n) : x(n), objective(std::numeric_limits::quiet_NaN()), lower_bound(-inf) + mip_solution_t(i_t n) + : x(n), + objective(std::numeric_limits::quiet_NaN()), + lower_bound(-inf), + has_incumbent(false) { } @@ -69,8 +73,9 @@ class mip_solution_t { void set_incumbent_solution(f_t primal_objective, const std::vector& primal_solution) { - x = primal_solution; - objective = primal_objective; + x = primal_solution; + objective = primal_objective; + has_incumbent = true; } // Primal solution vector @@ -79,6 +84,7 @@ class mip_solution_t { f_t lower_bound; i_t nodes_explored; i_t simplex_iterations; + bool has_incumbent; }; } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/utilities/pcg.hpp b/cpp/src/utilities/pcg.hpp new file mode 100644 index 0000000000..5d1737410d --- /dev/null +++ b/cpp/src/utilities/pcg.hpp @@ -0,0 +1,155 @@ +/* + * 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 + +// Copied from raft/PCGenerator (rng_device.cuh). +// It is based on the PCG code (https://www.pcg-random.org/). +namespace cuopt { +class PCG { + public: + static constexpr uint64_t default_seed = 0x853c49e6748fea9bULL; + static constexpr uint64_t default_stream = 0xda3e39cb94b95bdbULL; + + /** + * @brief ctor. Initializes the PCG + * @param rng_state is the generator state used for initializing the generator + * @param subsequence specifies the subsequence to be generated out of 2^64 possible subsequences + * In a parallel setting, like threads of a CUDA kernel, each thread is required to generate a + * unique set of random numbers. This can be achieved by initializing the generator with same + * rng_state for all the threads and diststreamt values for subsequence. + */ + PCG(const uint64_t seed = default_seed, + const uint64_t subsequence = default_stream, + uint64_t offset = 0) + { + set_seed(seed, subsequence, offset); + } + + // Set the seed, subsequence and offset of the PCG + void set_seed(uint64_t seed, const uint64_t subsequence = default_stream, uint64_t offset = 0) + { + state = uint64_t(0); + stream = (subsequence << 1u) | 1u; + uint32_t discard; + next(discard); + state += seed; + next(discard); + skipahead(offset); + } + + // Based on "Random Number Generation with Arbitrary Strides" F. B. Brown + // Link https://mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf + void skipahead(uint64_t offset) + { + uint64_t G = 1; + uint64_t h = 6364136223846793005ULL; + uint64_t C = 0; + uint64_t f = stream; + while (offset) { + if (offset & 1) { + G = G * h; + C = C * h + f; + } + f = f * (h + 1); + h = h * h; + offset >>= 1; + } + state = state * G + C; + } + + /** + * @defgroup NextRand Generate the next random number + * @brief This code is derived from PCG basic code + * @{ + */ + uint32_t next_u32() + { + uint32_t ret; + uint64_t oldstate = state; + state = oldstate * 6364136223846793005ULL + stream; + uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; + uint32_t rot = oldstate >> 59u; + ret = (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); + return ret; + } + + uint64_t next_u64() + { + uint64_t ret; + uint32_t a, b; + a = next_u32(); + b = next_u32(); + ret = uint64_t(a) | (uint64_t(b) << 32); + return ret; + } + + int32_t next_i32() + { + int32_t ret; + uint32_t val; + val = next_u32(); + ret = int32_t(val & 0x7fffffff); + return ret; + } + + int64_t next_i64() + { + int64_t ret; + uint64_t val; + val = next_u64(); + ret = int64_t(val & 0x7fffffffffffffff); + return ret; + } + + float next_float() { return static_cast((next_u32() >> 8) * 0x1.0p-24); } + + double next_double() { return static_cast((next_u64() >> 11) * 0x1.0p-53); } + + template + T next() + { + T val; + next(val); + return val; + } + + void next(uint32_t& ret) { ret = next_u32(); } + void next(uint64_t& ret) { ret = next_u64(); } + void next(int32_t& ret) { ret = next_i32(); } + void next(int64_t& ret) { ret = next_i64(); } + void next(float& ret) { ret = next_float(); } + void next(double& ret) { ret = next_double(); } + + /// Draws a sample from a uniform distribution. The samples are uniformly distributed over + /// the semi-closed interval `[low, high)`. This routine may have a **slight bias** toward + /// some numbers in the range (scaling by floating-point). + template + T uniform(T low, T high) + { + double val = next_double(); + T range = high - low; + return low + (val * range); + } + + private: + uint64_t state; + uint64_t stream; +}; +} // namespace cuopt