diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 2eb877b3e9..a7c341f555 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -48,7 +48,9 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump_kernels.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/fj_cpu.cu ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_cpufj.cu - ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_gpufj.cu) + ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_gpufj.cu + ${CMAKE_CURRENT_SOURCE_DIR}/structural/early_structural.cu + ${CMAKE_CURRENT_SOURCE_DIR}/structural/arc_flow.cu) # Choose which files to include based on build mode if(BUILD_LP_ONLY) diff --git a/cpp/src/mip_heuristics/early_heuristic.cuh b/cpp/src/mip_heuristics/early_heuristic.cuh index 6654470732..cb0be4200a 100644 --- a/cpp/src/mip_heuristics/early_heuristic.cuh +++ b/cpp/src/mip_heuristics/early_heuristic.cuh @@ -74,7 +74,9 @@ class early_heuristic_t { // NOT thread-safe. solver_obj is in solver-space (always minimization). // Uses a private CUDA stream to avoid racing with the FJ solver's stream. - void try_update_best(f_t solver_obj, const std::vector& assignment) + void try_update_best(f_t solver_obj, + const std::vector& assignment, + const char* heuristic_name = Derived::name()) { if (solver_obj >= best_objective_) { return; } best_objective_ = solver_obj; @@ -92,7 +94,7 @@ class early_heuristic_t { // Log and callback are deferred to the shared incumbent_callback_ which enforces // global monotonicity across all early heuristic instances. if (incumbent_callback_) { - incumbent_callback_(solver_obj, user_obj, user_assignment, Derived::name()); + incumbent_callback_(solver_obj, user_obj, user_assignment, heuristic_name); } } diff --git a/cpp/src/mip_heuristics/mip_constants.hpp b/cpp/src/mip_heuristics/mip_constants.hpp index f3fb68343a..58a29c5182 100644 --- a/cpp/src/mip_heuristics/mip_constants.hpp +++ b/cpp/src/mip_heuristics/mip_constants.hpp @@ -15,11 +15,13 @@ #define PDLP_INSTANTIATE_FLOAT 1 /* @brief Minimimum number of threads to enable each part of the MIP Solver */ -#define CUOPT_MIP_FJ_REQUIRED_THREAD_COUNT 8 -#define CUOPT_MIP_EARLY_GPUFJ_REQUIRED_THREAD_COUNT 3 -#define CUOPT_MIP_EARLY_CPUFJ_REQUIRED_THREAD_COUNT 2 -#define CUOPT_MIP_BATCH_PDLP_REQUIRED_THREAD_COUNT 3 -#define CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT 3 +#define CUOPT_MIP_FJ_REQUIRED_THREAD_COUNT 8 +#define CUOPT_MIP_EARLY_GPUFJ_REQUIRED_THREAD_COUNT 3 +#define CUOPT_MIP_EARLY_CPUFJ_REQUIRED_THREAD_COUNT 2 +#define CUOPT_MIP_EARLY_STRUCTURAL_REQUIRED_THREAD_COUNT 2 +#define CUOPT_MIP_ROOT_STRUCTURAL_REQUIRED_THREAD_COUNT 3 +#define CUOPT_MIP_BATCH_PDLP_REQUIRED_THREAD_COUNT 3 +#define CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT 3 // MIP-only gate: skip the concurrent barrier when fewer threads are available than this // (1 PDLP + 1 dual simplex + 1 barrier). Stand-alone LP always runs all three. diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 162a5ba291..378ed08db3 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -65,6 +66,7 @@ #include #include +#include #include namespace cuopt::mathematical_optimization { @@ -249,7 +251,13 @@ mip_solution_t run_mip_solver( // optimization_problem_t). Its solver-space differs from both the first-pass FJ (original // problem) and B&B (post-trivial- presolve), so initial_upper_bound (user-space) is converted // via problem.get_solver_obj_from_user_obj. + + // Must outlive early_cpufj/early_structural below, whose destructors join the tasks. + std::mutex papilo_callback_mutex; + f_t papilo_best_solver_obj = std::numeric_limits::infinity(); + std::unique_ptr> early_cpufj; + std::unique_ptr> early_structural; bool run_early_cpufj = problem.has_papilo_presolve_data() && settings.determinism_mode != CUOPT_MODE_DETERMINISTIC && problem.original_problem_ptr->get_n_integers() > 0; @@ -268,13 +276,21 @@ mip_solution_t run_mip_solver( semi_continuous_original_num_variables = mip_solver_settings_accessor::get_semi_continuous_original_num_variables( settings), - ctx_ptr = &solver.context, + ctx_ptr = &solver.context, + papilo_num_original_vars = problem.get_papilo_original_num_variables(), + &papilo_callback_mutex, + &papilo_best_solver_obj, early_fj_start](f_t solver_obj, f_t user_obj, const std::vector& assignment, const char* heuristic_name) { + std::lock_guard lock(papilo_callback_mutex); + if (solver_obj >= papilo_best_solver_obj) { return; } + papilo_best_solver_obj = solver_obj; + std::vector user_assignment; presolver_ptr->uncrush_primal_solution(assignment, user_assignment); + cuopt_assert(user_assignment.size() == (size_t)papilo_num_original_vars, "Size mismatch"); ctx_ptr->initial_incumbent_assignment = user_assignment; ctx_ptr->initial_upper_bound = user_obj; double elapsed = @@ -302,6 +318,17 @@ mip_solution_t run_mip_solver( early_cpufj->start(); solver.context.early_cpufj_ptr = early_cpufj.get(); CUOPT_LOG_DEBUG("Started early CPUFJ on papilo-presolved problem during cuOpt presolve"); + + early_structural = mip::early_structural_t::create( + *problem.original_problem_ptr, settings.get_tolerances(), incumbent_callback); + if (early_structural) { + if (std::isfinite(initial_upper_bound)) { + early_structural->set_best_objective( + problem.get_solver_obj_from_user_obj(initial_upper_bound)); + } + early_structural->start(); + solver.context.early_structural_ptr = early_structural.get(); + } } auto presolved_sol = solver.run_solver(); @@ -498,6 +525,7 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p std::unique_ptr> early_cpufj; std::unique_ptr> early_gpufj; + std::unique_ptr> early_structural; bool run_early_fj = run_presolve && settings.determinism_mode != CUOPT_MODE_DETERMINISTIC && op_problem.get_n_integers() > 0 && op_problem.get_n_constraints() > 0; @@ -556,6 +584,9 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p std::make_unique>(op_problem, settings, early_fj_callback); early_gpufj->start(); CUOPT_LOG_DEBUG("Started early GPUFJ during presolve"); + early_structural = mip::early_structural_t::create( + op_problem, settings.get_tolerances(), early_fj_callback); + if (early_structural) { early_structural->start(); } } auto constexpr const dual_postsolve = false; @@ -650,6 +681,17 @@ mip_solution_t solve_mip_helper(optimization_problem_t& op_p early_cpufj.reset(); } + if (early_structural) { + early_structural->stop(); + if (early_structural->solution_found()) { + CUOPT_LOG_DEBUG( + "Early structural heuristic (original) found incumbent with objective %.6e " + "during presolve", + early_structural->get_best_objective()); + } + early_structural.reset(); + } + // Add early-heuristic incumbents (original-space) to initial_solutions. // PaPILO crushing + validation happens downstream in add_user_given_solutions(). if (!early_incumbent_pool.empty()) { diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index f8eac0c4d8..22b4672496 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -25,6 +25,7 @@ #include #include +#include #include #include @@ -227,6 +228,16 @@ solution_t mip_solver_t::run_solver() } } + if (context.early_structural_ptr) { + context.early_structural_ptr->stop(); + if (context.early_structural_ptr->solution_found()) { + CUOPT_LOG_DEBUG( + "Early structural heuristic found incumbent with user-space objective %g " + "during presolve", + context.early_structural_ptr->get_best_user_objective()); + } + } + if (!presolve_success) { CUOPT_LOG_INFO("Problem proven infeasible in presolve"); sol.set_problem_fully_reduced(); @@ -480,6 +491,20 @@ solution_t mip_solver_t::run_solver() } } + std::unique_ptr> root_structural; + if (num_threads >= CUOPT_MIP_ROOT_STRUCTURAL_REQUIRED_THREAD_COUNT && + context.settings.determinism_mode != CUOPT_MODE_DETERMINISTIC && + !context.settings.heuristics_only) { + root_structural = std::make_unique>( + *context.problem_ptr, + context.settings.get_tolerances(), + context.preempt_heuristic_solver_, + [&dm](const std::vector& assignment, f_t objective) { + dm.population.add_external_solution(assignment, objective, solution_origin_t::EXTERNAL); + }); + if (!root_structural->recognized()) { root_structural.reset(); } + } + #pragma omp taskgroup { if (!context.settings.heuristics_only) { @@ -489,11 +514,25 @@ solution_t mip_solver_t::run_solver() } } + if (root_structural) { +#pragma omp task default(shared) priority(CUOPT_DEFAULT_TASK_PRIORITY) + { + root_structural->run(); + } + } + // Start the primal heuristics context.diversity_manager_ptr = &dm; sol = dm.run_solver(); } // implicit barrier for all tasks created in B&B and heuristics + dm.population.add_external_solutions_to_population(); + if (dm.population.is_feasible() && + (!sol.get_feasible() || + dm.population.best_feasible().get_objective() < sol.get_objective())) { + sol = solution_t(dm.population.best_feasible()); + } + if (!context.settings.heuristics_only && branch_and_bound->has_solver_space_incumbent()) { solution_t branch_and_bound_sol(*context.problem_ptr); branch_and_bound_sol.copy_new_assignment(branch_and_bound_solution.x); diff --git a/cpp/src/mip_heuristics/solver_context.cuh b/cpp/src/mip_heuristics/solver_context.cuh index f98386cbaf..f656c50b45 100644 --- a/cpp/src/mip_heuristics/solver_context.cuh +++ b/cpp/src/mip_heuristics/solver_context.cuh @@ -33,6 +33,9 @@ class diversity_manager_t; template class early_cpufj_t; +template +class early_structural_t; + // Aggregate structure containing the global context of the solving process for convenience: // The current problem, user settings, raft handle and statistics objects template @@ -66,6 +69,7 @@ struct mip_solver_context_t { work_unit_scheduler_t work_unit_scheduler_{5.0}; early_cpufj_t* early_cpufj_ptr{nullptr}; + early_structural_t* early_structural_ptr{nullptr}; // Best upper bound from early heuristics, in user-space. // Must be converted to the target solver-space before use: // - B&B: problem_ptr->get_solver_obj_from_user_obj(initial_upper_bound) diff --git a/cpp/src/mip_heuristics/structural/arc_flow.cu b/cpp/src/mip_heuristics/structural/arc_flow.cu new file mode 100644 index 0000000000..a5d0bc2941 --- /dev/null +++ b/cpp/src/mip_heuristics/structural/arc_flow.cu @@ -0,0 +1,980 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "arc_flow.cuh" + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +namespace { + +constexpr int arcflow_paths_supported = 2; +constexpr int arcflow_max_tokens = 20000; +constexpr int arcflow_max_col_entries = 3; +constexpr size_t arcflow_history_bytes_max = size_t{32} << 20; +constexpr size_t arcflow_candidate_bytes_max = size_t{32} << 20; + +enum class row_role_t : uint8_t { flow, cover }; + +template +struct arc_t { + i_t from{-1}; + i_t to{-1}; + i_t label{-1}; + i_t col{-1}; + f_t cost{0}; +}; + +template +struct arc_flow_model_t { + i_t n_nodes{0}; + i_t n_labels{0}; + + std::vector phi; + std::vector path_start; + std::vector demand; + std::vector displacement; + std::vector slope; + std::vector arc_offset; + std::vector> arcs; + + // A negative terminator_col denotes conservation-row slack. + std::vector terminator_col; + std::vector terminator_cost; + std::vector terminator_capacity; +}; + +template +struct frontier_t { + std::array node{}; + f_t cost{0}; +}; + +template +struct parent_t { + i_t prev{-1}; + i_t arc{-1}; +}; + +template +struct candidate_t { + frontier_t front; + parent_t parent; +}; + +template +struct arc_flow_result_t { + std::vector columns; + bool exact{true}; +}; + +template +bool is_known(f_t v) +{ + return !std::isnan(v); +} + +template +struct arcflow_profile_t { + arcflow_profile_t(i_t n_variables = 0, i_t n_constraints = 0) + : col_entries(n_variables, 0), + row_min_mag(n_constraints, std::numeric_limits::infinity()), + row_max_mag(n_constraints, 0) + { + } + + std::vector col_entries; + std::vector row_min_mag; + std::vector row_max_mag; +}; + +template +struct host_problem_t { + i_t n_variables{0}; + i_t n_constraints{0}; + std::vector csr_values; + std::vector csr_cols; + std::vector csr_offsets; + std::vector row_lb; + std::vector row_ub; + std::vector obj; + std::vector var_lb; + std::vector var_ub; + std::vector var_types; +}; + +bool arcflow_accepts_shape(int64_t n_variables, int64_t n_constraints, int64_t nnz) +{ + if (n_variables <= 0 || n_constraints <= 0) { return false; } + return nnz > 0 && nnz <= (int64_t)arcflow_max_col_entries * n_variables; +} + +// Every row is either two sided or bounded from below only, and both kinds must occur: the first +// become conservation rows and the second covering rows. +template +bool arcflow_accepts_bounds(const std::vector& row_lb, const std::vector& row_ub) +{ + int64_t n_flow_candidates = 0; + int64_t n_cover_candidates = 0; + for (size_t r = 0; r < row_lb.size(); ++r) { + const bool lo_fin = std::isfinite(row_lb[r]); + const bool hi_fin = std::isfinite(row_ub[r]); + if (lo_fin && hi_fin) { + n_flow_candidates++; + } else if (lo_fin) { + n_cover_candidates++; + } else { + return false; + } + } + return n_flow_candidates > 0 && n_cover_candidates > 0; +} + +// The row loop needs every magnitude in the row, so the acceptance pass cannot merge into the +// nonzero pass that builds the profile. +template +bool arcflow_accepts_profile( + const host_problem_t& h, + const typename mip_solver_settings_t::tolerances_t& tolerances, + arcflow_profile_t& p) +{ + p = arcflow_profile_t(h.n_variables, h.n_constraints); + for (i_t r = 0; r < h.n_constraints; ++r) { + for (i_t k = h.csr_offsets[r]; k < h.csr_offsets[r + 1]; ++k) { + const i_t col = h.csr_cols[k]; + cuopt_assert(col >= 0 && col < h.n_variables, "Column index out of range"); + if (++p.col_entries[col] > arcflow_max_col_entries) { return false; } + const f_t mag = std::abs(h.csr_values[k]); + p.row_min_mag[r] = std::min(p.row_min_mag[r], mag); + p.row_max_mag[r] = std::max(p.row_max_mag[r], mag); + } + } + + f_t cover_demand = 0; + for (i_t r = 0; r < h.n_constraints; ++r) { + if (p.row_max_mag[r] == 0 || p.row_min_mag[r] <= tolerances.absolute_tolerance) { + return false; + } + if (p.row_max_mag[r] - p.row_min_mag[r] > tolerances.absolute_tolerance) { return false; } + if (!std::isfinite(h.row_ub[r])) { + const f_t demand = h.row_lb[r] / p.row_max_mag[r]; + if (!is_integer(demand, tolerances.integrality_tolerance) || + demand < 1 - tolerances.absolute_tolerance) { + return false; + } + cover_demand += std::round(demand); + } + } + return cover_demand > 0 && cover_demand <= arcflow_max_tokens; +} + +template +struct row_info_t { + row_role_t role{row_role_t::cover}; + f_t scale{1}; + f_t lo{0}; + f_t hi{0}; +}; + +template +bool classify_rows(const host_problem_t& h, + const arcflow_profile_t& profile, + const typename mip_solver_settings_t::tolerances_t& tolerances, + std::vector>& rows) +{ + cuopt_assert((i_t)profile.row_min_mag.size() == h.n_constraints, "Size mismatch"); + cuopt_assert((i_t)profile.row_max_mag.size() == h.n_constraints, "Size mismatch"); + rows.assign(h.n_constraints, row_info_t{}); + + for (i_t r = 0; r < h.n_constraints; ++r) { + row_info_t info; + info.scale = profile.row_max_mag[r]; + const f_t lo = h.row_lb[r] / info.scale; + const f_t hi = h.row_ub[r] / info.scale; + const bool lo_fin = std::isfinite(lo); + const bool hi_fin = std::isfinite(hi); + + if (lo_fin && hi_fin) { + info.role = row_role_t::flow; + info.lo = lo; + info.hi = hi; + } else if (lo_fin) { + if (!is_integer(lo, tolerances.integrality_tolerance) || + lo < 1 - tolerances.absolute_tolerance) { + return false; + } + info.role = row_role_t::cover; + info.lo = std::round(lo); + info.hi = hi; + } else { + return false; + } + rows[r] = info; + } + return true; +} + +// The nonzero flow right-hand sides must all carry one sign. That sign says which incidence of an +// arc is its tail, and their absolute sum is the number of paths. Returns 0 when ambiguous. +template +f_t supply_orientation(const std::vector>& rows, + const typename mip_solver_settings_t::tolerances_t& tolerances, + f_t& total) +{ + const f_t absolute_tolerance = tolerances.absolute_tolerance; + f_t positive = 0; + f_t negative = 0; + for (const auto& info : rows) { + if (info.role != row_role_t::flow) { continue; } + if (info.lo > absolute_tolerance) { positive += info.lo; } + if (info.hi < -absolute_tolerance) { negative += -info.hi; } + } + if (positive > absolute_tolerance && negative > absolute_tolerance) { return 0; } + if (positive > absolute_tolerance) { + total = positive; + return 1; + } + if (negative > absolute_tolerance) { + total = negative; + return -1; + } + return 0; +} + +// One incidence per role at most, so a fourth column entry necessarily duplicates one of the three +// and is rejected without counting entries. +template +struct column_incidence_t { + i_t tail{-1}; + i_t head{-1}; + i_t label{-1}; +}; + +// A column carrying a single flow incidence is accepted only when that incidence is the arc's +// tail: sources are read from the right-hand side, so the mirror encoding of an explicit +// injection arc is out of scope here. +template +bool build_structure(const host_problem_t& h, + std::vector>& rows, + const typename mip_solver_settings_t::tolerances_t& tolerances, + arc_flow_model_t& model) +{ + f_t supply_total = 0; + const f_t sign = supply_orientation(rows, tolerances, supply_total); + if (sign == 0 || !is_integer(supply_total, tolerances.integrality_tolerance)) { + return false; + } + if (std::round(supply_total) != arcflow_paths_supported) { return false; } + + // Reorient so a normalized coefficient of +1 always means "the arc leaves this node". + if (sign < 0) { + for (auto& info : rows) { + if (info.role != row_role_t::flow) { continue; } + const f_t lo = info.lo; + info.lo = -info.hi; + info.hi = -lo; + } + } + + std::vector node_of_row(h.n_constraints, -1); + std::vector label_of_row(h.n_constraints, -1); + for (i_t r = 0; r < h.n_constraints; ++r) { + if (rows[r].role == row_role_t::flow) { + node_of_row[r] = model.n_nodes++; + } else { + label_of_row[r] = model.n_labels++; + } + } + if (model.n_nodes == 0 || model.n_labels == 0) { return false; } + + model.demand.assign(model.n_labels, 0); + model.terminator_col.assign(model.n_nodes, -1); + model.terminator_cost.assign(model.n_nodes, 0); + model.terminator_capacity.assign(model.n_nodes, 0); + + int64_t total_demand = 0; + for (i_t r = 0; r < h.n_constraints; ++r) { + const i_t l = label_of_row[r]; + if (l < 0) { continue; } + model.demand[l] = std::round(rows[r].lo); + total_demand += model.demand[l]; + } + if (total_demand <= 0 || total_demand > arcflow_max_tokens) { return false; } + + // Negative net outflow encodes path termination after singleton-column substitution. + for (i_t r = 0; r < h.n_constraints; ++r) { + const i_t v = node_of_row[r]; + if (v < 0) { continue; } + const auto& info = rows[r]; + if (info.lo > tolerances.absolute_tolerance) { + if (std::abs(info.lo - info.hi) > tolerances.absolute_tolerance || + !is_integer(info.lo, tolerances.integrality_tolerance)) { + return false; + } + model.path_start.insert(model.path_start.end(), (size_t)std::round(info.lo), v); + } else if (std::abs(info.hi) <= tolerances.absolute_tolerance) { + if (info.lo < -tolerances.absolute_tolerance) { + if (!is_integer(info.lo, tolerances.integrality_tolerance)) { return false; } + model.terminator_capacity[v] = std::min((f_t)arcflow_paths_supported, std::round(-info.lo)); + } + } else { + return false; + } + } + if ((i_t)model.path_start.size() != arcflow_paths_supported) { return false; } + + std::vector> columns(h.n_variables); + for (i_t r = 0; r < h.n_constraints; ++r) { + for (i_t k = h.csr_offsets[r]; k < h.csr_offsets[r + 1]; ++k) { + const i_t j = h.csr_cols[k]; + cuopt_assert(j >= 0 && j < h.n_variables, "Column index out of range"); + auto& column = columns[j]; + const f_t unit = h.csr_values[k] / rows[r].scale; + if (rows[r].role == row_role_t::flow) { + const f_t oriented = unit * sign; + if (std::abs(oriented - 1) <= tolerances.absolute_tolerance) { + if (column.tail >= 0) { return false; } + column.tail = node_of_row[r]; + } else if (std::abs(oriented + 1) <= tolerances.absolute_tolerance) { + if (column.head >= 0) { return false; } + column.head = node_of_row[r]; + } else { + return false; + } + } else { + // A covering incidence is positive irrespective of the flow orientation. + if (std::abs(unit - 1) > tolerances.absolute_tolerance) { return false; } + if (column.label >= 0) { return false; } + column.label = label_of_row[r]; + } + } + } + + for (i_t j = 0; j < h.n_variables; ++j) { + if (h.var_types[j] != var_t::INTEGER) { return false; } + if (!std::isfinite(h.obj[j])) { return false; } + if (std::abs(h.var_lb[j]) > tolerances.absolute_tolerance) { return false; } + const f_t ub = h.var_ub[j]; + if (!std::isfinite(ub) || ub < 1 - tolerances.absolute_tolerance) { return false; } + + const i_t tail = columns[j].tail; + const i_t head = columns[j].head; + const i_t label = columns[j].label; + if (tail >= 0 && head >= 0) { + if (label < 0) { return false; } + // With this bound the arc capacity can never bind: a use of the arc consumes one of the + // label's tokens, and there are exactly demand[label] of them. + if ((f_t)model.demand[label] > ub + tolerances.absolute_tolerance) { return false; } + model.arcs.push_back(arc_t{tail, head, label, j, h.obj[j]}); + } else if (tail >= 0 && label < 0) { + // Explicit loss arc: it leaves the node and never arrives, so it ends a path. A node that + // already absorbs paths through row slack would need the two capacities apportioned, which + // no model in this family does, so reject rather than guess. + if (model.terminator_col[tail] >= 0 || model.terminator_capacity[tail] > 0) { return false; } + model.terminator_col[tail] = j; + model.terminator_cost[tail] = h.obj[j]; + model.terminator_capacity[tail] = + std::min((f_t)arcflow_paths_supported, std::floor(ub + tolerances.absolute_tolerance)); + } else { + return false; + } + } + + if (model.arcs.empty()) { return false; } + + std::stable_sort( + model.arcs.begin(), model.arcs.end(), [](const arc_t& a, const arc_t& b) { + if (a.label != b.label) { return a.label < b.label; } + if (a.from != b.from) { return a.from < b.from; } + return a.cost < b.cost; + }); + model.arc_offset.assign(model.n_labels + 1, 0); + for (const auto& arc : model.arcs) { + model.arc_offset[arc.label + 1]++; + } + for (i_t l = 0; l < model.n_labels; ++l) { + model.arc_offset[l + 1] += model.arc_offset[l]; + if (model.arc_offset[l] == model.arc_offset[l + 1]) { return false; } + } + cuopt_assert(model.arc_offset.back() == (i_t)model.arcs.size(), + "arc CSR offsets must cover every arc"); + return true; +} + +// The potential is recovered from the objective rather than from any index: within a label the +// cost is affine in the potential of the arc's tail, so one label with enough distinct costs fixes +// the potential on every node it touches, and the remaining labels are fitted and extended from +// there. The result is an affine image of the true potential, which leaves the order that +// consumes it unchanged, since a common positive factor cancels out of every comparison. +template +bool derive_potential(arc_flow_model_t& model, + const typename mip_solver_settings_t::tolerances_t& tolerances, + const std::atomic& preemption_flag) +{ + const i_t n_labels = model.n_labels; + + i_t reference = -1; + size_t best_count = 0; + std::vector costs; + for (i_t l = 0; l < n_labels; ++l) { + costs.clear(); + for (i_t k = model.arc_offset[l]; k < model.arc_offset[l + 1]; ++k) { + costs.push_back(model.arcs[k].cost); + } + std::sort(costs.begin(), costs.end()); + const size_t distinct = (size_t)(std::unique(costs.begin(), costs.end()) - costs.begin()); + if (distinct > best_count) { + best_count = distinct; + reference = l; + } + } + if (reference < 0 || best_count < 2) { return false; } + + const f_t unknown = std::numeric_limits::quiet_NaN(); + model.phi.assign(model.n_nodes, unknown); + model.slope.assign(n_labels, unknown); + model.displacement.assign(n_labels, unknown); + std::vector intercept(n_labels, unknown); + + for (i_t k = model.arc_offset[reference]; k < model.arc_offset[reference + 1]; ++k) { + model.phi[model.arcs[k].from] = model.arcs[k].cost; + } + + // Each productive round resolves at least one potential, slope, or displacement. + const long max_rounds = 2L * n_labels + model.n_nodes + 2L; + long rounds = 0; + for (; rounds < max_rounds; ++rounds) { + if (preemption_flag.load()) { return false; } + bool progress = false; + for (i_t l = 0; l < n_labels; ++l) { + const i_t begin = model.arc_offset[l]; + const i_t end = model.arc_offset[l + 1]; + + if (!is_known(model.slope[l])) { + i_t lowest = -1; + i_t highest = -1; + for (i_t k = begin; k < end; ++k) { + const f_t p = model.phi[model.arcs[k].from]; + if (!is_known(p)) { continue; } + if (lowest < 0 || p < model.phi[model.arcs[lowest].from]) { lowest = k; } + if (highest < 0 || p > model.phi[model.arcs[highest].from]) { highest = k; } + } + if (lowest >= 0 && highest >= 0) { + const f_t lo_phi = model.phi[model.arcs[lowest].from]; + const f_t hi_phi = model.phi[model.arcs[highest].from]; + if (std::abs(lo_phi - hi_phi) > tolerances.absolute_tolerance) { + model.slope[l] = + (model.arcs[highest].cost - model.arcs[lowest].cost) / (hi_phi - lo_phi); + intercept[l] = model.arcs[lowest].cost - model.slope[l] * lo_phi; + progress = true; + } + } + } + if (is_known(model.slope[l]) && std::abs(model.slope[l]) > tolerances.absolute_tolerance) { + for (i_t k = begin; k < end; ++k) { + const i_t from = model.arcs[k].from; + if (is_known(model.phi[from])) { continue; } + model.phi[from] = (model.arcs[k].cost - intercept[l]) / model.slope[l]; + progress = true; + } + } + + if (!is_known(model.displacement[l])) { + for (i_t k = begin; k < end; ++k) { + const f_t from = model.phi[model.arcs[k].from]; + const f_t to = model.phi[model.arcs[k].to]; + if (is_known(from) && is_known(to)) { + model.displacement[l] = to - from; + progress = true; + break; + } + } + } + if (is_known(model.displacement[l])) { + for (i_t k = begin; k < end; ++k) { + const i_t from = model.arcs[k].from; + const i_t to = model.arcs[k].to; + if (is_known(model.phi[from]) && !is_known(model.phi[to])) { + model.phi[to] = model.phi[from] + model.displacement[l]; + progress = true; + } else if (is_known(model.phi[to]) && !is_known(model.phi[from])) { + model.phi[from] = model.phi[to] - model.displacement[l]; + progress = true; + } + } + } + } + if (!progress) { break; } + } + cuopt_assert(rounds < max_rounds, "propagation must reach a fixpoint within its progress bound"); + + for (f_t p : model.phi) { + if (!is_known(p)) { return false; } + } + for (f_t p : model.displacement) { + if (!is_known(p)) { return false; } + } + + f_t phi_scale = 0; + for (f_t p : model.phi) { + phi_scale = std::max(phi_scale, std::abs(p)); + } + if (phi_scale <= tolerances.absolute_tolerance) { return false; } + + // Orient the potential so displacements are positive. A consistent potential whose + // displacements are all strictly positive is exactly what makes the arc graph acyclic: a cycle + // would need its displacements to sum to zero. + f_t displacement_sum = 0; + for (f_t p : model.displacement) { + displacement_sum += p; + } + if (displacement_sum < 0) { + for (auto& p : model.phi) { + p = -p; + } + for (auto& p : model.displacement) { + p = -p; + } + for (auto& w : model.slope) { + w = -w; + } + } + for (f_t p : model.displacement) { + if (p <= tolerances.absolute_tolerance) { return false; } + } + + // Smith ordering assumes nonnegative job weights. + for (f_t w : model.slope) { + if (is_known(w) && w < -tolerances.absolute_tolerance) { return false; } + } + + // Verify against every arc, not just the two points each fit was built from. The cost residual + // matters as much as the potential: the order is read off the fitted slopes, so a label whose + // costs are not affine in the potential would be ordered on a meaningless quantity. + for (const auto& arc : model.arcs) { + const f_t displaced = model.phi[arc.from] + model.displacement[arc.label]; + if (std::abs(model.phi[arc.to] - displaced) > tolerances.absolute_tolerance) { return false; } + if (!is_known(model.slope[arc.label])) { continue; } + const f_t predicted = model.slope[arc.label] * model.phi[arc.from] + intercept[arc.label]; + if (std::abs(predicted - arc.cost) > tolerances.absolute_tolerance) { return false; } + } + return true; +} + +// Weighted shortest processing time orders labels by decreasing slope over displacement. +template +std::vector token_order( + const arc_flow_model_t& model, + const typename mip_solver_settings_t::tolerances_t& tolerances, + bool& all_slopes_identified) +{ + std::vector lowest_phi(model.n_labels, 0); + for (i_t l = 0; l < model.n_labels; ++l) { + f_t lowest = std::numeric_limits::infinity(); + for (i_t k = model.arc_offset[l]; k < model.arc_offset[l + 1]; ++k) { + lowest = std::min(lowest, model.phi[model.arcs[k].from]); + } + lowest_phi[l] = lowest; + } + + std::vector ordered; + std::vector unidentified; + ordered.reserve(model.n_labels); + for (i_t l = 0; l < model.n_labels; ++l) { + if (is_known(model.slope[l])) { + ordered.push_back(l); + } else { + unidentified.push_back(l); + } + } + + // Approximate equality is not transitive. Tolerance forms ratio classes before the final sort. + std::sort(ordered.begin(), ordered.end(), [&](i_t a, i_t b) { + const _Float128 lhs = (_Float128)model.slope[a] * (_Float128)model.displacement[b]; + const _Float128 rhs = (_Float128)model.slope[b] * (_Float128)model.displacement[a]; + if (lhs != rhs) { return lhs > rhs; } + return a < b; + }); + std::vector ratio_class(model.n_labels, 0); + for (size_t position = 1; position < ordered.size(); ++position) { + const i_t previous = ordered[position - 1]; + const i_t current = ordered[position]; + const _Float128 lhs = (_Float128)model.slope[previous] * (_Float128)model.displacement[current]; + const _Float128 rhs = (_Float128)model.slope[current] * (_Float128)model.displacement[previous]; + const _Float128 difference = lhs > rhs ? lhs - rhs : rhs - lhs; + const bool tied = difference <= (_Float128)tolerances.absolute_tolerance; + ratio_class[current] = ratio_class[previous] + (tied ? 0 : 1); + } + std::sort(ordered.begin(), ordered.end(), [&](i_t a, i_t b) { + if (ratio_class[a] != ratio_class[b]) { return ratio_class[a] < ratio_class[b]; } + if (lowest_phi[a] != lowest_phi[b]) { return lowest_phi[a] < lowest_phi[b]; } + if (model.displacement[a] != model.displacement[b]) { + return model.displacement[a] < model.displacement[b]; + } + return a < b; + }); + + int64_t total_demand = 0; + for (int64_t d : model.demand) { + total_demand += d; + } + std::vector tokens; + tokens.reserve((size_t)total_demand); + for (i_t l : ordered) { + tokens.insert(tokens.end(), (size_t)model.demand[l], l); + } + + // Labels without fitted slopes are placed at their first reachable potential. + all_slopes_identified = unidentified.empty(); + std::sort(unidentified.begin(), unidentified.end(), [&](i_t a, i_t b) { + if (lowest_phi[a] != lowest_phi[b]) { return lowest_phi[a] < lowest_phi[b]; } + return a < b; + }); + for (i_t l : unidentified) { + f_t consumed = 0; + size_t slot = 0; + while (slot < tokens.size() && consumed < lowest_phi[l]) { + consumed += model.displacement[tokens[slot]]; + ++slot; + } + tokens.insert(tokens.begin() + (ptrdiff_t)slot, (size_t)model.demand[l], l); + } + return tokens; +} + +template +std::optional> run_dp(const arc_flow_model_t& model, + const std::vector& tokens, + const std::atomic& preemption_flag) +{ + const i_t n_tokens = tokens.size(); + if (n_tokens == 0) { return std::nullopt; } + + arc_flow_result_t result; + + // The reconstruction budget is charged against retained states at each level. + size_t retained_bytes = 0; + + std::vector>> history; + history.reserve((size_t)n_tokens); + + frontier_t root; + for (i_t k = 0; k < arcflow_paths_supported; ++k) { + root.node[k] = model.path_start[k]; + } + std::sort(root.node.begin(), root.node.end()); + std::vector> current{root}; + + std::vector> candidates; + std::vector> next; + std::vector> parents; + for (i_t t = 0; t < n_tokens; ++t) { + if (preemption_flag.load()) { return std::nullopt; } + const i_t label = tokens[t]; + const auto arc_begin = model.arcs.begin() + model.arc_offset[label]; + const auto arc_end = model.arcs.begin() + model.arc_offset[label + 1]; + + const size_t candidate_limit = arcflow_candidate_bytes_max / sizeof(candidate_t); + size_t candidate_count = 0; + for (const auto& entry : current) { + for (i_t k = 0; k < arcflow_paths_supported; ++k) { + const i_t node = entry.node[k]; + const auto begin = std::lower_bound( + arc_begin, arc_end, node, [](const arc_t& a, i_t v) { return a.from < v; }); + const auto end = std::upper_bound( + begin, arc_end, node, [](i_t v, const arc_t& a) { return v < a.from; }); + const size_t added = end - begin; + if (added > candidate_limit - candidate_count) { return std::nullopt; } + candidate_count += added; + } + } + + candidates.clear(); + candidates.reserve(candidate_count); + for (i_t i = 0; i < (i_t)current.size(); ++i) { + const frontier_t& entry = current[i]; + for (i_t k = 0; k < arcflow_paths_supported; ++k) { + const i_t node = entry.node[k]; + const auto begin = std::lower_bound( + arc_begin, arc_end, node, [](const arc_t& a, i_t v) { return a.from < v; }); + const auto end = std::upper_bound( + begin, arc_end, node, [](i_t v, const arc_t& a) { return v < a.from; }); + for (auto it = begin; it != end; ++it) { + candidate_t candidate; + candidate.front = entry; + candidate.front.node[k] = it->to; + std::sort(candidate.front.node.begin(), candidate.front.node.end()); + candidate.front.cost = entry.cost + it->cost; + candidate.parent = parent_t{i, (i_t)(it - model.arcs.begin())}; + candidates.push_back(candidate); + } + } + } + if (candidates.empty()) { return std::nullopt; } + + // A total order makes representative selection independent of enumeration order. + std::sort(candidates.begin(), + candidates.end(), + [](const candidate_t& a, const candidate_t& b) { + if (a.front.node != b.front.node) { return a.front.node < b.front.node; } + if (a.front.cost != b.front.cost) { return a.front.cost < b.front.cost; } + if (a.parent.prev != b.parent.prev) { return a.parent.prev < b.parent.prev; } + return a.parent.arc < b.parent.arc; + }); + candidates.erase( + std::unique(candidates.begin(), + candidates.end(), + [](const candidate_t& a, const candidate_t& b) { + return a.front.node == b.front.node; + }), + candidates.end()); + + const size_t remaining = + arcflow_history_bytes_max > retained_bytes ? arcflow_history_bytes_max - retained_bytes : 0; + const size_t affordable = std::max(remaining / sizeof(parent_t), 1); + if (candidates.size() > affordable) { + std::stable_sort(candidates.begin(), + candidates.end(), + [](const candidate_t& a, const candidate_t& b) { + if (a.front.cost != b.front.cost) { return a.front.cost < b.front.cost; } + return a.front.node < b.front.node; + }); + candidates.resize(affordable); + std::sort(candidates.begin(), + candidates.end(), + [](const candidate_t& a, const candidate_t& b) { + return a.front.node < b.front.node; + }); + result.exact = false; + } + + next.clear(); + parents.clear(); + next.reserve(candidates.size()); + parents.reserve(candidates.size()); + for (const auto& candidate : candidates) { + next.push_back(candidate.front); + parents.push_back(candidate.parent); + } + retained_bytes += parents.size() * sizeof(parent_t); + history.push_back(std::move(parents)); + current.swap(next); + } + + i_t best_index = -1; + f_t best_total = std::numeric_limits::infinity(); + for (i_t i = 0; i < (i_t)current.size(); ++i) { + const frontier_t& entry = current[i]; + f_t total = entry.cost; + bool closable = true; + for (i_t k = 0; k < arcflow_paths_supported && closable; ++k) { + const i_t node = entry.node[k]; + int64_t sharing = 0; + for (i_t q = 0; q < arcflow_paths_supported; ++q) { + if (entry.node[q] == node) { sharing++; } + } + if (model.terminator_capacity[node] < sharing) { + closable = false; + } else { + total += model.terminator_cost[node]; + } + } + if (closable && total < best_total) { + best_total = total; + best_index = i; + } + } + if (best_index < 0) { return std::nullopt; } + + for (i_t k = 0; k < arcflow_paths_supported; ++k) { + const i_t col = model.terminator_col[current[best_index].node[k]]; + if (col >= 0) { result.columns.push_back(col); } + } + i_t index = best_index; + for (i_t t = n_tokens; t > 0; --t) { + const parent_t& step = history[t - 1][index]; + cuopt_assert(step.arc >= 0, "every level beyond the root records the arc it consumed"); + result.columns.push_back(model.arcs[step.arc].col); + index = step.prev; + } + cuopt_assert(index == 0, "reconstruction must terminate at the root state"); + return result; +} + +} // namespace + +template +struct arc_flow_t::host_state_t { + host_state_t(host_problem_t&& problem, arcflow_profile_t&& profile) + : h(std::move(problem)), profile(std::move(profile)) + { + } + + host_problem_t h; + arcflow_profile_t profile; +}; + +template +arc_flow_t::arc_flow_t() = default; + +template +arc_flow_t::~arc_flow_t() = default; + +template +bool arc_flow_t::recognize( + const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + const i_t n_variables = op_problem.get_n_variables(); + const i_t n_constraints = op_problem.get_n_constraints(); + if (!arcflow_accepts_shape(n_variables, n_constraints, op_problem.get_nnz())) { return false; } + if (op_problem.get_n_integers() != n_variables) { return false; } + if (!op_problem.get_variable_lower_bounds().is_empty() && + (i_t)op_problem.get_variable_lower_bounds().size() != n_variables) { + return false; + } + if ((i_t)op_problem.get_variable_upper_bounds().size() != n_variables) { return false; } + + auto stream = op_problem.get_handle_ptr()->get_stream(); + + host_problem_t h; + h.n_variables = n_variables; + h.n_constraints = n_constraints; + h.row_lb = cuopt::host_copy(op_problem.get_constraint_lower_bounds(), stream); + h.row_ub = cuopt::host_copy(op_problem.get_constraint_upper_bounds(), stream); + if (!arcflow_accepts_bounds(h.row_lb, h.row_ub)) { return false; } + + h.csr_values = cuopt::host_copy(op_problem.get_constraint_matrix_values(), stream); + h.csr_cols = cuopt::host_copy(op_problem.get_constraint_matrix_indices(), stream); + h.csr_offsets = cuopt::host_copy(op_problem.get_constraint_matrix_offsets(), stream); + arcflow_profile_t profile; + if (!arcflow_accepts_profile(h, tolerances, profile)) { return false; } + + h.obj = cuopt::host_copy(op_problem.get_objective_coefficients(), stream); + if (op_problem.get_sense()) { + for (auto& coefficient : h.obj) { + coefficient = -coefficient; + } + } + h.var_lb.assign(n_variables, f_t{0}); + if (!op_problem.get_variable_lower_bounds().is_empty()) { + h.var_lb = cuopt::host_copy(op_problem.get_variable_lower_bounds(), stream); + } + h.var_ub = cuopt::host_copy(op_problem.get_variable_upper_bounds(), stream); + h.var_types = cuopt::host_copy(op_problem.get_variable_types(), stream); + state_ = std::make_unique(std::move(h), std::move(profile)); + return true; +} + +template +bool arc_flow_t::recognize( + const problem_t& problem, + const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + const i_t n_variables = problem.n_variables; + const i_t n_constraints = problem.n_constraints; + if (!arcflow_accepts_shape(n_variables, n_constraints, problem.nnz)) { return false; } + if (problem.n_integer_vars != n_variables) { return false; } + + auto stream = problem.handle_ptr->get_stream(); + + host_problem_t h; + h.n_variables = n_variables; + h.n_constraints = n_constraints; + h.csr_values = cuopt::host_copy(problem.coefficients, stream); + h.csr_cols = cuopt::host_copy(problem.variables, stream); + h.csr_offsets = cuopt::host_copy(problem.offsets, stream); + h.row_lb = cuopt::host_copy(problem.constraint_lower_bounds, stream); + h.row_ub = cuopt::host_copy(problem.constraint_upper_bounds, stream); + if (!arcflow_accepts_bounds(h.row_lb, h.row_ub)) { return false; } + + arcflow_profile_t profile; + if (!arcflow_accepts_profile(h, tolerances, profile)) { return false; } + + h.obj = cuopt::host_copy(problem.objective_coefficients, stream); + h.var_types = cuopt::host_copy(problem.variable_types, stream); + std::tie(h.var_lb, h.var_ub) = + cuopt::extract_host_bounds(problem.variable_bounds, problem.handle_ptr); + state_ = std::make_unique(std::move(h), std::move(profile)); + return true; +} + +template +bool arc_flow_t::solve( + const typename mip_solver_settings_t::tolerances_t& tolerances, + std::atomic& preemption, + std::vector& assignment) +{ + cuopt_assert(state_ != nullptr, "solve called without a successful recognize"); + const auto& h = state_->h; + + std::vector> rows; + if (!classify_rows(h, state_->profile, tolerances, rows)) { + CUOPT_LOG_DEBUG("[ArcFlow] rejected: rows are not unit incidence after normalization"); + return false; + } + + arc_flow_model_t model; + if (!build_structure(h, rows, tolerances, model)) { + CUOPT_LOG_DEBUG("[ArcFlow] rejected: columns do not match the labelled arc pattern"); + return false; + } + if (preemption.load()) { return false; } + + if (!derive_potential(model, tolerances, preemption)) { + CUOPT_LOG_DEBUG("[ArcFlow] rejected: no consistent potential and affine cost model"); + return false; + } + if (preemption.load()) { return false; } + + bool all_slopes_identified = true; + const auto tokens = token_order(model, tolerances, all_slopes_identified); + CUOPT_LOG_DEBUG("[ArcFlow] detected %d nodes, %d labels, %d paths, %zu tokens, ordering %s", + (int)model.n_nodes, + (int)model.n_labels, + arcflow_paths_supported, + tokens.size(), + all_slopes_identified ? "identified" : "partly by reachability"); + + const auto result = run_dp(model, tokens, preemption); + if (!result.has_value()) { + CUOPT_LOG_DEBUG("[ArcFlow] no complete path set found in the ordered family"); + return false; + } + CUOPT_LOG_DEBUG("[ArcFlow] search %s", result->exact ? "exact" : "beamed by the history budget"); + + assignment.assign((size_t)h.n_variables, f_t{0}); + for (i_t col : result->columns) { + assignment[col] += f_t{1}; + } + return true; +} + +#if MIP_INSTANTIATE_FLOAT +template class arc_flow_t; +#endif + +#if MIP_INSTANTIATE_DOUBLE +template class arc_flow_t; +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/structural/arc_flow.cuh b/cpp/src/mip_heuristics/structural/arc_flow.cuh new file mode 100644 index 0000000000..9bb155637c --- /dev/null +++ b/cpp/src/mip_heuristics/structural/arc_flow.cuh @@ -0,0 +1,39 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +namespace cuopt::mathematical_optimization::mip { + +// Constructs paths through labelled arc-flow DAGs. Arc costs must be affine in tail potential; +// path termination may use explicit loss arcs or conservation-row slack. +template +class arc_flow_t : public structural_heuristic_t { + public: + arc_flow_t(); + ~arc_flow_t() override; + + const char* name() const override { return "ArcFlowDP"; } + + bool recognize(const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances) override; + bool recognize(const problem_t& problem, + const typename mip_solver_settings_t::tolerances_t& tolerances) override; + + bool solve(const typename mip_solver_settings_t::tolerances_t& tolerances, + std::atomic& preemption, + std::vector& assignment) override; + + private: + struct host_state_t; + + std::unique_ptr state_; +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/structural/early_structural.cu b/cpp/src/mip_heuristics/structural/early_structural.cu new file mode 100644 index 0000000000..0f64e49eff --- /dev/null +++ b/cpp/src/mip_heuristics/structural/early_structural.cu @@ -0,0 +1,217 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "early_structural.cuh" + +#include +#include +#include + +#include + +#include + +#include + +namespace cuopt::mathematical_optimization::mip { + +template +static bool validate(problem_t& problem, + const std::vector& assignment, + f_t& objective) +{ + if ((i_t)assignment.size() != problem.n_variables) { return false; } + solution_t solution(problem); + solution.copy_new_assignment(assignment); + if (has_variable_bounds_violation(problem.handle_ptr, solution.assignment, &problem) || + !solution.compute_feasibility()) { + return false; + } + objective = solution.get_objective(); + return true; +} + +template +static std::unique_ptr> make_structural_heuristic( + const model_t& model, const typename mip_solver_settings_t::tolerances_t& tolerances) +{ + auto heuristic = std::make_unique>(); + if (!heuristic->recognize(model, tolerances)) { return nullptr; } + return heuristic; +} + +template +std::unique_ptr> early_structural_t::create( + const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances, + early_incumbent_callback_t incumbent_callback) +{ + if (omp_get_num_threads() < CUOPT_MIP_EARLY_STRUCTURAL_REQUIRED_THREAD_COUNT) { return nullptr; } + auto active = make_structural_heuristic(op_problem, tolerances); + if (!active) { return nullptr; } + return std::unique_ptr(new early_structural_t( + op_problem, tolerances, std::move(incumbent_callback), std::move(active))); +} + +template +early_structural_t::early_structural_t( + const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances, + early_incumbent_callback_t incumbent_callback, + std::unique_ptr> active) + : early_heuristic_t>( + op_problem, tolerances, std::move(incumbent_callback)), + op_problem_(op_problem), + tolerances_(tolerances), + active_(std::move(active)) +{ + cuopt_assert(active_ != nullptr, "missing structural heuristic"); + CUOPT_LOG_DEBUG("[Early Structural] %s recognized the model", active_->name()); +} + +template +early_structural_t::~early_structural_t() +{ + stop(); +} + +template +void early_structural_t::start() +{ + if (task_launched_) { return; } + + preemption_flag_.store(false); + this->start_time_ = std::chrono::steady_clock::now(); + task_launched_ = true; + + // OpenMP depend clauses require a variable or array element. + auto* task_token = &preemption_flag_; + CUOPT_LOG_DEBUG("Launching early structural task for %s", active_->name()); +#pragma omp task priority(CUOPT_DEFAULT_TASK_PRIORITY) depend(out : *task_token) + this->run(); +} + +template +void early_structural_t::stop() +{ + if (!task_launched_) { return; } + + auto* task_token = &preemption_flag_; + preemption_flag_.store(true); +#pragma omp taskwait depend(in : *task_token) + task_launched_ = false; + + CUOPT_LOG_DEBUG("[Early Structural] Stopped, solution_found=%d", (int)this->solution_found_); +} + +template +bool early_structural_t::preprocessing_is_identity() const +{ + // Recognition produces assignments in the source problem's column space. + const auto& presolve_data = this->problem_ptr_->presolve_data; + if (this->problem_ptr_->n_variables != op_problem_.get_n_variables()) { return false; } + if ((i_t)presolve_data.variable_offsets.size() != this->problem_ptr_->n_variables) { + return false; + } + for (const f_t offset : presolve_data.variable_offsets) { + if (offset != f_t{0}) { return false; } + } + for (size_t j = 0; j < presolve_data.additional_var_used.size(); ++j) { + if (presolve_data.additional_var_used[j]) { return false; } + } + return true; +} + +template +void early_structural_t::run() +{ + cuopt_assert(active_ != nullptr, "task launched without a recognized structure"); + + std::vector assignment; + if (!active_->solve(tolerances_, preemption_flag_, assignment)) { + CUOPT_LOG_DEBUG("[Early Structural] %s constructed nothing", active_->name()); + return; + } + if (preemption_flag_.load()) { return; } + + if (!preprocessing_is_identity()) { + CUOPT_LOG_DEBUG( + "[Early Structural] %s constructed a point but preprocessing moved the columns, discarding", + active_->name()); + return; + } + + f_t objective{0}; + if (!validate(*this->problem_ptr_, assignment, objective)) { + CUOPT_LOG_DEBUG("[Early Structural] %s constructed a point that failed validation, discarding", + active_->name()); + return; + } + this->try_update_best(objective, assignment, active_->name()); +} + +template +root_structural_t::root_structural_t( + problem_t& problem, + const typename mip_solver_settings_t::tolerances_t& tolerances, + std::atomic& preemption, + structural_incumbent_callback_t incumbent_callback) + : tolerances_(tolerances), + preemption_(preemption), + incumbent_callback_(std::move(incumbent_callback)) +{ + RAFT_CUDA_TRY(cudaGetDevice(&device_id_)); + active_ = make_structural_heuristic(problem, tolerances); + if (!active_) { return; } + problem.handle_ptr->sync_stream(); + problem_ = std::make_unique>(problem, &handle_); + CUOPT_LOG_DEBUG("[Root Structural] %s recognized the model", active_->name()); +} + +template +root_structural_t::~root_structural_t() = default; + +template +void root_structural_t::run() +{ + if (!active_) { return; } + cuopt_assert(problem_ != nullptr, "missing structural problem"); + cuopt_assert(incumbent_callback_ != nullptr, "missing incumbent callback"); + + std::vector assignment; + if (!active_->solve(tolerances_, preemption_, assignment)) { + CUOPT_LOG_DEBUG("[Root Structural] %s constructed nothing", active_->name()); + return; + } + if (preemption_.load()) { return; } + + RAFT_CUDA_TRY(cudaSetDevice(device_id_)); + f_t objective{0}; + if (!validate(*problem_, assignment, objective)) { + CUOPT_LOG_DEBUG("[Root Structural] %s constructed a point that failed validation, discarding", + active_->name()); + return; + } + if (preemption_.load()) { return; } + + incumbent_callback_(assignment, objective); + CUOPT_LOG_DEBUG("[Root Structural] %s queued objective %+.6e", + active_->name(), + (double)problem_->get_user_obj_from_solver_obj(objective)); +} + +#if MIP_INSTANTIATE_FLOAT +template class early_structural_t; +template class root_structural_t; +#endif + +#if MIP_INSTANTIATE_DOUBLE +template class early_structural_t; +template class root_structural_t; +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/structural/early_structural.cuh b/cpp/src/mip_heuristics/structural/early_structural.cuh new file mode 100644 index 0000000000..883ec2e9fe --- /dev/null +++ b/cpp/src/mip_heuristics/structural/early_structural.cuh @@ -0,0 +1,100 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include + +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +template +class structural_heuristic_t { + public: + virtual ~structural_heuristic_t() = default; + + virtual const char* name() const = 0; + + virtual bool recognize( + const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances) = 0; + + virtual bool recognize(const problem_t&, + const typename mip_solver_settings_t::tolerances_t&) + { + return false; + } + + virtual bool solve(const typename mip_solver_settings_t::tolerances_t& tolerances, + std::atomic& preemption, + std::vector& assignment) = 0; +}; + +template +class early_structural_t : public early_heuristic_t> { + public: + static std::unique_ptr create( + const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances, + early_incumbent_callback_t incumbent_callback); + + ~early_structural_t(); + + static constexpr const char* name() { return "Structural"; } + + void start(); + void stop(); + + private: + early_structural_t(const optimization_problem_t& op_problem, + const typename mip_solver_settings_t::tolerances_t& tolerances, + early_incumbent_callback_t incumbent_callback, + std::unique_ptr> active); + + void run(); + + bool preprocessing_is_identity() const; + + const optimization_problem_t& op_problem_; + typename mip_solver_settings_t::tolerances_t tolerances_; + std::unique_ptr> active_; + std::atomic preemption_flag_{false}; + bool task_launched_{false}; +}; + +template +using structural_incumbent_callback_t = + std::function& assignment, f_t objective)>; + +template +class root_structural_t { + public: + root_structural_t(problem_t& problem, + const typename mip_solver_settings_t::tolerances_t& tolerances, + std::atomic& preemption, + structural_incumbent_callback_t incumbent_callback); + + ~root_structural_t(); + + bool recognized() const { return active_ != nullptr; } + + void run(); + + private: + int device_id_{0}; + raft::handle_t handle_; + std::unique_ptr> problem_; + typename mip_solver_settings_t::tolerances_t tolerances_; + std::atomic& preemption_; + structural_incumbent_callback_t incumbent_callback_; + std::unique_ptr> active_; +}; + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index b69fa08f46..9beef78df3 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -30,6 +30,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST ${CUOPT_TEST_DIR}/mip/block_bve_test.cu ${CUOPT_TEST_DIR}/mip/bhw_coeff_reduce_test.cpp ${CUOPT_TEST_DIR}/mip/gf2_presolve_test.cpp + ${CUOPT_TEST_DIR}/mip/arc_flow_test.cu ${CUOPT_TEST_DIR}/mip/single_lock_dual_aggregation_test.cpp ${CUOPT_TEST_DIR}/mip/termination_test.cu ${CUOPT_TEST_DIR}/mip/determinism_test.cu diff --git a/cpp/tests/mip/arc_flow_test.cu b/cpp/tests/mip/arc_flow_test.cu new file mode 100644 index 0000000000..ac0a37408f --- /dev/null +++ b/cpp/tests/mip/arc_flow_test.cu @@ -0,0 +1,412 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::test { + +namespace { + +struct built_model_t { + std::vector values; + std::vector indices; + std::vector offsets; + std::vector row_lb; + std::vector row_ub; + std::vector obj; + std::vector var_lb; + std::vector var_ub; + std::vector var_types; +}; + +struct build_options_t { + bool row_slack_terminators{false}; + bool permute{false}; + bool large_finite_capacities{false}; + double flow_row_factor{1.0}; + double cost_intercept{0.0}; + int perturbed_cost_label{-1}; + int negative_cost_label{-1}; + int single_arc_label{-1}; +}; + +enum state_t : int { t0, t1, t2, t3, t5, n_states }; + +struct arc_t { + state_t from; + state_t to; + int label; + double cost; +}; + +constexpr int n_labels = 3; +constexpr int n_paths = 2; +constexpr double expected_objective = 6.0; +constexpr std::array label_demand{2, 1, 1}; +constexpr int demanded_arcs = label_demand[0] + label_demand[1] + label_demand[2]; +constexpr std::array terminator_states{t2, t3, t5}; +constexpr std::array arcs{{{t0, t1, 0, 0.0}, + {t1, t2, 0, 3.0}, + {t2, t3, 0, 6.0}, + {t0, t3, 1, 0.0}, + {t2, t5, 1, 4.0}, + {t0, t2, 2, 0.0}, + {t1, t3, 2, 1.0}, + {t3, t5, 2, 3.0}}}; + +built_model_t build_arc_flow(const build_options_t& opts = {}) +{ + struct column_t { + std::vector> entries; + double cost; + double ub; + }; + + std::vector columns; + for (const auto& arc : arcs) { + if (arc.label == opts.single_arc_label && arc.from != t0) { continue; } + double cost = arc.cost + opts.cost_intercept; + if (arc.label == opts.negative_cost_label) { cost = -arc.cost; } + if (arc.label == opts.perturbed_cost_label && arc.from == t1) { cost += 1.0; } + const double ub = + opts.large_finite_capacities ? std::numeric_limits::max() : label_demand[arc.label]; + columns.push_back( + column_t{{{arc.from, 1.0}, {arc.to, -1.0}, {n_states + arc.label, 1.0}}, cost, ub}); + } + if (!opts.row_slack_terminators) { + const double ub = opts.large_finite_capacities ? std::numeric_limits::max() : 1.0; + for (const state_t state : terminator_states) { + columns.push_back(column_t{{{state, 1.0}}, 0.0, ub}); + } + } + + const int n_rows = n_states + n_labels; + std::vector row_lb(n_rows, 0.0); + std::vector row_ub(n_rows, 0.0); + row_lb[t0] = row_ub[t0] = n_paths; + if (opts.row_slack_terminators) { + for (const state_t state : terminator_states) { + row_lb[state] = -1.0; + row_ub[state] = 0.0; + } + } + for (int label = 0; label < n_labels; ++label) { + row_lb[n_states + label] = label_demand[label]; + row_ub[n_states + label] = std::numeric_limits::infinity(); + } + + const int n_cols = columns.size(); + std::vector row_perm(n_rows); + std::vector col_perm(n_cols); + std::iota(row_perm.begin(), row_perm.end(), 0); + std::iota(col_perm.begin(), col_perm.end(), 0); + if (opts.permute) { + std::reverse(row_perm.begin(), row_perm.end()); + for (int i = 0; i + 1 < n_cols; i += 2) { + std::swap(col_perm[i], col_perm[i + 1]); + } + } + + built_model_t model; + model.obj.assign(n_cols, 0.0); + model.var_lb.assign(n_cols, 0.0); + model.var_ub.assign(n_cols, 0.0); + model.var_types.assign(n_cols, var_t::INTEGER); + model.row_lb.assign(n_rows, 0.0); + model.row_ub.assign(n_rows, 0.0); + + for (int r = 0; r < n_rows; ++r) { + model.row_lb[row_perm[r]] = row_lb[r]; + model.row_ub[row_perm[r]] = row_ub[r]; + } + + std::vector>> by_row(n_rows); + for (int c = 0; c < n_cols; ++c) { + const auto& col = columns[c]; + const int mapped = col_perm[c]; + model.obj[mapped] = col.cost; + model.var_ub[mapped] = col.ub; + for (const auto& [row, value] : col.entries) { + by_row[row_perm[row]].emplace_back(mapped, value); + } + } + if (opts.flow_row_factor != 1.0) { + const int scaled = row_perm[1]; + for (auto& [col, value] : by_row[scaled]) { + value *= opts.flow_row_factor; + } + model.row_lb[scaled] *= opts.flow_row_factor; + model.row_ub[scaled] *= opts.flow_row_factor; + } + + model.offsets.push_back(0); + for (int r = 0; r < n_rows; ++r) { + std::sort(by_row[r].begin(), by_row[r].end()); + for (const auto& [col, value] : by_row[r]) { + model.indices.push_back(col); + model.values.push_back(value); + } + model.offsets.push_back(model.indices.size()); + } + return model; +} + +struct run_outcome_t { + bool prescreened{false}; + bool found{false}; + double objective{0.0}; + std::vector assignment; +}; + +struct input_options_t { + bool set_lower_bounds{true}; + bool set_upper_bounds{true}; + bool maximize{false}; + bool use_internal_problem{false}; +}; + +void expect_feasible(const built_model_t& model, const std::vector& assignment) +{ + ASSERT_EQ(assignment.size(), model.obj.size()); + for (size_t j = 0; j < assignment.size(); ++j) { + EXPECT_GE(assignment[j], model.var_lb[j]); + EXPECT_LE(assignment[j], model.var_ub[j]); + if (model.var_types[j] == var_t::INTEGER) { + EXPECT_DOUBLE_EQ(assignment[j], std::round(assignment[j])); + } + } + for (size_t r = 0; r < model.row_lb.size(); ++r) { + double activity = 0.0; + for (int k = model.offsets[r]; k < model.offsets[r + 1]; ++k) { + activity += model.values[k] * assignment[model.indices[k]]; + } + EXPECT_GE(activity, model.row_lb[r]); + EXPECT_LE(activity, model.row_ub[r]); + } +} + +run_outcome_t run_heuristic(const built_model_t& model, input_options_t options = {}) +{ + const raft::handle_t handle{}; + optimization_problem_t problem(&handle); + problem.set_csr_constraint_matrix(model.values.data(), + model.values.size(), + model.indices.data(), + model.indices.size(), + model.offsets.data(), + model.offsets.size()); + problem.set_objective_coefficients(model.obj.data(), model.obj.size()); + if (options.set_lower_bounds) { + problem.set_variable_lower_bounds(model.var_lb.data(), model.var_lb.size()); + } + if (options.set_upper_bounds) { + problem.set_variable_upper_bounds(model.var_ub.data(), model.var_ub.size()); + } + problem.set_variable_types(model.var_types.data(), model.var_types.size()); + problem.set_constraint_lower_bounds(model.row_lb.data(), model.row_lb.size()); + problem.set_constraint_upper_bounds(model.row_ub.data(), model.row_ub.size()); + problem.set_maximize(options.maximize); + + mip_solver_settings_t settings; + run_outcome_t outcome; + mip::arc_flow_t heuristic; + if (options.use_internal_problem) { + mip::problem_t internal_problem(problem, settings.get_tolerances(), false); + outcome.prescreened = heuristic.recognize(internal_problem, settings.get_tolerances()); + } else { + outcome.prescreened = heuristic.recognize(problem, settings.get_tolerances()); + } + if (!outcome.prescreened) { return outcome; } + + std::atomic preemption{false}; + outcome.found = heuristic.solve(settings.get_tolerances(), preemption, outcome.assignment); + if (outcome.found) { + expect_feasible(model, outcome.assignment); + outcome.objective = 0.0; + for (size_t j = 0; j < outcome.assignment.size(); ++j) { + outcome.objective += model.obj[j] * outcome.assignment[j]; + } + } + return outcome; +} + +} // namespace + +TEST(arc_flow, finds_exact_optimum_on_reduced_graph) +{ + const auto outcome = run_heuristic(build_arc_flow()); + ASSERT_TRUE(outcome.prescreened); + ASSERT_TRUE(outcome.found); + EXPECT_DOUBLE_EQ(outcome.objective, expected_objective); +} + +TEST(arc_flow, handles_maximization_in_both_recognizers) +{ + auto model = build_arc_flow(); + for (auto& coefficient : model.obj) { + coefficient = -coefficient; + } + const auto user_problem = run_heuristic(model, {.maximize = true}); + const auto internal_problem = + run_heuristic(model, {.maximize = true, .use_internal_problem = true}); + ASSERT_TRUE(user_problem.prescreened); + ASSERT_TRUE(user_problem.found); + ASSERT_TRUE(internal_problem.prescreened); + ASSERT_TRUE(internal_problem.found); + EXPECT_DOUBLE_EQ(user_problem.objective, -expected_objective); + EXPECT_DOUBLE_EQ(internal_problem.objective, -expected_objective); +} + +TEST(arc_flow, uses_solver_integrality_tolerance_for_scaled_demand) +{ + constexpr double test_integrality_tolerance = 1e-4; + mip_solver_settings_t settings; + auto tolerances = settings.get_tolerances(); + tolerances.integrality_tolerance = test_integrality_tolerance; + const double normalized_offset = tolerances.integrality_tolerance / 2.0; + + auto model = build_arc_flow(); + constexpr int cover_row = n_states; + constexpr int large_demand = 19998; + constexpr double row_scale = 0.1; + const double scaled_demand = (large_demand + normalized_offset) * row_scale; + model.row_lb[cover_row] = scaled_demand; + for (int k = model.offsets[cover_row]; k < model.offsets[cover_row + 1]; ++k) { + model.values[k] = row_scale; + } + std::fill(model.var_ub.begin(), model.var_ub.end(), large_demand); + + const double normalized_demand = scaled_demand / row_scale; + const double integrality_error = std::abs(normalized_demand - std::round(normalized_demand)); + EXPECT_GT(integrality_error, tolerances.absolute_tolerance); + EXPECT_LE(integrality_error, tolerances.integrality_tolerance); + + const raft::handle_t handle{}; + optimization_problem_t problem(&handle); + problem.set_csr_constraint_matrix(model.values.data(), + model.values.size(), + model.indices.data(), + model.indices.size(), + model.offsets.data(), + model.offsets.size()); + problem.set_objective_coefficients(model.obj.data(), model.obj.size()); + problem.set_variable_lower_bounds(model.var_lb.data(), model.var_lb.size()); + problem.set_variable_upper_bounds(model.var_ub.data(), model.var_ub.size()); + problem.set_variable_types(model.var_types.data(), model.var_types.size()); + problem.set_constraint_lower_bounds(model.row_lb.data(), model.row_lb.size()); + problem.set_constraint_upper_bounds(model.row_ub.data(), model.row_ub.size()); + + mip::arc_flow_t heuristic; + EXPECT_TRUE(heuristic.recognize(problem, tolerances)); +} + +TEST(arc_flow, single_arc_label_is_ordered_but_not_exact) +{ + const auto outcome = run_heuristic(build_arc_flow({.single_arc_label = 1})); + ASSERT_TRUE(outcome.prescreened); + ASSERT_TRUE(outcome.found); + EXPECT_DOUBLE_EQ(outcome.objective, expected_objective); +} + +TEST(arc_flow, accepts_supported_variants) +{ + struct test_case_t { + const char* name; + build_options_t build_options; + input_options_t input_options; + double objective; + }; + const std::array cases{{ + {.name = "permutation", .build_options = {.permute = true}, .objective = expected_objective}, + {.name = "row scaling", + .build_options = {.flow_row_factor = 4.0}, + .objective = expected_objective}, + {.name = "row-slack termination", + .build_options = {.row_slack_terminators = true}, + .objective = expected_objective}, + {.name = "implicit zero lower bounds", + .input_options = {.set_lower_bounds = false}, + .objective = expected_objective}, + {.name = "large finite capacities", + .build_options = {.large_finite_capacities = true}, + .objective = expected_objective}, + {.name = "affine cost intercept", + .build_options = {.cost_intercept = 7.0}, + .objective = expected_objective + 7.0 * demanded_arcs}, + }}; + for (const auto& test_case : cases) { + SCOPED_TRACE(test_case.name); + const auto outcome = + run_heuristic(build_arc_flow(test_case.build_options), test_case.input_options); + EXPECT_TRUE(outcome.prescreened); + EXPECT_TRUE(outcome.found); + if (outcome.found) { EXPECT_DOUBLE_EQ(outcome.objective, test_case.objective); } + } +} + +TEST(arc_flow, rejects_non_affine_costs) +{ + const auto outcome = run_heuristic(build_arc_flow({.perturbed_cost_label = 0})); + EXPECT_FALSE(outcome.found); +} + +TEST(arc_flow, rejects_negative_cost_slope) +{ + const auto outcome = run_heuristic(build_arc_flow({.negative_cost_label = 2})); + EXPECT_FALSE(outcome.found); +} + +TEST(arc_flow, rejects_model_without_unit_incidence) +{ + built_model_t knapsack; + knapsack.values = {2.0, 3.0}; + knapsack.indices = {0, 1}; + knapsack.offsets = {0, 2}; + knapsack.row_lb = {1.0}; + knapsack.row_ub = {std::numeric_limits::infinity()}; + knapsack.obj = {1.0, 1.0}; + knapsack.var_lb = {0.0, 0.0}; + knapsack.var_ub = {1.0, 1.0}; + knapsack.var_types = {var_t::INTEGER, var_t::INTEGER}; + const auto outcome = run_heuristic(knapsack); + EXPECT_FALSE(outcome.prescreened); + EXPECT_FALSE(outcome.found); +} + +TEST(arc_flow, rejects_implicit_infinite_upper_bounds) +{ + const auto outcome = run_heuristic(build_arc_flow(), {.set_upper_bounds = false}); + EXPECT_FALSE(outcome.prescreened); + EXPECT_FALSE(outcome.found); +} + +TEST(arc_flow, is_reproducible) +{ + const auto model = build_arc_flow(); + const auto first = run_heuristic(model); + const auto second = run_heuristic(model); + ASSERT_TRUE(first.found); + ASSERT_TRUE(second.found); + EXPECT_DOUBLE_EQ(first.objective, second.objective); + EXPECT_EQ(first.assignment, second.assignment); +} + +} // namespace cuopt::mathematical_optimization::test diff --git a/skills/cuopt-developer/SKILL.md b/skills/cuopt-developer/SKILL.md index ff1738bc5f..2aee1e7c75 100644 --- a/skills/cuopt-developer/SKILL.md +++ b/skills/cuopt-developer/SKILL.md @@ -170,6 +170,7 @@ cuopt/ ### CUDA/GPU Hygiene - Keep operations stream-ordered - Follow existing RAFT/RMM patterns +- `rmm::device_uvector` tests emptiness with `is_empty()`; it has no `empty()` member. - No raw `new`/`delete` - use RMM allocators - Prefer modern CCCL bit/math helpers in kernels (`cuda::bitfield_extract`, `cuda::bitmask`, pow2 utilities) over hand-rolled `%`/`/` by runtime powers of two — see [references/conventions.md](references/conventions.md) @@ -225,6 +226,8 @@ For pre-commit setup, DCO sign-off (`git commit -s`), the fork-based PR workflow ## Coding Conventions +Use `_Float128`, never `long double`, whenever extended precision arithmetic is required; `long double` is architecture-dependent and uses x87 on x86. + For C++ naming (`snake_case`, `d_`/`h_` prefixes, `_t` suffix), file extensions (`.hpp`/`.cpp`/`.cu`/`.cuh` and which compiler each uses), include order, Python style, error handling (`CUOPT_EXPECTS`, `RAFT_CUDA_TRY`), memory management (RMM patterns, no raw `new`/`delete`), CCCL bit/math helpers in device code, test-impact rules, volatile-comment rules (hardware names and self-referential issue/PR numbers in comments or skip messages go stale; issue links to a separate tracking issue are fine), **no large local lambdas** (extract named helpers instead), and **coarse work-estimate / time-limit gating** (phase/outer-loop only; no fine inner-loop or double checks), see [references/conventions.md](references/conventions.md). ## OpenMP task/runtime compatibility