From 0fcae7eaae16671d2c88f5285ff726752962fcb4 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 29 Aug 2025 06:52:42 -0700 Subject: [PATCH 01/14] First draft node presolve --- cpp/src/dual_simplex/branch_and_bound.cpp | 15 +- cpp/src/dual_simplex/presolve.cpp | 181 ++++++++++++++-------- cpp/src/dual_simplex/presolve.hpp | 6 + 3 files changed, 132 insertions(+), 70 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index ce92532044..f684b6d910 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -642,8 +642,19 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut std::vector leaf_edge_norms = edge_norms; // = node.steepest_edge_norms; simplex_solver_settings_t lp_settings = settings; lp_settings.set_log(false); - lp_settings.cut_off = upper_bound + settings.dual_tol; - lp_settings.inside_mip = 2; + lp_settings.cut_off = upper_bound + settings.dual_tol; + lp_settings.inside_mip = 2; + + // in B&B we only have equality constraints + std::vector row_sense(leaf_problem.num_rows, 'E'); + std::vector is_int(leaf_problem.num_cols, 0); + // FIXME:: populate is_int correctly + bool feasible = bound_strengthening(row_sense, lp_settings, leaf_problem, is_int); + + if (!feasible) { + // do something here + } + dual::status_t lp_status = dual_phase2(2, 0, lp_start_time, diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 457da9113d..725eca37c4 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -23,26 +23,55 @@ 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 -void bound_strengthening(const std::vector& row_sense, +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); +} + +template +bool bound_strengthening(const std::vector& row_sense, const simplex_solver_settings_t& settings, - lp_problem_t& problem) + lp_problem_t& problem, + const std::vector& is_int) { const i_t m = problem.num_rows; const i_t n = problem.num_cols; - std::vector constraint_lower(m); - std::vector num_lower_infinity(m); - std::vector num_upper_infinity(m); + std::vector min_activity(m); + std::vector max_activity(m); csc_matrix_t Arow(1, 1, 1); problem.A.transpose(Arow); - std::vector less_rows; - less_rows.reserve(m); + std::vector constraint_lb(m); + std::vector constraint_ub(m); for (i_t i = 0; i < m; ++i) { - if (row_sense[i] == 'L') { less_rows.push_back(i); } + if (row_sense[i] == 'L') { + constraint_ub[i] = problem.rhs[i]; + constraint_lb[i] = -inf; + } else if (row_sense[i] == 'G') { + constraint_lb[i] = problem.rhs[i]; + constraint_ub[i] = inf; + } else { + constraint_lb[i] = problem.rhs[i]; + constraint_ub[i] = problem.rhs[i]; + } } std::vector lower = problem.lower; @@ -55,73 +84,88 @@ void bound_strengthening(const std::vector& row_sense, i_t iter = 0; const i_t iter_limit = 10; i_t total_strengthened_variables = 0; - settings.log.printf("Less equal rows %d\n", less_rows.size()); - while (iter < iter_limit && less_rows.size() > 0) { + while (iter < iter_limit) { // Derive bounds on the constraints - settings.log.printf("Running bound strengthening on %d rows\n", - static_cast(less_rows.size())); - for (i_t i : less_rows) { - const i_t row_start = Arow.col_start[i]; - const i_t row_end = Arow.col_start[i + 1]; - num_lower_infinity[i] = 0; - num_upper_infinity[i] = 0; - - f_t lower_limit = 0.0; + settings.log.printf("Running bound strengthening on %d rows\n", static_cast(n)); + + // Compute the min and max activity of the constraints + // FIXME:: use changed constraints logic + for (i_t i = 0; i < m; ++i) { + 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]; if (a_ij > 0) { - lower_limit += a_ij * lower[j]; + min_a += a_ij * lower[j]; + max_a += a_ij * upper[j]; } else if (a_ij < 0) { - lower_limit += a_ij * upper[j]; - } - if (lower[j] == -inf && a_ij > 0) { - num_lower_infinity[i]++; - lower_limit = -inf; - } - if (upper[j] == inf && a_ij < 0) { - num_lower_infinity[i]++; - lower_limit = -inf; + 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; } } - constraint_lower[i] = lower_limit; + min_activity[i] = min_a; + max_activity[i] = max_a; } - // Use the constraint bounds to derive new bounds on the variables - for (i_t i : less_rows) { - if (std::isfinite(constraint_lower[i]) && num_lower_infinity[i] == 0) { - const i_t row_start = Arow.col_start[i]; - const i_t row_end = Arow.col_start[i + 1]; - for (i_t p = row_start; p < row_end; ++p) { - const i_t k = Arow.i[p]; - const f_t a_ik = Arow.x[p]; - if (a_ik > 0) { - const f_t new_upper = lower[k] + (problem.rhs[i] - constraint_lower[i]) / a_ik; - if (new_upper < upper[k]) { - upper[k] = new_upper; - if (lower[k] > upper[k]) { - settings.log.printf( - "\t INFEASIBLE!!!!!!!!!!!!!!!!! constraint_lower %e lower %e rhs %e\n", - constraint_lower[i], - lower[k], - problem.rhs[i]); - } - if (!updated_variables_mark[k]) { updated_variables_list.push_back(k); } - } - } else if (a_ik < 0) { - const f_t new_lower = upper[k] + (problem.rhs[i] - constraint_lower[i]) / a_ik; - if (new_lower > lower[k]) { - lower[k] = new_lower; - if (lower[k] > upper[k]) { - settings.log.printf("\t INFEASIBLE !!!!!!!!!!!!!!!!!!1\n"); - } - if (!updated_variables_mark[k]) { updated_variables_list.push_back(k); } - } - } + i_t num_bounds_changed = 0; + // FIXME:: Use the transpose of Arow here to update variable bounds + // Use the activities to derive new bounds on the variables + for (i_t i = 0; i < m; ++i) { + const i_t row_start = Arow.col_start[i]; + const i_t row_end = Arow.col_start[i + 1]; + for (i_t p = row_start; p < row_end; ++p) { + const i_t k = Arow.i[p]; + const f_t a_ik = Arow.x[p]; + + f_t min_a = min_activity[i]; + f_t max_a = max_activity[i]; + 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) { return false; } + + min_a -= (a_ik < 0) ? a_ik * cnst_ub : a_ik * cnst_lb; + max_a -= (a_ik > 0) ? a_ik * cnst_ub : a_ik * cnst_lb; + + f_t delta_min_act = cnst_ub - min_a; + f_t delta_max_act = cnst_lb - max_a; + f_t old_lb = lower[k]; + f_t old_ub = upper[k]; + f_t new_lb = update_lb(old_lb, a_ik, delta_min_act, delta_max_act); + f_t new_ub = update_ub(old_ub, a_ik, delta_min_act, delta_max_act); + + // Integer rounding + if (!is_int.empty() && is_int[k]) { + 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; + + if (lb_updated) { lower[k] = new_lb; } + if (ub_updated) { upper[k] = new_ub; } + + bool bounds_changed = lb_updated || ub_updated; + if (bounds_changed) { + num_bounds_changed++; + if (!updated_variables_mark[k]) { updated_variables_list.push_back(k); } } } } - less_rows.clear(); + + if (num_bounds_changed == 0) { break; } // Update the bounds on the constraints settings.log.printf("Round %d: Strengthend %d variables\n", @@ -130,12 +174,6 @@ void bound_strengthening(const std::vector& row_sense, total_strengthened_variables += updated_variables_list.size(); for (i_t j : updated_variables_list) { updated_variables_mark[j] = 0; - const i_t col_start = problem.A.col_start[j]; - const i_t col_end = problem.A.col_start[j + 1]; - for (i_t p = col_start; p < col_end; ++p) { - const i_t i = problem.A.i[p]; - less_rows.push_back(i); - } } updated_variables_list.clear(); iter++; @@ -143,6 +181,7 @@ void bound_strengthening(const std::vector& row_sense, settings.log.printf("Total strengthened variables %d\n", total_strengthened_variables); problem.lower = lower; problem.upper = upper; + return true; } template @@ -1077,6 +1116,12 @@ template void uncrush_solution(const presolve_info_t& const std::vector& crushed_z, std::vector& uncrushed_x, std::vector& uncrushed_z); + +template bool bound_strengthening( + const std::vector& row_sense, + const simplex_solver_settings_t& settings, + lp_problem_t& problem, + const std::vector& is_int); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 7a307e6f73..4d7ba2c901 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -131,4 +131,10 @@ void uncrush_solution(const presolve_info_t& presolve_info, std::vector& uncrushed_x, std::vector& uncrushed_z); +template +bool bound_strengthening(const std::vector& row_sense, + const simplex_solver_settings_t& settings, + lp_problem_t& problem, + const std::vector& is_int = {}); + } // namespace cuopt::linear_programming::dual_simplex From b5fd37e38ded83f5370d114d35f3ba2af07fdf00 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 29 Aug 2025 06:59:54 -0700 Subject: [PATCH 02/14] Don't run LP solve if infeasibility is detected in bounds strengthening --- cpp/src/dual_simplex/branch_and_bound.cpp | 35 ++++++++++++----------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index f684b6d910..dd0a78a1f3 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -651,24 +651,25 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // FIXME:: populate is_int correctly bool feasible = bound_strengthening(row_sense, lp_settings, leaf_problem, is_int); - if (!feasible) { - // do something here - } + dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; - dual::status_t lp_status = dual_phase2(2, - 0, - lp_start_time, - leaf_problem, - lp_settings, - leaf_vstatus, - leaf_solution, - node_iter, - leaf_edge_norms); - if (lp_status == dual::status_t::NUMERICAL) { - settings.log.printf("Numerical issue node %d. Resolving from scratch.\n", nodes_explored); - lp_status_t second_status = solve_linear_program_advanced( - leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); - lp_status = convert_lp_status_to_dual_status(second_status); + // If the problem is infeasible after bounds strengthening, we don't need to solve the LP + if (feasible) { + lp_status = dual_phase2(2, + 0, + lp_start_time, + leaf_problem, + lp_settings, + leaf_vstatus, + leaf_solution, + node_iter, + leaf_edge_norms); + if (lp_status == dual::status_t::NUMERICAL) { + settings.log.printf("Numerical issue node %d. Resolving from scratch.\n", nodes_explored); + lp_status_t second_status = solve_linear_program_advanced( + leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms); + lp_status = convert_lp_status_to_dual_status(second_status); + } } total_lp_solve_time += toc(lp_start_time); total_lp_iters += node_iter; From 9062867984ff61a1cb3e1a510bd264ea532f18e0 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 29 Aug 2025 07:02:24 -0700 Subject: [PATCH 03/14] Add a log message --- cpp/src/dual_simplex/branch_and_bound.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index dd0a78a1f3..3abcd8724a 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -676,6 +676,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut nodes_explored++; if (lp_status == dual::status_t::DUAL_UNBOUNDED) { + if (!feasible) { + settings.log.printf("Infeasible after bounds strengthening. Fathoming node %d.\n", + nodes_explored); + } node_ptr->lower_bound = inf; std::vector*> stack; node_ptr->set_status(node_status_t::INFEASIBLE, stack); From 4779b5d59f12a0d1b40d4af2471ad9ce773d70bd Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 29 Aug 2025 07:16:25 -0700 Subject: [PATCH 04/14] use variable types instead --- cpp/src/dual_simplex/branch_and_bound.cpp | 5 ++--- cpp/src/dual_simplex/presolve.cpp | 7 ++++--- cpp/src/dual_simplex/presolve.hpp | 3 ++- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 3abcd8724a..e0e1769b83 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -631,6 +631,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // Set the correct bounds for the leaf problem leaf_problem.lower = original_lp.lower; leaf_problem.upper = original_lp.upper; + // FIXME:: should we get the bounds from the node/parent instead of the original problem? node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper); std::vector& leaf_vstatus = node_ptr->vstatus; @@ -647,9 +648,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // in B&B we only have equality constraints std::vector row_sense(leaf_problem.num_rows, 'E'); - std::vector is_int(leaf_problem.num_cols, 0); - // FIXME:: populate is_int correctly - bool feasible = bound_strengthening(row_sense, lp_settings, leaf_problem, is_int); + bool feasible = bound_strengthening(row_sense, lp_settings, leaf_problem, var_types); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 725eca37c4..647065e400 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -47,7 +47,7 @@ template bool bound_strengthening(const std::vector& row_sense, const simplex_solver_settings_t& settings, lp_problem_t& problem, - const std::vector& is_int) + const std::vector& var_types) { const i_t m = problem.num_rows; const i_t n = problem.num_cols; @@ -146,7 +146,8 @@ bool bound_strengthening(const std::vector& row_sense, f_t new_ub = update_ub(old_ub, a_ik, delta_min_act, delta_max_act); // Integer rounding - if (!is_int.empty() && is_int[k]) { + 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); } @@ -1121,7 +1122,7 @@ template bool bound_strengthening( const std::vector& row_sense, const simplex_solver_settings_t& settings, lp_problem_t& problem, - const std::vector& is_int); + const std::vector& var_types); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 4d7ba2c901..c1ea5f9c74 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -131,10 +131,11 @@ void uncrush_solution(const presolve_info_t& presolve_info, std::vector& uncrushed_x, 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 std::vector& is_int = {}); + const std::vector& var_types = {}); } // namespace cuopt::linear_programming::dual_simplex From 1fcc491ca71ccd874a7d7da8277c42f302c0f95b Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 29 Aug 2025 11:08:00 -0700 Subject: [PATCH 05/14] Fix a bug --- cpp/src/dual_simplex/presolve.cpp | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 647065e400..a42549c31e 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -133,15 +133,21 @@ bool bound_strengthening(const std::vector& row_sense, bool is_infeasible = check_infeasibility(min_a, max_a, cnst_lb, cnst_ub, settings.primal_tol); - if (is_infeasible) { return false; } + if (is_infeasible) { + settings.log.printf("Infeasible constraint %d, cnst_lb %e, cnst_ub %e, min_a %e, max_a %e\n", i, cnst_lb, cnst_ub, min_a, max_a); + std::cout << "Iter:: " << iter << ", Infeasible constraint " << i << ", cnst_lb " << cnst_lb << ", cnst_ub " << cnst_ub << ", min_a " << min_a << ", max_a " << max_a << std::endl; + return false; + } + + f_t old_lb = lower[k]; + f_t old_ub = upper[k]; - min_a -= (a_ik < 0) ? a_ik * cnst_ub : a_ik * cnst_lb; - max_a -= (a_ik > 0) ? a_ik * cnst_ub : a_ik * cnst_lb; + min_a -= (a_ik < 0) ? a_ik * old_ub : a_ik * old_lb; + max_a -= (a_ik > 0) ? a_ik * old_ub : a_ik * old_lb; f_t delta_min_act = cnst_ub - min_a; f_t delta_max_act = cnst_lb - max_a; - f_t old_lb = lower[k]; - f_t old_ub = upper[k]; + f_t new_lb = update_lb(old_lb, a_ik, delta_min_act, delta_max_act); f_t new_ub = update_ub(old_ub, a_ik, delta_min_act, delta_max_act); From 01558854db4170145ce43d52dd51325a7fd8c2a5 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 29 Aug 2025 11:45:36 -0700 Subject: [PATCH 06/14] Fix compilation error --- cpp/src/dual_simplex/presolve.cpp | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index a42549c31e..c997be5d62 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -20,6 +20,7 @@ #include #include +#include namespace cuopt::linear_programming::dual_simplex { @@ -133,14 +134,22 @@ bool bound_strengthening(const std::vector& row_sense, bool is_infeasible = check_infeasibility(min_a, max_a, cnst_lb, cnst_ub, settings.primal_tol); - if (is_infeasible) { - settings.log.printf("Infeasible constraint %d, cnst_lb %e, cnst_ub %e, min_a %e, max_a %e\n", i, cnst_lb, cnst_ub, min_a, max_a); - std::cout << "Iter:: " << iter << ", Infeasible constraint " << i << ", cnst_lb " << cnst_lb << ", cnst_ub " << cnst_ub << ", min_a " << min_a << ", max_a " << max_a << std::endl; - return false; + if (is_infeasible) { + settings.log.printf( + "Infeasible constraint %d, cnst_lb %e, cnst_ub %e, min_a %e, max_a %e\n", + i, + cnst_lb, + cnst_ub, + min_a, + max_a); + std::cout << "Iter:: " << iter << ", Infeasible constraint " << i << ", cnst_lb " + << cnst_lb << ", cnst_ub " << cnst_ub << ", min_a " << min_a << ", max_a " + << max_a << std::endl; + return false; } - f_t old_lb = lower[k]; - f_t old_ub = upper[k]; + f_t old_lb = lower[k]; + f_t old_ub = upper[k]; min_a -= (a_ik < 0) ? a_ik * old_ub : a_ik * old_lb; max_a -= (a_ik > 0) ? a_ik * old_ub : a_ik * old_lb; @@ -148,8 +157,8 @@ bool bound_strengthening(const std::vector& row_sense, f_t delta_min_act = cnst_ub - min_a; f_t delta_max_act = cnst_lb - max_a; - f_t new_lb = update_lb(old_lb, a_ik, delta_min_act, delta_max_act); - f_t new_ub = update_ub(old_ub, a_ik, delta_min_act, delta_max_act); + f_t new_lb = update_lb(old_lb, a_ik, delta_min_act, delta_max_act); + f_t new_ub = update_ub(old_ub, a_ik, delta_min_act, delta_max_act); // Integer rounding if (!var_types.empty() && From e7f043e3e5eff8b5c2f94cde38f96483657b9f43 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Tue, 2 Sep 2025 14:13:10 -0700 Subject: [PATCH 07/14] Debug prints --- cpp/src/dual_simplex/presolve.cpp | 109 +++++++++++++++++++++++++++--- 1 file changed, 100 insertions(+), 9 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index c997be5d62..65444aeb17 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -44,6 +44,27 @@ static inline bool check_infeasibility(f_t min_a, f_t max_a, f_t cnst_lb, f_t cn 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, @@ -62,6 +83,9 @@ bool bound_strengthening(const std::vector& row_sense, std::vector constraint_lb(m); std::vector constraint_ub(m); + i_t num_finite_lower_bounds = 0; + i_t num_finite_upper_bounds = 0; + for (i_t i = 0; i < m; ++i) { if (row_sense[i] == 'L') { constraint_ub[i] = problem.rhs[i]; @@ -77,6 +101,14 @@ bool bound_strengthening(const std::vector& row_sense, std::vector lower = problem.lower; std::vector upper = problem.upper; + print_bounds_stats(lower, upper, settings, "Initial bounds"); + for (i_t i = 0; i < n; ++i) { + if (lower[i] > upper[i]) { + // std::cout << "Infeasible variable " << i << ", lower " << lower[i] << ", upper " << + // upper[i] << std::endl; + return false; + } + } std::vector updated_variables_list; updated_variables_list.reserve(n); @@ -86,8 +118,9 @@ bool bound_strengthening(const std::vector& row_sense, const i_t iter_limit = 10; i_t total_strengthened_variables = 0; while (iter < iter_limit) { + // std::cout << " XXXXXXXXXX Iter:: " << iter << std::endl; // Derive bounds on the constraints - settings.log.printf("Running bound strengthening on %d rows\n", static_cast(n)); + // settings.log.printf("Running bound strengthening on %d rows\n", static_cast(n)); // Compute the min and max activity of the constraints // FIXME:: use changed constraints logic @@ -113,6 +146,13 @@ bool bound_strengthening(const std::vector& row_sense, if (lower[j] == -inf && a_ij > 0) { min_a = -inf; } if (upper[j] == inf && a_ij < 0) { min_a = -inf; } } + + if (max_a < min_a) { + // std::cout << "Iter:: " << iter << ", Infeasible constraint " << i << ", min_a " << min_a + // << ", max_a " << max_a << std::endl; + return false; + } + min_activity[i] = min_a; max_activity[i] = max_a; } @@ -142,9 +182,6 @@ bool bound_strengthening(const std::vector& row_sense, cnst_ub, min_a, max_a); - std::cout << "Iter:: " << iter << ", Infeasible constraint " << i << ", cnst_lb " - << cnst_lb << ", cnst_ub " << cnst_ub << ", min_a " << min_a << ", max_a " - << max_a << std::endl; return false; } @@ -167,12 +204,22 @@ bool bound_strengthening(const std::vector& row_sense, new_ub = std::floor(new_ub + settings.integer_tol); } + if (std::isfinite(old_lb)) { assert(std::isfinite(new_lb)); } + if (std::isfinite(old_ub)) { assert(std::isfinite(new_ub)); } + bool lb_updated = abs(new_lb - old_lb) > 1e3 * settings.primal_tol; bool ub_updated = abs(new_ub - old_ub) > 1e3 * settings.primal_tol; if (lb_updated) { lower[k] = new_lb; } if (ub_updated) { upper[k] = new_ub; } + if (lower[k] > upper[k] + 1e-6) { + // std::cout << "Iter:: " << iter << ", Infeasible variable after update " << k << ", + // lower " << lower[k] << ", upper " << upper[k] << std::endl; + return false; + } + lower[k] = std::min(lower[k], upper[k]); + bool bounds_changed = lb_updated || ub_updated; if (bounds_changed) { num_bounds_changed++; @@ -181,12 +228,17 @@ bool bound_strengthening(const std::vector& row_sense, } } - if (num_bounds_changed == 0) { break; } + if (num_bounds_changed == 0) { + settings.log.printf("No bounds changed, iter %d\n", iter + 1); + break; + } else { + settings.log.printf("Num bounds changed %d, iter %d\n", num_bounds_changed, iter + 1); + } // Update the bounds on the constraints - settings.log.printf("Round %d: Strengthend %d variables\n", - iter, - static_cast(updated_variables_list.size())); + // settings.log.printf("Round %d: Strengthend %d variables\n", + // iter, + // static_cast(updated_variables_list.size())); total_strengthened_variables += updated_variables_list.size(); for (i_t j : updated_variables_list) { updated_variables_mark[j] = 0; @@ -194,9 +246,48 @@ bool bound_strengthening(const std::vector& row_sense, updated_variables_list.clear(); iter++; } - settings.log.printf("Total strengthened variables %d\n", total_strengthened_variables); + // 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; } From e5eb0e7c18f231dfcb007ea84e2a4a295c65cac1 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 3 Sep 2025 12:43:43 -0700 Subject: [PATCH 08/14] Implement changed constraints and variables logic --- cpp/src/dual_simplex/presolve.cpp | 180 ++++++++++++++---------------- 1 file changed, 84 insertions(+), 96 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 65444aeb17..5dfd90efae 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -18,6 +18,7 @@ #include #include +#include #include #include @@ -74,17 +75,18 @@ bool bound_strengthening(const std::vector& row_sense, const i_t m = problem.num_rows; const i_t n = problem.num_cols; - std::vector min_activity(m); - std::vector max_activity(m); - + // FIXME:: pre-compute transpose and store csc_matrix_t Arow(1, 1, 1); problem.A.transpose(Arow); + std::vector min_activity(m); + std::vector max_activity(m); std::vector constraint_lb(m); std::vector constraint_ub(m); - i_t num_finite_lower_bounds = 0; - i_t num_finite_upper_bounds = 0; + std::vector constraint_changed(m, true); + std::vector variable_changed(n, false); + std::vector constraint_changed_next(m, false); for (i_t i = 0; i < m; ++i) { if (row_sense[i] == 'L') { @@ -102,29 +104,12 @@ bool bound_strengthening(const std::vector& row_sense, std::vector lower = problem.lower; std::vector upper = problem.upper; print_bounds_stats(lower, upper, settings, "Initial bounds"); - for (i_t i = 0; i < n; ++i) { - if (lower[i] > upper[i]) { - // std::cout << "Infeasible variable " << i << ", lower " << lower[i] << ", upper " << - // upper[i] << std::endl; - return false; - } - } - std::vector updated_variables_list; - updated_variables_list.reserve(n); - std::vector updated_variables_mark(n, 0); - - i_t iter = 0; - const i_t iter_limit = 10; - i_t total_strengthened_variables = 0; + i_t iter = 0; + const i_t iter_limit = 10; while (iter < iter_limit) { - // std::cout << " XXXXXXXXXX Iter:: " << iter << std::endl; - // Derive bounds on the constraints - // settings.log.printf("Running bound strengthening on %d rows\n", static_cast(n)); - - // Compute the min and max activity of the constraints - // FIXME:: use changed constraints logic 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]; @@ -133,6 +118,8 @@ bool bound_strengthening(const std::vector& row_sense, 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]; @@ -148,8 +135,24 @@ bool bound_strengthening(const std::vector& row_sense, } if (max_a < min_a) { - // std::cout << "Iter:: " << iter << ", Infeasible constraint " << i << ", min_a " << min_a - // << ", max_a " << max_a << std::endl; + settings.log.printf( + "Iter:: %d, Infeasible constraint %d, min_a %e, max_a %e\n", iter, i, min_a, max_a); + return false; + } + + 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; } @@ -158,94 +161,79 @@ bool bound_strengthening(const std::vector& row_sense, } i_t num_bounds_changed = 0; - // FIXME:: Use the transpose of Arow here to update variable bounds - // Use the activities to derive new bounds on the variables - for (i_t i = 0; i < m; ++i) { - const i_t row_start = Arow.col_start[i]; - const i_t row_end = Arow.col_start[i + 1]; - for (i_t p = row_start; p < row_end; ++p) { - const i_t k = Arow.i[p]; - const f_t a_ik = Arow.x[p]; - - f_t min_a = min_activity[i]; - f_t max_a = max_activity[i]; - 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( - "Infeasible constraint %d, cnst_lb %e, cnst_ub %e, min_a %e, max_a %e\n", - i, - cnst_lb, - cnst_ub, - min_a, - max_a); - return false; - } - f_t old_lb = lower[k]; - f_t old_ub = upper[k]; + 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 min_a = min_activity[i]; + f_t max_a = max_activity[i]; + f_t cnst_lb = constraint_lb[i]; + f_t cnst_ub = constraint_ub[i]; min_a -= (a_ik < 0) ? a_ik * old_ub : a_ik * old_lb; max_a -= (a_ik > 0) ? a_ik * old_ub : a_ik * old_lb; f_t delta_min_act = cnst_ub - min_a; f_t delta_max_act = cnst_lb - max_a; - f_t new_lb = update_lb(old_lb, a_ik, delta_min_act, delta_max_act); - f_t 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); - } - - if (std::isfinite(old_lb)) { assert(std::isfinite(new_lb)); } - if (std::isfinite(old_ub)) { assert(std::isfinite(new_ub)); } + 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)); + } - bool lb_updated = abs(new_lb - old_lb) > 1e3 * settings.primal_tol; - bool ub_updated = abs(new_ub - old_ub) > 1e3 * settings.primal_tol; + // 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); + } - if (lb_updated) { lower[k] = new_lb; } - if (ub_updated) { upper[k] = new_ub; } + bool lb_updated = abs(new_lb - old_lb) > 1e3 * settings.primal_tol; + bool ub_updated = abs(new_ub - old_ub) > 1e3 * settings.primal_tol; - if (lower[k] > upper[k] + 1e-6) { - // std::cout << "Iter:: " << iter << ", Infeasible variable after update " << k << ", - // lower " << lower[k] << ", upper " << upper[k] << std::endl; - return false; - } - lower[k] = std::min(lower[k], upper[k]); + new_lb = std::max(new_lb, problem.lower[k]); + new_ub = std::min(new_ub, problem.upper[k]); - bool bounds_changed = lb_updated || ub_updated; - if (bounds_changed) { - num_bounds_changed++; - if (!updated_variables_mark[k]) { updated_variables_list.push_back(k); } + if (new_lb > new_ub + 1e-6) { + // std::cout << "Iter:: " << iter << ", Infeasible variable after update " << k << ", " << + // new_lb << " > " << new_ub << std::endl; + 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; } } - } - if (num_bounds_changed == 0) { - settings.log.printf("No bounds changed, iter %d\n", iter + 1); - break; - } else { - settings.log.printf("Num bounds changed %d, iter %d\n", num_bounds_changed, iter + 1); - } + lower[k] = std::min(new_lb, new_ub); + upper[k] = std::max(new_lb, new_ub); - // Update the bounds on the constraints - // settings.log.printf("Round %d: Strengthend %d variables\n", - // iter, - // static_cast(updated_variables_list.size())); - total_strengthened_variables += updated_variables_list.size(); - for (i_t j : updated_variables_list) { - updated_variables_mark[j] = 0; + bool bounds_changed = lb_updated || ub_updated; + if (bounds_changed) { num_bounds_changed++; } } - updated_variables_list.clear(); + + 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 From 8eb8d9727b1dcf090d97fd89cdcb3aad11c3e055 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 3 Sep 2025 12:44:53 -0700 Subject: [PATCH 09/14] style checks --- cpp/src/dual_simplex/presolve.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 5dfd90efae..e20cd20f31 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -84,6 +84,10 @@ bool bound_strengthening(const std::vector& row_sense, 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); From 08e51e3eafe5afef35bb6eaa037701517dc88a2b Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 3 Sep 2025 13:02:34 -0700 Subject: [PATCH 10/14] Add status to string --- cpp/src/dual_simplex/phase2.hpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp index ef52c0a190..39311e0607 100644 --- a/cpp/src/dual_simplex/phase2.hpp +++ b/cpp/src/dual_simplex/phase2.hpp @@ -38,7 +38,22 @@ enum class status_t { CONCURRENT_LIMIT = 6, UNSET = 7 }; + +static std::string status_to_string(status_t status) +{ + switch (status) { + case status_t::OPTIMAL: return "OPTIMAL"; + case status_t::DUAL_UNBOUNDED: return "DUAL_UNBOUNDED"; + case status_t::NUMERICAL: return "NUMERICAL"; + case status_t::CUTOFF: return "CUTOFF"; + case status_t::TIME_LIMIT: return "TIME_LIMIT"; + case status_t::ITERATION_LIMIT: return "ITERATION_LIMIT"; + case status_t::CONCURRENT_LIMIT: return "CONCURRENT_LIMIT"; + case status_t::UNSET: return "UNSET"; + } + return "UNKNOWN"; } +} // namespace dual template dual::status_t dual_phase2(i_t phase, From 11673e0e3f4d6559889482619bec5128a88f62b4 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Wed, 3 Sep 2025 13:08:31 -0700 Subject: [PATCH 11/14] cleanup --- cpp/src/dual_simplex/presolve.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index e20cd20f31..6d1998e598 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -138,12 +138,6 @@ bool bound_strengthening(const std::vector& row_sense, if (upper[j] == inf && a_ij < 0) { min_a = -inf; } } - if (max_a < min_a) { - settings.log.printf( - "Iter:: %d, Infeasible constraint %d, min_a %e, max_a %e\n", iter, i, min_a, max_a); - return false; - } - f_t cnst_lb = constraint_lb[i]; f_t cnst_ub = constraint_ub[i]; bool is_infeasible = From 8fdc8cbf4dddcde54f7bd7f1cfa04f4ac903a854 Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Fri, 5 Sep 2025 08:29:42 -0700 Subject: [PATCH 12/14] Allow specifying models in a file --- .../linear_programming/run_mps_files.sh | 53 +++++++++++++++++-- 1 file changed, 50 insertions(+), 3 deletions(-) diff --git a/benchmarks/linear_programming/run_mps_files.sh b/benchmarks/linear_programming/run_mps_files.sh index 61a8af0e04..5baf729b90 100755 --- a/benchmarks/linear_programming/run_mps_files.sh +++ b/benchmarks/linear_programming/run_mps_files.sh @@ -79,6 +79,7 @@ Optional Arguments: --batch-num Batch number --n-batches Number of batches --log-to-console Log to console + --model-list File containing a list of models to run -h, --help Show this help message and exit Examples: @@ -168,6 +169,11 @@ while [[ $# -gt 0 ]]; do LOG_TO_CONSOLE="$2" shift 2 ;; + --model-list) + echo "MODEL_LIST: $2" + MODEL_LIST="$2" + shift 2 + ;; *) echo "Unknown argument: $1" print_help @@ -194,7 +200,7 @@ PRESOLVE=${PRESOLVE:-true} BATCH_NUM=${BATCH_NUM:-0} N_BATCHES=${N_BATCHES:-1} LOG_TO_CONSOLE=${LOG_TO_CONSOLE:-true} - +MODEL_LIST=${MODEL_LIST:-} # Determine GPU list if [[ -n "$CUDA_VISIBLE_DEVICES" ]]; then IFS=',' read -ra GPU_LIST <<< "$CUDA_VISIBLE_DEVICES" @@ -206,8 +212,49 @@ else fi GPU_COUNT=${#GPU_LIST[@]} -# Gather all mps files into an array -mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps) +# Ensure all entries in MODEL_LIST have .mps extension +if [[ -n "$MODEL_LIST" && -f "$MODEL_LIST" ]]; then + # Create a temporary file to store the updated model list + TMP_MODEL_LIST=$(mktemp) + while IFS= read -r line || [[ -n "$line" ]]; do + # Skip empty lines + [[ -z "$line" ]] && continue + # If the line does not end with .mps, append it + if [[ "$line" != *.mps ]]; then + echo "${line}.mps" >> "$TMP_MODEL_LIST" + else + echo "$line" >> "$TMP_MODEL_LIST" + fi + done < "$MODEL_LIST" + # Overwrite the original MODEL_LIST with the updated one + mv "$TMP_MODEL_LIST" "$MODEL_LIST" +fi + + +# Gather all mps files into an array, either from the model list or from the directory +if [[ -n "$MODEL_LIST" ]]; then + if [[ ! -f "$MODEL_LIST" ]]; then + echo "Model list file not found: $MODEL_LIST" + exit 1 + fi + mapfile -t mps_files < <(grep -v '^\s*$' "$MODEL_LIST" | sed "s|^|$MPS_DIR/|") + # Optionally, check that all files exist + missing_files=() + for f in "${mps_files[@]}"; do + if [[ ! -f "$f" ]]; then + missing_files+=("$f") + fi + done + if (( ${#missing_files[@]} > 0 )); then + echo "The following files from the model list do not exist in $MPS_DIR:" + for f in "${missing_files[@]}"; do + echo " $f" + done + exit 1 + fi +else + mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps) +fi # Calculate batch size and start/end indices batch_size=$(( (${#mps_files[@]} + N_BATCHES - 1) / N_BATCHES )) From d6d4030a95d21708e0e6e571e4e3378abbb89ebb Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Mon, 8 Sep 2025 11:11:26 -0700 Subject: [PATCH 13/14] address PR comments --- cpp/src/dual_simplex/branch_and_bound.cpp | 4 +- cpp/src/dual_simplex/presolve.cpp | 54 ++++++++++++----------- 2 files changed, 30 insertions(+), 28 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index e0e1769b83..7730bbc0e3 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -646,8 +646,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut lp_settings.cut_off = upper_bound + settings.dual_tol; lp_settings.inside_mip = 2; - // in B&B we only have equality constraints - std::vector row_sense(leaf_problem.num_rows, 'E'); + // 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, var_types); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index 6d1998e598..b2c24b889e 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -79,29 +79,36 @@ bool bound_strengthening(const std::vector& row_sense, csc_matrix_t Arow(1, 1, 1); problem.A.transpose(Arow); - std::vector min_activity(m); - std::vector max_activity(m); + 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 + // 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 + // 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); - for (i_t i = 0; i < m; ++i) { - if (row_sense[i] == 'L') { - constraint_ub[i] = problem.rhs[i]; - constraint_lb[i] = -inf; - } else if (row_sense[i] == 'G') { - constraint_lb[i] = problem.rhs[i]; - constraint_ub[i] = inf; - } else { - constraint_lb[i] = problem.rhs[i]; - constraint_ub[i] = problem.rhs[i]; + 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; + } } } @@ -154,8 +161,8 @@ bool bound_strengthening(const std::vector& row_sense, return false; } - min_activity[i] = min_a; - max_activity[i] = max_a; + delta_min_activity[i] = cnst_ub - min_a; + delta_max_activity[i] = cnst_lb - max_a; } i_t num_bounds_changed = 0; @@ -175,15 +182,12 @@ bool bound_strengthening(const std::vector& row_sense, if (!constraint_changed[i]) { continue; } const f_t a_ik = problem.A.x[p]; - f_t min_a = min_activity[i]; - f_t max_a = max_activity[i]; - f_t cnst_lb = constraint_lb[i]; - f_t cnst_ub = constraint_ub[i]; - min_a -= (a_ik < 0) ? a_ik * old_ub : a_ik * old_lb; - max_a -= (a_ik > 0) ? a_ik * old_ub : a_ik * old_lb; - f_t delta_min_act = cnst_ub - min_a; - f_t delta_max_act = cnst_lb - max_a; + 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)); @@ -203,8 +207,6 @@ bool bound_strengthening(const std::vector& row_sense, new_ub = std::min(new_ub, problem.upper[k]); if (new_lb > new_ub + 1e-6) { - // std::cout << "Iter:: " << iter << ", Infeasible variable after update " << k << ", " << - // new_lb << " > " << new_ub << std::endl; settings.log.printf( "Iter:: %d, Infeasible variable after update %d, %e > %e\n", iter, k, new_lb, new_ub); return false; From 3a007551320a2419c53cc4a977dfe67e9b1812fd Mon Sep 17 00:00:00 2001 From: Rajesh Gandham Date: Mon, 8 Sep 2025 13:45:14 -0700 Subject: [PATCH 14/14] Optimize with changed variables --- cpp/src/dual_simplex/branch_and_bound.cpp | 12 ++++++--- cpp/src/dual_simplex/mip_node.hpp | 13 +++++++--- cpp/src/dual_simplex/presolve.cpp | 30 +++++++++++++++++------ cpp/src/dual_simplex/presolve.hpp | 4 ++- 4 files changed, 44 insertions(+), 15 deletions(-) diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 7730bbc0e3..898f0a85ef 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -529,6 +529,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut std::move(up_child)); // child pointers moved into the tree lp_problem_t leaf_problem = original_lp; // Make a copy of the original LP. We will modify its bounds at each leaf + csc_matrix_t Arow(1, 1, 1); + original_lp.A.transpose(Arow); f_t gap = get_upper_bound() - lower_bound; i_t nodes_explored = 0; settings.log.printf( @@ -631,8 +633,11 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // Set the correct bounds for the leaf problem leaf_problem.lower = original_lp.lower; leaf_problem.upper = original_lp.upper; - // FIXME:: should we get the bounds from the node/parent instead of the original problem? - node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper); + std::vector bounds_changed(original_lp.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); std::vector& leaf_vstatus = node_ptr->vstatus; lp_solution_t leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols); @@ -648,7 +653,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut // 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, var_types); + bool feasible = + bound_strengthening(row_sense, lp_settings, leaf_problem, Arow, var_types, bounds_changed); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index cae93a4dc7..c4d83a8333 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -79,22 +79,27 @@ class mip_node_t { children[1] = nullptr; } - void get_variable_bounds(std::vector& lower, std::vector& upper) const + 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; 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; - parent_ptr = parent_ptr->parent; + 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; } } diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index b2c24b889e..48d696bf96 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -70,15 +70,13 @@ template bool bound_strengthening(const std::vector& row_sense, const simplex_solver_settings_t& settings, lp_problem_t& problem, - const std::vector& var_types) + 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; - // FIXME:: pre-compute transpose and store - csc_matrix_t Arow(1, 1, 1); - problem.A.transpose(Arow); - std::vector delta_min_activity(m); std::vector delta_max_activity(m); std::vector constraint_lb(m); @@ -92,6 +90,20 @@ bool bound_strengthening(const std::vector& row_sense, 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()); @@ -825,7 +837,9 @@ void convert_user_problem(const user_problem_t& user_problem, constexpr bool run_bound_strengthening = false; if constexpr (run_bound_strengthening) { settings.log.printf("Running bound strengthening\n"); - bound_strengthening(row_sense, settings, problem); + csc_matrix_t Arow(1, 1, 1); + problem.A.transpose(Arow); + bound_strengthening(row_sense, settings, problem, Arow); } // The original problem may have a variable without a lower bound @@ -1216,7 +1230,9 @@ template bool bound_strengthening( const std::vector& row_sense, const simplex_solver_settings_t& settings, lp_problem_t& problem, - const std::vector& var_types); + 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 c1ea5f9c74..83b8c7d751 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -136,6 +136,8 @@ template bool bound_strengthening(const std::vector& row_sense, const simplex_solver_settings_t& settings, lp_problem_t& problem, - const std::vector& var_types = {}); + const csc_matrix_t& Arow, + const std::vector& var_types = {}, + const std::vector& bounds_changed = {}); } // namespace cuopt::linear_programming::dual_simplex