diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index cc051ab49..15082254e 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -57,6 +57,7 @@ #define CUOPT_MIP_HEURISTICS_ONLY "mip_heuristics_only" #define CUOPT_MIP_SCALING "mip_scaling" #define CUOPT_MIP_PRESOLVE "mip_presolve" +#define CUOPT_MIP_CUT_PASSES "mip_cut_passes" #define CUOPT_SOLUTION_FILE "solution_file" #define CUOPT_NUM_CPU_THREADS "num_cpu_threads" #define CUOPT_USER_PROBLEM_FILE "user_problem_file" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 2c62f1b44..72026d7d1 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -81,6 +81,7 @@ class mip_solver_settings_t { f_t time_limit = std::numeric_limits::infinity(); bool heuristics_only = false; i_t num_cpu_threads = -1; // -1 means use default number of threads in branch and bound + i_t max_cut_passes = 10; // number of cut passes to make bool log_to_console = true; std::string log_file; std::string sol_file; diff --git a/cpp/src/dual_simplex/CMakeLists.txt b/cpp/src/dual_simplex/CMakeLists.txt index a376ee23d..157a00a07 100644 --- a/cpp/src/dual_simplex/CMakeLists.txt +++ b/cpp/src/dual_simplex/CMakeLists.txt @@ -10,6 +10,7 @@ set(DUAL_SIMPLEX_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/basis_updates.cpp ${CMAKE_CURRENT_SOURCE_DIR}/bound_flipping_ratio_test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/cuts.cpp ${CMAKE_CURRENT_SOURCE_DIR}/crossover.cpp ${CMAKE_CURRENT_SOURCE_DIR}/folding.cpp ${CMAKE_CURRENT_SOURCE_DIR}/initial_basis.cpp @@ -33,7 +34,7 @@ set(DUAL_SIMPLEX_SRC_FILES ) # Uncomment to enable debug info -#set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g1") +set_source_files_properties(${DUAL_SIMPLEX_SRC_FILES} DIRECTORY ${CMAKE_SOURCE_DIR} PROPERTIES COMPILE_OPTIONS "-g1") set(CUOPT_SRC_FILES ${CUOPT_SRC_FILES} ${DUAL_SIMPLEX_SRC_FILES} PARENT_SCOPE) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 3e16411f4..fd70194e1 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -1108,6 +1108,212 @@ i_t basis_update_t::lower_triangular_multiply(const csc_matrix_t +i_t basis_update_mpf_t::append_cuts(const csr_matrix_t& cuts_basic) +{ + const i_t m = L0_.m; + + // Solve for U^T W^T = C_B^T + // We do this one row at a time of C_B + csc_matrix_t WT(m, cuts_basic.m, 0); + + i_t WT_nz = 0; + for (i_t k = 0; k < cuts_basic.m; k++) { + sparse_vector_t rhs(cuts_basic, k); + u_transpose_solve(rhs); + WT.col_start[k] = WT_nz; + for (i_t q = 0; q < rhs.i.size(); q++) { + WT.i.push_back(rhs.i[q]); + WT.x.push_back(rhs.x[q]); + WT_nz++; + } + } + WT.col_start[cuts_basic.m] = WT_nz; + + +#ifdef CHECK_W + { + for (i_t k = 0; k < cuts_basic.m; k++) { + std::vector WT_col(m, 0.0); + WT.load_a_column(k, WT_col); + std::vector CBT_col(m, 0.0); + matrix_transpose_vector_multiply(U0_, 1.0, WT_col, 0.0, CBT_col); + sparse_vector_t CBT_col_sparse(cuts_basic, k); + std::vector CBT_col_dense(m); + CBT_col_sparse.to_dense(CBT_col_dense); + for (i_t h = 0; h < m; h++) { + if (std::abs(CBT_col_dense[h] - CBT_col[h]) > 1e-6) { + printf("col %d CBT_col_dense[%d] = %e CBT_col[%d] = %e\n", k, h, CBT_col_dense[h], h, CBT_col[h]); + exit(1); + } + } + } + } +#endif + + csc_matrix_t V(cuts_basic.m, m, 0); + if (num_updates_ > 0) { + // W = V T_0 ... T_{num_updates_ - 1} + // or V = W T_{num_updates_ - 1}^{-1} ... T_0^{-1} + // or V^T = T_0^{-T} ... T_{num_updates_ - 1}^{-T} W^T + // We can compute V^T column by column so that we have + // V^T(:, h) = T_0^{-T} ... T_{num_updates_ - 1}^{-T} W^T(:, h) + // or + // V(h, :) = T_0^{-T} ... T_{num_updates_ - 1}^{-T} W^T(:, h) + // So we can form V row by row in CSR and then covert it to CSC + // for appending to L0 + + csr_matrix_t V_row(cuts_basic.m, m, 0); + i_t V_nz = 0; + const f_t zero_tol = 1e-13; + for (i_t h = 0; h < cuts_basic.m; h++) { + sparse_vector_t rhs(WT, h); + scatter_into_workspace(rhs); + i_t nz = rhs.i.size(); + for (i_t k = num_updates_ - 1; k >= 0; --k) { + // T_k^{-T} = ( I - v u^T/(1 + u^T v)) + // T_k^{-T} * b = b - v * (u^T * b) / (1 + u^T * v) = b - theta * v, theta = u^T b / mu + + const i_t u_col = 2 * k; + const i_t v_col = 2 * k + 1; + const f_t mu = mu_values_[k]; + + // dot = u^T * b + f_t dot = dot_product(u_col, xi_workspace_, x_workspace_); + const f_t theta = dot / mu; + if (std::abs(theta) > zero_tol) { + add_sparse_column(S_, v_col, -theta, xi_workspace_, nz, x_workspace_); + } + } + gather_into_sparse_vector(nz, rhs); + V_row.row_start[h] = V_nz; + for (i_t q = 0; q < rhs.i.size(); q++) { + V_row.j.push_back(rhs.i[q]); + V_row.x.push_back(rhs.x[q]); + V_nz++; + } + } + V_row.row_start[cuts_basic.m] = V_nz; + + V_row.to_compressed_col(V); + + +#ifdef CHECK_V + csc_matrix_t CB_col(cuts_basic.m, m, 0); + cuts_basic.to_compressed_col(CB_col); + for (i_t k = 0; k < m; k++) { + std::vector U_col(m, 0.0); + U0_.load_a_column(k, U_col); + for (i_t h = num_updates_ - 1; h >= 0; --h) { + // T_h = ( I + u_h v_h^T) + // T_h * x = x + u_h * v_h^T * x = x + theta * u_h + const i_t u_col = 2 * h; + const i_t v_col = 2 * h + 1; + f_t theta = dot_product(v_col, U_col); + const i_t col_start = S_.col_start[u_col]; + const i_t col_end = S_.col_start[u_col + 1]; + for (i_t p = col_start; p < col_end; ++p) { + const i_t i = S_.i[p]; + U_col[i] += theta * S_.x[p]; + } + } + std::vector CB_column(cuts_basic.m, 0.0); + matrix_vector_multiply(V, 1.0, U_col, 0.0, CB_column); + std::vector CB_col_dense(cuts_basic.m); + CB_col.load_a_column(k, CB_col_dense); + for (i_t l = 0; l < cuts_basic.m; l++) { + if (std::abs(CB_col_dense[l] - CB_column[l]) > 1e-6) { + printf("col %d CB_col_dense[%d] = %e CB_column[%d] = %e\n", k, l, CB_col_dense[l], l, CB_column[l]); + exit(1); + } + } + } +#endif + } else { + // W = V + WT.transpose(V); + } + + // Extend u_i, v_i for i = 0, ..., num_updates_ - 1 + S_.m += cuts_basic.m; + + // Adjust L and U + // L = [ L0 0 ] + // [ V I ] + + i_t V_nz = V.col_start[m]; + i_t L_nz = L0_.col_start[m]; + csc_matrix_t new_L(m + cuts_basic.m, m + cuts_basic.m, L_nz + V_nz + cuts_basic.m); + i_t predicted_nz = L_nz + V_nz + cuts_basic.m; + L_nz = 0; + for (i_t j = 0; j < m; ++j) { + new_L.col_start[j] = L_nz; + const i_t col_start = L0_.col_start[j]; + const i_t col_end = L0_.col_start[j + 1]; + for (i_t p = col_start; p < col_end; ++p) { + new_L.i[L_nz] = L0_.i[p]; + new_L.x[L_nz] = L0_.x[p]; + L_nz++; + } + const i_t V_col_start = V.col_start[j]; + const i_t V_col_end = V.col_start[j + 1]; + for (i_t p = V_col_start; p < V_col_end; ++p) { + new_L.i[L_nz] = V.i[p] + m; + new_L.x[L_nz] = V.x[p]; + L_nz++; + } + } + for (i_t j = m; j < m + cuts_basic.m; ++j) { + new_L.col_start[j] = L_nz; + new_L.i[L_nz] = j; + new_L.x[L_nz] = 1.0; + L_nz++; + } + new_L.col_start[m + cuts_basic.m] = L_nz; + if (L_nz != predicted_nz) { + printf("L_nz %d predicted_nz %d\n", L_nz, predicted_nz); + exit(1); + } + + L0_ = new_L; + + // Adjust U + // U = [ U0 0 ] + // [ 0 I ] + + i_t U_nz = U0_.col_start[m]; + U0_.col_start.resize(m + cuts_basic.m + 1); + U0_.i.resize(U_nz + cuts_basic.m); + U0_.x.resize(U_nz + cuts_basic.m); + for (i_t k = m; k < m + cuts_basic.m; ++k) { + U0_.col_start[k] = U_nz; + U0_.i[U_nz] = k; + U0_.x[U_nz] = 1.0; + U_nz++; + } + U0_.col_start[m + cuts_basic.m] = U_nz; + U0_.n = m + cuts_basic.m; + U0_.m = m + cuts_basic.m; + + compute_transposes(); + + // Adjust row_permutation_ and inverse_row_permutation_ + row_permutation_.resize(m + cuts_basic.m); + inverse_row_permutation_.resize(m + cuts_basic.m); + for (i_t k = m; k < m + cuts_basic.m; ++k) { + row_permutation_[k] = k; + } + inverse_permutation(row_permutation_, inverse_row_permutation_); + + // Adjust workspace sizes + xi_workspace_.resize(2 * (m + cuts_basic.m), 0); + x_workspace_.resize(m + cuts_basic.m, 0.0); + + return 0; +} + template void basis_update_mpf_t::gather_into_sparse_vector(i_t nz, sparse_vector_t& out) const @@ -2061,7 +2267,7 @@ int basis_update_mpf_t::refactor_basis( q, deficient, slacks_needed) == -1) { - settings.log.debug("Initial factorization failed\n"); + settings.log.printf("Initial factorization failed\n"); basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); #ifdef CHECK_BASIS_REPAIR diff --git a/cpp/src/dual_simplex/basis_updates.hpp b/cpp/src/dual_simplex/basis_updates.hpp index 078dfffeb..283a1513e 100644 --- a/cpp/src/dual_simplex/basis_updates.hpp +++ b/cpp/src/dual_simplex/basis_updates.hpp @@ -291,6 +291,8 @@ class basis_update_mpf_t { reset_stats(); } + i_t append_cuts(const csr_matrix_t& cuts_basic); + f_t estimate_solution_density(f_t rhs_nz, f_t sum, i_t& num_calls, bool& use_hypersparse) const { num_calls++; diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 9f207b6a6..6f6917fec 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -5,9 +5,9 @@ */ /* clang-format on */ -#include -#include #include +#include +#include #include #include #include @@ -19,6 +19,9 @@ #include #include +#include + +#include #include #include #include @@ -206,6 +209,10 @@ branch_and_bound_t::branch_and_bound_t( { stats_.start_time = tic(); dualize_info_t dualize_info; +#ifdef PRINT_A + settings_.log.printf("A"); + original_problem_.A.print_matrix(); +#endif convert_user_problem(original_problem_, settings_, original_lp_, new_slacks_, dualize_info); full_variable_types(original_problem_, original_lp_, var_types_); @@ -250,6 +257,7 @@ i_t branch_and_bound_t::get_heap_size() template void branch_and_bound_t::set_new_solution(const std::vector& solution) { + mutex_original_lp_.lock(); if (solution.size() != original_problem_.num_cols) { settings_.log.printf( "Solution size mismatch %ld %d\n", solution.size(), original_problem_.num_cols); @@ -258,16 +266,22 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu crush_primal_solution( original_problem_, original_lp_, solution, new_slacks_, crushed_solution); f_t obj = compute_objective(original_lp_, crushed_solution); + mutex_original_lp_.unlock(); bool is_feasible = false; bool attempt_repair = false; mutex_upper_.lock(); - if (obj < upper_bound_) { + f_t current_upper_bound = upper_bound_; + mutex_upper_.unlock(); + if (obj < current_upper_bound) { f_t primal_err; f_t bound_err; i_t num_fractional; + mutex_original_lp_.lock(); is_feasible = check_guess( original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); - if (is_feasible) { + mutex_original_lp_.unlock(); + mutex_upper_.lock(); + if (is_feasible && obj < upper_bound_) { upper_bound_ = obj; incumbent_.set_incumbent_solution(obj, crushed_solution); } else { @@ -282,8 +296,8 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu num_fractional); } } + mutex_upper_.unlock(); } - mutex_upper_.unlock(); if (is_feasible) { if (status_ == mip_exploration_status_t::RUNNING) { @@ -292,7 +306,7 @@ void branch_and_bound_t::set_new_solution(const std::vector& solu std::string gap = user_mip_gap(user_obj, user_lower); settings_.log.printf( - "H %+13.6e %+10.6e %s %9.2f\n", + "H %+13.6e %+10.6e %s %9.2f\n", user_obj, user_lower, gap.c_str(), @@ -405,7 +419,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(), @@ -512,12 +526,13 @@ void branch_and_bound_t::add_feasible_solution(f_t leaf_objective, f_t lower_bound = get_lower_bound(); f_t obj = compute_user_objective(original_lp_, upper_bound_); f_t lower = compute_user_objective(original_lp_, lower_bound); - settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf("%c%10d %10lu %+13.6e %+10.6e %6d %6d %7.1e %s %9.2f\n", thread_type, nodes_explored, nodes_unexplored, obj, lower, + 0, leaf_depth, nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, user_mip_gap(obj, lower).c_str(), @@ -660,7 +675,7 @@ node_status_t branch_and_bound_t::solve_node(search_tree_t& assert(leaf_vstatus.size() == leaf_problem.num_cols); search_tree.branch( - node_ptr, branch_var, leaf_solution.x[branch_var], leaf_vstatus, original_lp_, log); + node_ptr, branch_var, leaf_solution.x[branch_var], leaf_num_fractional, leaf_vstatus, original_lp_, log); node_ptr->status = node_status_t::HAS_CHILDREN; return node_status_t::HAS_CHILDREN; @@ -736,11 +751,12 @@ void branch_and_bound_t::exploration_ramp_up(search_tree_t* f_t user_lower = compute_user_objective(original_lp_, root_objective_); std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %6d %7.1e %s %9.2f\n", nodes_explored, nodes_unexplored, obj, user_lower, + node->integer_infeasible, node->depth, nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, gap_user.c_str(), @@ -836,11 +852,12 @@ void branch_and_bound_t::explore_subtree(i_t task_id, f_t obj = compute_user_objective(original_lp_, upper_bound); f_t user_lower = compute_user_objective(original_lp_, get_lower_bound()); std::string gap_user = user_mip_gap(obj, user_lower); - settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %7.1e %s %9.2f\n", + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %6d %7.1e %s %9.2f\n", nodes_explored, nodes_unexplored, obj, user_lower, + node_ptr->integer_infeasible, node_ptr->depth, nodes_explored > 0 ? stats_.total_lp_iters / nodes_explored : 0, gap_user.c_str(), @@ -1060,10 +1077,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); settings_.log.printf("Solving LP root relaxation\n"); + i_t original_rows = original_lp_.num_rows; simplex_solver_settings_t lp_settings = settings_; lp_settings.inside_mip = 1; - lp_status_t root_status = solve_linear_program_advanced( - original_lp_, stats_.start_time, lp_settings, root_relax_soln_, root_vstatus_, edge_norms_); + lp_settings.scale_columns = false; + std::vector basic_list(original_lp_.num_rows); + std::vector nonbasic_list; + basis_update_mpf_t basis_update(original_lp_.num_rows, settings_.refactor_frequency); + lp_status_t root_status = solve_linear_program_with_advanced_basis( + original_lp_, stats_.start_time, lp_settings, root_relax_soln_, basis_update, basic_list, nonbasic_list, 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) { @@ -1111,31 +1133,195 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } std::vector fractional; - const i_t num_fractional = + i_t num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); - if (num_fractional == 0) { - mutex_upper_.lock(); - incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); - upper_bound_ = root_objective_; - mutex_upper_.unlock(); - // We should be done here - uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); - solution.objective = incumbent_.objective; - solution.lower_bound = root_objective_; - solution.nodes_explored = 0; - solution.simplex_iterations = root_relax_soln_.iterations; - settings_.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", - compute_user_objective(original_lp_, root_objective_), - toc(stats_.start_time)); + csc_matrix_t Arow(1, 1, 1); + original_lp_.A.transpose(Arow); - if (settings_.solution_callback != nullptr) { - settings_.solution_callback(solution.x, solution.objective); - } - if (settings_.heuristic_preemption_callback != nullptr) { - settings_.heuristic_preemption_callback(); + status_ = mip_exploration_status_t::RUNNING; + lower_bound_ceiling_ = inf; + + if (num_fractional != 0) { + settings_.log.printf( + " | Explored | Unexplored | Objective | Bound | IntInf | Depth | Iter/Node | Gap " + "| Time |\n"); + } + + cut_pool_t cut_pool(original_lp_.num_cols, settings_); + cut_generation_t cut_generation(cut_pool); + + for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { + if (num_fractional == 0) { +#ifdef PRINT_SOLUTION + for (i_t j = 0; j < original_lp_.num_cols; j++) { + if (var_types_[j] == variable_type_t::INTEGER) { + settings_.log.printf("Variable %d type %d val %e\n", j, var_types_[j], root_relax_soln_.x[j]); + } + } +#endif + mutex_upper_.lock(); + incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x); + upper_bound_ = root_objective_; + mutex_upper_.unlock(); + // We should be done here + uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x); + solution.objective = incumbent_.objective; + solution.lower_bound = root_objective_; + solution.nodes_explored = 0; + solution.simplex_iterations = root_relax_soln_.iterations; + settings_.log.printf("Optimal solution found at root node. Objective %.16e. Time %.2f.\n", + compute_user_objective(original_lp_, root_objective_), + toc(stats_.start_time)); + + if (settings_.solution_callback != nullptr) { + settings_.solution_callback(solution.x, solution.objective); + } + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + return mip_status_t::OPTIMAL; + } else { +#ifdef PRINT_FRACTIONAL_INFO + settings_.log.printf("Found %d fractional variables on cut pass %d\n", num_fractional, cut_pass); + for (i_t j: fractional) { + settings_.log.printf("Fractional variable %d lower %e value %e upper %e\n", j, original_lp_.lower[j], root_relax_soln_.x[j], original_lp_.upper[j]); + } +#endif + + // Generate cuts and add them to the cut pool + cut_generation.generate_cuts(original_lp_, settings_, Arow, new_slacks_, var_types_, basis_update, root_relax_soln_.x, basic_list, nonbasic_list); + + // Score the cuts + cut_pool.score_cuts(root_relax_soln_.x); + // Get the best cuts from the cut pool + csr_matrix_t cuts_to_add(0, original_lp_.num_cols, 0); + std::vector cut_rhs; + i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs); + + cuts_to_add.check_matrix(); + +#ifdef PRINT_MIN_CUT_VIOLATION + f_t min_cut_violation = minimum_violation(cuts_to_add, cut_rhs, root_relax_soln_.x); + settings_.log.printf("Min cut violation %e\n", min_cut_violation); +#endif + + // Resolve the LP with the new cuts + settings_.log.debug("Solving LP with %d cuts (%d cut nonzeros). Cuts in pool %d. Total constraints %d\n", + num_cuts, + cuts_to_add.row_start[cuts_to_add.m], + cut_pool.pool_size(), + cuts_to_add.m + original_lp_.num_rows); + lp_settings.log.log = false; + + mutex_original_lp_.lock(); + i_t add_cuts_status = add_cuts(settings_, + cuts_to_add, + cut_rhs, + original_lp_, + new_slacks_, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + mutex_original_lp_.unlock(); + if (add_cuts_status != 0) { + settings_.log.printf("Failed to add cuts\n"); + exit(1); + } + + // Try to do bound strengthening + var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; +#ifdef CHECK_MATRICES + settings_.log.printf("Before A check\n"); + original_lp_.A.check_matrix(); +#endif + original_lp_.A.transpose(Arow); + bool feasible = + bound_strengthening(row_sense, settings_, original_lp_, Arow, var_types_, bounds_changed); + + if (!feasible) { + settings_.log.printf("Bound strengthening failed\n"); + exit(1); + } + + // Adjust the solution + root_relax_soln_.x.resize(original_lp_.num_cols, 0.0); + root_relax_soln_.y.resize(original_lp_.num_rows, 0.0); + root_relax_soln_.z.resize(original_lp_.num_cols, 0.0); + + // For now just clear the edge norms + edge_norms_.clear(); + i_t iter = 0; + bool initialize_basis = false; + dual::status_t cut_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + + settings_.log.debug("Cut LP iterations %d. A nz %d\n", + iter, + original_lp_.A.col_start[original_lp_.A.n]); + stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + + if (cut_status != dual::status_t::OPTIMAL) { + settings_.log.printf("Cut status %d\n", cut_status); + exit(1); + } + + local_lower_bounds_.assign(settings_.num_bfs_threads, root_objective_); + + remove_cuts(original_lp_, + settings_, + Arow, + new_slacks_, + original_rows, + var_types_, + root_vstatus_, + root_relax_soln_.x, + root_relax_soln_.y, + root_relax_soln_.z, + basic_list, + nonbasic_list, + basis_update); + + fractional.clear(); + num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + + // TODO: Get upper bound from heuristics + f_t upper_bound = get_upper_bound(); + f_t obj = num_fractional != 0 ? get_upper_bound() : compute_user_objective(original_lp_, root_objective_); + f_t user_obj = compute_user_objective(original_lp_, obj); + f_t user_lower = compute_user_objective(original_lp_, root_objective_); + std::string gap = num_fractional != 0 ? user_mip_gap(user_obj, user_lower) : "0.0%"; + + + settings_.log.printf(" %10d %10lu %+13.6e %+10.6e %6d %6d %7.1e %s %9.2f\n", + 0, + 0, + user_obj, + user_lower, + num_fractional, + 0, + stats_.total_lp_iters.load(), + gap.c_str(), + toc(stats_.start_time)); } - return mip_status_t::OPTIMAL; } pc_.resize(original_lp_.num_cols); @@ -1163,6 +1349,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut search_tree.branch(&search_tree.root, branch_var, root_relax_soln_.x[branch_var], + num_fractional, root_vstatus_, original_lp_, log); @@ -1173,12 +1360,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut settings_.num_diving_threads, settings_.num_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; @@ -1186,8 +1369,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut stats_.last_log = tic(); active_subtrees_ = 0; min_diving_queue_size_ = 4 * settings_.num_diving_threads; - status_ = mip_exploration_status_t::RUNNING; - lower_bound_ceiling_ = inf; + #pragma omp parallel num_threads(settings_.num_threads) { diff --git a/cpp/src/dual_simplex/branch_and_bound.hpp b/cpp/src/dual_simplex/branch_and_bound.hpp index 5b304addd..ccbad335a 100644 --- a/cpp/src/dual_simplex/branch_and_bound.hpp +++ b/cpp/src/dual_simplex/branch_and_bound.hpp @@ -149,6 +149,9 @@ class branch_and_bound_t { // Local lower bounds for each thread std::vector> local_lower_bounds_; + // Mutex for the original LP + omp_mutex_t mutex_original_lp_; + // Mutex for upper bound omp_mutex_t mutex_upper_; diff --git a/cpp/src/dual_simplex/cuts.cpp b/cpp/src/dual_simplex/cuts.cpp new file mode 100644 index 000000000..8fbf1a275 --- /dev/null +++ b/cpp/src/dual_simplex/cuts.cpp @@ -0,0 +1,1155 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + + +namespace cuopt::linear_programming::dual_simplex { + + +template +void cut_pool_t::add_cut(const sparse_vector_t& cut, f_t rhs) +{ + // TODO: Need to deduplicate cuts and only add if the cut is not already in the pool + + for (i_t p = 0; p < cut.i.size(); p++) { + const i_t j = cut.i[p]; + if (j >= original_vars_) { + settings_.log.printf( + "Cut has variable %d that is greater than original_vars_ %d\n", j, original_vars_); + return; + } + } + + cut_storage_.append_row(cut); + rhs_storage_.push_back(rhs); + cut_age_.push_back(0); +} + + +template +f_t cut_pool_t::cut_distance(i_t row, const std::vector& x, f_t& cut_violation, f_t &cut_norm) +{ + const i_t row_start = cut_storage_.row_start[row]; + const i_t row_end = cut_storage_.row_start[row + 1]; + f_t cut_x = 0.0; + f_t dot = 0.0; + for (i_t p = row_start; p < row_end; p++) { + const i_t j = cut_storage_.j[p]; + const f_t cut_coeff = cut_storage_.x[p]; + cut_x += cut_coeff * x[j]; + dot += cut_coeff * cut_coeff; + } + cut_violation = rhs_storage_[row] - cut_x; + cut_norm = std::sqrt(dot); + const f_t distance = cut_violation / cut_norm; + return distance; +} + +template +f_t cut_pool_t::cut_density(i_t row) +{ + const i_t row_start = cut_storage_.row_start[row]; + const i_t row_end = cut_storage_.row_start[row + 1]; + const i_t cut_nz = row_end - row_start; + const i_t original_vars = original_vars_; + return static_cast(cut_nz) / original_vars; +} + +template +f_t cut_pool_t::cut_orthogonality(i_t i, i_t j) +{ + const i_t i_start = cut_storage_.row_start[i]; + const i_t i_end = cut_storage_.row_start[i + 1]; + const i_t i_nz = i_end - i_start; + const i_t j_start = cut_storage_.row_start[j]; + const i_t j_end = cut_storage_.row_start[j + 1]; + const i_t j_nz = j_end - j_start; + + f_t dot = sparse_dot(cut_storage_.j.data() + i_start, cut_storage_.x.data() + i_start, i_nz, + cut_storage_.j.data() + j_start, cut_storage_.x.data() + j_start, j_nz); + + f_t norm_i = cut_norms_[i]; + f_t norm_j = cut_norms_[j]; + return 1.0 - std::abs(dot) / (norm_i * norm_j); +} + +template +void cut_pool_t::score_cuts(std::vector& x_relax) +{ + const f_t weight_distance = 1.0; + const f_t weight_orthogonality = 1.0; + cut_distances_.resize(cut_storage_.m, 0.0); + cut_norms_.resize(cut_storage_.m, 0.0); + cut_orthogonality_.resize(cut_storage_.m, 1); + cut_scores_.resize(cut_storage_.m, 0.0); + for (i_t i = 0; i < cut_storage_.m; i++) { + f_t violation; + cut_distances_[i] = cut_distance(i, x_relax, violation, cut_norms_[i]); + cut_scores_[i] = weight_distance * cut_distances_[i] + weight_orthogonality * cut_orthogonality_[i]; + //settings_.log.printf("Cut %d distance %e violation %e orthogonality %e score %e\n", i, cut_distances_[i], violation, cut_orthogonality_[i], cut_scores_[i]); + } + + std::vector sorted_indices(cut_storage_.m); + std::iota(sorted_indices.begin(), sorted_indices.end(), 0); + std::sort(sorted_indices.begin(), sorted_indices.end(), [&](i_t a, i_t b) { + return cut_scores_[a] > cut_scores_[b]; + }); + + std::vector indices; + indices.reserve(sorted_indices.size()); + + + const i_t max_cuts = 2000; + const f_t min_orthogonality = 0.5; + const f_t min_cut_distance = 1e-4; + best_cuts_.reserve(std::min(max_cuts, cut_storage_.m)); + + while (scored_cuts_ < max_cuts && !sorted_indices.empty()) { + const i_t i = sorted_indices[0]; + + if (cut_distances_[i] <= min_cut_distance) { + break; + } + + if (cut_age_[i] > 0) { + settings_.log.printf("Adding cut with age %d\n", cut_age_[i]); + } + //settings_.log.printf("Scored cuts %d. Adding cut %d score %e\n", scored_cuts_, i, cut_scores_[i]); + + best_cuts_.push_back(i); + scored_cuts_++; + + // Recompute the orthogonality for the remaining cuts + for (i_t k = 1; k < sorted_indices.size(); k++) { + const i_t j = sorted_indices[k]; + cut_orthogonality_[j] = std::min(cut_orthogonality_[j], cut_orthogonality(i, j)); + if (cut_orthogonality_[j] >= min_orthogonality) { + indices.push_back(j); + cut_scores_[j] = weight_distance * cut_distances_[j] + weight_orthogonality * cut_orthogonality_[j]; + //settings_.log.printf("Recomputed cut %d score %e\n", j, cut_scores_[j]); + } + } + + sorted_indices = indices; + indices.clear(); + //settings_.log.printf("Sorting %d cuts\n", sorted_indices.size()); + + std::sort(sorted_indices.begin(), sorted_indices.end(), [&](i_t a, i_t b) { + return cut_scores_[a] > cut_scores_[b]; + }); + } +} + +template +i_t cut_pool_t::get_best_cuts(csr_matrix_t& best_cuts, std::vector& best_rhs) +{ + best_cuts.m = 0; + best_cuts.n = original_vars_; + best_cuts.row_start.clear(); + best_cuts.j.clear(); + best_cuts.x.clear(); + best_cuts.row_start.reserve(scored_cuts_ + 1); + best_cuts.row_start.push_back(0); + + for (i_t i: best_cuts_) { + sparse_vector_t cut(cut_storage_, i); + cut.negate(); + best_cuts.append_row(cut); + //settings_.log.printf("Best cuts nz %d\n", best_cuts.row_start[best_cuts.m]); + best_rhs.push_back(-rhs_storage_[i]); + } + + return static_cast(best_cuts_.size()); +} + + +template +void cut_pool_t::age_cuts() +{ + for (i_t i = 0; i < cut_age_.size(); i++) { + cut_age_[i]++; + } +} + +template +void cut_pool_t::drop_cuts() +{ + // TODO: Implement this +} + +template +void cut_generation_t::generate_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + basis_update_mpf_t& basis_update, + const std::vector& xstar, + const std::vector& basic_list, + const std::vector& nonbasic_list) +{ + // Generate Gomory Cuts + generate_gomory_cuts( + lp, settings, Arow, new_slacks, var_types, basis_update, xstar, basic_list, nonbasic_list); + + + // Generate MIR cuts + // generate_mir_cuts(lp, settings, Arow, var_types, xstar); +} + +template +void cut_generation_t::generate_mir_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + const std::vector& xstar) +{ + mixed_integer_rounding_cut_t mir(lp.num_cols, settings); + mir.initialize(lp, new_slacks, xstar); + + for (i_t i = 0; i < lp.num_rows; i++) { + sparse_vector_t inequality(Arow, i); + f_t inequality_rhs = lp.rhs[i]; + + const i_t row_start = Arow.col_start[i]; + const i_t row_end = Arow.col_start[i + 1]; + i_t last_slack = -1; + for (i_t p = row_start; p < row_end; p++) { + const i_t j = Arow.i[p]; + const f_t a = Arow.x[p]; + if (var_types[j] == variable_type_t::CONTINUOUS && a == 1.0 && lp.lower[j] == 0.0) { + last_slack = j; + } + } + + if (last_slack != -1) { + // Remove the slack from the equality to get an inequality + for (i_t k = 0; k < inequality.i.size(); k++) { + const i_t j = inequality.i[k]; + if (j == last_slack) { + inequality.x[k] = 0.0; + } + } + + // inequaility'*x <= inequality_rhs + // But for MIR we need: inequality'*x >= inequality_rhs + inequality_rhs *= -1; + inequality.negate(); + + sparse_vector_t cut(lp.num_cols, 0); + f_t cut_rhs; + i_t mir_status = mir.generate_cut(inequality, inequality_rhs, lp.upper, lp.lower, var_types, cut, cut_rhs); + if (mir_status == 0) { + f_t dot = 0.0; + f_t cut_norm = 0.0; + for (i_t k = 0; k < cut.i.size(); k++) { + const i_t jj = cut.i[k]; + const f_t aj = cut.x[k]; + dot += aj * xstar[jj]; + cut_norm += aj * aj; + } + if (dot >= cut_rhs) { + continue; + } + } + + settings.log.printf("Adding MIR cut %d\n", i); + cut_pool_.add_cut(cut, cut_rhs); + } + } +} + + +template +void cut_generation_t::generate_gomory_cuts( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + basis_update_mpf_t& basis_update, + const std::vector& xstar, + const std::vector& basic_list, + const std::vector& nonbasic_list) +{ + mixed_integer_gomory_base_inequality_t gomory(lp, basis_update, nonbasic_list); + mixed_integer_rounding_cut_t mir(lp.num_cols, settings); + + mir.initialize(lp, new_slacks, xstar); + + for (i_t i = 0; i < lp.num_rows; i++) { + sparse_vector_t inequality(lp.num_cols, 0); + f_t inequality_rhs; + const i_t j = basic_list[i]; + if (var_types[j] != variable_type_t::INTEGER) { continue; } + const f_t x_j = xstar[j]; + if (std::abs(x_j - std::round(x_j)) < settings.integer_tol) { continue; } + i_t gomory_status = gomory.generate_base_inequality(lp, + settings, + Arow, + var_types, + basis_update, + xstar, + basic_list, + nonbasic_list, + i, + inequality, + inequality_rhs); + if (gomory_status == 0) { + // Given the base inequality, generate a MIR cut + sparse_vector_t cut_A(lp.num_cols, 0); + f_t cut_A_rhs; + i_t mir_status = + mir.generate_cut(inequality, inequality_rhs, lp.upper, lp.lower, var_types, cut_A, cut_A_rhs); + bool A_valid = false; + f_t cut_A_distance = 0.0; + if (mir_status == 0) { + mir.substitute_slacks(lp, Arow, cut_A, cut_A_rhs); + // Check that the cut is violated + f_t dot = cut_A.dot(xstar); + f_t cut_norm = cut_A.norm2_squared(); + if (dot >= cut_A_rhs) { + settings.log.printf("Cut %d is not violated. Skipping\n", i); + continue; + } + cut_A_distance = (cut_A_rhs - dot) / std::sqrt(cut_norm); + A_valid = true; + //cut_pool_.add_cut(lp.num_cols, cut, cut_rhs); + } + + // Negate the base inequality + inequality.negate(); + inequality_rhs *= -1; + + sparse_vector_t cut_B(lp.num_cols, 0); + f_t cut_B_rhs; + + mir_status = + mir.generate_cut(inequality, inequality_rhs, lp.upper, lp.lower, var_types, cut_B, cut_B_rhs); + bool B_valid = false; + f_t cut_B_distance = 0.0; + if (mir_status == 0) { + mir.substitute_slacks(lp, Arow, cut_B, cut_B_rhs); + // Check that the cut is violated + f_t dot = cut_B.dot(xstar); + f_t cut_norm = cut_B.norm2_squared(); + if (dot >= cut_B_rhs) { + settings.log.printf("Cut %d is not violated. Skipping\n", i); + continue; + } + cut_B_distance = (cut_B_rhs - dot) / std::sqrt(cut_norm); + B_valid = true; + // cut_pool_.add_cut(lp.num_cols, cut_B, cut_B_rhs); + } + + if ((cut_A_distance > cut_B_distance) && A_valid) { + cut_pool_.add_cut(cut_A, cut_A_rhs); + } else if (B_valid) { + cut_pool_.add_cut(cut_B, cut_B_rhs); + } + } + } +} + +template +i_t mixed_integer_gomory_base_inequality_t::generate_base_inequality( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& var_types, + basis_update_mpf_t& basis_update, + const std::vector& xstar, + const std::vector& basic_list, + const std::vector& nonbasic_list, + i_t i, + sparse_vector_t& inequality, + f_t& inequality_rhs) +{ + // Let's look for Gomory cuts + const i_t j = basic_list[i]; + if (var_types[j] != variable_type_t::INTEGER) { return -1; } + const f_t x_j = xstar[j]; + if (std::abs(x_j - std::round(x_j)) < settings.integer_tol) { return -1; } +#ifdef PRINT_CUT_INFO + settings_.log.printf("Generating cut for variable %d relaxed value %e row %d\n", j, x_j, i); +#endif +#ifdef PRINT_BASIS + for (i_t h = 0; h < basic_list.size(); h++) { + settings_.log.printf("basic_list[%d] = %d\n", h, basic_list[h]); + } +#endif + + // Solve B^T u_bar = e_i + sparse_vector_t e_i(lp.num_rows, 1); + e_i.i[0] = i; + e_i.x[0] = 1.0; + sparse_vector_t u_bar(lp.num_rows, 0); + basis_update.b_transpose_solve(e_i, u_bar); + + +#ifdef CHECK_B_TRANSPOSE_SOLVE + std::vector u_bar_dense(lp.num_rows); + u_bar.to_dense(u_bar_dense); + + std::vector BTu_bar(lp.num_rows); + b_transpose_multiply(lp, basic_list, u_bar_dense, BTu_bar); + for (i_t k = 0; k < lp.num_rows; k++) { + if (k == i) { + if (std::abs(BTu_bar[k] - 1.0) > 1e-6) { + settings_.log.printf("BTu_bar[%d] = %e i %d\n", k, BTu_bar[k], i); + exit(1); + } + } else { + if (std::abs(BTu_bar[k]) > 1e-6) { + settings_.log.printf("BTu_bar[%d] = %e i %d\n", k, BTu_bar[k], i); + exit(1); + } + } + } +#endif + + // Compute a_bar = N^T u_bar + // TODO: This is similar to a function in phase2 of dual simplex. See if it can be reused. + const i_t nz_ubar = u_bar.i.size(); + std::vector abar_indices; + abar_indices.reserve(nz_ubar); + for (i_t k = 0; k < nz_ubar; k++) { + const i_t ii = u_bar.i[k]; + const f_t u_bar_i = u_bar.x[k]; + const i_t row_start = Arow.col_start[ii]; + const i_t row_end = Arow.col_start[ii + 1]; + for (i_t p = row_start; p < row_end; p++) { + const i_t jj = Arow.i[p]; + if (nonbasic_mark_[jj] == 1) { + x_workspace_[jj] += u_bar_i * Arow.x[p]; + if (!x_mark_[jj]) { + x_mark_[jj] = 1; + abar_indices.push_back(jj); + } + } + } + } + + sparse_vector_t a_bar(lp.num_cols, abar_indices.size() + 1); + for (i_t k = 0; k < abar_indices.size(); k++) { + const i_t jj = abar_indices[k]; + a_bar.i[k] = jj; + a_bar.x[k] = x_workspace_[jj]; + } + + // Clear the workspace + for (i_t jj : abar_indices) { + x_workspace_[jj] = 0.0; + x_mark_[jj] = 0; + } + abar_indices.clear(); + + // We should now have the base inequality + // x_j + a_bar^T x_N >= b_bar_i + // We add x_j into a_bar so that everything is in a single sparse_vector_t + a_bar.i[a_bar.i.size() - 1] = j; + a_bar.x[a_bar.x.size() - 1] = 1.0; + +#ifdef CHECK_A_BAR_DENSE_DOT + std::vector a_bar_dense(lp.num_cols); + a_bar.to_dense(a_bar_dense); + + f_t a_bar_dense_dot = dot(a_bar_dense, xstar); + if (std::abs(a_bar_dense_dot - b_bar[i]) > 1e-6) { + settings_.log.printf("a_bar_dense_dot = %e b_bar[%d] = %e\n", a_bar_dense_dot, i, b_bar[i]); + settings_.log.printf("x_j %e b_bar_i %e\n", x_j, b_bar[i]); + exit(1); + } +#endif + + // We have that x_j + a_bar^T x_N == b_bar_i + // So x_j + a_bar^T x_N >= b_bar_i + // And x_j + a_bar^T x_N <= b_bar_i + // Or -x_j - a_bar^T x_N >= -b_bar_i + +#ifdef PRINT_CUT + { + settings_.log.printf("Cut %d\n", i); + for (i_t k = 0; k < a_bar.i.size(); k++) { + const i_t jj = a_bar.i[k]; + const f_t aj = a_bar.x[k]; + settings_.log.printf("(%d, %e) ", jj, aj); + } + settings_.log.printf("\nEnd cut %d b_bar[%d] = %e\n", i, b_bar[i]); + } +#endif + + // Skip cuts that are shallow + const f_t shallow_tol = 1e-2; + if (std::abs(x_j - std::round(x_j)) < shallow_tol) { + //settings_.log.printf("Skipping shallow cut %d. b_bar[%d] = %e x_j %e\n", i, i, b_bar[i], x_j); + return -1; + } + + const f_t f_val = b_bar_[i] - std::floor(b_bar_[i]); + if (f_val < 0.01 || f_val > 0.99) { + //settings_.log.printf("Skipping cut %d. b_bar[%d] = %e f_val %e\n", i, i, b_bar[i], f_val); + return -1; + } + +#ifdef PRINT_BASE_INEQUALITY + // Print out the base inequality + for (i_t k = 0; k < a_bar.i.size(); k++) { + const i_t jj = a_bar.i[k]; + const f_t aj = a_bar.x[k]; + settings_.log.printf("a_bar[%d] = %e\n", k, aj); + } + settings_.log.printf("b_bar[%d] = %e\n", i, b_bar[i]); +#endif + + inequality = a_bar; + inequality_rhs = b_bar_[i]; + + return 0; +} + +template +void mixed_integer_rounding_cut_t::initialize(const lp_problem_t& lp, + const std::vector& new_slacks, + const std::vector& xstar) +{ + + if (lp.num_cols != num_vars_) { + num_vars_ = lp.num_cols; + x_workspace_.resize(num_vars_, 0.0); + x_mark_.resize(num_vars_, 0); + has_lower_.resize(num_vars_, 0); + has_upper_.resize(num_vars_, 0); + } + + is_slack_.clear(); + is_slack_.resize(num_vars_, 0); + slack_rows_.clear(); + slack_rows_.resize(num_vars_, 0); + + for (i_t j : new_slacks) { + is_slack_[j] = 1; + const i_t col_start = lp.A.col_start[j]; + const i_t i = lp.A.i[col_start]; + slack_rows_[j] = i; + } + + needs_complement_ = false; + for (i_t j = 0; j < lp.num_cols; j++) { + if (lp.lower[j] < 0) { + settings_.log.printf("Variable %d has negative lower bound %e\n", j, lp.lower[j]); + exit(1); + } + const f_t uj = lp.upper[j]; + const f_t lj = lp.lower[j]; + if (uj != inf || lj != 0.0) { needs_complement_ = true; } + const f_t xstar_j = xstar[j]; + if (uj < inf) { + if (uj - xstar_j <= xstar_j - lj) { + has_upper_[j] = 1; + } else { + has_lower_[j] = 1; + } + continue; + } + + if (lj > -inf) { has_lower_[j] = 1; } + } +} + +template +i_t mixed_integer_rounding_cut_t::generate_cut( + const sparse_vector_t& a, + f_t beta, + const std::vector& upper_bounds, + const std::vector& lower_bounds, + const std::vector& var_types, + sparse_vector_t& cut, + f_t& cut_rhs) +{ + auto f = [](f_t q_1, f_t q_2) -> f_t { + f_t q_1_hat = q_1 - std::floor(q_1); + f_t q_2_hat = q_2 - std::floor(q_2); + return std::min(q_1_hat, q_2_hat) + q_2_hat * std::floor(q_1); + }; + + auto h = [](f_t q) -> f_t { return std::max(q, 0.0); }; + + std::vector cut_indices; + cut_indices.reserve(a.i.size()); + f_t R; + if (!needs_complement_) { + R = (beta - std::floor(beta)) * std::ceil(beta); + + for (i_t k = 0; k < a.i.size(); k++) { + const i_t jj = a.i[k]; + f_t aj = a.x[k]; + if (var_types[jj] == variable_type_t::INTEGER) { + x_workspace_[jj] += f(aj, beta); + if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { + x_mark_[jj] = 1; + cut_indices.push_back(jj); + } + } else { + x_workspace_[jj] += h(aj); + if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { + x_mark_[jj] = 1; + cut_indices.push_back(jj); + } + } + } + } else { + // Compute r + f_t r = beta; + for (i_t k = 0; k < a.i.size(); k++) { + const i_t jj = a.i[k]; + if (has_upper_[jj]) { + const f_t uj = upper_bounds[jj]; + r -= uj * a.x[k]; + continue; + } + if (has_lower_[jj]) { + const f_t lj = lower_bounds[jj]; + r -= lj * a.x[k]; + } + } + + // Compute R + R = std::ceil(r) * (r - std::floor(r)); + for (i_t k = 0; k < a.i.size(); k++) { + const i_t jj = a.i[k]; + const f_t aj = a.x[k]; + if (has_upper_[jj]) { + const f_t uj = upper_bounds[jj]; + if (var_types[jj] == variable_type_t::INTEGER) { + R -= f(-aj, r) * uj; + } else { + R -= h(-aj) * uj; + } + } else if (has_lower_[jj]) { + const f_t lj = lower_bounds[jj]; + if (var_types[jj] == variable_type_t::INTEGER) { + R += f(aj, r) * lj; + } else { + R += h(aj) * lj; + } + } + } + + // Compute the cut coefficients + for (i_t k = 0; k < a.i.size(); k++) { + const i_t jj = a.i[k]; + const f_t aj = a.x[k]; + if (has_upper_[jj]) { + if (var_types[jj] == variable_type_t::INTEGER) { + // Upper intersect I + x_workspace_[jj] -= f(-aj, r); + if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { + x_mark_[jj] = 1; + cut_indices.push_back(jj); + } + } else { + // Upper intersect C + f_t h_j = h(-aj); + if (h_j != 0.0) { + x_workspace_[jj] -= h_j; + if (!x_mark_[jj]) { + x_mark_[jj] = 1; + cut_indices.push_back(jj); + } + } + } + } else if (var_types[jj] == variable_type_t::INTEGER) { + // I \ Upper + x_workspace_[jj] += f(aj, r); + if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { + x_mark_[jj] = 1; + cut_indices.push_back(jj); + } + } else { + // C \ Upper + f_t h_j = h(aj); + if (h_j != 0.0) { + x_workspace_[jj] += h_j; + if (!x_mark_[jj]) { + x_mark_[jj] = 1; + cut_indices.push_back(jj); + } + } + } + } + } + + cut.i.reserve(cut_indices.size()); + cut.x.reserve(cut_indices.size()); + for (i_t k = 0; k < cut_indices.size(); k++) { + const i_t jj = cut_indices[k]; + + // Check for small coefficients + const f_t aj = x_workspace_[jj]; + if (std::abs(aj) < 1e-6) { + if (aj >= 0.0 && upper_bounds[jj] < inf) { + // Move this to the right-hand side + R -= aj * upper_bounds[jj]; + continue; + } else if (aj <= 0.0 && lower_bounds[jj] > -inf) { + R += aj * lower_bounds[jj]; + continue; + } else { + } + } + cut.i.push_back(jj); + cut.x.push_back(x_workspace_[jj]); + } + + // Clear the workspace + for (i_t jj : cut_indices) { + x_workspace_[jj] = 0.0; + x_mark_[jj] = 0; + } + + + // The new cut is: g'*x >= R + // But we want to have it in the form h'*x <= b + cut.sort(); + + cut_rhs = R; + return 0; +} + +template +void mixed_integer_rounding_cut_t::substitute_slacks(const lp_problem_t& lp, + csc_matrix_t& Arow, + sparse_vector_t& cut, + f_t& cut_rhs) +{ + // Remove slacks from the cut + // So that the cut is only over the original variables + bool found_slack = false; + i_t cut_nz = 0; + std::vector cut_indices; + cut_indices.reserve(cut.i.size()); + for (i_t k = 0; k < cut.i.size(); k++) { + const i_t j = cut.i[k]; + const f_t cj = cut.x[k]; + if (is_slack_[j]) { + found_slack = true; + // Do the substitution + // Slack variable s_j participates in row i of the constraint matrix + // Row i is of the form: + // sum_{k != j} A(i, k) * x_k + A(i, j) * s_j = rhs_i + /// So we have that + // s_j = rhs_i - sum_{k != j} A(i, k) * x_k + + // Our cut is of the form: + // sum_{k != j} C(k) * x_k + C(j) * s_j >= cut_rhs + // So the cut becomes + // sum_{k != j} C(k) * x_k + C(j) * (rhs_i - sum_{k != j} A(i, k) * x_k) >= cut_rhs + // This is equivalent to: + // sum_{k != j} C(k) * x_k + sum_{k != j} -C(k) * A(i, k) * x_k >= cut_rhs - C(j) * rhs_i + const i_t i = slack_rows_[j]; + cut_rhs -= cj * lp.rhs[i]; + const i_t row_start = Arow.col_start[i]; + const i_t row_end = Arow.col_start[i + 1]; + for (i_t q = row_start; q < row_end; q++) { + const i_t k = Arow.i[q]; + if (k != j) { + const f_t aik = Arow.x[q]; + x_workspace_[k] -= cj * aik; + if (!x_mark_[k]) { + x_mark_[k] = 1; + cut_indices.push_back(k); + cut_nz++; + } + } + } + + } else { + x_workspace_[j] += cj; + if (!x_mark_[j]) { + x_mark_[j] = 1; + cut_indices.push_back(j); + cut_nz++; + } + } + } + + if (found_slack) { + //printf("Found slack. Nz increased from %d to %d: %d\n", cut.i.size(), cut_nz, cut_nz - cut.i.size()); + cut.i.reserve(cut_nz); + cut.x.reserve(cut_nz); + cut.i.clear(); + cut.x.clear(); + + for (i_t k = 0; k < cut_nz; k++) { + const i_t j = cut_indices[k]; + cut.i.push_back(j); + cut.x.push_back(x_workspace_[j]); + } + // Sort the cut + cut.sort(); + + // Clear the workspace + for (i_t jj : cut_indices) { + x_workspace_[jj] = 0.0; + x_mark_[jj] = 0; + } + } +} + +template +i_t add_cuts(const simplex_solver_settings_t& settings, + const csr_matrix_t& cuts, + const std::vector& cut_rhs, + lp_problem_t& lp, + std::vector& new_slacks, + lp_solution_t& solution, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + std::vector& edge_norms) + +{ + // Given a set of cuts: C*x <= d that are currently violated + // by the current solution x* (i.e. C*x* > d), this function + // adds the cuts into the LP and solves again. + +#ifdef CHECK_BASIS + { + csc_matrix_t Btest(lp.num_rows, lp.num_rows, 1); + basis_update.multiply_lu(Btest); + csc_matrix_t B(lp.num_rows, lp.num_rows, 1); + form_b(lp.A, basic_list, B); + csc_matrix_t Diff(lp.num_rows, lp.num_rows, 1); + add(Btest, B, 1.0, -1.0, Diff); + const f_t err = Diff.norm1(); + settings.log.printf("Before || B - L*U || %e\n", err); + if (err > 1e-6) { exit(1); } + } +#endif + + const i_t p = cuts.m; + if (cut_rhs.size() != static_cast(p)) { + settings.log.printf("cut_rhs must have the same number of rows as cuts\n"); + return -1; + } + settings.log.debug("Number of cuts %d\n", p); + settings.log.debug("Original lp rows %d\n", lp.num_rows); + settings.log.debug("Original lp cols %d\n", lp.num_cols); + + csr_matrix_t new_A_row(lp.num_rows, lp.num_cols, 1); + lp.A.to_compressed_row(new_A_row); + + i_t append_status = new_A_row.append_rows(cuts); + if (append_status != 0) { + settings.log.printf("append_rows error: %d\n", append_status); + exit(1); + } + + csc_matrix_t new_A_col(lp.num_rows + p, lp.num_cols, 1); + new_A_row.to_compressed_col(new_A_col); + + // Add in slacks variables for the new rows + lp.lower.resize(lp.num_cols + p); + lp.upper.resize(lp.num_cols + p); + lp.objective.resize(lp.num_cols + p); + i_t nz = new_A_col.col_start[lp.num_cols]; + new_A_col.col_start.resize(lp.num_cols + p + 1); + new_A_col.i.resize(nz + p); + new_A_col.x.resize(nz + p); + i_t k = lp.num_rows; + for (i_t j = lp.num_cols; j < lp.num_cols + p; j++) { + new_A_col.col_start[j] = nz; + new_A_col.i[nz] = k++; + new_A_col.x[nz] = 1.0; + nz++; + lp.lower[j] = 0.0; + lp.upper[j] = inf; + lp.objective[j] = 0.0; + new_slacks.push_back(j); + } + settings.log.debug("Done adding slacks\n"); + new_A_col.col_start[lp.num_cols + p] = nz; + new_A_col.n = lp.num_cols + p; + + lp.A = new_A_col; + i_t old_rows = lp.num_rows; + lp.num_rows += p; + i_t old_cols = lp.num_cols; + lp.num_cols += p; + + lp.rhs.resize(lp.num_rows); + for (i_t k = old_rows; k < old_rows + p; k++) { + const i_t h = k - old_rows; + lp.rhs[k] = cut_rhs[h]; + } + settings.log.debug("Done adding rhs\n"); + + // Construct C_B = C(:, basic_list) + std::vector C_col_degree(lp.num_cols, 0); + i_t cuts_nz = cuts.row_start[p]; + for (i_t q = 0; q < cuts_nz; q++) { + const i_t j = cuts.j[q]; + if (j >= lp.num_cols) { + settings.log.printf("j %d is greater than p %d\n", j, p); + return -1; + } + C_col_degree[j]++; + } + settings.log.debug("Done computing C_col_degree\n"); + + std::vector in_basis(old_cols, -1); + const i_t num_basic = static_cast(basic_list.size()); + i_t C_B_nz = 0; + for (i_t k = 0; k < num_basic; k++) { + const i_t j = basic_list[k]; + if (j < 0 || j >= old_cols) { + settings.log.printf( + "basic_list[%d] = %d is out of bounds %d old_cols %d\n", k, j, j, old_cols); + return -1; + } + in_basis[j] = k; + if (j < cuts.n) { C_B_nz += C_col_degree[j]; } + } + settings.log.debug("Done estimating C_B_nz\n"); + + csr_matrix_t C_B(p, num_basic, C_B_nz); + nz = 0; + for (i_t i = 0; i < p; i++) { + C_B.row_start[i] = nz; + const i_t row_start = cuts.row_start[i]; + const i_t row_end = cuts.row_start[i + 1]; + for (i_t q = row_start; q < row_end; q++) { + const i_t j = cuts.j[q]; + const i_t j_basis = in_basis[j]; + if (j_basis == -1) { continue; } + C_B.j[nz] = j_basis; + C_B.x[nz] = cuts.x[q]; + nz++; + } + } + C_B.row_start[p] = nz; + + if (nz != C_B_nz) { + settings.log.printf("predicted nz %d actual nz %d\n", C_B_nz, nz); + return -1; + } + settings.log.debug("C_B rows %d cols %d nz %d\n", C_B.m, C_B.n, nz); + + // Adjust the basis update to include the new cuts + basis_update.append_cuts(C_B); + + basic_list.resize(lp.num_rows, 0); + i_t h = old_cols; + for (i_t j = old_rows; j < lp.num_rows; j++) { + basic_list[j] = h++; + } + +#ifdef CHECK_BASIS + // Check the basis update + csc_matrix_t Btest(lp.num_rows, lp.num_rows, 1); + basis_update.multiply_lu(Btest); + + csc_matrix_t B(lp.num_rows, lp.num_rows, 1); + form_b(lp.A, basic_list, B); + + csc_matrix_t Diff(lp.num_rows, lp.num_rows, 1); + add(Btest, B, 1.0, -1.0, Diff); + const f_t err = Diff.norm1(); + settings.log.printf("After || B - L*U || %e\n", err); + if (err > 1e-6) { + settings.log.printf("Diff matrix\n"); + // Diff.print_matrix(); + exit(1); + } +#endif + // Adjust the vstatus + vstatus.resize(lp.num_cols); + for (i_t j = old_cols; j < lp.num_cols; j++) { + vstatus[j] = variable_status_t::BASIC; + } + + return 0; +} + +template +void remove_cuts(lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + std::vector& new_slacks, + i_t original_rows, + std::vector& var_types, + std::vector& vstatus, + std::vector& x, + std::vector& y, + std::vector& z, + std::vector& basic_list, + std::vector& nonbasic_list, + basis_update_mpf_t& basis_update) +{ + std::vector cuts_to_remove; + cuts_to_remove.reserve(lp.num_rows - original_rows); + std::vector slacks_to_remove; + slacks_to_remove.reserve(lp.num_rows - original_rows); + const f_t dual_tol = 1e-10; + + std::vector is_slack(lp.num_cols, 0); + for (i_t j : new_slacks) { + is_slack[j] = 1; + } + + for (i_t k = original_rows; k < lp.num_rows; k++) { + if (std::abs(y[k]) < dual_tol) { + const i_t row_start = Arow.col_start[k]; + const i_t row_end = Arow.col_start[k + 1]; + i_t last_slack = -1; + const f_t slack_tol = 1e-3; + for (i_t p = row_start; p < row_end; p++) { + const i_t j = Arow.i[p]; + if (is_slack[j]) { + if (vstatus[j] == variable_status_t::BASIC && x[j] > slack_tol) { last_slack = j; } + } + } + if (last_slack != -1) { + cuts_to_remove.push_back(k); + slacks_to_remove.push_back(last_slack); + } + } + } + + if (cuts_to_remove.size() > 0) { + settings.log.printf("Removing %d cuts\n", cuts_to_remove.size()); + std::vector marked_rows(lp.num_rows, 0); + for (i_t i : cuts_to_remove) { + marked_rows[i] = 1; + } + std::vector marked_cols(lp.num_cols, 0); + for (i_t j : slacks_to_remove) { + marked_cols[j] = 1; + } + + std::vector new_rhs(lp.num_rows - cuts_to_remove.size()); + std::vector new_solution_y(lp.num_rows - cuts_to_remove.size()); + i_t h = 0; + for (i_t i = 0; i < lp.num_rows; i++) { + if (!marked_rows[i]) { + new_rhs[h] = lp.rhs[i]; + new_solution_y[h] = y[i]; + h++; + } + } + + Arow.remove_columns(marked_rows); + Arow.transpose(lp.A); + + std::vector new_objective(lp.num_cols - slacks_to_remove.size()); + std::vector new_lower(lp.num_cols - slacks_to_remove.size()); + std::vector new_upper(lp.num_cols - slacks_to_remove.size()); + std::vector new_var_types(lp.num_cols - slacks_to_remove.size()); + std::vector new_vstatus(lp.num_cols - slacks_to_remove.size()); + std::vector new_basic_list; + new_basic_list.reserve(lp.num_rows - slacks_to_remove.size()); + std::vector new_nonbasic_list; + new_nonbasic_list.reserve(nonbasic_list.size()); + std::vector new_solution_x(lp.num_cols - slacks_to_remove.size()); + std::vector new_solution_z(lp.num_cols - slacks_to_remove.size()); + std::vector new_is_slacks(lp.num_cols - slacks_to_remove.size(), 0); + h = 0; + for (i_t k = 0; k < lp.num_cols; k++) { + if (!marked_cols[k]) { + new_objective[h] = lp.objective[k]; + new_lower[h] = lp.lower[k]; + new_upper[h] = lp.upper[k]; + new_var_types[h] = var_types[k]; + new_vstatus[h] = vstatus[k]; + new_solution_x[h] = x[k]; + new_solution_z[h] = z[k]; + new_is_slacks[h] = is_slack[k]; + if (new_vstatus[h] != variable_status_t::BASIC) { + new_nonbasic_list.push_back(h); + } else { + new_basic_list.push_back(h); + } + h++; + } + } + lp.A.remove_columns(marked_cols); + lp.A.transpose(Arow); + lp.objective = new_objective; + lp.lower = new_lower; + lp.upper = new_upper; + lp.rhs = new_rhs; + var_types = new_var_types; + lp.num_cols = lp.A.n; + lp.num_rows = lp.A.m; + + new_slacks.clear(); + new_slacks.resize(lp.num_cols); + for (i_t j = 0; j < lp.num_cols; j++) { + if (new_is_slacks[j]) { + new_slacks.push_back(j); + } + } + basic_list = new_basic_list; + nonbasic_list = new_nonbasic_list; + vstatus = new_vstatus; + x = new_solution_x; + y = new_solution_y; + z = new_solution_z; + + settings.log.printf("After removal %d rows %d columns %d nonzeros\n", + lp.num_rows, + lp.num_cols, + lp.A.col_start[lp.A.n]); + + basis_update.resize(lp.num_rows); + basis_update.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus); + } +} + + +#ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE +template class cut_pool_t; +template class cut_generation_t; +template class mixed_integer_gomory_base_inequality_t; +template class mixed_integer_rounding_cut_t; + +template +int add_cuts(const simplex_solver_settings_t& settings, + const csr_matrix_t& cuts, + const std::vector& cut_rhs, + lp_problem_t& lp, + std::vector& new_slacks, + lp_solution_t& solution, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + std::vector& edge_norms); + +template +void remove_cuts(lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + std::vector& new_slacks, + int original_rows, + std::vector& var_types, + std::vector& vstatus, + std::vector& x, + std::vector& y, + std::vector& z, + std::vector& basic_list, + std::vector& nonbasic_list, + basis_update_mpf_t& basis_update); +#endif + +} // namespace cuopt::linear_programming::dual_simplex + + diff --git a/cpp/src/dual_simplex/cuts.hpp b/cpp/src/dual_simplex/cuts.hpp new file mode 100644 index 000000000..9113b926e --- /dev/null +++ b/cpp/src/dual_simplex/cuts.hpp @@ -0,0 +1,235 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ +#pragma once + +#include +#include +#include +#include +#include +#include + + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +f_t minimum_violation(const csr_matrix_t& C, + const std::vector& cut_rhs, + const std::vector& x) +{ + // Check to see that this is a cut i.e C*x > d + std::vector Cx(C.m); + csc_matrix_t C_col(C.m, C.n, 0); + C.to_compressed_col(C_col); + matrix_vector_multiply(C_col, 1.0, x, 0.0, Cx); + f_t min_cut_violation = inf; + for (i_t k = 0; k < Cx.size(); k++) { + if (Cx[k] <= cut_rhs[k]) { + printf("C*x <= d for cut %d. C*x %e rhs %e\n", k, Cx[k], cut_rhs[k]); + } + min_cut_violation = std::min(min_cut_violation, Cx[k] - cut_rhs[k]); + } + return min_cut_violation; +} + +template +class cut_pool_t { + public: + cut_pool_t(i_t original_vars, const simplex_solver_settings_t& settings) + : original_vars_(original_vars), + settings_(settings), + cut_storage_(0, original_vars, 0), + rhs_storage_(0), + cut_age_(0), + scored_cuts_(0) + { + } + + // Add a cut in the form: cut'*x >= rhs. + // We expect that the cut is violated by the current relaxation + // cut'*xstart < rhs + void add_cut(const sparse_vector_t& cut, f_t rhs); + + void score_cuts(std::vector& x_relax); + + // We return the cuts in the form best_cuts*x <= best_rhs + i_t get_best_cuts(csr_matrix_t& best_cuts, std::vector& best_rhs); + + void age_cuts(); + + void drop_cuts(); + + i_t pool_size() const { return cut_storage_.m; } + + private: + f_t cut_distance(i_t row, const std::vector& x, f_t& cut_violation, f_t &cut_norm); + f_t cut_density(i_t row); + f_t cut_orthogonality(i_t i, i_t j); + + i_t original_vars_; + const simplex_solver_settings_t& settings_; + + csr_matrix_t cut_storage_; + std::vector rhs_storage_; + std::vector cut_age_; + + i_t scored_cuts_; + std::vector cut_distances_; + std::vector cut_norms_; + std::vector cut_orthogonality_; + std::vector cut_scores_; + std::vector best_cuts_; +}; + +template +class cut_generation_t { + public: + cut_generation_t(cut_pool_t& cut_pool) : cut_pool_(cut_pool) {} + + + void generate_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + basis_update_mpf_t& basis_update, + const std::vector& xstar, + const std::vector& basic_list, + const std::vector& nonbasic_list); + private: + + void generate_gomory_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + basis_update_mpf_t& basis_update, + const std::vector& xstar, + const std::vector& basic_list, + const std::vector& nonbasic_list); + + void generate_mir_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + const std::vector& xstar); + cut_pool_t& cut_pool_; +}; + +template +class mixed_integer_gomory_base_inequality_t { + public: + mixed_integer_gomory_base_inequality_t(const lp_problem_t& lp, + basis_update_mpf_t& basis_update, + const std::vector nonbasic_list) + : b_bar_(lp.num_rows, 0.0), + nonbasic_mark_(lp.num_cols, 0), + x_workspace_(lp.num_cols, 0.0), + x_mark_(lp.num_cols, 0) + { + basis_update.b_solve(lp.rhs, b_bar_); + for (i_t j : nonbasic_list) { + nonbasic_mark_[j] = 1; + } + } + + // Generates the base inequalities: C*x == d that will be turned into cuts + i_t generate_base_inequality(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + const std::vector& var_types, + basis_update_mpf_t& basis_update, + const std::vector& xstar, + const std::vector& basic_list, + const std::vector& nonbasic_list, + i_t i, + sparse_vector_t& inequality, + f_t& inequality_rhs); + + private: + std::vector b_bar_; + std::vector nonbasic_mark_; + std::vector x_workspace_; + std::vector x_mark_; +}; + +template +class mixed_integer_rounding_cut_t { + public: + mixed_integer_rounding_cut_t(i_t num_vars, const simplex_solver_settings_t& settings) + : num_vars_(num_vars), + settings_(settings), + x_workspace_(num_vars, 0.0), + x_mark_(num_vars, 0), + has_lower_(num_vars, 0), + has_upper_(num_vars, 0), + needs_complement_(false) + { + } + + void initialize(const lp_problem_t& lp, + const std::vector& new_slacks, + const std::vector& xstar); + + i_t generate_cut(const sparse_vector_t& a, + f_t beta, + const std::vector& upper_bounds, + const std::vector& lower_bounds, + const std::vector& var_types, + sparse_vector_t& cut, + f_t& cut_rhs); + + void substitute_slacks(const lp_problem_t& lp, + csc_matrix_t& Arow, + sparse_vector_t& cut, + f_t& cut_rhs); + + private: + i_t num_vars_; + const simplex_solver_settings_t& settings_; + std::vector x_workspace_; + std::vector x_mark_; + std::vector has_lower_; + std::vector has_upper_; + std::vector is_slack_; + std::vector slack_rows_; + bool needs_complement_; +}; + +template +i_t add_cuts(const simplex_solver_settings_t& settings, + const csr_matrix_t& cuts, + const std::vector& cut_rhs, + lp_problem_t& lp, + std::vector& new_slacks, + lp_solution_t& solution, + basis_update_mpf_t& basis_update, + std::vector& basic_list, + std::vector& nonbasic_list, + std::vector& vstatus, + std::vector& edge_norms); + +template +void remove_cuts(lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csc_matrix_t& Arow, + std::vector& new_slacks, + i_t original_rows, + std::vector& var_types, + std::vector& vstatus, + std::vector& x, + std::vector& y, + std::vector& z, + std::vector& basic_list, + std::vector& nonbasic_list, + basis_update_mpf_t& basis_update); + +} + diff --git a/cpp/src/dual_simplex/mip_node.hpp b/cpp/src/dual_simplex/mip_node.hpp index 9034bfa22..f18ae0072 100644 --- a/cpp/src/dual_simplex/mip_node.hpp +++ b/cpp/src/dual_simplex/mip_node.hpp @@ -41,6 +41,7 @@ class mip_node_t { node_id(0), branch_var(-1), branch_dir(-1), + integer_infeasible(-1), vstatus(basis) { children[0] = nullptr; @@ -53,6 +54,7 @@ class mip_node_t { i_t branch_variable, i_t branch_direction, f_t branch_var_value, + i_t integer_inf, const std::vector& basis) : status(node_status_t::ACTIVE), lower_bound(parent_node->lower_bound), @@ -62,8 +64,8 @@ class mip_node_t { branch_var(branch_variable), branch_dir(branch_direction), fractional_val(branch_var_value), + integer_infeasible(integer_inf), vstatus(basis) - { branch_var_lower = branch_direction == 0 ? problem.lower[branch_var] : std::ceil(branch_var_value); @@ -217,6 +219,7 @@ class mip_node_t { f_t branch_var_lower; f_t branch_var_upper; f_t fractional_val; + i_t integer_infeasible; mip_node_t* parent; std::unique_ptr children[2]; @@ -272,6 +275,7 @@ class search_tree_t { void branch(mip_node_t* parent_node, const i_t branch_var, const f_t fractional_val, + const i_t integer_infeasible, const std::vector& parent_vstatus, const lp_problem_t& original_lp, logger_t& log) @@ -280,13 +284,13 @@ class search_tree_t { // down child auto down_child = std::make_unique>( - original_lp, parent_node, ++id, branch_var, 0, fractional_val, parent_vstatus); + original_lp, parent_node, ++id, branch_var, 0, fractional_val, integer_infeasible, 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); + original_lp, parent_node, ++id, branch_var, 1, fractional_val, integer_infeasible, parent_vstatus); graphviz_edge(log, parent_node, up_child.get(), branch_var, 1, std::ceil(fractional_val)); diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 39ea9b465..2ff075c15 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -2397,6 +2397,39 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, } timers.pricing_time += timers.stop_timer(); if (leaving_index == -1) { + + +#ifdef CHECK_BASIS_UPDATE + for (i_t k = 0; k < basic_list.size(); k++) { + const i_t jj = basic_list[k]; + sparse_vector_t ei_sparse(m, 1); + ei_sparse.i[0] = k; + ei_sparse.x[0] = 1.0; + sparse_vector_t ubar_sparse(m, 0); + ft.b_transpose_solve(ei_sparse, ubar_sparse); + std::vector ubar_dense(m); + ubar_sparse.to_dense(ubar_dense); + std::vector BTu_dense(m); + b_transpose_multiply(lp, basic_list, ubar_dense, BTu_dense); + for (i_t l = 0; l < m; l++) { + if (l != k) { + settings.log.printf("BTu_dense[%d] = %e i %d\n", l, BTu_dense[l], k); + } else { + settings.log.printf("BTu_dense[%d] = %e != 1.0 i %d\n", l, BTu_dense[l], k); + } + } + for (i_t h = 0; h < m; h++) { + settings.log.printf("i %d ubar_dense[%d] = %.16e\n", k, h, ubar_dense[h]); + } + } + settings.log.printf("ft.num_updates() %d\n", ft.num_updates()); + for (i_t h = 0; h < m; h++) { + settings.log.printf("basic_list[%d] = %d\n", h, basic_list[h]); + } + +#endif + + phase2::prepare_optimality(lp, settings, ft, diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 8e54c40bb..6fffbfe1f 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -69,6 +69,7 @@ struct simplex_solver_settings_t { num_threads(omp_get_max_threads() - 1), num_bfs_threads(std::min(num_threads / 4, 1)), num_diving_threads(std::min(num_threads - num_bfs_threads, 1)), + max_cut_passes(10), random_seed(0), inside_mip(0), solution_callback(nullptr), @@ -134,6 +135,7 @@ struct simplex_solver_settings_t { i_t random_seed; // random seed i_t num_bfs_threads; // number of threads dedicated to the best-first search i_t num_diving_threads; // number of threads dedicated to diving + i_t max_cut_passes; // number of cut passes to make i_t inside_mip; // 0 if outside MIP, 1 if inside MIP at root node, 2 if inside MIP at leaf node std::function&, f_t)> solution_callback; std::function&, f_t)> node_processed_callback; diff --git a/cpp/src/dual_simplex/solve.cpp b/cpp/src/dual_simplex/solve.cpp index 5c5f9e165..951945092 100644 --- a/cpp/src/dual_simplex/solve.cpp +++ b/cpp/src/dual_simplex/solve.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #include #include diff --git a/cpp/src/dual_simplex/sparse_matrix.cpp b/cpp/src/dual_simplex/sparse_matrix.cpp index cdd45f720..6160cf1b4 100644 --- a/cpp/src/dual_simplex/sparse_matrix.cpp +++ b/cpp/src/dual_simplex/sparse_matrix.cpp @@ -357,6 +357,75 @@ i_t csc_matrix_t::remove_row(i_t row) return 0; } +template +i_t csr_matrix_t::append_rows(const csr_matrix_t& C) +{ + const i_t old_m = this->m; + const i_t n = this->n; + const i_t old_nz = this->row_start[old_m]; + const i_t C_row = C.m; + if (C.n > n) { + printf("append_rows error: C.n %d n %d\n", C.n, n); + return -1; + } + const i_t C_nz = C.row_start[C_row]; + const i_t new_nz = old_nz + C_nz; + const i_t new_m = old_m + C_row; + + this->j.resize(new_nz); + this->x.resize(new_nz); + this->row_start.resize(new_m + 1); + + i_t nz = old_nz; + for (i_t i = old_m; i < new_m; i++) { + const i_t k = i - old_m; + const i_t nz_row = C.row_start[k + 1] - C.row_start[k]; + this->row_start[i] = nz; + nz += nz_row; + } + this->row_start[new_m] = nz; + + for (i_t p = old_nz; p < new_nz; p++) { + const i_t q = p - old_nz; + this->j[p] = C.j[q]; + } + + for (i_t p = old_nz; p < new_nz; p++) { + const i_t q = p - old_nz; + this->x[p] = C.x[q]; + } + + this->m = new_m; + this->nz_max = new_nz; + return 0; +} + +template +i_t csr_matrix_t::append_row(const sparse_vector_t& c) +{ + const i_t old_m = this->m; + const i_t old_nz = this->row_start[old_m]; + const i_t c_nz = c.i.size(); + const i_t new_nz = old_nz + c_nz; + const i_t new_m = old_m + 1; + + this->j.resize(new_nz); + this->x.resize(new_nz); + this->row_start.resize(new_m + 1); + this->row_start[new_m] = new_nz; + + i_t nz = old_nz; + for (i_t k = 0; k < c_nz; k++) { + this->j[nz] = c.i[k]; + this->x[nz] = c.x[k]; + nz++; + } + + this->m = new_m; + this->nz_max = new_nz; + return 0; +} + template void csc_matrix_t::print_matrix(FILE* fid) const { @@ -498,6 +567,10 @@ i_t csc_matrix_t::check_matrix() const { std::vector row_marker(this->m, -1); for (i_t j = 0; j < this->n; ++j) { + if (j >= col_start.size()) { + printf("Col start too small size %ld n %d\n", col_start.size(), this->n); + return -1; + } const i_t col_start = this->col_start[j]; const i_t col_end = this->col_start[j + 1]; if (col_start > col_end || col_start > this->col_start[this->n]) { @@ -556,6 +629,7 @@ void csr_matrix_t::check_matrix() const const i_t row_end = this->row_start[i + 1]; for (i_t p = row_start; p < row_end; ++p) { const i_t j = this->j[p]; + if (j < 0 || j >= this->n) { printf("CSR Error: column index %d not in range [0, %d)\n", j, this->n); } if (col_marker[j] == i) { printf("CSR Error: repeated column index %d in row %d\n", j, i); } col_marker[j] = i; } diff --git a/cpp/src/dual_simplex/sparse_matrix.hpp b/cpp/src/dual_simplex/sparse_matrix.hpp index c14e6d0f1..49c5c185a 100644 --- a/cpp/src/dual_simplex/sparse_matrix.hpp +++ b/cpp/src/dual_simplex/sparse_matrix.hpp @@ -136,6 +136,12 @@ class csr_matrix_t { // Create a new matrix with the marked rows removed i_t remove_rows(std::vector& row_marker, csr_matrix_t& Aout) const; + // Append rows from another CSR matrix + i_t append_rows(const csr_matrix_t& C); + + // Append a row from a sparse vector + i_t append_row(const sparse_vector_t& c); + // Ensures no repeated column indices within a row void check_matrix() const; diff --git a/cpp/src/dual_simplex/sparse_vector.cpp b/cpp/src/dual_simplex/sparse_vector.cpp index 2d4745650..3ba981539 100644 --- a/cpp/src/dual_simplex/sparse_vector.cpp +++ b/cpp/src/dual_simplex/sparse_vector.cpp @@ -28,6 +28,21 @@ sparse_vector_t::sparse_vector_t(const csc_matrix_t& A, i_t } } +template +sparse_vector_t::sparse_vector_t(const csr_matrix_t& A, i_t row) +{ + const i_t row_start = A.row_start[row]; + const i_t row_end = A.row_start[row + 1]; + const i_t nz = row_end - row_start; + n = A.n; + i.reserve(nz); + x.reserve(nz); + for (i_t k = row_start; k < row_end; ++k) { + i.push_back(A.j[k]); + x.push_back(A.x[k]); + } +} + template void sparse_vector_t::from_dense(const std::vector& in) { @@ -106,6 +121,17 @@ void sparse_vector_t::inverse_permute_vector(const std::vector& p y.i = i_perm; } +template +f_t sparse_vector_t::dot(const std::vector& x_dense) const +{ + const i_t nz = i.size(); + f_t dot = 0.0; + for (i_t k = 0; k < nz; ++k) { + dot += x[k] * x_dense[i[k]]; + } + return dot; +} + template f_t sparse_vector_t::sparse_dot(const csc_matrix_t& Y, i_t y_col) const { diff --git a/cpp/src/dual_simplex/sparse_vector.hpp b/cpp/src/dual_simplex/sparse_vector.hpp index 7acfdc8b5..3badeed12 100644 --- a/cpp/src/dual_simplex/sparse_vector.hpp +++ b/cpp/src/dual_simplex/sparse_vector.hpp @@ -25,6 +25,8 @@ class sparse_vector_t { sparse_vector_t(const std::vector& in) { from_dense(in); } // Construct a sparse vector from a column of a CSC matrix sparse_vector_t(const csc_matrix_t& A, i_t col); + // Construct a sparse vector from a row of a CSR matrix + sparse_vector_t(const csr_matrix_t& A, i_t row); // gather a dense vector into a sparse vector void from_dense(const std::vector& in); // convert a sparse vector into a CSC matrix with a single column @@ -38,6 +40,8 @@ class sparse_vector_t { void inverse_permute_vector(const std::vector& p); // inverse permute a sparse vector into another sparse vector void inverse_permute_vector(const std::vector& p, sparse_vector_t& y) const; + // compute the dot product of a sparse vector with a dense vector + f_t dot(const std::vector& x) const; // compute the dot product of a sparse vector with a column of a CSC matrix f_t sparse_dot(const csc_matrix_t& Y, i_t y_col) const; // ensure the coefficients in the sparse vectory are sorted in terms of increasing index diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index b5da4f095..0d3874321 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -86,7 +86,8 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_FOLDING, &pdlp_settings.folding, -1, 1, -1}, {CUOPT_DUALIZE, &pdlp_settings.dualize, -1, 1, -1}, {CUOPT_ORDERING, &pdlp_settings.ordering, -1, 1, -1}, - {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1} + {CUOPT_BARRIER_DUAL_INITIAL_POINT, &pdlp_settings.barrier_dual_initial_point, -1, 1, -1}, + {CUOPT_MIP_CUT_PASSES, &mip_settings.max_cut_passes, -1, std::numeric_limits::max(), 10} }; // Bool parameters diff --git a/cpp/src/mip/diversity/diversity_manager.cu b/cpp/src/mip/diversity/diversity_manager.cu index 78f2b9fa0..cbc59b2ab 100644 --- a/cpp/src/mip/diversity/diversity_manager.cu +++ b/cpp/src/mip/diversity/diversity_manager.cu @@ -409,7 +409,7 @@ solution_t diversity_manager_t::run_solver() run_fj_alone(sol); return sol; } - rins.enable(); + //rins.enable(); generate_solution(timer.remaining_time(), false); if (timer.check_time_limit()) { diff --git a/cpp/src/mip/diversity/lns/rins.cu b/cpp/src/mip/diversity/lns/rins.cu index 1efc971b2..23522f07d 100644 --- a/cpp/src/mip/diversity/lns/rins.cu +++ b/cpp/src/mip/diversity/lns/rins.cu @@ -245,6 +245,7 @@ void rins_t::run_rins() branch_and_bound_settings.num_bfs_threads = 1; branch_and_bound_settings.num_diving_threads = 1; branch_and_bound_settings.log.log_prefix = "[RINS] "; + branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.solution_callback = [this, &rins_solution_queue]( std::vector& solution, f_t objective) { rins_solution_queue.push_back(solution); diff --git a/cpp/src/mip/solver.cu b/cpp/src/mip/solver.cu index 28659cccd..cc9f9f6c5 100644 --- a/cpp/src/mip/solver.cu +++ b/cpp/src/mip/solver.cu @@ -109,7 +109,8 @@ solution_t mip_solver_t::run_solver() diversity_manager_t dm(context); dm.timer = timer_; - bool presolve_success = dm.run_presolve(timer_.remaining_time()); + //bool presolve_success = dm.run_presolve(timer_.remaining_time()); + bool presolve_success = true; if (!presolve_success) { CUOPT_LOG_INFO("Problem proven infeasible in presolve"); solution_t sol(*context.problem_ptr); @@ -117,7 +118,7 @@ solution_t mip_solver_t::run_solver() context.problem_ptr->post_process_solution(sol); return sol; } - if (context.problem_ptr->empty) { + if (0 && context.problem_ptr->empty) { CUOPT_LOG_INFO("Problem full reduced in presolve"); solution_t sol(*context.problem_ptr); sol.set_problem_fully_reduced(); @@ -126,7 +127,7 @@ solution_t mip_solver_t::run_solver() } // if the problem was reduced to a LP: run concurrent LP - if (context.problem_ptr->n_integer_vars == 0) { + if (0 && context.problem_ptr->n_integer_vars == 0) { CUOPT_LOG_INFO("Problem reduced to a LP, running concurrent LP"); pdlp_solver_settings_t settings{}; settings.time_limit = timer_.remaining_time(); @@ -167,6 +168,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.absolute_mip_gap_tol = context.settings.tolerances.absolute_mip_gap; branch_and_bound_settings.relative_mip_gap_tol = context.settings.tolerances.relative_mip_gap; branch_and_bound_settings.integer_tol = context.settings.tolerances.integrality_tolerance; + branch_and_bound_settings.max_cut_passes = context.settings.max_cut_passes; if (context.settings.num_cpu_threads < 0) { branch_and_bound_settings.num_threads = omp_get_max_threads() - 1; @@ -224,6 +226,9 @@ solution_t mip_solver_t::run_solver() std::ref(branch_and_bound_solution)); } + //auto bb_status = branch_and_bound_status_future.get(); + //CUOPT_LOG_INFO("BB status: %d", bb_status); + // Start the primal heuristics auto sol = dm.run_solver(); if (!context.settings.heuristics_only) { diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 7aed72fe0..66a2347d1 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/cpp/tests/dual_simplex/unit_tests/solve.cpp @@ -326,4 +326,155 @@ TEST(dual_simplex, dual_variable_greater_than) EXPECT_NEAR(solution.z[1], 0.0, 1e-6); } +#if 0 +TEST(dual_simplex, simple_cuts) +{ + // minimize x + y + 2 z + // subject to x + y + z == 1 + // x, y, z >= 0 + + raft::handle_t handle{}; + cuopt::linear_programming::dual_simplex::user_problem_t user_problem(&handle); + constexpr int m = 1; + constexpr int n = 3; + constexpr int nz = 3; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective.resize(n); + user_problem.objective[0] = 1.0; + user_problem.objective[1] = 1.0; + user_problem.objective[2] = 2.0; + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start.resize(n + 1); + user_problem.A.col_start[0] = 0; + user_problem.A.col_start[1] = 1; + user_problem.A.col_start[2] = 2; + user_problem.A.col_start[3] = 3; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = 1.0; + user_problem.A.i[2] = 0; + user_problem.A.x[2] = 1.0; + user_problem.lower.resize(n, 0.0); + user_problem.upper.resize(n, dual_simplex::inf); + user_problem.num_range_rows = 0; + user_problem.problem_name = "simple_cuts"; + user_problem.obj_scale = 1.0; + user_problem.obj_constant = 0.0; + user_problem.rhs.resize(m, 1.0); + user_problem.row_sense.resize(m, 'E'); + user_problem.var_types.resize( + n, cuopt::linear_programming::dual_simplex::variable_type_t::CONTINUOUS); + + cuopt::init_logger_t logger("", true); + + cuopt::linear_programming::dual_simplex::lp_problem_t lp( + user_problem.handle_ptr, 1, 1, 1); + cuopt::linear_programming::dual_simplex::simplex_solver_settings_t settings; + settings.barrier = false; + settings.barrier_presolve = false; + settings.log.log = true; + settings.log.log_to_console = true; + settings.log.printf("Test print\n"); + std::vector new_slacks; + cuopt::linear_programming::dual_simplex::dualize_info_t dualize_info; + cuopt::linear_programming::dual_simplex::convert_user_problem( + user_problem, settings, lp, new_slacks, dualize_info); + cuopt::linear_programming::dual_simplex::lp_solution_t solution(lp.num_rows, + lp.num_cols); + std::vector vstatus; + std::vector edge_norms; + std::vector basic_list(lp.num_rows); + std::vector nonbasic_list; + cuopt::linear_programming::dual_simplex::basis_update_mpf_t basis_update( + lp.num_cols, settings.refactor_frequency); + double start_time = dual_simplex::tic(); + printf("Calling solve linear program with advanced basis\n"); + EXPECT_EQ((cuopt::linear_programming::dual_simplex::solve_linear_program_with_advanced_basis( + lp, + start_time, + settings, + solution, + basis_update, + basic_list, + nonbasic_list, + vstatus, + edge_norms)), + cuopt::linear_programming::dual_simplex::lp_status_t::OPTIMAL); + printf("Solution objective: %e\n", solution.objective); + printf("Solution x: %e %e %e\n", solution.x[0], solution.x[1], solution.x[2]); + printf("Solution y: %e\n", solution.y[0]); + printf("Solution z: %e %e %e\n", solution.z[0], solution.z[1], solution.z[2]); + EXPECT_NEAR(solution.objective, 1.0, 1e-6); + EXPECT_NEAR(solution.x[0], 1.0, 1e-6); + + // Add a cut z >= 1/3. Needs to be in the form C*x <= d + csr_matrix_t cuts(1, n, 1); + cuts.row_start[0] = 0; + cuts.j[0] = 2; + cuts.x[0] = -1.0; + cuts.row_start[1] = 1; + printf("cuts m %d n %d\n", cuts.m, cuts.n); + std::vector cut_rhs(1); + cut_rhs[0] = -1.0 / 3.0; + + std::vector var_types; + EXPECT_EQ(cuopt::linear_programming::dual_simplex::solve_linear_program_with_cuts(start_time, + settings, + cuts, + cut_rhs, + lp, + solution, + basis_update, + basic_list, + nonbasic_list, + vstatus, + edge_norms, + var_types), + cuopt::linear_programming::dual_simplex::lp_status_t::OPTIMAL); + printf("Solution objective: %e\n", solution.objective); + printf("Solution x: %e %e %e\n", solution.x[0], solution.x[1], solution.x[2]); + EXPECT_NEAR(solution.objective, 4.0 / 3.0, 1e-6); + + cuts.row_start.resize(3); + cuts.j.resize(2); + cuts.x.resize(2); + // Add cut y >= 1/3 + cuts.j[0] = 1; + cuts.row_start[2] = 2; + // Add cut x <= 0.0 + cuts.j[1] = 0; + cuts.x[1] = 1.0; + cuts.m = 2; + cut_rhs.resize(2); + cut_rhs[1] = 0.0; + + EXPECT_EQ(cuopt::linear_programming::dual_simplex::solve_linear_program_with_cuts(start_time, + settings, + cuts, + cut_rhs, + lp, + solution, + basis_update, + basic_list, + nonbasic_list, + vstatus, + edge_norms, + var_types), + cuopt::linear_programming::dual_simplex::lp_status_t::OPTIMAL); + printf("Solution objective: %e\n", solution.objective); + printf("Solution x: %e %e %e\n", solution.x[0], solution.x[1], solution.x[2]); + EXPECT_NEAR(solution.objective, 4.0 / 3.0, 1e-6); + EXPECT_NEAR(solution.x[0], 0.0, 1e-6); + EXPECT_NEAR(solution.x[1], 2.0 / 3.0, 1e-6); + EXPECT_NEAR(solution.x[2], 1.0 / 3.0, 1e-6); + +} +#endif + } // namespace cuopt::linear_programming::dual_simplex::test