diff --git a/benchmarks/linear_programming/cuopt/run_mip.cpp b/benchmarks/linear_programming/cuopt/run_mip.cpp index 9b79dff8af..e01e533a65 100644 --- a/benchmarks/linear_programming/cuopt/run_mip.cpp +++ b/benchmarks/linear_programming/cuopt/run_mip.cpp @@ -183,6 +183,10 @@ int run_single_file(std::string file_path, CUOPT_LOG_ERROR("Parsing MPS failed exiting!"); return -1; } + // Use the benchmark filename for downstream instance-level reporting. + // This keeps per-instance metrics aligned with the run list even if the MPS NAME card differs. + mps_data_model.set_problem_name(base_filename); + if (initial_solution_dir.has_value()) { auto initial_solutions = read_solution_from_dir( initial_solution_dir.value(), base_filename, mps_data_model.get_variable_names()); @@ -209,6 +213,7 @@ int run_single_file(std::string file_path, settings.tolerances.absolute_tolerance = 1e-6; settings.presolver = cuopt::linear_programming::presolver_t::Default; settings.reliability_branching = reliability_branching; + settings.clique_cuts = -1; settings.seed = 42; cuopt::linear_programming::benchmark_info_t benchmark_info; settings.benchmark_info_ptr = &benchmark_info; @@ -413,7 +418,16 @@ int main(int argc, char* argv[]) int reliability_branching = program.get("--reliability-branching"); bool deterministic = program.get("--determinism"); - if (num_cpu_threads < 0) { num_cpu_threads = omp_get_max_threads() / n_gpus; } + if (num_cpu_threads < 0) { + num_cpu_threads = omp_get_max_threads() / n_gpus; + // std::ifstream smt_file("/sys/devices/system/cpu/smt/active"); + // if (smt_file.is_open()) { + // int smt_active = 0; + // smt_file >> smt_active; + // if (smt_active) { num_cpu_threads /= 2; } + // } + num_cpu_threads = std::max(num_cpu_threads, 1); + } if (program.is_used("--out-dir")) { out_dir = program.get("--out-dir"); diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index d9dfbce16d..86becfe06d 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -64,6 +64,7 @@ #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" #define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS "mip_mixed_integer_gomory_cuts" #define CUOPT_MIP_KNAPSACK_CUTS "mip_knapsack_cuts" +#define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" #define CUOPT_MIP_CUT_CHANGE_THRESHOLD "mip_cut_change_threshold" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 6d32cd5ed9..a9e404d14e 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -93,6 +93,7 @@ class mip_solver_settings_t { i_t mir_cuts = -1; i_t mixed_integer_gomory_cuts = -1; i_t knapsack_cuts = -1; + i_t clique_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; f_t cut_change_threshold = 1e-3; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 6ce9a4f4d0..41d23bc0ff 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -241,9 +242,11 @@ template branch_and_bound_t::branch_and_bound_t( const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, - f_t start_time) + f_t start_time, + std::shared_ptr> clique_table) : original_problem_(user_problem), settings_(solver_settings), + clique_table_(std::move(clique_table)), original_lp_(user_problem.handle_ptr, 1, 1, 1), Arow_(1, 1, 0), incumbent_(1), @@ -1967,6 +1970,31 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols); + if (settings_.clique_cuts != 0 && clique_table_ == nullptr) { + signal_extend_cliques_.store(false, std::memory_order_release); + typename ::cuopt::linear_programming::mip_solver_settings_t::tolerances_t + tolerances_for_clique{}; + tolerances_for_clique.presolve_absolute_tolerance = settings_.primal_tol; + tolerances_for_clique.absolute_tolerance = settings_.primal_tol; + tolerances_for_clique.relative_tolerance = settings_.zero_tol; + tolerances_for_clique.integrality_tolerance = settings_.integer_tol; + tolerances_for_clique.absolute_mip_gap = settings_.absolute_mip_gap_tol; + tolerances_for_clique.relative_mip_gap = settings_.relative_mip_gap_tol; + auto* signal_ptr = &signal_extend_cliques_; + clique_table_future_ = + std::async(std::launch::async, + [this, + tolerances_for_clique, + signal_ptr]() -> std::shared_ptr> { + user_problem_t problem_copy = original_problem_; + cuopt::timer_t timer(std::numeric_limits::infinity()); + std::shared_ptr> table; + detail::find_initial_cliques( + problem_copy, tolerances_for_clique, &table, timer, false, signal_ptr); + return table; + }); + } + i_t original_rows = original_lp_.num_rows; simplex_solver_settings_t lp_settings = settings_; lp_settings.inside_mip = 1; @@ -2002,14 +2030,16 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.total_lp_iters = root_relax_soln_.iterations; exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time); + auto finish_clique_thread = [this]() { + if (clique_table_future_.valid()) { + signal_extend_cliques_.store(true, std::memory_order_release); + clique_table_ = clique_table_future_.get(); + } + }; + if (root_status == lp_status_t::INFEASIBLE) { settings_.log.printf("MIP Infeasible\n"); - // FIXME: rarely dual simplex detects infeasible whereas it is feasible. - // to add a small safety net, check if there is a primal solution already. - // Uncomment this if the issue with cost266-UUE is resolved - // if (settings.heuristic_preemption_callback != nullptr) { - // settings.heuristic_preemption_callback(); - // } + finish_clique_thread(); return mip_status_t::INFEASIBLE; } if (root_status == lp_status_t::UNBOUNDED) { @@ -2017,23 +2047,27 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (settings_.heuristic_preemption_callback != nullptr) { settings_.heuristic_preemption_callback(); } + finish_clique_thread(); return mip_status_t::UNBOUNDED; } if (root_status == lp_status_t::TIME_LIMIT) { solver_status_ = mip_status_t::TIME_LIMIT; set_final_solution(solution, -inf); + finish_clique_thread(); return solver_status_; } if (root_status == lp_status_t::WORK_LIMIT) { solver_status_ = mip_status_t::WORK_LIMIT; set_final_solution(solution, -inf); + finish_clique_thread(); return solver_status_; } if (root_status == lp_status_t::NUMERICAL_ISSUES) { solver_status_ = mip_status_t::NUMERICAL; set_final_solution(solution, -inf); + finish_clique_thread(); return solver_status_; } @@ -2064,6 +2098,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (num_fractional == 0) { set_solution_at_root(solution, cut_info); + finish_clique_thread(); return mip_status_t::OPTIMAL; } @@ -2078,8 +2113,16 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } cut_pool_t cut_pool(original_lp_.num_cols, settings_); - cut_generation_t cut_generation( - cut_pool, original_lp_, settings_, Arow_, new_slacks_, var_types_); + cut_generation_t cut_generation(cut_pool, + original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + original_problem_, + clique_table_, + &clique_table_future_, + &signal_extend_cliques_); std::vector saved_solution; #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION @@ -2090,7 +2133,8 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t last_objective = root_objective_; f_t root_relax_objective = root_objective_; - i_t cut_pool_size = 0; + f_t cut_generation_start_time = tic(); + i_t cut_pool_size = 0; for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) { if (num_fractional == 0) { set_solution_at_root(solution, cut_info); @@ -2109,16 +2153,25 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut #endif // Generate cuts and add them to the cut pool - f_t cut_start_time = tic(); - cut_generation.generate_cuts(original_lp_, - settings_, - Arow_, - new_slacks_, - var_types_, - basis_update, - root_relax_soln_.x, - basic_list, - nonbasic_list); + f_t cut_start_time = tic(); + bool problem_feasible = cut_generation.generate_cuts(original_lp_, + settings_, + Arow_, + new_slacks_, + var_types_, + basis_update, + root_relax_soln_.x, + root_relax_soln_.z, + basic_list, + nonbasic_list, + exploration_stats_.start_time); + if (!problem_feasible) { + if (settings_.heuristic_preemption_callback != nullptr) { + settings_.heuristic_preemption_callback(); + } + finish_clique_thread(); + return mip_status_t::INFEASIBLE; + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); @@ -2339,8 +2392,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } print_cut_info(settings_, cut_info); - + f_t cut_generation_time = toc(cut_generation_start_time); if (cut_info.has_cuts()) { + settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time); settings_.log.printf("Cut pool size : %d\n", cut_pool_size); settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n", original_lp_.num_rows, diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index a13d5cedcf..f44142f37f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -32,9 +32,17 @@ #include +#include #include +#include +#include #include +namespace cuopt::linear_programming::detail { +template +struct clique_table_t; +} + namespace cuopt::linear_programming::dual_simplex { enum class mip_status_t { @@ -68,7 +76,8 @@ class branch_and_bound_t { public: branch_and_bound_t(const user_problem_t& user_problem, const simplex_solver_settings_t& solver_settings, - f_t start_time); + f_t start_time, + std::shared_ptr> clique_table = nullptr); // Set an initial guess based on the user_problem. This should be called before solve. void set_initial_guess(const std::vector& user_guess) { guess_ = user_guess; } @@ -106,8 +115,6 @@ class branch_and_bound_t { void set_concurrent_lp_root_solve(bool enable) { enable_concurrent_lp_root_solve_ = enable; } - bool stop_for_time_limit(mip_solution_t& solution); - // Repair a low-quality solution from the heuristics. bool repair_solution(const std::vector& leaf_edge_norms, const std::vector& potential_solution, @@ -141,6 +148,9 @@ class branch_and_bound_t { private: const user_problem_t& original_problem_; const simplex_solver_settings_t settings_; + std::shared_ptr> clique_table_; + std::future>> clique_table_future_; + std::atomic signal_extend_cliques_{false}; work_limit_context_t work_unit_context_{"B&B"}; diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index cc2555611a..f4b3106112 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -9,11 +9,486 @@ #include #include +#include +#include + +#include +#include +#include +#include #include namespace cuopt::linear_programming::dual_simplex { +namespace { + +#define DEBUG_CLIQUE_CUTS 0 +#define CHECK_WORKSPACE 0 + +enum class clique_cut_build_status_t : int8_t { NO_CUT = 0, CUT_ADDED = 1, INFEASIBLE = 2 }; + +#if DEBUG_CLIQUE_CUTS +#define CLIQUE_CUTS_DEBUG(...) \ + do { \ + std::fprintf(stderr, "[DEBUG_CLIQUE_CUTS] "); \ + std::fprintf(stderr, __VA_ARGS__); \ + std::fprintf(stderr, "\n"); \ + } while (0) +#else +#define CLIQUE_CUTS_DEBUG(...) \ + do { \ + } while (0) +#endif + +template +clique_cut_build_status_t build_clique_cut(const std::vector& clique_vertices, + i_t num_vars, + const std::vector& var_types, + const std::vector& lower_bounds, + const std::vector& upper_bounds, + const std::vector& xstar, + f_t bound_tol, + f_t min_violation, + sparse_vector_t& cut, + f_t& cut_rhs, + f_t* work_estimate, + f_t max_work_estimate) +{ + if (clique_vertices.size() < 2) { return clique_cut_build_status_t::NO_CUT; } + const f_t clique_size = static_cast(clique_vertices.size()); + CLIQUE_CUTS_DEBUG("build_clique_cut start clique_size=%lld", + static_cast(clique_vertices.size())); + const f_t sort_work = clique_size > 0.0 ? 2.0 * clique_size * std::log2(clique_size + 1.0) : 0.0; + const f_t dot_work = 2.0 * clique_size; + const f_t estimated_work = 9.0 * clique_size + sort_work + dot_work; + if (add_work_estimate(estimated_work, work_estimate, max_work_estimate)) { + CLIQUE_CUTS_DEBUG("build_clique_cut skip work_limit clique_size=%lld work=%g limit=%g", + static_cast(clique_vertices.size()), + work_estimate == nullptr ? -1.0 : static_cast(*work_estimate), + static_cast(max_work_estimate)); + return clique_cut_build_status_t::NO_CUT; + } + + cuopt_assert(num_vars > 0, "Clique cut num_vars must be positive"); + cuopt_assert(static_cast(num_vars) <= lower_bounds.size(), + "Clique cut lower bounds size mismatch"); + cuopt_assert(static_cast(num_vars) <= xstar.size(), "Clique cut xstar size mismatch"); + + cut.i.clear(); + cut.x.clear(); + i_t num_complements = 0; + std::unordered_set seen_original; + std::unordered_set seen_complement; + seen_original.reserve(clique_vertices.size()); + seen_complement.reserve(clique_vertices.size()); + for (const auto vertex_idx : clique_vertices) { + cuopt_assert(vertex_idx >= 0 && vertex_idx < 2 * num_vars, "Clique vertex out of range"); + const i_t var_idx = vertex_idx % num_vars; + const bool complement = vertex_idx >= num_vars; + const f_t lower_bound = lower_bounds[var_idx]; + const f_t upper_bound = upper_bounds[var_idx]; + + cuopt_assert(var_types[var_idx] != variable_type_t::CONTINUOUS, + "Clique contains continuous variable"); + cuopt_assert(lower_bound >= -bound_tol, "Clique variable lower bound below zero"); + cuopt_assert(upper_bound <= 1 + bound_tol, "Clique variable upper bound above one"); + + // we store the cut in the form of >= 1, for easy violation check with dot product + // that's why compelements have 1 as coeff and normal vars have -1 + if (complement) { + if (seen_original.count(var_idx) > 0) { + // FIXME: this is temporary, fix all the vars of all other vars in the clique + return clique_cut_build_status_t::NO_CUT; + CLIQUE_CUTS_DEBUG("build_clique_cut infeasible var=%lld appears as variable and complement", + static_cast(var_idx)); + return clique_cut_build_status_t::INFEASIBLE; + } + cuopt_assert(seen_complement.count(var_idx) == 0, "Duplicate complement in clique"); + seen_complement.insert(var_idx); + num_complements++; + cut.i.push_back(var_idx); + cut.x.push_back(1.0); + } else { + if (seen_complement.count(var_idx) > 0) { + // FIXME: this is temporary, fix all the vars of all other vars in the clique + return clique_cut_build_status_t::NO_CUT; + CLIQUE_CUTS_DEBUG("build_clique_cut infeasible var=%lld appears as variable and complement", + static_cast(var_idx)); + return clique_cut_build_status_t::INFEASIBLE; + } + cuopt_assert(seen_original.count(var_idx) == 0, "Duplicate variable in clique"); + seen_original.insert(var_idx); + cut.i.push_back(var_idx); + cut.x.push_back(-1.0); + } + } + + if (cut.i.empty()) { + CLIQUE_CUTS_DEBUG("build_clique_cut no_cut empty support"); + return clique_cut_build_status_t::NO_CUT; + } + + cut_rhs = static_cast(num_complements - 1); + cut.sort(); + + const f_t dot = cut.dot(xstar); + const f_t violation = cut_rhs - dot; + if (violation > min_violation) { + CLIQUE_CUTS_DEBUG( + "build_clique_cut accepted nz=%lld rhs=%g dot=%g violation=%g threshold=%g complements=%lld", + static_cast(cut.i.size()), + static_cast(cut_rhs), + static_cast(dot), + static_cast(violation), + static_cast(min_violation), + static_cast(num_complements)); + return clique_cut_build_status_t::CUT_ADDED; + } + CLIQUE_CUTS_DEBUG( + "build_clique_cut rejected nz=%lld rhs=%g dot=%g violation=%g threshold=%g complements=%lld", + static_cast(cut.i.size()), + static_cast(cut_rhs), + static_cast(dot), + static_cast(violation), + static_cast(min_violation), + static_cast(num_complements)); + return clique_cut_build_status_t::NO_CUT; +} + +template +struct bk_bitset_context_t { + const std::vector>& adj; + const std::vector& weights; + f_t min_weight; + i_t max_calls; + f_t start_time; + f_t time_limit; + size_t words; + f_t* work_estimate; + f_t max_work_estimate; + i_t num_calls{0}; + bool work_limit_reached{false}; + bool call_limit_reached{false}; + std::vector> cliques; + + bool add_work(f_t accesses) + { + return add_work_estimate(accesses, work_estimate, max_work_estimate, &work_limit_reached); + } + + bool over_work_limit() const + { + if (work_limit_reached) { return true; } + if (work_estimate == nullptr) { return false; } + return *work_estimate > max_work_estimate; + } + + bool over_call_limit() const { return call_limit_reached || num_calls >= max_calls; } +}; + +inline size_t bitset_words(size_t n) { return (n + 63) / 64; } + +inline bool bitset_any(const std::vector& bs) +{ + for (auto word : bs) { + if (word != 0) { return true; } + } + return false; +} + +inline void bitset_set(std::vector& bs, size_t idx) +{ + bs[idx >> 6] |= (uint64_t(1) << (idx & 63)); +} + +inline void bitset_clear(std::vector& bs, size_t idx) +{ + bs[idx >> 6] &= ~(uint64_t(1) << (idx & 63)); +} + +template +f_t sum_weights_bitset(const std::vector& bs, const std::vector& weights) +{ + f_t sum = 0.0; + for (size_t w = 0; w < bs.size(); ++w) { + uint64_t word = bs[w]; + while (word) { + const int bit = __builtin_ctzll(word); + const size_t idx = w * 64 + static_cast(bit); + sum += weights[idx]; + word &= (word - 1); + } + } + return sum; +} + +template +void bron_kerbosch(bk_bitset_context_t& ctx, + std::vector& R, // current clique + std::vector& P, // potential candidates + std::vector& X, // already in the clique + f_t weight_R) +{ + if (ctx.over_work_limit() || ctx.over_call_limit()) { return; } + if (toc(ctx.start_time) >= ctx.time_limit) { return; } + ctx.num_calls++; + // stop the recursion, for perf reasons + if (ctx.num_calls > ctx.max_calls) { + ctx.call_limit_reached = true; + return; + } + if (ctx.add_work(static_cast(4 * ctx.words))) { return; } + + // if P and X are empty, we are at maximal clique + if (!bitset_any(P) && !bitset_any(X)) { + // if the weight is enough, add and exit + if (weight_R >= ctx.min_weight) { + ctx.add_work(static_cast(R.size())); + ctx.cliques.push_back(R); + } + return; + } + + const f_t sumP = sum_weights_bitset(P, ctx.weights); + // check if all P is added to clique, would we exceed the weight? + if (weight_R + sumP < ctx.min_weight) { return; } + + i_t pivot = -1; + i_t max_deg = -1; + i_t pivot_vertices_examined = 0; + // pivoting rule according to the highest degree vertex + // TODO try other pivoting strategies, we can also implement some online learning like MAB + for (size_t w = 0; w < ctx.words; ++w) { + // union of P and X + uint64_t word = P[w] | X[w]; + while (word) { + pivot_vertices_examined++; + // least significant set bit idnex + const int bit = __builtin_ctzll(word); + // overall vertex index + const i_t v = static_cast(w * 64 + static_cast(bit)); + // clear the least significant set bit (v) + word &= (word - 1); + i_t count = 0; + // count the number of neighbors of v in P + for (size_t k = 0; k < ctx.words; ++k) { + count += __builtin_popcountll(P[k] & ctx.adj[v][k]); + } + // chose the highest degree v as the pivot + // we choose the highest degree as the pivot to reduce the recursion size + // later in this function we recurse on the candidate P / N(v) + // so it is good to maximize P n N(v) + if (count > max_deg) { + max_deg = count; + pivot = v; + } + } + } + ctx.add_work(static_cast(2 * ctx.words) + + static_cast(pivot_vertices_examined) * static_cast(2 * ctx.words)); + + std::vector candidates; + candidates.reserve(ctx.weights.size()); + cuopt_assert(pivot >= 0, "Pivot must be valid when P or X is non-empty"); + for (size_t w = 0; w < ctx.words; ++w) { + // P / N(pivot) + uint64_t word = P[w] & ~ctx.adj[pivot][w]; + while (word) { + const int bit = __builtin_ctzll(word); + const i_t v = static_cast(w * 64 + static_cast(bit)); + word &= (word - 1); + candidates.push_back(v); + } + } + const i_t num_candidates = static_cast(candidates.size()); + ctx.add_work(static_cast(2 * ctx.words + num_candidates)); + ctx.add_work(static_cast(num_candidates) * static_cast(7 * ctx.words + 6)); + // note that candidates will include pivot if it is in P + for (auto v : candidates) { + if (ctx.over_call_limit()) { + ctx.call_limit_reached = true; + return; + } + if (toc(ctx.start_time) >= ctx.time_limit) { return; } + + R.push_back(v); + std::vector P_next(ctx.words, 0); + std::vector X_next(ctx.words, 0); + for (size_t k = 0; k < ctx.words; ++k) { + P_next[k] = P[k] & ctx.adj[v][k]; + X_next[k] = X[k] & ctx.adj[v][k]; + } + + bron_kerbosch(ctx, R, P_next, X_next, weight_R + ctx.weights[v]); + if (ctx.over_work_limit()) { return; } + if (ctx.over_call_limit()) { + ctx.call_limit_reached = true; + return; + } + R.pop_back(); + bitset_clear(P, static_cast(v)); + bitset_set(X, static_cast(v)); + } +} + +template +void extend_clique_vertices(std::vector& clique_vertices, + detail::clique_table_t& graph, + const std::vector& xstar, + const std::vector& reduced_costs, + i_t num_vars, + f_t integer_tol, + f_t start_time, + f_t time_limit, + f_t* work_estimate, + f_t max_work_estimate) +{ + if (toc(start_time) >= time_limit) { return; } + if (clique_vertices.empty()) { return; } +#if DEBUG_CLIQUE_CUTS + const size_t initial_clique_vertices = clique_vertices.size(); +#endif + CLIQUE_CUTS_DEBUG("extend_clique_vertices start size=%lld", + static_cast(clique_vertices.size())); + const f_t initial_clique_size = static_cast(clique_vertices.size()); + + i_t smallest_degree = std::numeric_limits::max(); + i_t smallest_degree_var = -1; + for (auto v : clique_vertices) { + if (toc(start_time) >= time_limit) { return; } + i_t degree = graph.get_degree_of_var(v); + if (degree < smallest_degree) { + smallest_degree = degree; + smallest_degree_var = v; + } + } + + auto adj_set = graph.get_adj_set_of_var(smallest_degree_var); + std::unordered_set clique_members(clique_vertices.begin(), clique_vertices.end()); + std::vector candidates; + candidates.reserve(adj_set.size()); + // the candidate list if only the integer valued vertices + for (const auto& candidate : adj_set) { + if (toc(start_time) >= time_limit) { return; } + if (clique_members.count(candidate) != 0) { continue; } + i_t var_idx = candidate % num_vars; + f_t value = candidate >= num_vars ? (1.0 - xstar[var_idx]) : xstar[var_idx]; + if (std::abs(value - std::round(value)) <= integer_tol) { candidates.push_back(candidate); } + } + CLIQUE_CUTS_DEBUG( + "extend_clique_vertices anchor=%lld degree=%lld adj_size=%lld integer_candidates=%lld", + static_cast(smallest_degree_var), + static_cast(smallest_degree), + static_cast(adj_set.size()), + static_cast(candidates.size())); + const f_t candidate_size = static_cast(candidates.size()); + const f_t sort_work = + candidate_size > 0.0 ? 2.0 * candidate_size * std::log2(candidate_size + 1.0) : 0.0; + const f_t adj_set_build_cost = 2.0 * static_cast(adj_set.size()); + const f_t adj_check_cost = 5.0; + const f_t estimated_preloop_work = 2.0 * initial_clique_size + adj_set_build_cost + + 3.0 * static_cast(adj_set.size()) + sort_work + + 2.0 * candidate_size; + if (add_work_estimate(estimated_preloop_work, work_estimate, max_work_estimate)) { + CLIQUE_CUTS_DEBUG("extend_clique_vertices skip work_limit work=%g limit=%g", + work_estimate == nullptr ? -1.0 : static_cast(*work_estimate), + static_cast(max_work_estimate)); + return; + } + + // sort the candidates by reduced cost. + // smaller reduce cost disturbs dual simplex less + // less refactors and less iterations after resolve. + // it also increases the cut's effectiveness by keeping xstar not disturbed much + // if it is disturbed too much, the cut might become non-binding + auto reduced_cost = [&](i_t vertex_idx) -> f_t { + i_t var_idx = vertex_idx % num_vars; + cuopt_assert(var_idx >= 0 && var_idx < static_cast(reduced_costs.size()), + "Variable index out of range"); + f_t rc = reduced_costs[var_idx]; + if (!std::isfinite(rc)) { rc = 0.0; } + return vertex_idx >= num_vars ? -rc : rc; + }; + + std::sort(candidates.begin(), candidates.end(), [&](i_t a, i_t b) { + return reduced_cost(a) < reduced_cost(b); + }); + + for (const auto candidate : candidates) { + bool add = true; + i_t checks = 0; + for (const auto v : clique_vertices) { + checks++; + if (!graph.check_adjacency(candidate, v)) { + add = false; + break; + } + } + if (add_work_estimate( + adj_check_cost * static_cast(checks), work_estimate, max_work_estimate)) { + break; + } + if (add) { + clique_vertices.push_back(candidate); + clique_members.insert(candidate); + } + } + CLIQUE_CUTS_DEBUG("extend_clique_vertices done start=%lld final=%lld added=%lld", + static_cast(initial_clique_vertices), + static_cast(clique_vertices.size()), + static_cast(clique_vertices.size() - initial_clique_vertices)); +} + +} // namespace + +// This function is only used in tests +std::vector> find_maximal_cliques_for_test( + const std::vector>& adjacency_list, + const std::vector& weights, + double min_weight, + int max_calls, + double time_limit) +{ + const size_t n_vertices = adjacency_list.size(); + if (n_vertices == 0) { return {}; } + cuopt_assert(weights.size() == n_vertices, "Weights size mismatch in clique test helper"); + cuopt_assert(max_calls > 0, "max_calls must be positive in clique test helper"); + + const size_t words = bitset_words(n_vertices); + std::vector> adj_bitset(n_vertices, std::vector(words, 0)); + for (size_t v = 0; v < n_vertices; ++v) { + for (const auto& nbr : adjacency_list[v]) { + cuopt_assert(nbr >= 0 && static_cast(nbr) < n_vertices, + "Neighbor index out of range in clique test helper"); + bitset_set(adj_bitset[v], static_cast(nbr)); + } + } + + double work_estimate = 0.0; + const double max_work_estimate = std::numeric_limits::infinity(); + const double start_time = tic(); + + bk_bitset_context_t ctx{adj_bitset, + weights, + min_weight, + max_calls, + start_time, + time_limit, + words, + &work_estimate, + max_work_estimate}; + + std::vector R; + std::vector P(words, 0); + std::vector X(words, 0); + for (size_t idx = 0; idx < n_vertices; ++idx) { + bitset_set(P, idx); + } + bron_kerbosch(ctx, R, P, X, 0.0); + return ctx.cliques; +} + template void cut_pool_t::add_cut(cut_type_t cut_type, const sparse_vector_t& cut, @@ -553,15 +1028,17 @@ f_t knapsack_generation_t::solve_knapsack_problem(const std::vector -void cut_generation_t::generate_cuts(const lp_problem_t& lp, +bool cut_generation_t::generate_cuts(const lp_problem_t& lp, const simplex_solver_settings_t& settings, csr_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& reduced_costs, const std::vector& basic_list, - const std::vector& nonbasic_list) + const std::vector& nonbasic_list, + f_t start_time) { // Generate Gomory and CG Cuts if (settings.mixed_integer_gomory_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { @@ -593,6 +1070,21 @@ void cut_generation_t::generate_cuts(const lp_problem_t& lp, settings.log.debug("MIR and CG cut generation time %.2f seconds\n", cut_generation_time); } } + + // Generate Clique cuts (last to give background clique table generation maximum time) + if (settings.clique_cuts != 0) { + f_t cut_start_time = tic(); + bool feasible = generate_clique_cuts(lp, settings, var_types, xstar, reduced_costs, start_time); + if (!feasible) { + settings.log.printf("Clique cuts proved infeasible\n"); + return false; + } + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Clique cut generation time %.2f seconds\n", cut_generation_time); + } + } + return true; } template @@ -615,6 +1107,279 @@ void cut_generation_t::generate_knapsack_cuts( } } +template +bool cut_generation_t::generate_clique_cuts( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + const std::vector& reduced_costs, + f_t start_time) +{ + if (settings.clique_cuts == 0) { return true; } + if (toc(start_time) >= settings.time_limit) { return true; } + + const i_t num_vars = user_problem_.num_cols; + CLIQUE_CUTS_DEBUG("generate_clique_cuts start num_vars=%lld time_limit=%g elapsed=%g", + static_cast(num_vars), + static_cast(settings.time_limit), + static_cast(toc(start_time))); + + if (clique_table_ == nullptr && clique_table_future_ != nullptr && + clique_table_future_->valid()) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts signaling background thread and waiting"); + if (signal_extend_) { signal_extend_->store(true, std::memory_order_release); } + clique_table_ = clique_table_future_->get(); + clique_table_future_ = nullptr; + if (clique_table_) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts received clique table first=%lld addtl=%lld", + static_cast(clique_table_->first.size()), + static_cast(clique_table_->addtl_cliques.size())); + } + } + + if (clique_table_ == nullptr) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts no clique table available, skipping"); + return true; + } + CLIQUE_CUTS_DEBUG("generate_clique_cuts using clique table first=%lld addtl=%lld", + static_cast(clique_table_->first.size()), + static_cast(clique_table_->addtl_cliques.size())); + + if (clique_table_->first.empty() && clique_table_->addtl_cliques.empty()) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts empty clique table, nothing to separate"); + return true; + } + + cuopt_assert(clique_table_->n_variables == num_vars, "Clique table variable count mismatch"); + cuopt_assert(static_cast(num_vars) <= xstar.size(), "Clique cut xstar size mismatch"); + + const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); + const f_t bound_tol = settings.primal_tol; + const f_t min_weight = 1.0 + min_violation; + // TODO this can be problem dependent + const i_t max_calls = 100000; + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e8; + + cuopt_assert(user_problem_.var_types.size() == static_cast(num_vars), + "User problem var_types size mismatch"); + + std::vector vertices; + std::vector weights; + vertices.reserve(num_vars * 2); + weights.reserve(num_vars * 2); + + // create the sub graph induced by fractional binary variables + for (i_t j = 0; j < num_vars; ++j) { + if (user_problem_.var_types[j] == variable_type_t::CONTINUOUS) { continue; } + const f_t lower_bound = user_problem_.lower[j]; + const f_t upper_bound = user_problem_.upper[j]; + if (lower_bound < -bound_tol || upper_bound > 1 + bound_tol) { continue; } + const f_t xj = xstar[j]; + if (std::abs(xj - std::round(xj)) <= settings.integer_tol) { continue; } + vertices.push_back(j); + weights.push_back(xj); + vertices.push_back(j + num_vars); + weights.push_back(1.0 - xj); + } + // Coarse loop estimate: variable scans + selected vertex/weight writes + work_estimate += 4.0 * static_cast(num_vars) + 2.0 * static_cast(vertices.size()); + if (work_estimate > max_work_estimate) { return true; } + + if (vertices.empty()) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts no fractional binary vertices"); + return true; + } + CLIQUE_CUTS_DEBUG("generate_clique_cuts fractional subgraph vertices=%lld (literals=%lld)", + static_cast(vertices.size() / 2), + static_cast(vertices.size())); + + std::vector vertex_to_local(2 * num_vars, -1); + std::vector in_subgraph(2 * num_vars, 0); + for (size_t idx = 0; idx < vertices.size(); ++idx) { + if (toc(start_time) >= settings.time_limit) { return true; } + const i_t vertex_idx = vertices[idx]; + vertex_to_local[vertex_idx] = static_cast(idx); + in_subgraph[vertex_idx] = 1; + } + work_estimate += 3.0 * static_cast(vertices.size()); + if (work_estimate > max_work_estimate) { return true; } + + std::vector> adj_local(vertices.size()); + size_t total_adj_entries = 0; + size_t kept_adj_entries = 0; + for (size_t idx = 0; idx < vertices.size(); ++idx) { + if (toc(start_time) >= settings.time_limit) { return true; } + i_t vertex_idx = vertices[idx]; + // returns the complement as well + auto adj_set = clique_table_->get_adj_set_of_var(vertex_idx); + total_adj_entries += adj_set.size(); + auto& adj = adj_local[idx]; + adj.reserve(adj_set.size()); + for (const auto neighbor : adj_set) { + if (toc(start_time) >= settings.time_limit) { return true; } + cuopt_assert(neighbor >= 0 && neighbor < 2 * num_vars, "Neighbor out of range"); + if (!in_subgraph[neighbor]) { continue; } + i_t local_neighbor = vertex_to_local[neighbor]; + cuopt_assert(local_neighbor >= 0, "Local neighbor out of range"); + adj.push_back(local_neighbor); + } + kept_adj_entries += adj.size(); +#ifdef ASSERT_MODE + { + std::unordered_set adj_global; + adj_global.reserve(adj.size()); + for (const auto neighbor : adj) { + i_t v = vertices[neighbor]; + cuopt_assert(adj_global.insert(v).second, "Duplicate neighbor in adjacency list"); + i_t complement = (v >= num_vars) ? (v - num_vars) : (v + num_vars); + cuopt_assert(adj_global.find(complement) == adj_global.end(), + "Adjacency list contains complementing variable"); + } + } +#endif + } + work_estimate += static_cast(vertices.size()) + static_cast(total_adj_entries) + + 2.0 * static_cast(kept_adj_entries); + if (work_estimate > max_work_estimate) { return true; } + CLIQUE_CUTS_DEBUG("generate_clique_cuts adjacency raw_entries=%lld kept_entries=%lld", + static_cast(total_adj_entries), + static_cast(kept_adj_entries)); + + const size_t words = bitset_words(vertices.size()); + std::vector> adj_bitset(vertices.size(), std::vector(words, 0)); + size_t local_adj_entries = 0; + for (size_t v = 0; v < adj_local.size(); ++v) { + local_adj_entries += adj_local[v].size(); + for (const auto neighbor : adj_local[v]) { + bitset_set(adj_bitset[v], static_cast(neighbor)); + } + } + work_estimate += static_cast(adj_local.size()) + 3.0 * static_cast(local_adj_entries); + if (work_estimate > max_work_estimate) { return true; } + CLIQUE_CUTS_DEBUG("generate_clique_cuts bitset graph words=%lld local_entries=%lld", + static_cast(words), + static_cast(local_adj_entries)); + + bk_bitset_context_t ctx{adj_bitset, + weights, + min_weight, + max_calls, + start_time, + settings.time_limit, + words, + &work_estimate, + max_work_estimate}; + std::vector R; + std::vector P(words, 0); + std::vector X(words, 0); + for (size_t idx = 0; idx < vertices.size(); ++idx) { + bitset_set(P, idx); + } + work_estimate += 2.0 * static_cast(vertices.size()); + if (work_estimate > max_work_estimate) { return true; } + bron_kerbosch(ctx, R, P, X, 0.0); + CLIQUE_CUTS_DEBUG( + "generate_clique_cuts maximal cliques found=%lld bk_calls=%lld work=%g work_limit=%d " + "call_limit=%d", + static_cast(ctx.cliques.size()), + static_cast(ctx.num_calls), + static_cast(work_estimate), + ctx.over_work_limit() ? 1 : 0, + ctx.over_call_limit() ? 1 : 0); + if (ctx.over_call_limit()) { return true; } + if (ctx.over_work_limit()) { return true; } + if (toc(start_time) >= settings.time_limit) { return true; } + if (work_estimate > max_work_estimate) { return true; } + + sparse_vector_t cut(lp.num_cols, 0); + f_t cut_rhs = 0.0; +#if DEBUG_CLIQUE_CUTS + size_t candidate_cliques = 0; + size_t added_cuts = 0; + size_t rejected_cliques = 0; + size_t extension_gain = 0; +#endif + for (auto& clique_local : ctx.cliques) { + if (toc(start_time) >= settings.time_limit) { return true; } +#if DEBUG_CLIQUE_CUTS + candidate_cliques++; +#endif + std::vector clique_vertices; + clique_vertices.reserve(clique_local.size()); + for (auto local_idx : clique_local) { + clique_vertices.push_back(vertices[local_idx]); + } + work_estimate += 3.0 * static_cast(clique_local.size()); + if (work_estimate > max_work_estimate) { return true; } +#if DEBUG_CLIQUE_CUTS + const size_t size_before_extension = clique_vertices.size(); +#endif + extend_clique_vertices(clique_vertices, + *clique_table_, + xstar, + reduced_costs, + num_vars, + settings.integer_tol, + start_time, + settings.time_limit, + &work_estimate, + max_work_estimate); +#if DEBUG_CLIQUE_CUTS + extension_gain += clique_vertices.size() - size_before_extension; +#endif + if (work_estimate > max_work_estimate) { return true; } + if (toc(start_time) >= settings.time_limit) { return true; } + const auto build_status = build_clique_cut(clique_vertices, + num_vars, + var_types, + user_problem_.lower, + user_problem_.upper, + xstar, + bound_tol, + min_violation, + cut, + cut_rhs, + &work_estimate, + max_work_estimate); + if (work_estimate > max_work_estimate) { return true; } + if (build_status == clique_cut_build_status_t::INFEASIBLE) { + settings.log.debug("Detected contradictory variable/complement clique\n"); + CLIQUE_CUTS_DEBUG( + "generate_clique_cuts infeasible clique detected after processing=%lld cliques", + static_cast(candidate_cliques)); + return false; + } + if (build_status == clique_cut_build_status_t::CUT_ADDED) { + cut_pool_.add_cut(cut_type_t::CLIQUE, cut, cut_rhs); +#if DEBUG_CLIQUE_CUTS + added_cuts++; + CLIQUE_CUTS_DEBUG("generate_clique_cuts added cut nz=%lld rhs=%g clique_size=%lld", + static_cast(cut.i.size()), + static_cast(cut_rhs), + static_cast(clique_vertices.size())); +#endif + } +#if DEBUG_CLIQUE_CUTS + else { + rejected_cliques++; + } +#endif + } +#if DEBUG_CLIQUE_CUTS + CLIQUE_CUTS_DEBUG( + "generate_clique_cuts done candidate_cliques=%lld added=%lld rejected=%lld extension_gain=%lld " + "final_work=%g", + static_cast(candidate_cliques), + static_cast(added_cuts), + static_cast(rejected_cliques), + static_cast(extension_gain), + static_cast(work_estimate)); +#endif + return true; +} + template void cut_generation_t::generate_mir_cuts( const lp_problem_t& lp, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 4f55e96e4d..fffa4b10fe 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -15,6 +15,9 @@ #include #include +#include +#include +#include #include #include #include @@ -22,6 +25,11 @@ #include #include +namespace cuopt::linear_programming::detail { +template +struct clique_table_t; +} + namespace cuopt::linear_programming::dual_simplex { enum cut_type_t : int8_t { @@ -29,9 +37,31 @@ enum cut_type_t : int8_t { MIXED_INTEGER_ROUNDING = 1, KNAPSACK = 2, CHVATAL_GOMORY = 3, - MAX_CUT_TYPE = 4 + CLIQUE = 4, + MAX_CUT_TYPE = 5 +}; + +template +struct cut_gap_closure_t { + f_t initial_gap{0.0}; + f_t final_gap{0.0}; + f_t gap_closed{0.0}; + f_t gap_closed_ratio{0.0}; }; +template +cut_gap_closure_t compute_cut_gap_closure(f_t objective_reference, + f_t objective_before_cuts, + f_t objective_after_cuts) +{ + const f_t initial_gap = std::abs(objective_reference - objective_before_cuts); + const f_t final_gap = std::abs(objective_reference - objective_after_cuts); + const f_t gap_closed = initial_gap - final_gap; + constexpr f_t eps = static_cast(1e-12); + const f_t gap_closed_ratio = initial_gap > eps ? gap_closed / initial_gap : static_cast(0.0); + return {initial_gap, final_gap, gap_closed, gap_closed_ratio}; +} + template struct cut_info_t { bool has_cuts() const @@ -48,8 +78,9 @@ struct cut_info_t { num_cuts[static_cast(cut_type)]++; } } - const char* cut_type_names[MAX_CUT_TYPE] = {"Gomory ", "MIR ", "Knapsack ", "Strong CG"}; - std::array num_cuts = {0}; + const char* cut_type_names[MAX_CUT_TYPE] = { + "Gomory ", "MIR ", "Knapsack ", "Strong CG", "Clique "}; + std::array num_cuts = {0}; }; template @@ -84,6 +115,19 @@ f_t fractional_part(f_t a) return a - std::floor(a); } +template +bool add_work_estimate(f_t accesses, + f_t* work_estimate, + f_t max_work_estimate, + bool* work_limit_reached = nullptr) +{ + if (work_estimate == nullptr) { return false; } + *work_estimate += accesses; + const bool over_work_limit = *work_estimate > max_work_estimate; + if (over_work_limit && work_limit_reached != nullptr) { *work_limit_reached = true; } + return over_work_limit; +} + // Computes a permutation of a score vector that puts the highest scores first template void best_score_first_permutation(std::vector& scores, std::vector& permutation) @@ -119,6 +163,15 @@ void verify_cuts_against_saved_solution(const csr_matrix_t& cuts, const std::vector& cut_rhs, const std::vector& saved_solution); +// Test-only helper to run the production maximal-clique algorithm used by clique cuts. +// adjacency_list must contain local vertex indices in [0, n_vertices). +std::vector> find_maximal_cliques_for_test( + const std::vector>& adjacency_list, + const std::vector& weights, + double min_weight, + int max_calls, + double time_limit); + template class cut_pool_t { public: @@ -221,25 +274,37 @@ class mixed_integer_rounding_cut_t; template class cut_generation_t { public: - cut_generation_t(cut_pool_t& cut_pool, - const lp_problem_t& lp, - const simplex_solver_settings_t& settings, - csr_matrix_t& Arow, - const std::vector& new_slacks, - const std::vector& var_types) - : cut_pool_(cut_pool), knapsack_generation_(lp, settings, Arow, new_slacks, var_types) + cut_generation_t( + cut_pool_t& cut_pool, + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + csr_matrix_t& Arow, + const std::vector& new_slacks, + const std::vector& var_types, + const user_problem_t& user_problem, + std::shared_ptr> clique_table = nullptr, + std::future>>* clique_table_future = nullptr, + std::atomic* signal_extend = nullptr) + : cut_pool_(cut_pool), + knapsack_generation_(lp, settings, Arow, new_slacks, var_types), + user_problem_(user_problem), + clique_table_(std::move(clique_table)), + clique_table_future_(clique_table_future), + signal_extend_(signal_extend) { } - void generate_cuts(const lp_problem_t& lp, + bool generate_cuts(const lp_problem_t& lp, const simplex_solver_settings_t& settings, csr_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& reduced_costs, const std::vector& basic_list, - const std::vector& nonbasic_list); + const std::vector& nonbasic_list, + f_t start_time); private: // Generate all mixed integer gomory cuts @@ -269,8 +334,20 @@ class cut_generation_t { const std::vector& var_types, const std::vector& xstar); + // Generate clique cuts from conflict graph cliques + bool generate_clique_cuts(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + const std::vector& reduced_costs, + f_t start_time); + cut_pool_t& cut_pool_; knapsack_generation_t knapsack_generation_; + const user_problem_t& user_problem_; + std::shared_ptr> clique_table_; + std::future>>* clique_table_future_{nullptr}; + std::atomic* signal_extend_{nullptr}; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 815e229232..a8369a1858 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -100,6 +100,7 @@ struct simplex_solver_settings_t { mir_cuts(-1), mixed_integer_gomory_cuts(-1), knapsack_cuts(-1), + clique_cuts(-1), strong_chvatal_gomory_cuts(-1), reduced_cost_strengthening(-1), cut_change_threshold(1e-3), @@ -178,6 +179,7 @@ struct simplex_solver_settings_t { i_t mixed_integer_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable mixed integer Gomory // cuts i_t knapsack_cuts; // -1 automatic, 0 to disable, >0 to enable knapsack cuts + i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory // cuts i_t reduced_cost_strengthening; // -1 automatic, 0 to disable, >0 to enable reduced cost diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 7435bb37fa..f7d0408aba 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -94,6 +94,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS, &mip_settings.mir_cuts, -1, 1, -1}, {CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS, &mip_settings.mixed_integer_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_KNAPSACK_CUTS, &mip_settings.knapsack_cuts, -1, 1, -1}, + {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, {CUOPT_NUM_GPUS, &pdlp_settings.num_gpus, 1, 2, 1}, diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index ed165fe610..0ded8337d8 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -9,6 +9,7 @@ #include "diversity_manager.cuh" #include + #include #include #include @@ -78,9 +79,8 @@ diversity_manager_t::diversity_manager_t(mip_solver_context_t::n_of_arms, cuopt::seed_generator::get_seed(), ls_alpha, "ls"), ls_hash_map(*context.problem_ptr) { - // Read configuration ID from environment variable - int max_config = -1; - // Read max configuration value from environment variable + int max_config = -1; + int env_config_id = -1; const char* env_max_config = std::getenv("CUOPT_MAX_CONFIG"); if (env_max_config != nullptr) { try { @@ -90,17 +90,21 @@ diversity_manager_t::diversity_manager_t(mip_solver_context_t 1) { - [[maybe_unused]] int config_id = -1; // Default value - const char* env_config_id = std::getenv("CUOPT_CONFIG_ID"); - if (env_config_id != nullptr) { - try { - config_id = std::stoi(env_config_id); - CUOPT_LOG_INFO("Using configuration ID from environment: %d", config_id); - } catch (const std::exception& e) { - CUOPT_LOG_WARN("Failed to parse CUOPT_CONFIG_ID environment variable: %s", e.what()); - } - } + + const char* env_config_id_raw = std::getenv("CUOPT_CONFIG_ID"); + if (env_config_id_raw == nullptr) { return; } + + try { + env_config_id = std::stoi(env_config_id_raw); + } catch (const std::exception& e) { + CUOPT_LOG_WARN("Failed to parse CUOPT_CONFIG_ID environment variable: %s", e.what()); + return; + } + + if (max_config > 0 && env_config_id >= max_config) { + CUOPT_LOG_WARN( + "CUOPT_CONFIG_ID=%d is outside [0, %d). Ignoring cut override.", env_config_id, max_config); + return; } } @@ -205,28 +209,32 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ const bool remap_cache_ids = true; if (!global_timer.check_time_limit()) { trivial_presolve(*problem_ptr, remap_cache_ids); } if (!problem_ptr->empty && !check_bounds_sanity(*problem_ptr)) { return false; } - if (!presolve_timer.check_time_limit() && !context.settings.heuristics_only && - !problem_ptr->empty) { - f_t time_limit_for_clique_table = std::min(3., presolve_timer.remaining_time() / 5); - timer_t clique_timer(time_limit_for_clique_table); - dual_simplex::user_problem_t host_problem(problem_ptr->handle_ptr); - problem_ptr->get_host_user_problem(host_problem); - std::shared_ptr> clique_table; - constexpr bool modify_problem_with_cliques = false; - find_initial_cliques( - host_problem, context.settings.tolerances, clique_timer, modify_problem_with_cliques); - if (modify_problem_with_cliques) { - problem_ptr->set_constraints_from_host_user_problem(host_problem); - cuopt_assert(host_problem.lower.size() == static_cast(problem_ptr->n_variables), - "host lower bound size mismatch"); - cuopt_assert(host_problem.upper.size() == static_cast(problem_ptr->n_variables), - "host upper bound size mismatch"); - std::vector all_var_indices(problem_ptr->n_variables); - std::iota(all_var_indices.begin(), all_var_indices.end(), 0); - problem_ptr->update_variable_bounds(all_var_indices, host_problem.lower, host_problem.upper); - trivial_presolve(*problem_ptr, remap_cache_ids); - } - } + // if (!presolve_timer.check_time_limit() && !context.settings.heuristics_only && + // !problem_ptr->empty) { + // f_t time_limit_for_clique_table = std::min(3., presolve_timer.remaining_time() / 5); + // timer_t clique_timer(time_limit_for_clique_table); + // dual_simplex::user_problem_t host_problem(problem_ptr->handle_ptr); + // problem_ptr->get_host_user_problem(host_problem); + // std::shared_ptr> clique_table; + // constexpr bool modify_problem_with_cliques = false; + // find_initial_cliques(host_problem, + // context.settings.tolerances, + // &clique_table, + // clique_timer, + // modify_problem_with_cliques, + // nullptr); + // if (modify_problem_with_cliques) { + // problem_ptr->set_constraints_from_host_user_problem(host_problem); + // cuopt_assert(host_problem.lower.size() == static_cast(problem_ptr->n_variables), + // "host lower bound size mismatch"); + // cuopt_assert(host_problem.upper.size() == static_cast(problem_ptr->n_variables), + // "host upper bound size mismatch"); + // std::vector all_var_indices(problem_ptr->n_variables); + // std::iota(all_var_indices.begin(), all_var_indices.end(), 0); + // problem_ptr->update_variable_bounds(all_var_indices, host_problem.lower, + // host_problem.upper); trivial_presolve(*problem_ptr, remap_cache_ids); + // } + // } // May overconstrain if Papilo presolve has been run before if (context.settings.presolver == presolver_t::None) { if (!problem_ptr->empty) { diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cuh b/cpp/src/mip_heuristics/diversity/diversity_manager.cuh index d4e24bdeaf..4f86192db8 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cuh +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cuh @@ -27,6 +27,8 @@ #include #include +#include + namespace cuopt::linear_programming::detail { template @@ -69,7 +71,6 @@ class diversity_manager_t { void set_simplex_solution(const std::vector& solution, const std::vector& dual_solution, f_t objective); - mip_solver_context_t& context; dual_simplex::branch_and_bound_t* branch_and_bound_ptr; problem_t* problem_ptr; diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index 2fbe79ba34..d7d7601014 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -266,6 +266,7 @@ void rins_t::run_rins() branch_and_bound_settings.num_threads = 1; branch_and_bound_settings.reliability_branching = 0; branch_and_bound_settings.max_cut_passes = 0; + branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.sub_mip = 1; branch_and_bound_settings.log.log = false; branch_and_bound_settings.log.log_prefix = "[RINS] "; diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index b2f7f80066..a867141d0a 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -108,6 +108,7 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.num_threads = 1; branch_and_bound_settings.reliability_branching = 0; branch_and_bound_settings.max_cut_passes = 0; + branch_and_bound_settings.clique_cuts = 0; branch_and_bound_settings.sub_mip = 1; branch_and_bound_settings.solution_callback = [this](std::vector& solution, f_t objective) { diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu index e21af7f69f..70855267df 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu @@ -15,7 +15,7 @@ * limitations under the License. */ -#define DEBUG_KNAPSACK_CONSTRAINTS 1 +#define DEBUG_KNAPSACK_CONSTRAINTS 0 #include "clique_table.cuh" @@ -34,11 +34,13 @@ namespace cuopt::linear_programming::detail { // do constraints with only binary variables. template void find_cliques_from_constraint(const knapsack_constraint_t& kc, - clique_table_t& clique_table) + clique_table_t& clique_table, + cuopt::timer_t& timer) { i_t size = kc.entries.size(); cuopt_assert(size > 1, "Constraint has not enough variables"); if (kc.entries[size - 1].val + kc.entries[size - 2].val <= kc.rhs) { return; } + std::vector clique; i_t k = size - 1; // find the first clique, which is the largest @@ -55,6 +57,7 @@ void find_cliques_from_constraint(const knapsack_constraint_t& kc, // find the additional cliques k--; while (k >= 0) { + if (timer.check_time_limit()) { return; } f_t curr_val = kc.entries[k].val; i_t curr_col = kc.entries[k].col; // do a binary search in the clique coefficients to find f, such that coeff_k + coeff_f > rhs @@ -210,13 +213,14 @@ void fill_knapsack_constraints(const dual_simplex::user_problem_t& pro } template -void remove_small_cliques(clique_table_t& clique_table) +void remove_small_cliques(clique_table_t& clique_table, cuopt::timer_t& timer) { i_t num_removed_first = 0; i_t num_removed_addtl = 0; std::vector to_delete(clique_table.first.size(), false); // if a clique is small, we remove it from the cliques and add it to adjlist for (size_t clique_idx = 0; clique_idx < clique_table.first.size(); clique_idx++) { + if (timer.check_time_limit()) { return; } const auto& clique = clique_table.first[clique_idx]; if (clique.size() <= (size_t)clique_table.min_clique_size) { for (size_t i = 0; i < clique.size(); i++) { @@ -316,15 +320,16 @@ std::unordered_set clique_table_t::get_adj_set_of_var(i_t var_idx addtl_cliques[addtl_clique_idx].start_pos_on_clique, first[addtl_cliques[addtl_clique_idx].clique_idx].end()); } - // Memory-neutral reverse lookup for additional cliques: + // Reverse lookup for additional cliques using position map: // if var_idx is in first[clique_idx][start_pos_on_clique:], it is adjacent to vertex_idx. for (const auto& addtl : addtl_cliques) { if (addtl.vertex_idx == var_idx) { continue; } - const auto& clique = first[addtl.clique_idx]; - size_t start_pos = static_cast(addtl.start_pos_on_clique); - if (start_pos < clique.size() && - std::find(clique.begin() + start_pos, clique.end(), var_idx) != clique.end()) { - adj_set.insert(addtl.vertex_idx); + if (static_cast(addtl.clique_idx) < first_var_positions.size()) { + const auto& pos_map = first_var_positions[addtl.clique_idx]; + auto it = pos_map.find(var_idx); + if (it != pos_map.end() && it->second >= addtl.start_pos_on_clique) { + adj_set.insert(addtl.vertex_idx); + } } } @@ -349,44 +354,39 @@ i_t clique_table_t::get_degree_of_var(i_t var_idx) template bool clique_table_t::check_adjacency(i_t var_idx1, i_t var_idx2) { - // if passed same variable if (var_idx1 == var_idx2) { return false; } - // in case they are complements of each other if (var_idx1 % n_variables == var_idx2 % n_variables) { return true; } - if (adj_list_small_cliques[var_idx1].count(var_idx2) > 0) { return true; } - // Check first cliques: var_clique_map_first stores clique indices - for (const auto& clique_idx : var_clique_map_first[var_idx1]) { - const auto& clique = first[clique_idx]; - // TODO: we can also keep a set of the clique if the memory allows, instead of doing linear - // search - if (std::find(clique.begin(), clique.end(), var_idx2) != clique.end()) { return true; } + + { + auto it = adj_list_small_cliques.find(var_idx1); + if (it != adj_list_small_cliques.end() && it->second.count(var_idx2) > 0) { return true; } } - // Check additional cliques: var_clique_map_addtl stores indices into addtl_cliques - for (const auto& addtl_idx : var_clique_map_addtl[var_idx1]) { - const auto& addtl = addtl_cliques[addtl_idx]; - const auto& clique = first[addtl.clique_idx]; - // addtl clique is: vertex_idx + first[clique_idx][start_pos_on_clique:] - if (addtl.vertex_idx == var_idx2) { return true; } - if (addtl.start_pos_on_clique < static_cast(clique.size())) { - if (std::find(clique.begin() + addtl.start_pos_on_clique, clique.end(), var_idx2) != - clique.end()) { - return true; - } + // Iterate whichever variable belongs to fewer first-cliques + { + i_t probe_var = var_idx1; + i_t target_var = var_idx2; + if (var_clique_map_first[var_idx1].size() > var_clique_map_first[var_idx2].size()) { + probe_var = var_idx2; + target_var = var_idx1; + } + for (const auto& clique_idx : var_clique_map_first[probe_var]) { + if (first_var_positions[clique_idx].count(target_var) > 0) { return true; } } } - // var_clique_map_addtl is keyed by addtl.vertex_idx, so also check the reverse direction. + for (const auto& addtl_idx : var_clique_map_addtl[var_idx1]) { + const auto& addtl = addtl_cliques[addtl_idx]; + const auto& pos_map = first_var_positions[addtl.clique_idx]; + auto it = pos_map.find(var_idx2); + if (it != pos_map.end() && it->second >= addtl.start_pos_on_clique) { return true; } + } + for (const auto& addtl_idx : var_clique_map_addtl[var_idx2]) { - const auto& addtl = addtl_cliques[addtl_idx]; - const auto& clique = first[addtl.clique_idx]; - if (addtl.vertex_idx == var_idx1) { return true; } - if (addtl.start_pos_on_clique < static_cast(clique.size())) { - if (std::find(clique.begin() + addtl.start_pos_on_clique, clique.end(), var_idx1) != - clique.end()) { - return true; - } - } + const auto& addtl = addtl_cliques[addtl_idx]; + const auto& pos_map = first_var_positions[addtl.clique_idx]; + auto it = pos_map.find(var_idx1); + if (it != pos_map.end() && it->second >= addtl.start_pos_on_clique) { return true; } } return false; @@ -548,6 +548,232 @@ bool extend_clique(const std::vector& clique, return new_clique.size() > clique.size(); } +template +struct clique_sig_t { + i_t knapsack_idx; + i_t size; + long long signature; +}; + +template +struct extension_candidate_t { + i_t knapsack_idx; + i_t estimated_gain; + i_t clique_size; +}; + +template +bool compare_clique_sig(const clique_sig_t& a, const clique_sig_t& b) +{ + if (a.signature != b.signature) { return a.signature < b.signature; } + return a.size < b.size; +} + +template +bool compare_signature_value(long long value, const clique_sig_t& a) +{ + return value < a.signature; +} + +template +bool compare_extension_candidate(const extension_candidate_t& a, + const extension_candidate_t& b) +{ + if (a.estimated_gain != b.estimated_gain) { return a.estimated_gain > b.estimated_gain; } + if (a.clique_size != b.clique_size) { return a.clique_size < b.clique_size; } + return a.knapsack_idx < b.knapsack_idx; +} + +template +bool is_sorted_subset(const std::vector& a, const std::vector& b) +{ + size_t i = 0; + size_t j = 0; + while (i < a.size() && j < b.size()) { + if (a[i] == b[j]) { + i++; + j++; + } else if (a[i] > b[j]) { + j++; + } else { + return false; + } + } + return i == a.size(); +} + +template +void fix_difference(const std::vector& superset, + const std::vector& subset, + dual_simplex::user_problem_t& problem) +{ + cuopt_assert(std::is_sorted(subset.begin(), subset.end()), + "subset vector passed to fix_difference is not sorted"); + for (auto var_idx : superset) { + if (std::binary_search(subset.begin(), subset.end(), var_idx)) { continue; } + if (var_idx >= problem.num_cols) { + i_t orig_idx = var_idx - problem.num_cols; + CUOPT_LOG_DEBUG("Fixing variable %d", orig_idx); + cuopt_assert(problem.lower[orig_idx] != 0 || problem.upper[orig_idx] != 0, + "Variable is fixed to other side"); + problem.lower[orig_idx] = 1; + problem.upper[orig_idx] = 1; + } else { + CUOPT_LOG_DEBUG("Fixing variable %d", var_idx); + cuopt_assert(problem.lower[var_idx] != 1 || problem.upper[var_idx] != 1, + "Variable is fixed to other side"); + problem.lower[var_idx] = 0; + problem.upper[var_idx] = 0; + } + } +} + +template +void remove_marked_elements(std::vector& vec, const std::vector& removal_marker) +{ + size_t write_idx = 0; + for (size_t i = 0; i < vec.size(); i++) { + if (!removal_marker[i]) { + if (write_idx != i) { vec[write_idx] = std::move(vec[i]); } + write_idx++; + } + } + vec.resize(write_idx); +} + +template +void remove_dominated_cliques_in_problem_for_single_extended_clique( + const std::vector& curr_clique, + f_t coeff_scale, + i_t remaining_rows_budget, + i_t remaining_nnz_budget, + i_t& inserted_row_nnz, + const std::vector>& sp_sigs, + const std::vector>& cstr_vars, + const std::vector>& knapsack_constraints, + std::vector& original_to_current_row_idx, + dual_simplex::user_problem_t& problem, + dual_simplex::csr_matrix_t& A, + cuopt::timer_t& timer) +{ + inserted_row_nnz = 0; + if (curr_clique.empty() || sp_sigs.empty()) { return; } + std::vector curr_clique_vars(curr_clique.begin(), curr_clique.end()); + std::sort(curr_clique_vars.begin(), curr_clique_vars.end()); + curr_clique_vars.erase(std::unique(curr_clique_vars.begin(), curr_clique_vars.end()), + curr_clique_vars.end()); + long long signature = 0; + for (auto v : curr_clique_vars) { + signature += static_cast(v); + } + constexpr size_t dominance_window = 20000; + auto end_it = + std::upper_bound(sp_sigs.begin(), sp_sigs.end(), signature, compare_signature_value); + size_t end = static_cast(std::distance(sp_sigs.begin(), end_it)); + size_t start = (end > dominance_window) ? (end - dominance_window) : 0; + std::vector rows_to_remove; + bool covering_clique_implied_by_partitioning = false; + for (size_t idx = end; idx > start; idx--) { + if (timer.check_time_limit()) { break; } + const auto& sp = sp_sigs[idx - 1]; + const auto& vars_sp = cstr_vars[sp.knapsack_idx]; + if (vars_sp.size() > curr_clique_vars.size()) { continue; } + cuopt_assert(std::is_sorted(vars_sp.begin(), vars_sp.end()), + "vars_sp vector passed to is_sorted_subset is not sorted"); + if (!is_sorted_subset(vars_sp, curr_clique_vars)) { continue; } + if (knapsack_constraints[sp.knapsack_idx].is_set_partitioning) { + if (vars_sp.size() != curr_clique_vars.size()) { + fix_difference(curr_clique_vars, vars_sp, problem); + covering_clique_implied_by_partitioning = true; + } + continue; + } + i_t original_row_idx = knapsack_constraints[sp.knapsack_idx].cstr_idx; + if (original_row_idx < 0) { continue; } + cuopt_assert(original_row_idx < static_cast(original_to_current_row_idx.size()), + "Invalid original row index in knapsack constraint"); + i_t current_row_idx = original_to_current_row_idx[original_row_idx]; + if (current_row_idx < 0) { continue; } + cuopt_assert(current_row_idx < static_cast(problem.row_sense.size()), + "Invalid current row index in row mapping"); + rows_to_remove.push_back(current_row_idx); + } + if (rows_to_remove.empty()) { return; } + std::sort(rows_to_remove.begin(), rows_to_remove.end()); + rows_to_remove.erase(std::unique(rows_to_remove.begin(), rows_to_remove.end()), + rows_to_remove.end()); + if (!covering_clique_implied_by_partitioning) { + if (remaining_rows_budget <= 0 || + remaining_nnz_budget < static_cast(curr_clique_vars.size())) { + return; + } + insert_clique_into_problem(curr_clique_vars, problem, A, coeff_scale); + inserted_row_nnz = static_cast(curr_clique_vars.size()); + } + std::vector removal_marker(problem.row_sense.size(), 0); + for (auto row_idx : rows_to_remove) { + cuopt_assert(row_idx >= 0 && row_idx < static_cast(removal_marker.size()), + "Invalid dominated row index"); + CUOPT_LOG_DEBUG("Removing dominated row %d", row_idx); + removal_marker[row_idx] = true; + } + dual_simplex::csr_matrix_t A_removed(0, 0, 0); + A.remove_rows(removal_marker, A_removed); + A = std::move(A_removed); + problem.num_rows = A.m; + remove_marked_elements(problem.row_sense, removal_marker); + remove_marked_elements(problem.rhs, removal_marker); + remove_marked_elements(problem.row_names, removal_marker); + cuopt_assert(problem.rhs.size() == problem.row_sense.size(), "rhs and row sense size mismatch"); + cuopt_assert(problem.row_names.size() == problem.rhs.size(), "row names and rhs size mismatch"); + cuopt_assert(problem.num_rows == static_cast(problem.rhs.size()), + "matrix and num rows mismatch after removal"); + if (!problem.range_rows.empty()) { + std::vector old_to_new_indices; + old_to_new_indices.reserve(removal_marker.size()); + i_t new_idx = 0; + for (size_t i = 0; i < removal_marker.size(); ++i) { + if (!removal_marker[i]) { + old_to_new_indices.push_back(new_idx++); + } else { + old_to_new_indices.push_back(-1); + } + } + std::vector new_range_rows; + std::vector new_range_values; + for (size_t i = 0; i < problem.range_rows.size(); ++i) { + i_t old_row = problem.range_rows[i]; + cuopt_assert(old_row >= 0 && old_row < static_cast(removal_marker.size()), + "Invalid row index in range_rows"); + if (!removal_marker[old_row]) { + i_t new_row = old_to_new_indices[old_row]; + cuopt_assert(new_row != -1, "Invalid new row index for ranged row renumbering"); + new_range_rows.push_back(new_row); + new_range_values.push_back(problem.range_value[i]); + } + } + problem.range_rows = std::move(new_range_rows); + problem.range_value = std::move(new_range_values); + } + problem.num_range_rows = static_cast(problem.range_rows.size()); + std::vector removed_prefix(removal_marker.size() + 1, 0); + for (size_t row_idx = 0; row_idx < removal_marker.size(); row_idx++) { + removed_prefix[row_idx + 1] = + removed_prefix[row_idx] + static_cast(removal_marker[row_idx]); + } + for (i_t row_idx = 0; row_idx < static_cast(original_to_current_row_idx.size()); row_idx++) { + i_t current_row_idx = original_to_current_row_idx[row_idx]; + if (current_row_idx < 0) { continue; } + cuopt_assert(current_row_idx < static_cast(removal_marker.size()), + "Row index map is out of bounds"); + if (removal_marker[current_row_idx]) { + original_to_current_row_idx[row_idx] = -1; + } else { + original_to_current_row_idx[row_idx] = current_row_idx - removed_prefix[current_row_idx]; + } + } +} + // Also known as clique merging. Infer larger clique constraints which allows inclusion of vars from // other constraints. This only extends the original cliques in the formulation for now. // TODO: consider a heuristic on how much of the cliques derived from knapsacks to include here @@ -558,12 +784,17 @@ i_t extend_cliques(const std::vector>& knapsack_ dual_simplex::user_problem_t& problem, dual_simplex::csr_matrix_t& A, bool modify_problem, - cuopt::timer_t& timer) + cuopt::timer_t& timer, + double* work_estimate_out, + double max_work_estimate) { constexpr i_t min_extension_gain = 2; constexpr i_t extension_yield_window = 64; constexpr i_t min_successes_per_window = 1; + double local_work = 0.0; + double& work = work_estimate_out ? *work_estimate_out : local_work; + i_t base_rows = A.m; i_t base_nnz = A.row_start[A.m]; i_t max_added_rows = std::max(8, base_rows / 50); @@ -579,12 +810,7 @@ i_t extend_cliques(const std::vector>& knapsack_ max_added_rows, max_added_nnz); std::vector> cstr_vars(knapsack_constraints.size()); - struct clique_sig_t { - i_t knapsack_idx; - i_t size; - long long signature; - }; - std::vector sp_sigs; + std::vector> sp_sigs; sp_sigs.reserve(set_packing_constraints.size()); for (const auto knapsack_idx : set_packing_constraints) { cuopt_assert(knapsack_idx >= 0 && knapsack_idx < static_cast(knapsack_constraints.size()), @@ -603,207 +829,20 @@ i_t extend_cliques(const std::vector>& knapsack_ signature += static_cast(v); } sp_sigs.push_back({knapsack_idx, static_cast(cstr_vars[knapsack_idx].size()), signature}); + work += cstr_vars[knapsack_idx].size(); } - std::sort(sp_sigs.begin(), sp_sigs.end(), [](const auto& a, const auto& b) { - if (a.signature != b.signature) { return a.signature < b.signature; } - return a.size < b.size; - }); + if (work > max_work_estimate) { return 0; } + std::sort(sp_sigs.begin(), sp_sigs.end(), compare_clique_sig); std::vector original_to_current_row_idx(problem.row_sense.size(), -1); for (i_t row_idx = 0; row_idx < static_cast(original_to_current_row_idx.size()); row_idx++) { original_to_current_row_idx[row_idx] = row_idx; } - auto is_subset = [](const std::vector& a, const std::vector& b) { - size_t i = 0; - size_t j = 0; - while (i < a.size() && j < b.size()) { - if (a[i] == b[j]) { - i++; - j++; - } else if (a[i] > b[j]) { - j++; - } else { - return false; - } - } - return i == a.size(); - }; - auto fix_difference = [&](const std::vector& superset, const std::vector& subset) { - if (!modify_problem) { return; } - cuopt_assert(std::is_sorted(subset.begin(), subset.end()), - "subset vector passed to fix_difference is not sorted"); - for (auto var_idx : superset) { - if (std::binary_search(subset.begin(), subset.end(), var_idx)) { continue; } - if (var_idx >= problem.num_cols) { - i_t orig_idx = var_idx - problem.num_cols; - CUOPT_LOG_DEBUG("Fixing variable %d", orig_idx); - cuopt_assert(problem.lower[orig_idx] != 0 || problem.upper[orig_idx] != 0, - "Variable is fixed to other side"); - problem.lower[orig_idx] = 1; - problem.upper[orig_idx] = 1; - } else { - CUOPT_LOG_DEBUG("Fixing variable %d", var_idx); - cuopt_assert(problem.lower[var_idx] != 1 || problem.upper[var_idx] != 1, - "Variable is fixed to other side"); - problem.lower[var_idx] = 0; - problem.upper[var_idx] = 0; - } - } - }; - auto remove_dominated_cliques_in_problem_for_single_extended_clique = - [&](const std::vector& curr_clique, - f_t coeff_scale, - i_t remaining_rows_budget, - i_t remaining_nnz_budget, - i_t& inserted_row_nnz) { - inserted_row_nnz = 0; - if (curr_clique.empty() || sp_sigs.empty()) { return; } - std::vector curr_clique_vars(curr_clique.begin(), curr_clique.end()); - std::sort(curr_clique_vars.begin(), curr_clique_vars.end()); - curr_clique_vars.erase(std::unique(curr_clique_vars.begin(), curr_clique_vars.end()), - curr_clique_vars.end()); - long long signature = 0; - for (auto v : curr_clique_vars) { - signature += static_cast(v); - } - constexpr size_t dominance_window = 20000; - auto end_it = std::upper_bound( - sp_sigs.begin(), sp_sigs.end(), signature, [](long long value, const auto& a) { - return value < a.signature; - }); - size_t end = static_cast(std::distance(sp_sigs.begin(), end_it)); - size_t start = (end > dominance_window) ? (end - dominance_window) : 0; - std::vector rows_to_remove; - bool covering_clique_implied_by_partitioning = false; - for (size_t idx = end; idx > start; idx--) { - if (timer.check_time_limit()) { break; } - const auto& sp = sp_sigs[idx - 1]; - const auto& vars_sp = cstr_vars[sp.knapsack_idx]; - if (vars_sp.size() > curr_clique_vars.size()) { continue; } - cuopt_assert(std::is_sorted(vars_sp.begin(), vars_sp.end()), - "vars_sp vector passed to is_subset is not sorted"); - if (!is_subset(vars_sp, curr_clique_vars)) { continue; } - if (knapsack_constraints[sp.knapsack_idx].is_set_partitioning) { - if (vars_sp.size() != curr_clique_vars.size()) { - fix_difference(curr_clique_vars, vars_sp); - covering_clique_implied_by_partitioning = true; - } - continue; - } - i_t original_row_idx = knapsack_constraints[sp.knapsack_idx].cstr_idx; - if (original_row_idx < 0) { continue; } - cuopt_assert(original_row_idx < static_cast(original_to_current_row_idx.size()), - "Invalid original row index in knapsack constraint"); - i_t current_row_idx = original_to_current_row_idx[original_row_idx]; - if (current_row_idx < 0) { continue; } - cuopt_assert(current_row_idx < static_cast(problem.row_sense.size()), - "Invalid current row index in row mapping"); - rows_to_remove.push_back(current_row_idx); - } - if (rows_to_remove.empty()) { return; } - if (!modify_problem) { return; } - std::sort(rows_to_remove.begin(), rows_to_remove.end()); - rows_to_remove.erase(std::unique(rows_to_remove.begin(), rows_to_remove.end()), - rows_to_remove.end()); - if (!covering_clique_implied_by_partitioning) { - if (remaining_rows_budget <= 0 || - remaining_nnz_budget < static_cast(curr_clique_vars.size())) { - return; - } - // Replace dominated rows with this stronger clique row. - insert_clique_into_problem(curr_clique_vars, problem, A, coeff_scale); - inserted_row_nnz = static_cast(curr_clique_vars.size()); - } - std::vector removal_marker(problem.row_sense.size(), 0); - for (auto row_idx : rows_to_remove) { - cuopt_assert(row_idx >= 0 && row_idx < static_cast(removal_marker.size()), - "Invalid dominated row index"); - CUOPT_LOG_DEBUG("Removing dominated row %d", row_idx); - removal_marker[row_idx] = true; - } - dual_simplex::csr_matrix_t A_removed(0, 0, 0); - A.remove_rows(removal_marker, A_removed); - A = std::move(A_removed); - problem.num_rows = A.m; - i_t n = 0; - auto new_end = std::remove_if( - problem.row_sense.begin(), problem.row_sense.end(), [&removal_marker, &n](char) mutable { - return removal_marker[n++]; - }); - problem.row_sense.erase(new_end, problem.row_sense.end()); - n = 0; - auto new_end_rhs = - std::remove_if(problem.rhs.begin(), problem.rhs.end(), [&removal_marker, &n](f_t) mutable { - return removal_marker[n++]; - }); - problem.rhs.erase(new_end_rhs, problem.rhs.end()); - n = 0; - auto new_end_row_names = std::remove_if( - problem.row_names.begin(), - problem.row_names.end(), - [&removal_marker, &n](const std::string&) mutable { return removal_marker[n++]; }); - problem.row_names.erase(new_end_row_names, problem.row_names.end()); - cuopt_assert(problem.rhs.size() == problem.row_sense.size(), - "rhs and row sense size mismatch"); - cuopt_assert(problem.row_names.size() == problem.rhs.size(), - "row names and rhs size mismatch"); - cuopt_assert(problem.num_rows == static_cast(problem.rhs.size()), - "matrix and num rows mismatch after removal"); - if (!problem.range_rows.empty()) { - std::vector old_to_new_indices; - old_to_new_indices.reserve(removal_marker.size()); - i_t new_idx = 0; - for (size_t i = 0; i < removal_marker.size(); ++i) { - if (!removal_marker[i]) { - old_to_new_indices.push_back(new_idx++); - } else { - old_to_new_indices.push_back(-1); - } - } - std::vector new_range_rows; - std::vector new_range_values; - for (size_t i = 0; i < problem.range_rows.size(); ++i) { - i_t old_row = problem.range_rows[i]; - cuopt_assert(old_row >= 0 && old_row < static_cast(removal_marker.size()), - "Invalid row index in range_rows"); - if (!removal_marker[old_row]) { - i_t new_row = old_to_new_indices[old_row]; - cuopt_assert(new_row != -1, "Invalid new row index for ranged row renumbering"); - new_range_rows.push_back(new_row); - new_range_values.push_back(problem.range_value[i]); - } - } - problem.range_rows = std::move(new_range_rows); - problem.range_value = std::move(new_range_values); - } - problem.num_range_rows = static_cast(problem.range_rows.size()); - std::vector removed_prefix(removal_marker.size() + 1, 0); - for (size_t row_idx = 0; row_idx < removal_marker.size(); row_idx++) { - removed_prefix[row_idx + 1] = - removed_prefix[row_idx] + static_cast(removal_marker[row_idx]); - } - for (i_t row_idx = 0; row_idx < static_cast(original_to_current_row_idx.size()); - row_idx++) { - i_t current_row_idx = original_to_current_row_idx[row_idx]; - if (current_row_idx < 0) { continue; } - cuopt_assert(current_row_idx < static_cast(removal_marker.size()), - "Row index map is out of bounds"); - if (removal_marker[current_row_idx]) { - original_to_current_row_idx[row_idx] = -1; - } else { - original_to_current_row_idx[row_idx] = current_row_idx - removed_prefix[current_row_idx]; - } - } - }; - struct extension_candidate_t { - i_t knapsack_idx; - i_t estimated_gain; - i_t clique_size; - }; - std::vector extension_worklist; + std::vector> extension_worklist; extension_worklist.reserve(knapsack_constraints.size()); for (i_t knapsack_idx = 0; knapsack_idx < static_cast(knapsack_constraints.size()); knapsack_idx++) { if (timer.check_time_limit()) { break; } + if (work > max_work_estimate) { break; } const auto& knapsack_constraint = knapsack_constraints[knapsack_idx]; if (!knapsack_constraint.is_set_packing) { continue; } i_t clique_size = static_cast(knapsack_constraint.entries.size()); @@ -812,26 +851,19 @@ i_t extend_cliques(const std::vector>& knapsack_ for (const auto& entry : knapsack_constraint.entries) { smallest_degree = std::min(smallest_degree, clique_table.get_degree_of_var(entry.col)); } - // The smallest-degree vertex upper-bounds how many new literals can be added. i_t estimated_gain = std::max(0, smallest_degree - (clique_size - 1)); if (estimated_gain < min_extension_gain) { continue; } extension_worklist.push_back({knapsack_idx, estimated_gain, clique_size}); + work += knapsack_constraint.entries.size(); } - std::stable_sort(extension_worklist.begin(), - extension_worklist.end(), - [](const extension_candidate_t& a, const extension_candidate_t& b) { - if (a.estimated_gain != b.estimated_gain) { - return a.estimated_gain > b.estimated_gain; - } - if (a.clique_size != b.clique_size) { return a.clique_size < b.clique_size; } - return a.knapsack_idx < b.knapsack_idx; - }); + std::stable_sort( + extension_worklist.begin(), extension_worklist.end(), compare_extension_candidate); CUOPT_LOG_DEBUG("Clique extension candidates after scoring: %zu", extension_worklist.size()); i_t n_extended_cliques = 0; - // Try highest estimated gain candidates first so budget is spent on promising rows. for (const auto& candidate : extension_worklist) { if (timer.check_time_limit()) { break; } + if (work > max_work_estimate) { break; } if (added_rows >= max_added_rows || added_nnz >= max_added_nnz) { CUOPT_LOG_DEBUG( "Stopping clique extension: budget reached (rows=%d nnz=%d)", added_rows, added_nnz); @@ -855,14 +887,24 @@ i_t extend_cliques(const std::vector>& knapsack_ max_added_rows - added_rows, max_added_nnz - added_nnz, inserted_row_nnz); + work += clique.size() * clique.size(); if (extended_clique) { n_extended_cliques++; i_t replacement_row_nnz = 0; - remove_dominated_cliques_in_problem_for_single_extended_clique(clique_table.first.back(), - coeff_scale, - max_added_rows - added_rows, - max_added_nnz - added_nnz, - replacement_row_nnz); + if (modify_problem) { + remove_dominated_cliques_in_problem_for_single_extended_clique(clique_table.first.back(), + coeff_scale, + max_added_rows - added_rows, + max_added_nnz - added_nnz, + replacement_row_nnz, + sp_sigs, + cstr_vars, + knapsack_constraints, + original_to_current_row_idx, + problem, + A, + timer); + } if (replacement_row_nnz > 0) { window_successes++; added_rows++; @@ -890,11 +932,15 @@ i_t extend_cliques(const std::vector>& knapsack_ template void fill_var_clique_maps(clique_table_t& clique_table) { + clique_table.first_var_positions.resize(clique_table.first.size()); for (size_t clique_idx = 0; clique_idx < clique_table.first.size(); clique_idx++) { const auto& clique = clique_table.first[clique_idx]; + auto& pos_map = clique_table.first_var_positions[clique_idx]; + pos_map.reserve(clique.size()); for (size_t idx = 0; idx < clique.size(); idx++) { i_t var_idx = clique[idx]; clique_table.var_clique_map_first[var_idx].insert(clique_idx); + pos_map[var_idx] = static_cast(idx); } } for (size_t addtl_c = 0; addtl_c < clique_table.addtl_cliques.size(); addtl_c++) { @@ -903,6 +949,37 @@ void fill_var_clique_maps(clique_table_t& clique_table) } } +template +void build_clique_table(const dual_simplex::user_problem_t& problem, + clique_table_t& clique_table, + typename mip_solver_settings_t::tolerances_t tolerances, + bool remove_small_cliques_flag, + bool fill_var_clique_maps_flag, + cuopt::timer_t& timer) +{ + if (timer.check_time_limit()) { return; } + cuopt_assert(clique_table.n_variables == problem.num_cols, "Clique table size mismatch"); + cuopt_assert(problem.var_types.size() == static_cast(problem.num_cols), + "Problem variable types size mismatch"); + std::vector> knapsack_constraints; + std::unordered_set set_packing_constraints; + dual_simplex::csr_matrix_t A(problem.num_rows, problem.num_cols, 0); + problem.A.to_compressed_row(A); + fill_knapsack_constraints(problem, knapsack_constraints, A); + make_coeff_positive_knapsack_constraint( + problem, knapsack_constraints, set_packing_constraints, tolerances); + sort_csr_by_constraint_coefficients(knapsack_constraints); + clique_table.tolerances = tolerances; + for (const auto& knapsack_constraint : knapsack_constraints) { + if (timer.check_time_limit()) { return; } + find_cliques_from_constraint(knapsack_constraint, clique_table, timer); + } + if (timer.check_time_limit()) { return; } + if (remove_small_cliques_flag) { remove_small_cliques(clique_table, timer); } + if (timer.check_time_limit()) { return; } + if (fill_var_clique_maps_flag) { fill_var_clique_maps(clique_table); } +} + template void print_knapsack_constraints( const std::vector>& knapsack_constraints, @@ -946,8 +1023,10 @@ void print_clique_table(const clique_table_t& clique_table) template void find_initial_cliques(dual_simplex::user_problem_t& problem, typename mip_solver_settings_t::tolerances_t tolerances, + std::shared_ptr>* clique_table_out, cuopt::timer_t& timer, - bool modify_problem) + bool modify_problem, + std::atomic* signal_extend) { cuopt::timer_t stage_timer(std::numeric_limits::infinity()); #ifdef DEBUG_CLIQUE_TABLE @@ -977,43 +1056,61 @@ void find_initial_cliques(dual_simplex::user_problem_t& problem, #ifdef DEBUG_CLIQUE_TABLE t_sort = stage_timer.elapsed_time(); #endif - // print_knapsack_constraints(knapsack_constraints); - // TODO think about getting min_clique_size according to some problem property clique_config_t clique_config; - clique_table_t clique_table(2 * problem.num_cols, - clique_config.min_clique_size, - clique_config.max_clique_size_for_extension); - clique_table.tolerances = tolerances; + std::shared_ptr> clique_table_shared; + clique_table_t clique_table_local(2 * problem.num_cols, + clique_config.min_clique_size, + clique_config.max_clique_size_for_extension); + clique_table_t* clique_table_ptr = &clique_table_local; + if (clique_table_out != nullptr) { + clique_table_shared = + std::make_shared>(2 * problem.num_cols, + clique_config.min_clique_size, + clique_config.max_clique_size_for_extension); + clique_table_ptr = clique_table_shared.get(); + } + clique_table_ptr->tolerances = tolerances; + double time_limit_for_additional_cliques = timer.remaining_time() / 2; + cuopt::timer_t additional_cliques_timer(time_limit_for_additional_cliques); + double find_work_estimate = 0.0; for (const auto& knapsack_constraint : knapsack_constraints) { if (timer.check_time_limit()) { break; } - find_cliques_from_constraint(knapsack_constraint, clique_table); + if (signal_extend && signal_extend->load(std::memory_order_acquire)) { break; } + find_cliques_from_constraint(knapsack_constraint, *clique_table_ptr, additional_cliques_timer); + find_work_estimate += knapsack_constraint.entries.size(); } - if (timer.check_time_limit()) { return; } #ifdef DEBUG_CLIQUE_TABLE t_find = stage_timer.elapsed_time(); #endif - CUOPT_LOG_DEBUG("Number of cliques: %d, additional cliques: %d", - clique_table.first.size(), - clique_table.addtl_cliques.size()); - // print_clique_table(clique_table); - // remove small cliques and add them to adj_list - remove_small_cliques(clique_table); + CUOPT_LOG_DEBUG("Number of cliques: %d, additional cliques: %d, find_work=%.0f", + clique_table_ptr->first.size(), + clique_table_ptr->addtl_cliques.size(), + find_work_estimate); + remove_small_cliques(*clique_table_ptr, timer); #ifdef DEBUG_CLIQUE_TABLE t_small = stage_timer.elapsed_time(); #endif - // fill var clique maps - fill_var_clique_maps(clique_table); + fill_var_clique_maps(*clique_table_ptr); #ifdef DEBUG_CLIQUE_TABLE t_maps = stage_timer.elapsed_time(); #endif - extend_cliques( - knapsack_constraints, set_packing_constraints, clique_table, problem, A, modify_problem, timer); + if (clique_table_out != nullptr) { *clique_table_out = std::move(clique_table_shared); } + double extend_work = 0.0; + constexpr double max_extend_work = 2e9; + i_t n_extended_cliques = extend_cliques(knapsack_constraints, + set_packing_constraints, + *clique_table_ptr, + problem, + A, + modify_problem, + timer, + &extend_work, + max_extend_work); #ifdef DEBUG_CLIQUE_TABLE t_extend = stage_timer.elapsed_time(); - t_remove = t_extend; CUOPT_LOG_DEBUG( "Clique table timing (s): fill=%.6f coeff=%.6f sort=%.6f find=%.6f small=%.6f maps=%.6f " - "extend=%.6f remove=%.6f total=%.6f", + "extend=%.6f total=%.6f find_work=%.0f extend_work=%.0f", t_fill, t_coeff - t_fill, t_sort - t_coeff, @@ -1021,8 +1118,9 @@ void find_initial_cliques(dual_simplex::user_problem_t& problem, t_small - t_find, t_maps - t_small, t_extend - t_maps, - t_remove - t_extend, - t_remove); + t_extend, + find_work_estimate, + extend_work); #endif } @@ -1030,8 +1128,18 @@ void find_initial_cliques(dual_simplex::user_problem_t& problem, template void find_initial_cliques( \ dual_simplex::user_problem_t & problem, \ typename mip_solver_settings_t::tolerances_t tolerances, \ + std::shared_ptr> * clique_table_out, \ cuopt::timer_t & timer, \ - bool modify_problem); + bool modify_problem, \ + std::atomic* signal_extend); \ + template void build_clique_table( \ + const dual_simplex::user_problem_t& problem, \ + clique_table_t& clique_table, \ + typename mip_solver_settings_t::tolerances_t tolerances, \ + bool remove_small_cliques_flag, \ + bool fill_var_clique_maps_flag, \ + cuopt::timer_t& timer); \ + template class clique_table_t; #if MIP_INSTANTIATE_FLOAT INSTANTIATE(float) diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh index a8261685ab..944241b4f0 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh @@ -20,8 +20,10 @@ #include #include +#include #include +#include #include #include #include @@ -83,6 +85,8 @@ struct clique_table_t { std::vector> var_clique_map_first; // keeps the indices of additional cliques that contain variable x std::vector> var_clique_map_addtl; + // var_idx -> position mapping for each first clique, enabling O(1) membership/position checks + std::vector> first_var_positions; // adjacency list to keep small cliques, this basically keeps the vars share a small clique // constraint std::unordered_map> adj_list_small_cliques; @@ -98,8 +102,18 @@ struct clique_table_t { template void find_initial_cliques(dual_simplex::user_problem_t& problem, typename mip_solver_settings_t::tolerances_t tolerances, + std::shared_ptr>* clique_table_out, cuopt::timer_t& timer, - bool modify_problem); + bool modify_problem, + std::atomic* signal_extend = nullptr); + +template +void build_clique_table(const dual_simplex::user_problem_t& problem, + clique_table_t& clique_table, + typename mip_solver_settings_t::tolerances_t tolerances, + bool remove_small_cliques, + bool fill_var_clique_maps, + cuopt::timer_t& timer); } // namespace cuopt::linear_programming::detail diff --git a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp index 20a586f6fb..2fa10f421b 100644 --- a/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp +++ b/cpp/src/mip_heuristics/presolve/third_party_presolve.cpp @@ -590,7 +590,7 @@ std::optional> third_party_presolve_t( pslp_presolver_, op_problem.get_handle_ptr(), maximize_, original_obj_offset); - + opt_problem.set_problem_name(op_problem.get_problem_name()); return std::make_optional(third_party_presolve_result_t{opt_problem, {}}); } else { cuopt_expects( @@ -669,6 +669,7 @@ std::optional> third_party_presolve_t( papilo_problem, op_problem.get_handle_ptr(), category, maximize_); + opt_problem.set_problem_name(op_problem.get_problem_name()); auto col_flags = papilo_problem.getColFlags(); std::vector implied_integer_indices; for (size_t i = 0; i < col_flags.size(); i++) { diff --git a/cpp/src/mip_heuristics/problem/problem.cu b/cpp/src/mip_heuristics/problem/problem.cu index fadc00a850..90d80f5948 100644 --- a/cpp/src/mip_heuristics/problem/problem.cu +++ b/cpp/src/mip_heuristics/problem/problem.cu @@ -143,6 +143,7 @@ problem_t::problem_t( objective_offset(problem_.get_objective_offset()), lp_state(*this, problem_.get_handle_ptr()->get_stream()), fixing_helpers(n_constraints, n_variables, handle_ptr), + clique_table(nullptr), Q_offsets(problem_.get_quadratic_objective_offsets()), Q_indices(problem_.get_quadratic_objective_indices()), Q_values(problem_.get_quadratic_objective_values()) @@ -199,6 +200,7 @@ problem_t::problem_t(const problem_t& problem_) objective_is_integral(problem_.objective_is_integral), lp_state(problem_.lp_state), fixing_helpers(problem_.fixing_helpers, handle_ptr), + clique_table(problem_.clique_table), vars_with_objective_coeffs(problem_.vars_with_objective_coeffs), expensive_to_fix_vars(problem_.expensive_to_fix_vars), Q_offsets(problem_.Q_offsets), @@ -255,6 +257,7 @@ problem_t::problem_t(const problem_t& problem_, objective_is_integral(problem_.objective_is_integral), lp_state(problem_.lp_state, handle_ptr), fixing_helpers(problem_.fixing_helpers, handle_ptr), + clique_table(problem_.clique_table), vars_with_objective_coeffs(problem_.vars_with_objective_coeffs), expensive_to_fix_vars(problem_.expensive_to_fix_vars), Q_offsets(problem_.Q_offsets), @@ -279,6 +282,7 @@ problem_t::problem_t(const problem_t& problem_, bool no_deep maximize(problem_.maximize), empty(problem_.empty), is_binary_pb(problem_.is_binary_pb), + clique_table(problem_.clique_table), // Copy constructor used by PDLP and MIP // PDLP uses the version with no_deep_copy = false which deep copy some fields but doesn't // allocate others that are not needed in PDLP diff --git a/cpp/src/mip_heuristics/problem/problem.cuh b/cpp/src/mip_heuristics/problem/problem.cuh index b9ca420820..9771bab568 100644 --- a/cpp/src/mip_heuristics/problem/problem.cuh +++ b/cpp/src/mip_heuristics/problem/problem.cuh @@ -24,6 +24,7 @@ #include +#include #include #include #include @@ -37,6 +38,9 @@ namespace cuopt { namespace linear_programming::detail { +template +struct clique_table_t; + template class solution_t; @@ -119,6 +123,8 @@ class problem_t { bool is_integer(f_t val) const; bool integer_equal(f_t val1, f_t val2) const; + std::shared_ptr> clique_table; + void get_host_user_problem( cuopt::linear_programming::dual_simplex::user_problem_t& user_problem) const; void set_constraints_from_host_user_problem( diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index e6f6d50b62..da6c7c0b81 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -220,6 +220,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.mixed_integer_gomory_cuts = context.settings.mixed_integer_gomory_cuts; branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; + branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; branch_and_bound_settings.reduced_cost_strengthening = @@ -261,7 +262,10 @@ solution_t mip_solver_t::run_solver() // Create the branch and bound object branch_and_bound = std::make_unique>( - branch_and_bound_problem, branch_and_bound_settings, timer_.get_tic_start()); + branch_and_bound_problem, + branch_and_bound_settings, + timer_.get_tic_start(), + context.problem_ptr->clique_table); context.branch_and_bound_ptr = branch_and_bound.get(); auto* stats_ptr = &context.stats; branch_and_bound->set_user_bound_callback( diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 22fff31906..2fc9ec08d5 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -1571,9 +1571,8 @@ cuopt::linear_programming::optimization_problem_t mps_data_model_to_op if (data_model.get_objective_name().size() != 0) { op_problem.set_objective_name(data_model.get_objective_name()); } - if (data_model.get_problem_name().size() != 0) { - op_problem.set_problem_name(data_model.get_problem_name().data()); - } + auto problem_name = data_model.get_problem_name(); + op_problem.set_problem_name(problem_name); if (data_model.get_variable_names().size() != 0) { op_problem.set_variable_names(data_model.get_variable_names()); } diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 1a360b41eb..67baab5f12 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -8,25 +8,829 @@ #include "../linear_programming/utilities/pdlp_test_utilities.cuh" #include "mip_utils.cuh" +#include +#include #include +#include +#include +#include #include #include +#include #include +#include #include #include #include +#include #include #include #include +#include +#include +#include #include #include +#include #include namespace cuopt::linear_programming::test { +namespace { + +constexpr double kCliqueTestTol = 1e-6; + +mps_parser::mps_data_model_t create_pairwise_triangle_set_packing_problem() +{ + // Maximize x0 + x1 + x2 via minimizing -x0 - x1 - x2. + // Pairwise conflicts: + // x0 + x1 <= 1 + // x1 + x2 <= 1 + // x0 + x2 <= 1 + mps_parser::mps_data_model_t problem; + std::vector offsets = {0, 2, 4, 6}; + std::vector indices = {0, 1, 1, 2, 0, 2}; + std::vector coefficients = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + std::vector lower_bounds = {-std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity()}; + std::vector upper_bounds = {1.0, 1.0, 1.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + std::vector var_lower_bounds = {0.0, 0.0, 0.0}; + std::vector var_upper_bounds = {1.0, 1.0, 1.0}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + std::vector objective_coefficients = {-1.0, -1.0, -1.0}; + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + std::vector variable_types = {'I', 'I', 'I'}; + problem.set_variable_types(variable_types); + problem.set_maximize(false); + return problem; +} + +mps_parser::mps_data_model_t create_pairwise_triangle_with_isolated_variable_problem() +{ + // Same triangle conflicts as create_pairwise_triangle_set_packing_problem(), + // plus an isolated binary variable x3 with no conflict rows. + mps_parser::mps_data_model_t problem; + std::vector offsets = {0, 2, 4, 6}; + std::vector indices = {0, 1, 1, 2, 0, 2}; + std::vector coefficients = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + std::vector lower_bounds = {-std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), + -std::numeric_limits::infinity()}; + std::vector upper_bounds = {1.0, 1.0, 1.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + std::vector var_lower_bounds = {0.0, 0.0, 0.0, 0.0}; + std::vector var_upper_bounds = {1.0, 1.0, 1.0, 1.0}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + std::vector objective_coefficients = {-1.0, -1.0, -1.0, 0.0}; + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + std::vector variable_types = {'I', 'I', 'I', 'I'}; + problem.set_variable_types(variable_types); + problem.set_maximize(false); + return problem; +} + +mps_parser::mps_data_model_t create_binary_continuous_mixed_conflict_problem() +{ + // x0 + y1 <= 1 (must be ignored for clique graph because y1 is continuous) + // x0 + x2 <= 1 (must generate a conflict edge) + mps_parser::mps_data_model_t problem; + std::vector offsets = {0, 2, 4}; + std::vector indices = {0, 1, 0, 2}; + std::vector coefficients = {1.0, 1.0, 1.0, 1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + std::vector lower_bounds = {-std::numeric_limits::infinity(), + -std::numeric_limits::infinity()}; + std::vector upper_bounds = {1.0, 1.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + std::vector var_lower_bounds = {0.0, 0.0, 0.0}; + std::vector var_upper_bounds = {1.0, 1.0, 1.0}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + std::vector objective_coefficients = {0.0, 0.0, 0.0}; + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + std::vector variable_types = {'I', 'C', 'I'}; + problem.set_variable_types(variable_types); + problem.set_maximize(false); + return problem; +} + +mps_parser::mps_data_model_t create_near_binary_bound_conflict_problem() +{ + // x0 + x1 <= 1 but x1 has upper bound 0.9999999, so this row should not be + // treated as a binary conflict row. + mps_parser::mps_data_model_t problem; + std::vector offsets = {0, 2}; + std::vector indices = {0, 1}; + std::vector coefficients = {1.0, 1.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + std::vector lower_bounds = {-std::numeric_limits::infinity()}; + std::vector upper_bounds = {1.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + std::vector var_lower_bounds = {0.0, 0.0}; + std::vector var_upper_bounds = {1.0, 0.9999999}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + std::vector objective_coefficients = {0.0, 0.0}; + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + std::vector variable_types = {'I', 'I'}; + problem.set_variable_types(variable_types); + problem.set_maximize(false); + return problem; +} + +mps_parser::mps_data_model_t create_weighted_addtl_conflict_problem() +{ + // One weighted binary knapsack row: + // 1*x0 + 2*x1 + 3*x2 + 4*x3 <= 5 + // This creates base clique {x2, x3} and additional clique inducing conflict {x1, x3}. + mps_parser::mps_data_model_t problem; + std::vector offsets = {0, 4}; + std::vector indices = {0, 1, 2, 3}; + std::vector coefficients = {1.0, 2.0, 3.0, 4.0}; + problem.set_csr_constraint_matrix(coefficients.data(), + coefficients.size(), + indices.data(), + indices.size(), + offsets.data(), + offsets.size()); + std::vector lower_bounds = {-std::numeric_limits::infinity()}; + std::vector upper_bounds = {5.0}; + problem.set_constraint_lower_bounds(lower_bounds.data(), lower_bounds.size()); + problem.set_constraint_upper_bounds(upper_bounds.data(), upper_bounds.size()); + std::vector var_lower_bounds = {0.0, 0.0, 0.0, 0.0}; + std::vector var_upper_bounds = {1.0, 1.0, 1.0, 1.0}; + problem.set_variable_lower_bounds(var_lower_bounds.data(), var_lower_bounds.size()); + problem.set_variable_upper_bounds(var_upper_bounds.data(), var_upper_bounds.size()); + std::vector objective_coefficients = {0.0, 0.0, 0.0, 0.0}; + problem.set_objective_coefficients(objective_coefficients.data(), objective_coefficients.size()); + std::vector variable_types = {'I', 'I', 'I', 'I'}; + problem.set_variable_types(variable_types); + problem.set_maximize(false); + return problem; +} + +detail::clique_table_t build_clique_table_for_model_with_min_size( + const raft::handle_t& handle, + const mps_parser::mps_data_model_t& model, + int min_clique_size) +{ + auto op_problem = mps_data_model_to_optimization_problem(&handle, model); + detail::problem_t mip_problem(op_problem); + dual_simplex::user_problem_t host_problem(op_problem.get_handle_ptr()); + mip_problem.get_host_user_problem(host_problem); + + detail::clique_config_t clique_config; + clique_config.min_clique_size = min_clique_size; + detail::clique_table_t clique_table(2 * host_problem.num_cols, + clique_config.min_clique_size, + clique_config.max_clique_size_for_extension); + + mip_solver_settings_t settings; + cuopt::timer_t timer(std::numeric_limits::infinity()); + detail::build_clique_table(host_problem, clique_table, settings.tolerances, true, true, timer); + return clique_table; +} + +detail::clique_table_t build_clique_table_for_model( + const raft::handle_t& handle, const mps_parser::mps_data_model_t& model) +{ + return build_clique_table_for_model_with_min_size(handle, model, 1); +} + +mps_parser::mps_data_model_t& get_neos8_model_cached() +{ + static std::once_flag init_flag; + static std::unique_ptr> model_ptr; + std::call_once(init_flag, []() { + const auto neos8_path = make_path_absolute("mip/neos8.mps"); + auto neos8_model = cuopt::mps_parser::parse_mps(neos8_path, false); + model_ptr = std::make_unique>(std::move(neos8_model)); + }); + cuopt_assert(model_ptr != nullptr, "Failed to initialize cached neos8 model"); + return *model_ptr; +} + +detail::clique_table_t& get_neos8_clique_table_cached() +{ + static std::once_flag init_flag; + static std::unique_ptr> clique_table_ptr; + std::call_once(init_flag, []() { + const raft::handle_t handle{}; + auto& neos8_model = get_neos8_model_cached(); + auto clique_table = build_clique_table_for_model(handle, neos8_model); + clique_table_ptr = + std::make_unique>(std::move(clique_table)); + }); + cuopt_assert(clique_table_ptr != nullptr, "Failed to initialize cached neos8 clique table"); + return *clique_table_ptr; +} + +std::vector> build_original_adjacency_matrix( + detail::clique_table_t& clique_table, int num_vars) +{ + std::vector> adj(num_vars, std::vector(num_vars, 0)); + for (int i = 0; i < num_vars; ++i) { + for (int j = i + 1; j < num_vars; ++j) { + if (clique_table.check_adjacency(i, j)) { + adj[i][j] = 1; + adj[j][i] = 1; + } + } + } + return adj; +} + +std::vector> maximal_cliques_bruteforce(const std::vector>& adj) +{ + const int n = static_cast(adj.size()); + if (n <= 0 || n > 20) { return {}; } + const uint64_t total_masks = (uint64_t{1} << n); + std::vector> maximal_cliques; + + auto is_mask_clique = [&](uint64_t mask) { + for (int i = 0; i < n; ++i) { + if ((mask & (uint64_t{1} << i)) == 0) { continue; } + for (int j = i + 1; j < n; ++j) { + if ((mask & (uint64_t{1} << j)) == 0) { continue; } + if (!adj[i][j]) { return false; } + } + } + return true; + }; + + for (uint64_t mask = 1; mask < total_masks; ++mask) { + if (!is_mask_clique(mask)) { continue; } + bool is_maximal = true; + for (int v = 0; v < n && is_maximal; ++v) { + if (mask & (uint64_t{1} << v)) { continue; } + bool can_extend = true; + for (int u = 0; u < n; ++u) { + if ((mask & (uint64_t{1} << u)) == 0) { continue; } + if (!adj[v][u]) { + can_extend = false; + break; + } + } + if (can_extend) { is_maximal = false; } + } + if (!is_maximal) { continue; } + std::vector clique; + for (int u = 0; u < n; ++u) { + if (mask & (uint64_t{1} << u)) { clique.push_back(u); } + } + maximal_cliques.push_back(std::move(clique)); + } + return maximal_cliques; +} + +std::vector> canonicalize_cliques(std::vector> cliques) +{ + for (auto& clique : cliques) { + std::sort(clique.begin(), clique.end()); + } + std::sort(cliques.begin(), cliques.end(), [](const auto& a, const auto& b) { + if (a.size() != b.size()) { return a.size() < b.size(); } + return a < b; + }); + cliques.erase(std::unique(cliques.begin(), cliques.end()), cliques.end()); + return cliques; +} + +std::vector> adjacency_matrix_to_list(const std::vector>& adj) +{ + const int n = static_cast(adj.size()); + std::vector> adj_list(n); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + if (adj[i][j]) { adj_list[i].push_back(j); } + } + } + return adj_list; +} + +std::vector> maximal_cliques_from_production_algorithm( + const std::vector>& adj) +{ + const auto adj_list = adjacency_matrix_to_list(adj); + std::vector weights(adj_list.size(), 1.0); + auto cliques = dual_simplex::find_maximal_cliques_for_test( + adj_list, weights, 0.0, 100000, std::numeric_limits::infinity()); + return canonicalize_cliques(std::move(cliques)); +} + +double original_clique_sum(const std::vector& clique_vars, + const std::vector& assignment) +{ + double lhs = 0.0; + for (const auto var : clique_vars) { + lhs += assignment[var]; + } + return lhs; +} + +std::string format_phase2_panic_dump(const mps_parser::mps_data_model_t& problem, + const std::vector& clique_vars, + const std::vector& x_star) +{ + std::ostringstream out; + const auto& var_lb = problem.get_variable_lower_bounds(); + const auto& var_ub = problem.get_variable_upper_bounds(); + out << "\nClique vars:"; + for (auto v : clique_vars) { + out << " x" << v << "(value=" << x_star[v] << ", lb=" << var_lb[v] << ", ub=" << var_ub[v] + << ")"; + } + + std::unordered_set clique_var_set(clique_vars.begin(), clique_vars.end()); + const auto& values = problem.get_constraint_matrix_values(); + const auto& cols = problem.get_constraint_matrix_indices(); + const auto& rows = problem.get_constraint_matrix_offsets(); + const auto& clb = problem.get_constraint_lower_bounds(); + const auto& cub = problem.get_constraint_upper_bounds(); + + out << "\nRelated constraints:"; + for (size_t row = 0; row + 1 < rows.size(); ++row) { + bool touches_clique = false; + for (int p = rows[row]; p < rows[row + 1]; ++p) { + if (clique_var_set.count(cols[p]) > 0) { + touches_clique = true; + break; + } + } + if (!touches_clique) { continue; } + out << "\n row " << row << ": "; + for (int p = rows[row]; p < rows[row + 1]; ++p) { + if (p > rows[row]) { out << " + "; } + out << values[p] << "*x" << cols[p]; + } + out << " in [" << clb[row] << ", " << cub[row] << "]"; + } + return out.str(); +} + +void disable_non_clique_cuts(mip_solver_settings_t& settings) +{ + settings.clique_cuts = 1; + settings.max_cut_passes = 10; + settings.mixed_integer_gomory_cuts = 0; + settings.knapsack_cuts = 0; + settings.mir_cuts = 0; + settings.strong_chvatal_gomory_cuts = 0; +} + +void disable_all_cuts(mip_solver_settings_t& settings) +{ + settings.max_cut_passes = 0; + settings.clique_cuts = 0; + settings.mixed_integer_gomory_cuts = 0; + settings.knapsack_cuts = 0; + settings.mir_cuts = 0; + settings.strong_chvatal_gomory_cuts = 0; +} + +bool cut_is_invalid_for_incumbent(const std::vector& cut_vars, + const std::vector& incumbent, + double tol) +{ + return original_clique_sum(cut_vars, incumbent) > 1.0 + tol; +} + +bool prefix_has_invalid_cut(const std::vector>& dumped_cuts, + size_t prefix_end_exclusive, + const std::vector& incumbent, + double tol) +{ + for (size_t i = 0; i < prefix_end_exclusive; ++i) { + if (cut_is_invalid_for_incumbent(dumped_cuts[i], incumbent, tol)) { return true; } + } + return false; +} + +std::optional isolate_first_invalid_cut_by_bisection( + const std::vector>& dumped_cuts, + const std::vector& incumbent, + double tol) +{ + if (!prefix_has_invalid_cut(dumped_cuts, dumped_cuts.size(), incumbent, tol)) { + return std::nullopt; + } + size_t lo = 0; + size_t hi = dumped_cuts.size() - 1; + while (lo < hi) { + const size_t mid = lo + (hi - lo) / 2; + if (prefix_has_invalid_cut(dumped_cuts, mid + 1, incumbent, tol)) { + hi = mid; + } else { + lo = mid + 1; + } + } + return lo; +} + +struct neos8_mip_solution_cache_t { + mip_termination_status_t status; + std::vector primal; + double objective; +}; + +struct neos8_lp_solution_cache_t { + pdlp_termination_status_t status; + std::vector primal; +}; + +neos8_mip_solution_cache_t& get_neos8_optimal_solution_no_cuts_cached() +{ + static std::once_flag init_flag; + static std::unique_ptr solution_ptr; + std::call_once(init_flag, []() { + const raft::handle_t handle{}; + auto& neos8_model = get_neos8_model_cached(); + mip_solver_settings_t settings; + settings.time_limit = 120.0; + settings.presolver = presolver_t::None; + disable_all_cuts(settings); + + auto mip_solution = solve_mip(&handle, neos8_model, settings); + auto cache = std::make_unique(); + cache->status = mip_solution.get_termination_status(); + cache->objective = mip_solution.get_objective_value(); + cache->primal = cuopt::host_copy(mip_solution.get_solution(), handle.get_stream()); + solution_ptr = std::move(cache); + }); + cuopt_assert(solution_ptr != nullptr, "Failed to initialize cached neos8 no-cut MIP solution"); + return *solution_ptr; +} + +neos8_lp_solution_cache_t& get_neos8_lp_relaxation_solution_cached() +{ + static std::once_flag init_flag; + static std::unique_ptr solution_ptr; + std::call_once(init_flag, []() { + const raft::handle_t handle{}; + auto lp_relaxation = get_neos8_model_cached(); + std::vector all_continuous(lp_relaxation.get_n_variables(), 'C'); + lp_relaxation.set_variable_types(all_continuous); + + pdlp_solver_settings_t lp_settings{}; + lp_settings.time_limit = 120.0; + lp_settings.presolver = presolver_t::None; + lp_settings.set_optimality_tolerance(1e-8); + + auto lp_solution = solve_lp(&handle, lp_relaxation, lp_settings); + auto cache = std::make_unique(); + cache->status = lp_solution.get_termination_status(); + cache->primal = cuopt::host_copy(lp_solution.get_primal_solution(), handle.get_stream()); + solution_ptr = std::move(cache); + }); + cuopt_assert(solution_ptr != nullptr, "Failed to initialize cached neos8 LP relaxation solution"); + return *solution_ptr; +} + +bool is_binary_var_for_clique_literals(const mps_parser::mps_data_model_t& problem, + int var_idx, + double bound_tol) +{ + const auto& var_types = problem.get_variable_types(); + const auto& var_lb = problem.get_variable_lower_bounds(); + const auto& var_ub = problem.get_variable_upper_bounds(); + return var_types[var_idx] != 'C' && var_lb[var_idx] >= -bound_tol && + var_ub[var_idx] <= 1.0 + bound_tol; +} + +std::vector> build_fractional_literal_cliques_for_assignment( + const mps_parser::mps_data_model_t& problem, + detail::clique_table_t& clique_table, + const std::vector& assignment, + double integer_tol, + double bound_tol, + int max_calls) +{ + const int num_vars = problem.get_n_variables(); + cuopt_assert(static_cast(assignment.size()) >= num_vars, + "Assignment size mismatch in fractional literal clique builder"); + + std::vector vertices; + std::vector weights; + vertices.reserve(2 * num_vars); + weights.reserve(2 * num_vars); + for (int j = 0; j < num_vars; ++j) { + if (!is_binary_var_for_clique_literals(problem, j, bound_tol)) { continue; } + const double xj = assignment[j]; + if (std::abs(xj - std::round(xj)) <= integer_tol) { continue; } + vertices.push_back(j); + weights.push_back(xj); + vertices.push_back(j + num_vars); + weights.push_back(1.0 - xj); + } + if (vertices.empty()) { return {}; } + + std::vector vertex_to_local(2 * num_vars, -1); + std::vector in_subgraph(2 * num_vars, 0); + for (size_t idx = 0; idx < vertices.size(); ++idx) { + vertex_to_local[vertices[idx]] = static_cast(idx); + in_subgraph[vertices[idx]] = 1; + } + + std::vector> adj_local(vertices.size()); + for (size_t idx = 0; idx < vertices.size(); ++idx) { + const auto vertex_idx = vertices[idx]; + auto adj_set = clique_table.get_adj_set_of_var(vertex_idx); + auto& adj = adj_local[idx]; + adj.reserve(adj_set.size()); + for (const auto neighbor : adj_set) { + cuopt_assert(neighbor >= 0 && neighbor < 2 * num_vars, + "Neighbor out of range in fractional literal clique builder"); + if (!in_subgraph[neighbor]) { continue; } + const auto local_neighbor = vertex_to_local[neighbor]; + if (local_neighbor >= 0) { adj.push_back(local_neighbor); } + } + } + + auto cliques_local = dual_simplex::find_maximal_cliques_for_test( + adj_local, weights, 1.0 + kCliqueTestTol, max_calls, std::numeric_limits::infinity()); + std::vector> cliques_global; + cliques_global.reserve(cliques_local.size()); + for (auto& local_clique : cliques_local) { + std::vector global_clique; + global_clique.reserve(local_clique.size()); + for (const auto local_idx : local_clique) { + cuopt_assert(local_idx >= 0 && static_cast(local_idx) < vertices.size(), + "Local clique index out of range"); + global_clique.push_back(vertices[local_idx]); + } + cliques_global.push_back(std::move(global_clique)); + } + return canonicalize_cliques(std::move(cliques_global)); +} + +std::vector>& get_neos8_fractional_literal_cliques_cached() +{ + static std::once_flag init_flag; + static std::unique_ptr>> cliques_ptr; + std::call_once(init_flag, []() { + auto& neos8_model = get_neos8_model_cached(); + auto& clique_table = get_neos8_clique_table_cached(); + auto& lp_relaxation = get_neos8_lp_relaxation_solution_cached(); + auto cliques = build_fractional_literal_cliques_for_assignment( + neos8_model, clique_table, lp_relaxation.primal, kCliqueTestTol, kCliqueTestTol, 100000); + cliques_ptr = std::make_unique>>(std::move(cliques)); + }); + cuopt_assert(cliques_ptr != nullptr, "Failed to initialize cached neos8 dumped literal cliques"); + return *cliques_ptr; +} + +double literal_clique_cut_violation(const std::vector& literal_clique, + const std::vector& assignment, + int num_vars) +{ + cuopt_assert(static_cast(assignment.size()) >= num_vars, + "Assignment size mismatch in literal clique violation"); + double dot = 0.0; + int num_complement_vars = 0; + for (const auto literal : literal_clique) { + cuopt_assert(literal >= 0 && literal < 2 * num_vars, "Literal out of range"); + const int var_idx = literal % num_vars; + const bool is_complement = literal >= num_vars; + if (is_complement) { + num_complement_vars++; + dot += assignment[var_idx]; + } else { + dot -= assignment[var_idx]; + } + } + const double rhs = static_cast(num_complement_vars - 1); + return rhs - dot; +} + +std::string format_phase2_literal_panic_dump(const std::vector& literal_clique, + const std::vector& incumbent, + int num_vars) +{ + std::ostringstream out; + out << "\nLiteral clique:"; + for (const auto literal : literal_clique) { + const bool is_complement = literal >= num_vars; + const int var_idx = literal % num_vars; + out << " " << (is_complement ? "~x" : "x") << var_idx << "(value=" << incumbent[var_idx] << ")"; + } + out << "\nViolation: " << literal_clique_cut_violation(literal_clique, incumbent, num_vars); + return out.str(); +} + +bool literal_cut_is_invalid_for_incumbent(const std::vector& literal_clique, + const std::vector& incumbent, + int num_vars, + double tol) +{ + return literal_clique_cut_violation(literal_clique, incumbent, num_vars) > tol; +} + +bool prefix_has_invalid_literal_cut(const std::vector>& dumped_cuts, + size_t prefix_end_exclusive, + const std::vector& incumbent, + int num_vars, + double tol) +{ + for (size_t i = 0; i < prefix_end_exclusive; ++i) { + if (literal_cut_is_invalid_for_incumbent(dumped_cuts[i], incumbent, num_vars, tol)) { + return true; + } + } + return false; +} + +std::optional isolate_first_invalid_literal_cut_by_bisection( + const std::vector>& dumped_cuts, + const std::vector& incumbent, + int num_vars, + double tol) +{ + if (!prefix_has_invalid_literal_cut(dumped_cuts, dumped_cuts.size(), incumbent, num_vars, tol)) { + return std::nullopt; + } + size_t lo = 0; + size_t hi = dumped_cuts.size() - 1; + while (lo < hi) { + const size_t mid = lo + (hi - lo) / 2; + if (prefix_has_invalid_literal_cut(dumped_cuts, mid + 1, incumbent, num_vars, tol)) { + hi = mid; + } else { + lo = mid + 1; + } + } + return lo; +} + +mps_parser::mps_data_model_t& get_neos8_lp_relaxation_model_cached() +{ + static std::once_flag init_flag; + static std::unique_ptr> model_ptr; + std::call_once(init_flag, []() { + auto lp_relaxation = get_neos8_model_cached(); + std::vector all_continuous(lp_relaxation.get_n_variables(), 'C'); + lp_relaxation.set_variable_types(all_continuous); + model_ptr = + std::make_unique>(std::move(lp_relaxation)); + }); + cuopt_assert(model_ptr != nullptr, "Failed to initialize cached neos8 LP relaxation model"); + return *model_ptr; +} + +mps_parser::mps_data_model_t append_literal_cut_prefix_to_lp_model( + const mps_parser::mps_data_model_t& base_lp_model, + const std::vector>& dumped_cuts, + size_t prefix_end_exclusive, + int num_vars) +{ + auto model_with_cuts = base_lp_model; + if (prefix_end_exclusive == 0) { return model_with_cuts; } + + std::vector matrix_values = base_lp_model.get_constraint_matrix_values(); + std::vector matrix_indices = base_lp_model.get_constraint_matrix_indices(); + std::vector matrix_offsets = base_lp_model.get_constraint_matrix_offsets(); + std::vector constraint_lbs = base_lp_model.get_constraint_lower_bounds(); + std::vector constraint_ubs = base_lp_model.get_constraint_upper_bounds(); + std::vector row_names = base_lp_model.get_row_names(); + if (matrix_offsets.empty()) { matrix_offsets.push_back(0); } + + const size_t cuts_to_apply = std::min(prefix_end_exclusive, dumped_cuts.size()); + for (size_t cut_idx = 0; cut_idx < cuts_to_apply; ++cut_idx) { + const auto& literal_cut = dumped_cuts[cut_idx]; + + std::vector row_vars; + std::vector row_coeffs; + row_vars.reserve(literal_cut.size()); + row_coeffs.reserve(literal_cut.size()); + + int num_complements = 0; + for (const auto literal : literal_cut) { + cuopt_assert(literal >= 0 && literal < 2 * num_vars, + "Literal out of range for LP cut append"); + const int var_idx = literal % num_vars; + const bool is_complement = literal >= num_vars; + if (is_complement) { num_complements++; } + const double coeff = is_complement ? 1.0 : -1.0; + + bool found = false; + for (size_t t = 0; t < row_vars.size(); ++t) { + if (row_vars[t] == var_idx) { + row_coeffs[t] += coeff; + found = true; + break; + } + } + if (!found) { + row_vars.push_back(var_idx); + row_coeffs.push_back(coeff); + } + } + + std::vector order(row_vars.size()); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](int a, int b) { return row_vars[a] < row_vars[b]; }); + for (const auto pos : order) { + const double coeff = row_coeffs[pos]; + if (std::abs(coeff) <= 1e-12) { continue; } + matrix_indices.push_back(row_vars[pos]); + matrix_values.push_back(coeff); + } + matrix_offsets.push_back(static_cast(matrix_indices.size())); + constraint_lbs.push_back(static_cast(num_complements - 1)); + constraint_ubs.push_back(std::numeric_limits::infinity()); + row_names.push_back("literal_cut_" + std::to_string(cut_idx)); + } + + model_with_cuts.set_csr_constraint_matrix(matrix_values.data(), + matrix_values.size(), + matrix_indices.data(), + matrix_indices.size(), + matrix_offsets.data(), + matrix_offsets.size()); + model_with_cuts.set_constraint_lower_bounds(constraint_lbs.data(), constraint_lbs.size()); + model_with_cuts.set_constraint_upper_bounds(constraint_ubs.data(), constraint_ubs.size()); + model_with_cuts.set_row_names(row_names); + return model_with_cuts; +} + +pdlp_termination_status_t solve_lp_with_literal_cut_prefix( + const std::vector>& dumped_cuts, size_t prefix_end_exclusive, int num_vars) +{ + const raft::handle_t handle{}; + auto& base_lp_model = get_neos8_lp_relaxation_model_cached(); + auto model_with_cuts = append_literal_cut_prefix_to_lp_model( + base_lp_model, dumped_cuts, prefix_end_exclusive, num_vars); + + pdlp_solver_settings_t lp_settings{}; + lp_settings.time_limit = 120.0; + lp_settings.presolver = presolver_t::None; + lp_settings.set_optimality_tolerance(1e-8); + + auto lp_solution = solve_lp(&handle, model_with_cuts, lp_settings); + return lp_solution.get_termination_status(); +} + +bool prefix_makes_lp_relaxation_infeasible(const std::vector>& dumped_cuts, + size_t prefix_end_exclusive, + int num_vars) +{ + const auto status = solve_lp_with_literal_cut_prefix(dumped_cuts, prefix_end_exclusive, num_vars); + return status == pdlp_termination_status_t::PrimalInfeasible; +} + +std::optional isolate_first_lp_infeasible_literal_cut_by_bisection( + const std::vector>& dumped_cuts, int num_vars) +{ + if (!prefix_makes_lp_relaxation_infeasible(dumped_cuts, dumped_cuts.size(), num_vars)) { + return std::nullopt; + } + size_t lo = 0; + size_t hi = dumped_cuts.size() - 1; + while (lo < hi) { + const size_t mid = lo + (hi - lo) / 2; + if (prefix_makes_lp_relaxation_infeasible(dumped_cuts, mid + 1, num_vars)) { + hi = mid; + } else { + lo = mid + 1; + } + } + return lo; +} + +} // namespace + // Problem data for the mixed integer linear programming problem mps_parser::mps_data_model_t create_cuts_problem_1() { @@ -165,4 +969,416 @@ TEST(cuts, test_cuts_2) EXPECT_EQ(solution.get_num_nodes(), 0); } +TEST(cuts, clique_phase1_smoke_conflict_graph_edges) +{ + const raft::handle_t handle{}; + auto problem = create_pairwise_triangle_with_isolated_variable_problem(); + auto clique_table = build_clique_table_for_model(handle, problem); + + // Positive edges from triangle. + EXPECT_TRUE(clique_table.check_adjacency(0, 1)); + EXPECT_TRUE(clique_table.check_adjacency(1, 0)); + EXPECT_TRUE(clique_table.check_adjacency(1, 2)); + EXPECT_TRUE(clique_table.check_adjacency(2, 1)); + EXPECT_TRUE(clique_table.check_adjacency(0, 2)); + EXPECT_TRUE(clique_table.check_adjacency(2, 0)); + + // Negative edges to isolated x3. + EXPECT_FALSE(clique_table.check_adjacency(0, 3)); + EXPECT_FALSE(clique_table.check_adjacency(3, 0)); + EXPECT_FALSE(clique_table.check_adjacency(1, 3)); + EXPECT_FALSE(clique_table.check_adjacency(3, 1)); + EXPECT_FALSE(clique_table.check_adjacency(2, 3)); + EXPECT_FALSE(clique_table.check_adjacency(3, 2)); + + // Self is never an edge. + EXPECT_FALSE(clique_table.check_adjacency(3, 3)); +} + +TEST(cuts, clique_phase1_unit_maximal_clique_finder_hardcoded_adj) +{ + // Hardcoded graph: + // triangle (0,1,2) and an extra edge (2,3) + std::vector> adj = { + {0, 1, 1, 0}, + {1, 0, 1, 0}, + {1, 1, 0, 1}, + {0, 0, 1, 0}, + }; + + auto maximal_bruteforce = canonicalize_cliques(maximal_cliques_bruteforce(adj)); + auto maximal_internal = maximal_cliques_from_production_algorithm(adj); + EXPECT_EQ(maximal_internal, maximal_bruteforce); + bool found_triangle = false; + for (const auto& clique : maximal_internal) { + if (clique.size() == 3 && clique[0] == 0 && clique[1] == 1 && clique[2] == 2) { + found_triangle = true; + break; + } + } + EXPECT_TRUE(found_triangle); +} + +TEST(cuts, clique_phase1_addtl_conflict_symmetry_and_reverse_lookup) +{ + const raft::handle_t handle{}; + auto problem = create_weighted_addtl_conflict_problem(); + auto clique_table = build_clique_table_for_model_with_min_size(handle, problem, 1); + + ASSERT_FALSE(clique_table.addtl_cliques.empty()); + + // Conflict introduced through additional clique path must be symmetric. + EXPECT_TRUE(clique_table.check_adjacency(1, 3)); + EXPECT_TRUE(clique_table.check_adjacency(3, 1)); + + // get_adj_set_of_var() must also include reverse lookup for addtl membership. + auto adj_of_1 = clique_table.get_adj_set_of_var(1); + auto adj_of_3 = clique_table.get_adj_set_of_var(3); + EXPECT_TRUE(adj_of_1.count(3) > 0); + EXPECT_TRUE(adj_of_3.count(1) > 0); +} + +TEST(cuts, clique_phase1_remove_small_cliques_preserves_addtl_conflicts) +{ + const raft::handle_t handle{}; + auto problem = create_weighted_addtl_conflict_problem(); + // Force base clique {x2,x3} to be considered "small" and removed. + auto clique_table = build_clique_table_for_model_with_min_size(handle, problem, 2); + + EXPECT_TRUE(clique_table.first.empty()); + EXPECT_TRUE(clique_table.addtl_cliques.empty()); + + // Conflicts must remain materialized in adj_list_small_cliques after removals. + EXPECT_TRUE(clique_table.check_adjacency(1, 3)); + EXPECT_TRUE(clique_table.check_adjacency(3, 1)); + EXPECT_TRUE(clique_table.check_adjacency(2, 3)); + EXPECT_TRUE(clique_table.check_adjacency(3, 2)); + EXPECT_FALSE(clique_table.check_adjacency(0, 3)); +} + +TEST(cuts, clique_phase2_no_cut_off_optimal_solution_validation) +{ + const raft::handle_t handle{}; + auto problem = create_pairwise_triangle_set_packing_problem(); + + mip_solver_settings_t settings; + settings.time_limit = 10.0; + settings.presolver = presolver_t::None; + disable_all_cuts(settings); + + auto mip_solution = solve_mip(&handle, problem, settings); + ASSERT_EQ(mip_solution.get_termination_status(), mip_termination_status_t::Optimal); + auto x_star = cuopt::host_copy(mip_solution.get_solution(), handle.get_stream()); + + auto clique_table = build_clique_table_for_model(handle, problem); + auto adj = build_original_adjacency_matrix(clique_table, problem.get_n_variables()); + auto maximal = maximal_cliques_bruteforce(adj); + ASSERT_FALSE(maximal.empty()); + + for (const auto& clique_vars : maximal) { + if (clique_vars.size() < 2) { continue; } + const double lhs = original_clique_sum(clique_vars, x_star); + ASSERT_LE(lhs, 1.0 + kCliqueTestTol) << format_phase2_panic_dump(problem, clique_vars, x_star); + } +} + +TEST(cuts, clique_phase3_fractional_separation_must_cut_off) +{ + const raft::handle_t handle{}; + auto mip_problem = create_pairwise_triangle_set_packing_problem(); + + auto lp_relaxation = mip_problem; + std::vector all_continuous(lp_relaxation.get_n_variables(), 'C'); + lp_relaxation.set_variable_types(all_continuous); + + pdlp_solver_settings_t lp_settings{}; + lp_settings.time_limit = 10.0; + lp_settings.presolver = presolver_t::None; + lp_settings.set_optimality_tolerance(1e-8); + + auto lp_solution = solve_lp(&handle, lp_relaxation, lp_settings); + ASSERT_EQ(lp_solution.get_termination_status(), pdlp_termination_status_t::Optimal); + auto x_bar = cuopt::host_copy(lp_solution.get_primal_solution(), handle.get_stream()); + + auto clique_table = build_clique_table_for_model(handle, mip_problem); + auto adj = build_original_adjacency_matrix(clique_table, mip_problem.get_n_variables()); + auto maximal = maximal_cliques_from_production_algorithm(adj); + + bool found_separating_clique = false; + for (const auto& clique_vars : maximal) { + if (clique_vars.size() < 2) { continue; } + const double lhs = original_clique_sum(clique_vars, x_bar); + if (lhs > 1.0 + kCliqueTestTol) { + found_separating_clique = true; + break; + } + } + EXPECT_TRUE(found_separating_clique); +} + +TEST(cuts, clique_phase4_fault_isolation_binary_search) +{ + // Simulated incumbent x* and dumped cuts. + // First invalid cut is at index 2: {0,1} gives 2 > 1. + const std::vector incumbent = {1.0, 1.0, 0.0, 0.0}; + const std::vector> dumped_cuts = { + {0, 2}, // valid + {1, 3}, // valid + {0, 1}, // invalid + {2, 3}, // valid + }; + + auto first_invalid = + isolate_first_invalid_cut_by_bisection(dumped_cuts, incumbent, kCliqueTestTol); + ASSERT_TRUE(first_invalid.has_value()); + EXPECT_EQ(first_invalid.value(), 2); +} + +TEST(cuts, clique_phase4_tree_depth_limit_smoke) +{ + const raft::handle_t handle{}; + auto problem = create_pairwise_triangle_set_packing_problem(); + + mip_solver_settings_t root_only_settings; + root_only_settings.time_limit = 10.0; + root_only_settings.presolver = presolver_t::None; + root_only_settings.node_limit = 0; + disable_non_clique_cuts(root_only_settings); + + mip_solver_settings_t deeper_settings = root_only_settings; + deeper_settings.node_limit = 100; + + auto root_only_solution = solve_mip(&handle, problem, root_only_settings); + auto deeper_solution = solve_mip(&handle, problem, deeper_settings); + + EXPECT_EQ(deeper_solution.get_termination_status(), mip_termination_status_t::Optimal); + EXPECT_NE(root_only_solution.get_termination_status(), mip_termination_status_t::Infeasible); + if (root_only_solution.get_termination_status() == mip_termination_status_t::Optimal) { + EXPECT_NEAR( + root_only_solution.get_objective_value(), deeper_solution.get_objective_value(), 1e-6); + } +} + +TEST(cuts, clique_phase5_ignores_non_binary_variables) +{ + const raft::handle_t handle{}; + auto problem = create_binary_continuous_mixed_conflict_problem(); + auto clique_table = build_clique_table_for_model(handle, problem); + + EXPECT_TRUE(clique_table.check_adjacency(0, 2)); + EXPECT_FALSE(clique_table.check_adjacency(0, 1)); + EXPECT_FALSE(clique_table.check_adjacency(1, 2)); +} + +TEST(cuts, clique_phase5_ignores_fractional_binary_bounds) +{ + const raft::handle_t handle{}; + auto problem = create_near_binary_bound_conflict_problem(); + auto clique_table = build_clique_table_for_model(handle, problem); + + EXPECT_FALSE(clique_table.check_adjacency(0, 1)); +} + +TEST(cuts, clique_neos8_phase1_addtl_indices_and_nonempty_graph) +{ + auto& clique_table = get_neos8_clique_table_cached(); + EXPECT_TRUE(!clique_table.first.empty() || !clique_table.addtl_cliques.empty()); + + const size_t max_addtl_to_check = std::min(clique_table.addtl_cliques.size(), 400); + for (size_t k = 0; k < max_addtl_to_check; ++k) { + const auto& addtl = clique_table.addtl_cliques[k]; + ASSERT_GE(addtl.clique_idx, 0); + ASSERT_LT(static_cast(addtl.clique_idx), clique_table.first.size()); + const auto& base = clique_table.first[addtl.clique_idx]; + ASSERT_GE(addtl.start_pos_on_clique, 0); + ASSERT_LE(static_cast(addtl.start_pos_on_clique), base.size()); + } +} + +TEST(cuts, clique_neos8_phase1_addtl_suffix_conflicts_materialized) +{ + auto& clique_table = get_neos8_clique_table_cached(); + if (clique_table.addtl_cliques.empty()) { + GTEST_SKIP() << "neos8 produced no additional cliques in this configuration"; + } + + size_t checked_addtl = 0; + const size_t max_addtl_to_check = std::min(clique_table.addtl_cliques.size(), 200); + for (size_t k = 0; k < max_addtl_to_check; ++k) { + const auto& addtl = clique_table.addtl_cliques[k]; + if (addtl.clique_idx < 0 || + static_cast(addtl.clique_idx) >= clique_table.first.size()) { + continue; + } + const auto& base = clique_table.first[addtl.clique_idx]; + const size_t start_at = static_cast(addtl.start_pos_on_clique); + if (start_at >= base.size()) { continue; } + + const size_t end_at = std::min(base.size(), start_at + 8); + for (size_t p = start_at; p < end_at; ++p) { + EXPECT_TRUE(clique_table.check_adjacency(addtl.vertex_idx, base[p])); + EXPECT_TRUE(clique_table.check_adjacency(base[p], addtl.vertex_idx)); + } + checked_addtl++; + } + EXPECT_GT(checked_addtl, 0); +} + +TEST(cuts, clique_neos8_phase1_symmetry_and_degree_cache_consistency) +{ + auto& clique_table = get_neos8_clique_table_cached(); + const int n_vertices = static_cast(clique_table.var_clique_map_first.size()); + ASSERT_GT(n_vertices, 0); + + const int sample_size = std::min(n_vertices, 24); + const int stride = std::max(1, n_vertices / sample_size); + std::vector sampled_vertices(sample_size); + for (int i = 0; i < sample_size; ++i) { + sampled_vertices[i] = (i * stride) % n_vertices; + } + + for (const auto v : sampled_vertices) { + const auto deg_cached = clique_table.get_degree_of_var(v); + const auto adj_set = clique_table.get_adj_set_of_var(v); + EXPECT_EQ(deg_cached, static_cast(adj_set.size())); + EXPECT_EQ(deg_cached, clique_table.get_degree_of_var(v)); + } + + for (int i = 0; i < sample_size; ++i) { + for (int j = i + 1; j < sample_size; ++j) { + const auto v1 = sampled_vertices[i]; + const auto v2 = sampled_vertices[j]; + EXPECT_EQ(clique_table.check_adjacency(v1, v2), clique_table.check_adjacency(v2, v1)); + } + } +} + +TEST(cuts, clique_neos8_phase2_no_cut_off_optimal_solution_validation) +{ + auto& no_cut_mip = get_neos8_optimal_solution_no_cuts_cached(); + ASSERT_EQ(no_cut_mip.status, mip_termination_status_t::Optimal); + + auto& lp_relaxation = get_neos8_lp_relaxation_solution_cached(); + ASSERT_EQ(lp_relaxation.status, pdlp_termination_status_t::Optimal); + + auto& dumped_literal_cuts = get_neos8_fractional_literal_cliques_cached(); + if (dumped_literal_cuts.empty()) { + GTEST_SKIP() << "neos8 produced no candidate literal cliques from LP relaxation"; + } + + const int num_vars = get_neos8_model_cached().get_n_variables(); + for (size_t i = 0; i < dumped_literal_cuts.size(); ++i) { + const double violation = + literal_clique_cut_violation(dumped_literal_cuts[i], no_cut_mip.primal, num_vars); + ASSERT_LE(violation, kCliqueTestTol) + << "Invalid clique cut at index " << i + << format_phase2_literal_panic_dump(dumped_literal_cuts[i], no_cut_mip.primal, num_vars); + } +} + +TEST(cuts, clique_neos8_phase3_fractional_separation_must_cut_off) +{ + auto& lp_relaxation = get_neos8_lp_relaxation_solution_cached(); + ASSERT_EQ(lp_relaxation.status, pdlp_termination_status_t::Optimal); + + auto& dumped_literal_cuts = get_neos8_fractional_literal_cliques_cached(); + if (dumped_literal_cuts.empty()) { + GTEST_SKIP() << "neos8 produced no candidate literal cliques from LP relaxation"; + } + + const int num_vars = get_neos8_model_cached().get_n_variables(); + for (size_t i = 0; i < dumped_literal_cuts.size(); ++i) { + const double violation = + literal_clique_cut_violation(dumped_literal_cuts[i], lp_relaxation.primal, num_vars); + ASSERT_GT(violation, kCliqueTestTol) + << "Non-separating clique cut at index " << i + << format_phase2_literal_panic_dump(dumped_literal_cuts[i], lp_relaxation.primal, num_vars); + } +} + +TEST(cuts, clique_neos8_phase4_fault_isolation_binary_search) +{ + auto& no_cut_mip = get_neos8_optimal_solution_no_cuts_cached(); + ASSERT_EQ(no_cut_mip.status, mip_termination_status_t::Optimal); + + auto& dumped_literal_cuts = get_neos8_fractional_literal_cliques_cached(); + if (dumped_literal_cuts.empty()) { + GTEST_SKIP() << "neos8 produced no candidate literal cliques from LP relaxation"; + } + + const auto& model = get_neos8_model_cached(); + const int num_vars = model.get_n_variables(); + + // Real dumped cuts should not invalidate the no-cut incumbent. + EXPECT_FALSE(prefix_has_invalid_literal_cut( + dumped_literal_cuts, dumped_literal_cuts.size(), no_cut_mip.primal, num_vars, kCliqueTestTol)); + + // Inject a known-invalid cut and verify bisection isolates it. + std::vector incumbent_ones; + incumbent_ones.reserve(2); + for (int j = 0; j < num_vars && incumbent_ones.size() < 2; ++j) { + if (!is_binary_var_for_clique_literals(model, j, kCliqueTestTol)) { continue; } + if (no_cut_mip.primal[j] >= 1.0 - kCliqueTestTol) { incumbent_ones.push_back(j); } + } + if (incumbent_ones.size() < 2) { + GTEST_SKIP() << "Could not find two binary variables fixed to one in neos8 incumbent"; + } + + auto cuts_with_injected_bug = dumped_literal_cuts; + const size_t injected_index = cuts_with_injected_bug.size(); + cuts_with_injected_bug.push_back({incumbent_ones[0], incumbent_ones[1]}); + + auto first_invalid = isolate_first_invalid_literal_cut_by_bisection( + cuts_with_injected_bug, no_cut_mip.primal, num_vars, kCliqueTestTol); + ASSERT_TRUE(first_invalid.has_value()); + EXPECT_EQ(first_invalid.value(), injected_index); +} + +TEST(cuts, clique_neos8_phase4_lp_infeasibility_binary_search) +{ + auto& dumped_literal_cuts = get_neos8_fractional_literal_cliques_cached(); + if (dumped_literal_cuts.empty()) { + GTEST_SKIP() << "neos8 produced no candidate literal cliques from LP relaxation"; + } + + const auto& model = get_neos8_model_cached(); + const int num_vars = model.get_n_variables(); + + std::vector> cuts_for_lp_search; + const size_t max_real_cuts = std::min(dumped_literal_cuts.size(), 64); + cuts_for_lp_search.insert(cuts_for_lp_search.end(), + dumped_literal_cuts.begin(), + dumped_literal_cuts.begin() + max_real_cuts); + + int inject_var = -1; + for (int j = 0; j < num_vars; ++j) { + if (is_binary_var_for_clique_literals(model, j, kCliqueTestTol)) { + inject_var = j; + break; + } + } + if (inject_var < 0) { + GTEST_SKIP() << "Could not find a binary variable for LP infeasibility injection"; + } + + const size_t injected_index = cuts_for_lp_search.size(); + cuts_for_lp_search.push_back( + {inject_var, inject_var, inject_var + num_vars, inject_var + num_vars}); + + // Prefix before injected cut should remain LP-feasible. + const auto status_before_injection = + solve_lp_with_literal_cut_prefix(cuts_for_lp_search, injected_index, num_vars); + EXPECT_NE(status_before_injection, pdlp_termination_status_t::PrimalInfeasible); + + // Full prefix should be LP-infeasible due to injected contradictory cut. + const auto status_with_injection = + solve_lp_with_literal_cut_prefix(cuts_for_lp_search, cuts_for_lp_search.size(), num_vars); + EXPECT_EQ(status_with_injection, pdlp_termination_status_t::PrimalInfeasible); + + auto first_infeasible = + isolate_first_lp_infeasible_literal_cut_by_bisection(cuts_for_lp_search, num_vars); + ASSERT_TRUE(first_infeasible.has_value()); + EXPECT_EQ(first_infeasible.value(), injected_index); +} + } // namespace cuopt::linear_programming::test diff --git a/datasets/mip/download_miplib_test_dataset.sh b/datasets/mip/download_miplib_test_dataset.sh index dc2dd79662..3040f0f543 100755 --- a/datasets/mip/download_miplib_test_dataset.sh +++ b/datasets/mip/download_miplib_test_dataset.sh @@ -20,6 +20,7 @@ INSTANCES=( "thor50dday" "stein9inf" "neos5" + "neos8" "swath1" "enlight_hard" "enlight11"