diff --git a/ci/validate_wheel.sh b/ci/validate_wheel.sh index a603c69098..9d53998a42 100755 --- a/ci/validate_wheel.sh +++ b/ci/validate_wheel.sh @@ -22,7 +22,7 @@ PYDISTCHECK_ARGS=( if [[ "${package_dir}" == "python/libcuopt" ]]; then if [[ "${RAPIDS_CUDA_MAJOR}" == "12" ]]; then PYDISTCHECK_ARGS+=( - --max-allowed-size-compressed '690Mi' + --max-allowed-size-compressed '695Mi' ) else PYDISTCHECK_ARGS+=( diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 86ed6965c8..a5eb4ed09b 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -144,6 +144,9 @@ #define CUOPT_MIP_HYPER_SUBMIP_ITERATION_LIMIT_RATIO "mip_hyper_submip_iteration_limit_ratio" #define CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ "mip_hyper_submip_enable_cpufj" +/* @brief Block bounded-variable-elimination step of cuOpt's internal MIP presolve */ +#define CUOPT_MIP_HYPER_BLOCK_BVE "mip_hyper_block_bve" + /* @brief QCQP (barrier) scaling hyper-parameters */ #define CUOPT_QCQP_HYPER_RUIZ_EQUILIBRATION "qcqp_hyper_ruiz_equilibration" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index 8ddcdbbb8a..7daee3c682 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -160,6 +160,15 @@ class mip_solver_settings_t { * When this is `false`, probing is skipped even if presolve is otherwise on. */ bool probing{true}; + /** + * @brief Enable the block bounded-variable-elimination step of cuOpt's MIP presolve. + * + * Runs after trivial_presolve and eliminates blocks of functionally-determined binary auxiliary + * variables discovered via the probing-cache implication closure, re-encoding each block's + * projected relation as certified prime-implicate clauses. Requires the probing-cache step; a + * no-op when no certified reduction exists. + */ + bool block_bve{true}; /** * @brief Determinism mode for MIP solver. * diff --git a/cpp/src/io/mps_writer.cpp b/cpp/src/io/mps_writer.cpp index d269d6ec8a..9275fd7685 100644 --- a/cpp/src/io/mps_writer.cpp +++ b/cpp/src/io/mps_writer.cpp @@ -228,8 +228,8 @@ void mps_writer_t::write(const std::string& mps_file_path) // save coefficients with full precision mps_file << std::setprecision(std::numeric_limits::max_digits10); - // NAME section - mps_file << "NAME " << problem_.get_problem_name() << "\n"; + const std::string& pname = problem_.get_problem_name(); + mps_file << "NAME " << (pname.empty() ? "cuopt" : pname) << "\n"; if (problem_.get_sense()) { mps_file << "OBJSENSE\n MAXIMIZE\n"; } diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 11f88b36d7..09289e05db 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -211,6 +211,9 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_HYPER_DIVING_SHOW_TYPE, &mip_settings.diving_params.show_type, false, "log diving heuristic type when it finds a new incumbent"}, // Recursive sub-MIP (RINS) hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_SUBMIP_ENABLE_CPUFJ, &mip_settings.submip_params.enable_cpufj, true, "run CPU FJ over the sub-MIP"}, + // Kept a hyper-parameter while block-BVE bakes in: settable so a run can be bisected against it, + // but not yet a documented knob (no constant in the proto / server surfaces). + {CUOPT_MIP_HYPER_BLOCK_BVE, &mip_settings.block_bve, true, "eliminate blocks of binaries in cuOpt's MIP presolve (needs " CUOPT_MIP_PROBING ")"}, }; // String parameters string_parameters = { diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 7705465512..6ad1009d84 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -33,6 +33,7 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/simple_rounding.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/feasibility_pump/feasibility_pump.cu ${CMAKE_CURRENT_SOURCE_DIR}/local_search/line_segment_search/line_segment_search.cu + ${CMAKE_CURRENT_SOURCE_DIR}/presolve/block_bve.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_presolve.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/bounds_update_data.cu ${CMAKE_CURRENT_SOURCE_DIR}/presolve/semi_continuous.cu diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index 61c90944f4..b65e392efa 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -11,6 +11,7 @@ #include #include +#include #include #include #include @@ -18,6 +19,10 @@ #include +#include +#include +#include +#include #include #include @@ -39,6 +44,59 @@ size_t sub_mip_recombiner_config_t::max_n_of_vars_from_other = template std::vector recombiner_t::enabled_recombiners; +template +static cuopt::mathematical_optimization::io::mps_data_model_t problem_to_mps_data_model( + const problem_t& problem) +{ + auto stream = problem.handle_ptr->get_stream(); + auto h_off = cuopt::host_copy(problem.offsets, stream); + auto h_ind = cuopt::host_copy(problem.variables, stream); + auto h_val = cuopt::host_copy(problem.coefficients, stream); + auto h_clb = cuopt::host_copy(problem.constraint_lower_bounds, stream); + auto h_cub = cuopt::host_copy(problem.constraint_upper_bounds, stream); + auto h_obj = cuopt::host_copy(problem.objective_coefficients, stream); + auto h_vb = cuopt::host_copy(problem.variable_bounds, stream); + auto h_vt = cuopt::host_copy(problem.variable_types, stream); + problem.handle_ptr->sync_stream(); + + const i_t n_vars = problem.n_variables; + std::vector var_lower(n_vars), var_upper(n_vars); + for (i_t v = 0; v < n_vars; ++v) { + var_lower[v] = get_lower(h_vb[v]); + var_upper[v] = get_upper(h_vb[v]); + } + std::vector var_types(n_vars); + for (i_t v = 0; v < n_vars; ++v) + var_types[v] = var_type_to_char(h_vt[v]); + + cuopt::mathematical_optimization::io::mps_data_model_t model; + model.set_maximize(false); + if (!h_off.empty()) { + model.set_csr_constraint_matrix(std::span{h_val.data(), h_val.size()}, + std::span{h_ind.data(), h_ind.size()}, + std::span{h_off.data(), h_off.size()}); + } + if (problem.n_constraints != 0) { + model.set_constraint_lower_bounds(std::span{h_clb.data(), h_clb.size()}); + model.set_constraint_upper_bounds(std::span{h_cub.data(), h_cub.size()}); + } + if (n_vars != 0) { + model.set_objective_coefficients(std::span{h_obj.data(), h_obj.size()}); + model.set_variable_lower_bounds(std::span{var_lower.data(), var_lower.size()}); + model.set_variable_upper_bounds(std::span{var_upper.data(), var_upper.size()}); + model.set_variable_types(var_types); + } + model.set_objective_scaling_factor(problem.presolve_data.objective_scaling_factor); + model.set_objective_offset(problem.presolve_data.objective_offset); + if (problem.original_problem_ptr != nullptr && + !problem.original_problem_ptr->get_problem_name().empty()) { + model.set_problem_name(problem.original_problem_ptr->get_problem_name()); + } else { + model.set_problem_name("cuopt"); + } + return model; +} + template diversity_manager_t::diversity_manager_t(mip_solver_context_t& context_) : context(context_), @@ -279,6 +337,42 @@ void diversity_manager_t::add_user_given_solutions( } } +// Pin variables that a BVE projection table showed to have a single admissible value. Ids arrive in +// the original frame and may repeat across blocks and rounds. Returns false when two blocks +// disagree on a variable, which proves infeasibility since each fixing is a consequence of its +// block alone. +template +static bool apply_bve_fixings(problem_t& problem, + const std::vector>& fixings, + i_t& n_applied) +{ + n_applied = 0; + if (fixings.empty()) { return true; } + std::vector> sorted(fixings); + std::sort(sorted.begin(), sorted.end()); + + const std::vector& reverse_original_ids = problem.reverse_original_ids; + std::vector var_indices; + std::vector lb_values; + std::vector ub_values; + for (size_t k = 0; k < sorted.size(); ++k) { + const auto [original_id, value] = sorted[k]; + if (k > 0 && original_id == sorted[k - 1].first) { + if (value != sorted[k - 1].second) { return false; } + continue; + } + if (original_id < 0 || original_id >= (i_t)reverse_original_ids.size()) { continue; } + const i_t column = reverse_original_ids[original_id]; + if (column < 0 || column >= problem.n_variables) { continue; } // already eliminated + var_indices.push_back(column); + lb_values.push_back(value ? f_t(1) : f_t(0)); + ub_values.push_back(value ? f_t(1) : f_t(0)); + } + n_applied = (i_t)var_indices.size(); + problem.update_variable_bounds(var_indices, lb_values, ub_values); + return true; +} + template bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_timer) { @@ -304,19 +398,102 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ CUOPT_LOG_INFO("Probing-cache step disabled via %s=false", CUOPT_MIP_PROBING); run_probing_cache = false; } - if (run_probing_cache) { - // Run probing cache before trivial presolve to discover variable implications + const bool remap_cache_ids = true; + problem_ptr->related_vars_time_limit = context.settings.heuristic_params.related_vars_time_limit; + + if (run_probing_cache && !global_timer.check_time_limit() && !presolve_timer.check_time_limit()) { const f_t max_time_on_probing = diversity_config.max_time_on_probing; - f_t time_for_probing_cache = std::min(max_time_on_probing, time_limit); + f_t time_for_probing_cache = + std::min(max_time_on_probing, std::min(time_limit, (f_t)presolve_timer.remaining_time())); timer_t probing_timer{time_for_probing_cache}; - // this function computes probing cache, finds singletons, substitutions and changes the problem bool problem_is_infeasible = compute_probing_cache(ls.constraint_prop.bounds_update, *problem_ptr, probing_timer); if (problem_is_infeasible) { return false; } } - const bool remap_cache_ids = true; - problem_ptr->related_vars_time_limit = context.settings.heuristic_params.related_vars_time_limit; + if (!global_timer.check_time_limit()) { trivial_presolve(*problem_ptr, remap_cache_ids); } + + i_t max_bve_rounds = 3; + const i_t n_vars_before_bve = problem_ptr->n_variables; + const i_t n_rows_before_bve = problem_ptr->n_constraints; + + if (!run_probing_cache) max_bve_rounds = 0; + // Implications read off the projection tables, accumulated across rounds. They feed the next + // round's adjacency (pairs the cache never held) and are folded back into the cache afterwards. + probe_findings_t bve_findings; + for (i_t bve_round = 0; bve_round < max_bve_rounds; ++bve_round) { + if (!context.settings.block_bve || problem_ptr->empty || global_timer.check_time_limit() || + presolve_timer.check_time_limit()) { + break; + } + + const i_t n_vars_before = problem_ptr->n_variables; + const i_t n_rows_before = problem_ptr->n_constraints; + auto impl_adj = bve_build_impl_adj(ls.constraint_prop.bounds_update.probing_cache, + problem_ptr->reverse_original_ids, + problem_ptr->n_variables, + &bve_findings); + double bve_work_units = 0.0; + timer_t bve_timer(global_timer.clamp_remaining_time(presolve_timer.remaining_time())); + const bool reduced = + block_bve_presolve(*problem_ptr, impl_adj, bve_timer, bve_work_units, &bve_findings); + CUOPT_LOG_DEBUG("Block-BVE outer round %d/%d: reduced=%d vars %d->%d rows %d->%d", + bve_round + 1, + max_bve_rounds, + (int)reduced, + n_vars_before, + problem_ptr->n_variables, + n_rows_before, + problem_ptr->n_constraints); + if (!reduced) { break; } + if (problem_ptr->n_variables >= n_vars_before) { break; } + } + + // Harvest the projections: tighten the cache in place, pin the variables the blocks left with a + // single value, then propagate. + ls.constraint_prop.bounds_update.probing_cache.merge_forcings(bve_findings.forcings, + bve_findings.fixings); + i_t n_bve_fixings = 0; + if (!global_timer.check_time_limit()) { + if (!apply_bve_fixings(*problem_ptr, bve_findings.fixings, n_bve_fixings)) { + stats.presolve_time = timer.elapsed_time(); + return false; + } + if (n_bve_fixings > 0) { trivial_presolve(*problem_ptr, remap_cache_ids); } + } + const bool bve_changed_model = problem_ptr->n_variables != n_vars_before_bve || + problem_ptr->n_constraints != n_rows_before_bve || + n_bve_fixings > 0; + if (bve_changed_model) { + CUOPT_LOG_DEBUG("Block-BVE projections fixed %d variables", n_bve_fixings); + // propagate fixings if any + if (!problem_ptr->empty && !global_timer.check_time_limit()) { + ls.constraint_prop.bounds_update.resize(*problem_ptr); + auto bve_term_crit = ls.constraint_prop.bounds_update.solve(*problem_ptr); + if (ls.constraint_prop.bounds_update.infeas_constraints_count > 0) { + stats.presolve_time = timer.elapsed_time(); + return false; + } + if (termination_criterion_t::NO_UPDATE != bve_term_crit) { + ls.constraint_prop.bounds_update.set_updated_bounds(*problem_ptr); + } + } + } + + if (const char* export_flag = std::getenv("CUOPT_EXPORT_GPU_PRESOLVED_PROBLEM"); + export_flag != nullptr && std::atoi(export_flag) != 0) { + const std::string instance_name = + (problem_ptr->original_problem_ptr != nullptr && + !problem_ptr->original_problem_ptr->get_problem_name().empty()) + ? problem_ptr->original_problem_ptr->get_problem_name() + : std::string("cuopt"); + const std::string mps_path = instance_name + "_gpupresolved.mps"; + CUOPT_LOG_DEBUG("Exporting GPU-presolved problem to %s", mps_path.c_str()); + auto model = problem_to_mps_data_model(*problem_ptr); + cuopt::mathematical_optimization::io::mps_writer_t writer(model); + writer.write(mps_path); + exit(0); + } if (!problem_ptr->empty && !check_bounds_sanity(*problem_ptr)) { return false; } // if (!presolve_timer.check_time_limit() && !context.settings.heuristics_only && // !problem_ptr->empty) { diff --git a/cpp/src/mip_heuristics/diversity/population.cu b/cpp/src/mip_heuristics/diversity/population.cu index 553e5d6e93..e6fffa97b2 100644 --- a/cpp/src/mip_heuristics/diversity/population.cu +++ b/cpp/src/mip_heuristics/diversity/population.cu @@ -233,6 +233,12 @@ std::vector> population_t::get_external_solutions sol.compute_number_of_integers(), problem_ptr->n_integer_vars); } + if (std::abs(sol.get_objective() - h_entry.objective) > OBJECTIVE_EPSILON) { + CUOPT_LOG_DEBUG( + "External solution objective mismatch: sol.get_objective() = %g, h_entry.objective = %g", + sol.get_objective(), + h_entry.objective); + } sol.handle_ptr->sync_stream(); return_vector.emplace_back(std::move(sol)); counter++; diff --git a/cpp/src/mip_heuristics/presolve/block_bve.cu b/cpp/src/mip_heuristics/presolve/block_bve.cu new file mode 100644 index 0000000000..708456b078 --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/block_bve.cu @@ -0,0 +1,1373 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "block_bve.cuh" +#include "trivial_presolve.cuh" + +#include +#include +#include + +#include // find_scaling_rational (exact row integerization) + +#include // raft::warpReduce +#include + +#include + +#include // cuda::bitfield_extract + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +// Block caps the header does not need to expose (its struct layouts and default arguments pin the +// rest). +static constexpr int BVE_MAX_INTERIOR = BVE_MAX_SCOPE - 1; +// Cap closure probes over high-degree implication neighborhoods. +static constexpr int BVE_MAX_GROWTH_NBRS = 256; +// Cap peak device allocation for each projection chunk. +static constexpr size_t BVE_PROJECT_DEVICE_BUDGET = 64ull << 20; // 64 MiB + +// Largest per-row rational multiplier / denominator we will apply. A row that would need a larger +// multiplier to become integer is treated as not exactly representable (passed to +// find_scaling_rational as its maxdnom/maxfinal caps). +static constexpr int64_t BVE_INT_SCALE_MAX = 1000000; // 1e6 + +// Closed-form part of the commit_projected work estimate: prime-cube enumeration in +// bve_greedy_prime_cover is Θ(nb · 3^nb); sanity check is Θ(2^nb · #clauses) with #clauses bounded +// by the growth gate (n_rows + clause_growth_margin). The cover build and the greedy selection +// scale with the prime count, which is only known after enumeration, so they are metered from the +// inside and reported through the commit_projected ops out-param instead. +static double bve_commit_wall_ops(int nb, int clause_budget) +{ + cuopt_assert(nb >= 0 && nb <= BVE_MAX_BOUNDARY, "nb out of BVE range"); + double three_nb = 1.0; + for (int i = 0; i < nb; ++i) + three_nb *= 3.0; + return double(nb) * three_nb + double(1 << nb) * double(clause_budget + 1); +} + +bool bve_sanity_check(const uint8_t* feas, int nb, const bve_clause_t* clauses, int n_clauses) +{ + const uint32_t full_mask = (1u << nb) - 1u; + for (int i = 0; i < n_clauses; ++i) + if (clauses[i].lit_mask & ~full_mask) return false; // literals must be on the boundary + for (uint32_t m = 0; m <= full_mask; ++m) { + bool crel = true; // CNF value: AND over clauses of (clause satisfied by pattern m) + for (int i = 0; i < n_clauses && crel; ++i) { + const uint32_t lit = clauses[i].lit_mask; + const uint32_t bit = clauses[i].bit_mask; + // clause satisfied iff some literal position differs from its forbidden bit under m + const bool satisfied = ((m ^ bit) & lit) != 0u; + if (!satisfied) crel = false; + } + const bool feasible = feas[m] != 0; + if (crel != feasible) return false; + } + return true; +} + +// =========================================================================================== +// Installed CNF: all prime forbidden cubes covered by max-gain greedy +// =========================================================================================== +// +// Two-level logic minimization in the shape of Quine, "The Problem of Simplifying Truth Functions" +// (Amer. Math. Monthly 1952) and McCluskey, "Minimization of Boolean Functions" (Bell System Tech. +// J. 1956): enumerate the prime implicants, then cover every minterm with a subset of them. Taking +// the primes of the infeasible patterns rather than the feasible ones makes each one a forbidden +// cube whose complement is a clause, so the cover comes out as a CNF instead of the usual DNF. +// +// The covering step is the greedy max-gain heuristic of Johnson, "Approximation Algorithms for +// Combinatorial Problems" (JCSS 1974), Lovász (Discrete Math. 1975) and Chvátal (Math. of OR 1979): +// repeatedly take the cube covering the most still-uncovered patterns, which lands within a factor +// 1 + ln m of the minimum cover for m infeasible patterns. An exact minimum cover (Petrick / unate +// covering) was measured against it: the greedy already hit the optimum on 90% of blocks and took +// 61% of the clauses an exact cover would have saved, which did not pay for owning a +// branch-and-bound with a node cap and a fallback path. + +static size_t bve_mask_words(int nb) { return size_t(((1u << nb) + 63u) / 64u); } + +static int bve_mask_size(const bve_mask_t& m) +{ + int n = 0; + for (uint64_t w : m) + n += std::popcount(w); + return n; +} + +static void bve_mask_set(bve_mask_t& m, uint32_t pattern) +{ + cuopt_assert(size_t(pattern >> 6) < m.size(), "pattern outside mask width"); + m[pattern >> 6] |= uint64_t{1} << (pattern & 63); +} + +static void bve_mask_subtract(bve_mask_t& m, const bve_mask_t& other) +{ + cuopt_assert(m.size() == other.size(), "mask width mismatch"); + for (size_t w = 0; w < m.size(); ++w) + m[w] &= ~other[w]; +} + +static int bve_mask_overlap(const bve_mask_t& a, const bve_mask_t& b) +{ + cuopt_assert(a.size() == b.size(), "mask width mismatch"); + int n = 0; + for (size_t w = 0; w < a.size(); ++w) + n += std::popcount(a[w] & b[w]); + return n; +} + +// valid(lit, bit): every boundary pattern matching cube (lit, bit) is infeasible, so the +// complementary clause excludes no feasible pattern. Adding a literal SHRINKS the cube, so the +// table is filled from the minterms (lit == full_mask) downward in literal count: +// valid(lit, bit) = valid(lit|j, bit) AND valid(lit|j, bit|j) for any j not in lit +// A cube is already the bve_clause_t (lit_mask, bit_mask) encoding, so no separate ternary cube +// code is needed. The dense table is 4^nb bytes (16 MiB at nb = 12), so `valid` is caller-owned and +// grown once rather than reallocated per block. It is never re-initialized: only cells with +// bit subset of lit are ever addressed, the minterm seeding plus the recurrence below write every +// such cell, and each pass reads only cells an earlier pass already wrote. +static void bve_enumerate_prime_cubes(const uint8_t* feas, + int nb, + std::vector& valid, + std::vector& primes) +{ + cuopt_assert(nb >= 1 && nb <= BVE_MAX_BOUNDARY, "nb out of BVE range"); + const uint32_t full_mask = (1u << nb) - 1u; + const size_t stride = size_t(full_mask) + 1; + if (valid.size() < stride * stride) valid.resize(stride * stride); + const auto at = [&](uint32_t lit, uint32_t bit) -> uint8_t& { + cuopt_assert((bit & ~lit) == 0u, "cube bit_mask outside its lit_mask"); + return valid[size_t(lit) * stride + bit]; + }; + + for (uint32_t m = 0; m <= full_mask; ++m) + at(full_mask, m) = feas[m] ? 0 : 1; + + for (int n_lits = nb - 1; n_lits >= 0; --n_lits) + for (uint32_t lit = 0; lit <= full_mask; ++lit) { + if (std::popcount(lit) != n_lits) continue; + const int j = std::countr_zero(~lit & full_mask); + const uint32_t child = lit | (1u << j); + for (uint32_t bit = lit;; bit = (bit - 1u) & lit) { + at(lit, bit) = at(child, bit) & at(child, bit | (1u << j)); + if (bit == 0u) break; + } + } + + primes.clear(); + for (uint32_t lit = 0; lit <= full_mask; ++lit) + for (uint32_t bit = lit;; bit = (bit - 1u) & lit) { + if (at(lit, bit)) { + bool prime = true; + for (int j = 0; j < nb && prime; ++j) + if ((lit & (1u << j)) != 0u && at(lit ^ (1u << j), bit & ~(1u << j))) prime = false; + if (prime) primes.push_back(bve_clause_t{lit, bit}); + } + if (bit == 0u) break; + } +} + +// Boundary patterns matching the cube; every one of them is infeasible when the cube is valid. +static void bve_cube_cover( + uint32_t lit, uint32_t bit, uint32_t full_mask, size_t n_words, bve_mask_t& cover) +{ + cover.assign(n_words, 0u); + const uint32_t free_positions = full_mask & ~lit; + for (uint32_t s = free_positions;; s = (s - 1u) & free_positions) { + bve_mask_set(cover, bit | s); + if (s == 0u) break; + } +} + +int bve_greedy_prime_cover(const uint8_t* feas, + int nb, + bve_clause_t* out, + int cap, + bve_cover_scratch_t& scratch, + int64_t* ops_out) +{ + cuopt_assert(nb >= 1 && nb <= BVE_MAX_BOUNDARY, "nb out of BVE range"); + cuopt_assert(cap >= 1, "clause cap leaves no room for a cover"); + int64_t ops = 0; + auto ops_guard = cuopt::scope_guard([&]() { + if (ops_out != nullptr) *ops_out += ops; + }); + + const uint32_t full_mask = (1u << nb) - 1u; + const uint32_t n_patterns = 1u << nb; + const size_t n_words = bve_mask_words(nb); + + bve_enumerate_prime_cubes(feas, nb, scratch.valid, scratch.primes); + const std::vector& primes = scratch.primes; + + bve_mask_t& uncovered = scratch.uncovered; + uncovered.assign(n_words, 0u); + for (uint32_t m = 0; m < n_patterns; ++m) + if (!feas[m]) bve_mask_set(uncovered, m); + if (bve_mask_size(uncovered) == 0) return 0; // nothing to forbid + cuopt_assert(!primes.empty(), "infeasible patterns exist but no prime cube was enumerated"); + + scratch.cover.resize(primes.size()); + for (size_t q = 0; q < primes.size(); ++q) { + bve_cube_cover(primes[q].lit_mask, primes[q].bit_mask, full_mask, n_words, scratch.cover[q]); + // Zeroing the words, then one set-bit per pattern the cube matches. + ops += (int64_t)n_words + (int64_t{1} << (nb - std::popcount(primes[q].lit_mask))); + } + + int n = 0; + while (bve_mask_size(uncovered) > 0) { + // Per pick: the size test above, one bve_mask_overlap per prime, then the subtract below. + ops += (int64_t)((primes.size() + 2) * n_words); + int best_q = -1; + int best_gain = 0; + for (size_t q = 0; q < primes.size(); ++q) { + const int gain = bve_mask_overlap(uncovered, scratch.cover[q]); + if (gain > best_gain) { + best_gain = gain; + best_q = (int)q; + } + } + cuopt_assert(best_q >= 0, "prime cubes do not cover the infeasible patterns"); + if (n >= cap) return -1; + out[n++] = primes[best_q]; + bve_mask_subtract(uncovered, scratch.cover[best_q]); + } + ops += (int64_t)n_words; // the size test that ended the loop + cuopt_assert(n >= 1, "non-empty infeasible set covered by zero clauses"); + return n; +} + +// Committed elimination in commit order. `witness[pattern]` packs interior values for the boundary +// pattern; reductions are replayed in reverse order during postsolve. +template +struct bve_reduction_t { + std::vector interior; + std::vector boundary; + std::vector witness; // size 2^boundary.size() +}; + +// A surviving clause row to append to problem_t (a set-covering no-good over boundary columns). +// No upper bound: every one of these is a >= no-good, so the upper side is always +inf. +template +struct bve_added_row_t { + std::vector> terms; + f_t lower; +}; + +template +struct bve_plan_t { + std::vector> reductions; // commit order + std::vector removed_rows; // original row ids to drop + std::vector> added_rows; // surviving clause rows +}; + +// Working model and accumulated reduction plan. Candidates are staged without mutation and +// committed only after projection and clause validation. +template +struct bve_reducer_t { + struct work_row_t { + std::vector> terms; + f_t lo, up; + bool active; + }; + + i_t n_vars, n_rows_orig; + f_t tol; + i_t boundary_cap, scope_cap, clause_growth_margin; + std::vector rows; + std::vector> col2rows; + std::vector is_bin, obj_nz, done; + bve_plan_t plan; + bve_cover_scratch_t cover_scratch; + + bve_reducer_t(i_t n_vars_, + i_t n_rows_orig_, + const std::vector& offsets, + const std::vector& variables, + const std::vector& coefficients, + const std::vector& row_lower, + const std::vector& row_upper, + const std::vector& col_lower, + const std::vector& col_upper, + const std::vector& is_integer, + const std::vector& obj, + f_t tol_, + i_t boundary_cap_, + i_t scope_cap_, + i_t clause_growth_margin_); + + // Rows spanned by `interior` and the boundary columns of those rows, both unsorted, with op + // accounting. Single traversal behind both the growth probe (which needs only the boundary size) + // and stage(); outputs are overwritten, so a caller in a loop can reuse them. + void scope_of(const std::vector& interior, + std::vector& rows_out, + std::vector& boundary_out, + int64_t& ops) const; + + // Gather and pack a candidate without projecting or mutating the working model. + bool stage(const std::vector& interior_in, + bve_candidate_t& out, + int64_t* ops_out = nullptr); + + // Validate and commit an already-projected candidate; return true iff reduced. `ops_out` receives + // the CNF construction cost that bve_commit_wall_ops cannot predict. + bool commit_projected(const bve_candidate_t& cand, int64_t* ops_out = nullptr); + + bve_plan_t finalize(); +}; + +template +bve_reducer_t::bve_reducer_t(i_t n_vars_, + i_t n_rows_orig_, + const std::vector& offsets, + const std::vector& variables, + const std::vector& coefficients, + const std::vector& row_lower, + const std::vector& row_upper, + const std::vector& col_lower, + const std::vector& col_upper, + const std::vector& is_integer, + const std::vector& obj, + f_t tol_, + i_t boundary_cap_, + i_t scope_cap_, + i_t clause_growth_margin_) + : n_vars(n_vars_), + n_rows_orig(n_rows_orig_), + tol(tol_), + boundary_cap(boundary_cap_), + scope_cap(scope_cap_), + clause_growth_margin(clause_growth_margin_), + col2rows(n_vars_), + is_bin(n_vars_), + obj_nz(n_vars_), + done(n_vars_, 0) +{ + const f_t INF = std::numeric_limits::infinity(); + for (i_t c = 0; c < n_vars; ++c) { + is_bin[c] = + (is_integer[c] && std::abs(col_lower[c]) < tol && std::abs(col_upper[c] - f_t(1)) < tol) ? 1 + : 0; + obj_nz[c] = (obj[c] != f_t(0)) ? 1 : 0; + } + rows.reserve(n_rows_orig * 2); + for (i_t r = 0; r < n_rows_orig; ++r) { + work_row_t R; + R.active = true; + R.lo = scaling_bound_finite(row_lower[r]) ? row_lower[r] : -INF; + R.up = scaling_bound_finite(row_upper[r]) ? row_upper[r] : INF; + for (i_t k = offsets[r]; k < offsets[r + 1]; ++k) + R.terms.emplace_back(variables[k], coefficients[k]); + i_t id = rows.size(); + rows.push_back(std::move(R)); + for (auto& p : rows[id].terms) + col2rows[p.first].insert(id); + } +} + +template +void bve_reducer_t::scope_of(const std::vector& interior, + std::vector& rows_out, + std::vector& boundary_out, + int64_t& ops) const +{ + ops += (int64_t)interior.size(); + std::unordered_set interior_set(interior.begin(), interior.end()); + std::unordered_set affected_rows; + for (i_t a : interior) + for (i_t r : col2rows[a]) { + ++ops; + affected_rows.insert(r); + } + std::unordered_set b; + for (i_t r : affected_rows) + for (const auto& p : rows[r].terms) { + ++ops; + if (!interior_set.count(p.first)) b.insert(p.first); + } + rows_out.assign(affected_rows.begin(), affected_rows.end()); + boundary_out.assign(b.begin(), b.end()); +} + +// This is where the block leaves floating point behind: rescale every row to integer coefficients +// and bounds so the projection can run at tolerance 0. Returns false if any row does not scale to +// bounded integers, which rejects the whole block rather than risk a tolerance-sensitive +// feasibility misclassification on large or non-rational coefficients. Only the projection's +// private copy is scaled -- the block rows are dropped from the model and the appended no-goods are +// scale-independent +/-1 clauses, so this never perturbs the installed model. +template +static bool integerize_projection_rows(bve_block_t& block) +{ + for (int rr = 0; rr < block.n_rows; ++rr) { + const int rb = block.row_off[rr]; + const int re = block.row_off[rr + 1]; + const double s = row_int_scale(block.row_coef + rb, + re - rb, + block.row_lo[rr], + block.row_up[rr], + BVE_MAX_ROW_LEN, + BVE_INT_SCALE_MAX); + if (s == 0.0) return false; + for (int k = rb; k < re; ++k) + block.row_coef[k] = (f_t)std::llround((double)block.row_coef[k] * s); + if (scaling_bound_finite(block.row_lo[rr])) + block.row_lo[rr] = (f_t)std::llround((double)block.row_lo[rr] * s); + if (scaling_bound_finite(block.row_up[rr])) + block.row_up[rr] = (f_t)std::llround((double)block.row_up[rr] * s); + } + return true; +} + +template +bool bve_reducer_t::stage(const std::vector& interior_in, + bve_candidate_t& out, + int64_t* ops_out) +{ + int64_t ops = 0; + auto ops_guard = cuopt::scope_guard([&]() { + if (ops_out != nullptr) *ops_out += ops; + }); + + std::vector interior(interior_in.begin(), interior_in.end()); + std::sort(interior.begin(), interior.end()); + std::vector affected_rows, boundary; + scope_of(interior, affected_rows, boundary, ops); + // row order is result-invariant; sorting improves GPU shape-binning + std::sort(affected_rows.begin(), affected_rows.end()); + ops += (int64_t)affected_rows.size(); + std::sort(boundary.begin(), boundary.end()); + ops += (int64_t)boundary.size(); + + const i_t nb = boundary.size(); + const i_t na = interior.size(); + if (nb == 0 || nb > boundary_cap || na + nb > scope_cap) return false; + for (i_t v : boundary) + if (!is_bin[v]) return false; + if (na > BVE_MAX_INTERIOR || nb > BVE_MAX_BOUNDARY || na + nb > BVE_MAX_SCOPE) return false; + if (affected_rows.size() > BVE_MAX_ROWS) return false; + + bve_block_t& blk = out.blk; + blk.na = na; + blk.nb = nb; + blk.n_rows = affected_rows.size(); + std::unordered_map local; + for (i_t j = 0; j < na; ++j) + local[interior[j]] = j; + for (i_t j = 0; j < nb; ++j) + local[boundary[j]] = na + j; + ops += (int64_t)(na + nb); + i_t nzc = 0; + bool row_overflow = false; + for (i_t rr = 0; rr < blk.n_rows && !row_overflow; ++rr) { + const i_t r = affected_rows[rr]; + blk.row_off[rr] = nzc; + if (rows[r].terms.size() > BVE_MAX_ROW_LEN || nzc + rows[r].terms.size() > BVE_MAX_NNZ) { + row_overflow = true; + break; + } + for (auto& p : rows[r].terms) { + blk.row_var[nzc] = local[p.first]; + blk.row_coef[nzc] = p.second; + ++nzc; + ++ops; + } + blk.row_lo[rr] = rows[r].lo; + blk.row_up[rr] = rows[r].up; + } + if (row_overflow) return false; + blk.row_off[blk.n_rows] = nzc; + + if (!integerize_projection_rows(blk)) return false; + + out.interior = std::move(interior); + out.boundary = std::move(boundary); + out.rows = std::move(affected_rows); + out.projection.feasible.assign(size_t(1) << nb, 0); + out.projection.witness.assign(size_t(1) << nb, 0u); + ops += (int64_t)(1 << nb); + return true; +} + +template +bool bve_reducer_t::commit_projected(const bve_candidate_t& cand, + int64_t* ops_out) +{ + const int nb = cand.blk.nb; + const uint8_t* feasible = cand.projection.feasible.data(); + cuopt_assert(cand.projection.feasible.size() == (size_t(1) << nb), "projection table unsized"); + bve_clause_t clauses[BVE_MAX_CLAUSES]; + const int n_clauses = + bve_greedy_prime_cover(feasible, nb, clauses, BVE_MAX_CLAUSES, cover_scratch, ops_out); + if (n_clauses < 0) return false; // clause explosion past cap + if (n_clauses > cand.blk.n_rows + clause_growth_margin) return false; // growth gate + if (!bve_sanity_check(feasible, nb, clauses, n_clauses)) + return false; // sanity check failed => keep block + + bve_reduction_t red; + red.interior = cand.interior; + red.boundary = cand.boundary; + red.witness = cand.projection.witness; + plan.reductions.push_back(std::move(red)); + + for (i_t r : cand.rows) { + for (auto& p : rows[r].terms) + col2rows[p.first].erase(r); + rows[r].active = false; + rows[r].terms.clear(); + } + const f_t INF = std::numeric_limits::infinity(); + for (i_t ci = 0; ci < n_clauses; ++ci) { + const uint32_t lit = clauses[ci].lit_mask; + const uint32_t bit = clauses[ci].bit_mask; + work_row_t R; + R.active = true; + R.up = INF; + i_t n1 = 0; + for (i_t j = 0; j < nb; ++j) + if (lit & (1u << j)) { + const i_t b = (bit >> j) & 1u; + R.terms.emplace_back(cand.boundary[j], b ? f_t(-1) : f_t(1)); + n1 += b; + } + R.lo = f_t(1 - n1); + i_t id = rows.size(); + rows.push_back(std::move(R)); + for (auto& p : rows[id].terms) + col2rows[p.first].insert(id); + } + for (i_t a : cand.interior) { + col2rows[a].clear(); + done[a] = 1; + } + return true; +} + +template +bve_plan_t bve_reducer_t::finalize() +{ + for (i_t r = 0; r < n_rows_orig; ++r) + if (!rows[r].active) plan.removed_rows.push_back(r); + for (size_t r = n_rows_orig; r < rows.size(); ++r) + if (rows[r].active) { + cuopt_assert(rows[r].up == std::numeric_limits::infinity(), + "clause rows carry no upper bound"); + bve_added_row_t ar; + ar.terms = std::move(rows[r].terms); + ar.lower = rows[r].lo; + plan.added_rows.push_back(std::move(ar)); + } + return plan; +} + +// =========================================================================================== +// GPU enumeration projection kernel +// =========================================================================================== + +// Exact-enumeration projection kernel, laid out to fill the GPU: +// +// grid : one CTA per assignment (block, boundary pattern m, interior pattern am), +// grid-strided over CTAs ( for assignment = blockIdx.x; ...; += gridDim.x ) +// CTA : one warp per row ( blockDim.x == min(nrows,32)*32; warps loop if nrows > 32 ) +// warp : reduces sum = Σ coeff * value over the row's entries, tests sum in [lower, upper] +// +// The CTA ANDs the per-row satisfied bits into a single "assignment feasible" bit. For each +// boundary pattern m, feasibility is the OR over its interior patterns am and the witness is the +// first feasible am; both are encoded by a single atomicMin into `out_witness` (sentinel 0xFFFFFFFF +// = no feasible interior), so downstream: +// feasible[block][m] == (out_witness[block][m] != 0xFFFFFFFF) +// witness [block][m] == out_witness[block][m] // the smallest feasible interior +// `out_witness` must be initialized to 0xFFFFFFFF by the caller before launch. +// +// Shape (nb, na, nrows, and the row layout) is passed at RUNTIME, not as template parameters: it +// would otherwise need one instantiation per distinct shape. All blocks in a single launch share +// the shape (they are pre-binned), so every CTA still runs the identical loop structure. +// `row_start` and `local_var_of_entry` describe that shared layout; `nnz == row_start[nrows]`. +// `row_satisfied` uses dynamic shared memory of `nrows` bytes. +template +__global__ void bve_enumerate_kernel( + i_t num_blocks, + i_t nb, + i_t na, + i_t nrows, + f_t tolerance, + const f_t* block_coeffs, // [num_blocks * nnz] + const i_t* local_var_of_entry, // [nnz] (shared by the bin) + const i_t* row_start, // [nrows + 1] (shared by the bin) + const f_t* block_row_lower, // [num_blocks * nrows] + const f_t* block_row_upper, // [num_blocks * nrows] + uint32_t* out_witness) // [num_blocks * (1<> (na + nb)); + + const f_t* coeffs = block_coeffs + block * nnz; + const f_t* lower = block_row_lower + block * nrows; + const f_t* upper = block_row_upper + block * nrows; + + // one warp per row (a warp loops over multiple rows when nrows > num_warps) + for (i_t row = warp_id; row < nrows; row += num_warps) { + f_t partial = 0; + for (i_t entry = row_start[row] + lane_id; entry < row_start[row + 1]; entry += 32) { + const i_t var = local_var_of_entry[entry]; + const f_t value = (var < na) ? (f_t)((interior_pattern >> var) & 1) + : (f_t)((boundary_pattern >> (var - na)) & 1); + partial += coeffs[entry] * value; + } + // Lane 0 holds the result for both XOR-butterfly and typical down-sweep warp reduces. + const f_t sum = raft::warpReduce(partial); + if (lane_id == 0) { + row_satisfied[row] = + (sum <= upper[row] + tolerance && sum >= lower[row] - tolerance) ? 1 : 0; + } + } + __syncthreads(); + + // AND the per-row bits; if this assignment is feasible, offer its interior as a witness + if (threadIdx.x == 0) { + uint8_t feasible = 1; + for (i_t row = 0; row < nrows; ++row) { + feasible &= row_satisfied[row]; + } + if (feasible) { + atomicMin(&out_witness[block * num_patterns + boundary_pattern], + (uint32_t)interior_pattern); + } + } + __syncthreads(); // guard row_satisfied before the next assignment overwrites it + } +} + +// ---- GPU batch projection: one enumeration-kernel launch per shape-bin ---- +// Returns raw work for the enumerations (sum over bins of assignments · nnz). +template +double bve_project_batch_gpu(const raft::handle_t& handle, + std::vector>& cands, + f_t tol) +{ + if (cands.empty()) return 0.0; + auto stream = handle.get_stream(); + double work_units = 0.0; + + // Bin candidates by identical shape so every CTA in a launch runs the same loop structure. The + // key is (na, nb, n_rows, nnz, row_off[...], row_var[...]) — everything the kernel reads as + // shared; only the coefficients and row bounds differ per block. Hash map avoids O(key_len · + // log n_bins) tree compares on long keys (up to ~1605 ints at the BVE caps). + struct shape_key_hash { + size_t operator()(const std::vector& key) const + { + size_t h = 0; + for (i_t x : key) { + h ^= std::hash{}(x) + 0x9e3779b9 + (h << 6) + (h >> 2); + } + return h; + } + }; + std::unordered_map, std::vector, shape_key_hash> bins; + for (size_t i = 0; i < cands.size(); ++i) { + const auto& blk = cands[i].blk; + const i_t nnz = blk.row_off[blk.n_rows]; + std::vector key; + key.reserve(4 + (blk.n_rows + 1) + nnz); + key.push_back(blk.na); + key.push_back(blk.nb); + key.push_back(blk.n_rows); + key.push_back(nnz); + for (i_t r = 0; r <= blk.n_rows; ++r) + key.push_back(blk.row_off[r]); + for (i_t k = 0; k < nnz; ++k) + key.push_back(blk.row_var[k]); + bins[std::move(key)].push_back(i); + } + + for (const auto& kv : bins) { + const std::vector& idxs = kv.second; + const auto& proto = cands[idxs[0]].blk; + const i_t na = proto.na; + const i_t nb = proto.nb; + const i_t nrows = proto.n_rows; + const i_t nnz = proto.row_off[nrows]; + const i_t patterns = i_t(1) << nb; + + // Shared layout is O(nnz) and identical for every candidate in the bin. + std::vector h_row_start(proto.row_off, proto.row_off + nrows + 1); + std::vector h_local_var(proto.row_var, proto.row_var + nnz); + rmm::device_uvector d_row_start(h_row_start.size(), stream); + rmm::device_uvector d_local_var(h_local_var.size(), stream); + raft::copy(d_row_start.data(), h_row_start.data(), h_row_start.size(), stream); + raft::copy(d_local_var.data(), h_local_var.data(), h_local_var.size(), stream); + + // Per-block device cost: coeffs + row bounds + witness table. + const size_t bytes_per_block = size_t(nnz) * sizeof(f_t) + 2 * size_t(nrows) * sizeof(f_t) + + size_t(patterns) * sizeof(uint32_t); + // Also clamp to i_t range: the kernel takes num_blocks as i_t. + const size_t chunk = + std::max(1, + std::min(size_t(std::numeric_limits::max()), + BVE_PROJECT_DEVICE_BUDGET / std::max(1, bytes_per_block))); + + // Launch dims are CUDA `int` by API. + const int num_warps = std::min(nrows, 32); + const int cta_dim = num_warps * 32; + const size_t shmem = size_t(nrows) * sizeof(uint8_t); + + for (size_t offset = 0; offset < idxs.size(); offset += chunk) { + const size_t num_sz = std::min(chunk, idxs.size() - offset); + const i_t num = num_sz; + + std::vector h_coeffs(num_sz * size_t(nnz)); + std::vector h_lower(num_sz * size_t(nrows)); + std::vector h_upper(num_sz * size_t(nrows)); + for (size_t g = 0; g < num_sz; ++g) { + const auto& blk = cands[idxs[offset + g]].blk; + std::copy(blk.row_coef, blk.row_coef + nnz, h_coeffs.begin() + g * nnz); + std::copy(blk.row_lo, blk.row_lo + nrows, h_lower.begin() + g * nrows); + std::copy(blk.row_up, blk.row_up + nrows, h_upper.begin() + g * nrows); + } + + rmm::device_uvector d_coeffs(h_coeffs.size(), stream); + rmm::device_uvector d_lower(h_lower.size(), stream); + rmm::device_uvector d_upper(h_upper.size(), stream); + rmm::device_uvector d_witness(num_sz * size_t(patterns), stream); + raft::copy(d_coeffs.data(), h_coeffs.data(), h_coeffs.size(), stream); + raft::copy(d_lower.data(), h_lower.data(), h_lower.size(), stream); + raft::copy(d_upper.data(), h_upper.data(), h_upper.size(), stream); + // sentinel 0xFFFFFFFF (every byte 0xFF) marks a boundary pattern with no feasible interior + // yet + RAFT_CUDA_TRY( + cudaMemsetAsync(d_witness.data(), 0xFF, d_witness.size() * sizeof(uint32_t), stream)); + + // one warp per row, one CTA per (block, m, am) assignment, grid-strided + const int64_t total = (int64_t)num * (int64_t)patterns * ((int64_t)1 << na); + const int grid = std::min(total, int64_t{65535}); + bve_enumerate_kernel<<>>(num, + nb, + na, + nrows, + tol, + d_coeffs.data(), + d_local_var.data(), + d_row_start.data(), + d_lower.data(), + d_upper.data(), + d_witness.data()); + RAFT_CUDA_TRY(cudaGetLastError()); + + // Unscaled op counts: host pack/unpack touches + one coeff read per assignment. + work_units += double(num_sz) * double(nnz + 2 * nrows + patterns); + work_units += double(total) * double(nnz); + + std::vector h_witness(num_sz * size_t(patterns)); + raft::copy(h_witness.data(), d_witness.data(), h_witness.size(), stream); + handle.sync_stream(); + for (size_t g = 0; g < num_sz; ++g) { + auto& cand = cands[idxs[offset + g]]; + // No-op for anything stage() produced; sizes a caller that assembled `blk` by hand. + cand.projection.feasible.resize(patterns); + cand.projection.witness.resize(patterns); + for (i_t m = 0; m < patterns; ++m) { + const uint32_t w = h_witness[g * patterns + m]; + const bool feasible = (w != 0xFFFFFFFFu); + cand.projection.feasible[m] = feasible ? 1 : 0; + cand.projection.witness[m] = feasible ? w : 0u; + } + } + } + } + return work_units; +} + +// ---- harvest unary-conditioned implications from an exactly projected block ---- +// +// `feas` is the block's exact existential projection onto its nb boundary columns, so for boundary +// position j and value a the feasible patterns agreeing with (j == a) describe every completion the +// block admits. Intersecting them (AND) gives the positions forced to 1 and the complement of their +// union (OR) gives those forced to 0; the same reasoning with no condition gives unconditional +// fixings. This is complete for the block's rows, where the probing cache only holds what bound +// propagation could prove, so these forcings can be strictly stronger. It holds whether or not the +// block is eventually eliminated, hence the call site harvests before the growth gate can reject. +// +// Ids are emitted in the current-problem frame; the caller maps them to original ids. +template +static void bve_extract_forcings(const bve_candidate_t& cand, probe_findings_t& out) +{ + const i_t nb = cand.blk.nb; + cuopt_assert(nb > 0 && nb <= BVE_MAX_BOUNDARY, "boundary width out of range"); + cuopt_assert((i_t)cand.boundary.size() == nb, "boundary id count disagrees with block width"); + const uint32_t n_patterns = 1u << nb; + + // Accumulators for condition s = 2*j + a; slot 2*nb holds the unconditional case. + constexpr i_t n_slots = 2 * BVE_MAX_BOUNDARY + 1; + const i_t unconditional = 2 * nb; + const uint32_t all_ones = n_patterns - 1u; + uint32_t and_acc[n_slots]; + uint32_t or_acc[n_slots]; + std::fill_n(and_acc, unconditional + 1, all_ones); + std::fill_n(or_acc, unconditional + 1, 0u); + + uint32_t n_feasible = 0; + for (uint32_t m = 0; m < n_patterns; ++m) { + if (!cand.projection.feasible[m]) continue; + ++n_feasible; + and_acc[unconditional] &= m; + or_acc[unconditional] |= m; + for (i_t j = 0; j < nb; ++j) { + const i_t s = 2 * j + ((m >> j) & 1u); + and_acc[s] &= m; + or_acc[s] |= m; + } + } + // Vacuous accumulators would otherwise read as "every position forced to 1". + if (n_feasible == 0u) return; // block alone is infeasible; left to the bound presolve + + // Positions the block fixes outright. Conditioning on anything would re-derive these, so they are + // recorded once here and skipped below. A position outside the mask takes both values among the + // feasible patterns, so both of its condition slots are non-empty. + const uint32_t fixed_mask = and_acc[unconditional] | (~or_acc[unconditional] & all_ones); + for (i_t j = 0; j < nb; ++j) { + if (!(fixed_mask & (1u << j))) continue; + out.fixings.emplace_back(cand.boundary[j], ((and_acc[unconditional] >> j) & 1u) != 0u); + } + + for (i_t j = 0; j < nb; ++j) { + if (fixed_mask & (1u << j)) continue; // condition never binds + for (i_t a = 0; a < 2; ++a) { + const i_t s = 2 * j + a; + for (i_t k = 0; k < nb; ++k) { + const uint32_t bit = 1u << k; + if (k == j || (fixed_mask & bit)) continue; + if (and_acc[s] & bit) { + out.forcings.push_back({cand.boundary[j], cand.boundary[k], a != 0, true}); + } else if (!(or_acc[s] & bit)) { + out.forcings.push_back({cand.boundary[j], cand.boundary[k], a != 0, false}); + } + } + } + } +} + +template +struct bve_growth_result_t { + std::vector interior; // sorted current-problem column ids, always contains the seed + int64_t ops = 0; // work performed, for the deterministic wall estimate +}; + +// Grows one seed into a block interior: starting from {seed}, repeatedly absorb the eligible +// implication-neighbor that shrinks the boundary the most, stopping when no neighbor strictly +// improves it or a cap is hit. Read-only on `reducer`, which is what lets the round run this across +// seeds under OpenMP against a frozen model. +template +static bve_growth_result_t grow_seed_interior( + i_t seed, + const bve_reducer_t& reducer, + const std::vector>& implication_adjacency) +{ + auto has_adj = [&](i_t v) { + return v >= 0 && v < (i_t)implication_adjacency.size() && !implication_adjacency[v].empty(); + }; + auto eligible = [&](i_t w) { + return reducer.is_bin[w] && !reducer.obj_nz[w] && !reducer.done[w] && + !reducer.col2rows[w].empty(); + }; + + bve_growth_result_t result; + std::unordered_set interior_set = {seed}; + std::vector probe_rows, probe_bnd; // scope_of scratch, reused across probes + for (;;) { + // Hub fast-path: raw implication degree upper-bounds |cands_w|. Skip boundary walks + // and adj materialization when the neighborhood is past the probe cap. + if (interior_set.size() == 1) { + const i_t s = *interior_set.begin(); + const i_t deg = has_adj(s) ? (i_t)implication_adjacency[s].size() : 0; + if (deg > BVE_MAX_GROWTH_NBRS) break; + } + std::vector candidate_interior(interior_set.begin(), interior_set.end()); + reducer.scope_of(candidate_interior, probe_rows, probe_bnd, result.ops); + const i_t cur = probe_bnd.size(); + // Implication-neighbors of the interior that are still eligible to enter it. + std::unordered_set cands_w; + bool gated = false; + for (i_t a : interior_set) { + if (!has_adj(a)) continue; + for (i_t w : implication_adjacency[a]) { + ++result.ops; + if (interior_set.count(w) || !eligible(w)) continue; + cands_w.insert(w); + if ((i_t)cands_w.size() > BVE_MAX_GROWTH_NBRS) { + gated = true; + break; + } + } + if (gated) break; + } + // Hub neighborhoods: full probe is Θ(|cands_w|) boundary walks and rarely absorbs. + if (gated) break; + // Pick the neighbor with the smallest boundary; stop when none strictly improves. + i_t best = -1; + i_t best_nb = cur; + for (i_t w : cands_w) { + candidate_interior.push_back(w); // probe interior ∪ {w}; the pop below restores it + const i_t na = candidate_interior.size(); + reducer.scope_of(candidate_interior, probe_rows, probe_bnd, result.ops); + const i_t nb = probe_bnd.size(); + candidate_interior.pop_back(); + if (nb < best_nb && na + nb <= reducer.scope_cap && na <= BVE_MAX_INTERIOR) { + best_nb = nb; + best = w; + } + } + if (best < 0) break; + interior_set.insert(best); + } + result.interior.assign(interior_set.begin(), interior_set.end()); + return result; +} + +// ---- production detector: round-based, scope-disjoint, one GPU projection launch per round ---- +// +// Implication-closure block growth over the probing-cache adjacency: each seed absorbs the +// implication-neighbor that most shrinks its boundary (subject to enum/interior caps) until no +// such neighbor remains. Restructured so many candidate blocks are projected in ONE GPU launch. +// Within a round the working model is FROZEN — every seed grows its interior against the same +// model. Because that growth is read-only on the model, it runs in an OpenMP parallel-for across +// the round's seeds; the results are deterministic per seed and acceptance is then applied +// serially in seed order, so the committed plan is identical to a serial run of the same frozen +// growth. Candidates are staged and only mutually SCOPE-DISJOINT ones (no shared interior or +// boundary column, which also forbids a shared row) are accepted into the batch. The batch is +// projected on the device (bve_project_batch_gpu), then committed on the host; because the accepted +// candidates touch disjoint columns/rows, commit order is irrelevant and each block's staged +// projection is still valid at commit time. Candidates deferred for overlap are retried in later +// rounds; the loop stops when a round accepts nothing or commits nothing (each committing round +// retires >= 1 column => terminates). The scope-disjoint rule is deliberately conservative (it also +// rejects candidates that merely share a boundary column, which would be safe); relax it if +// per-round batch sizes prove too small. TU-local (only the pass uses it). +template +static bve_plan_t bve_detect_closure_batched( + const raft::handle_t& handle, + bve_reducer_t& reducer, + const std::vector>& impl_adj, + timer_t& timer, + double& work_units, + probe_findings_t* findings) +{ + auto has_adj = [&](i_t v) { return v >= 0 && v < (i_t)impl_adj.size() && !impl_adj[v].empty(); }; + std::vector order; + for (i_t c = 0; c < reducer.n_vars; ++c) + if (reducer.is_bin[c] && !reducer.obj_nz[c] && !reducer.col2rows[c].empty() && has_adj(c)) + order.push_back(c); + std::sort(order.begin(), order.end(), [&](i_t a, i_t b) { + return reducer.col2rows[a].size() < reducer.col2rows[b].size(); + }); + + std::vector attempted(reducer.n_vars, + 0); // a seed is attempted once (whether or not it commits) + // Grow each seed at most once; overlap-deferred seeds only re-stage from the cached interior. + // Re-growing hubs every round dominated wall; retiring them on first overlap killed reductions. + std::vector growth_done(reducer.n_vars, 0); + std::vector> growth_interior(reducer.n_vars); + for (;;) { + if (timer.check_time_limit()) break; + + // This round's live seeds, in the deterministic growth order. + std::vector round_seeds; + for (i_t seed : order) + if (!attempted[seed] && !reducer.done[seed] && !reducer.col2rows[seed].empty()) + round_seeds.push_back(seed); + if (round_seeds.empty()) break; + + // Grow each seed against the frozen model (read-only on reducer → OMP-safe). Acceptance below + // is serial in round_seeds order, so the plan matches a serial frozen-growth run. + std::vector> interiors(round_seeds.size()); + std::vector growth_ops(round_seeds.size(), 0); +#pragma omp taskloop default(shared) priority(CUOPT_DEFAULT_TASK_PRIORITY) + for (i_t k = 0; k < (i_t)round_seeds.size(); ++k) { + const i_t seed = round_seeds[k]; + if (growth_done[seed]) { + interiors[k] = growth_interior[seed]; + continue; + } + bve_growth_result_t grown = grow_seed_interior(seed, reducer, impl_adj); + growth_ops[k] = grown.ops; + interiors[k] = std::move(grown.interior); + growth_interior[seed] = interiors[k]; + growth_done[seed] = 1; + } + // OMP growth: wall ≈ critical-path seed (max), not sum across threads. + int64_t max_growth_ops = 0; + for (int64_t ops : growth_ops) + max_growth_ops = std::max(max_growth_ops, ops); + work_units += double(max_growth_ops); + + if (timer.check_time_limit()) break; + + // Serial: stage each grown interior and greedily accept mutually SCOPE-DISJOINT candidates, in + // round_seeds order. Nothing mutates the model until commit, so this stays serial. + std::vector> cands; + std::unordered_set claimed; // interior+boundary columns of already-accepted candidates + for (size_t k = 0; k < round_seeds.size(); ++k) { + if (timer.check_time_limit()) break; + const i_t seed = round_seeds[k]; + bve_candidate_t cand; + int64_t stage_ops = 0; + if (!reducer.stage(interiors[k], cand, &stage_ops)) { + work_units += double(stage_ops); + attempted[seed] = + 1; // failed the caps against this model; treat as one touch, like sequential + continue; + } + work_units += double(stage_ops); + bool overlap = false; + for (i_t c : cand.interior) + if (claimed.count(c)) { + overlap = true; + break; + } + if (!overlap) + for (i_t c : cand.boundary) + if (claimed.count(c)) { + overlap = true; + break; + } + if (overlap) continue; // scope collides; retry stage later from cached interior + + attempted[seed] = 1; + for (i_t c : cand.interior) + claimed.insert(c); + for (i_t c : cand.boundary) + claimed.insert(c); + cands.push_back(std::move(cand)); + } + + if (cands.empty() || timer.check_time_limit()) break; + // Staged blocks are integerized (integerize_projection_rows), so the subset-sum feasibility + // test is exact: project with tolerance 0 rather than reducer.tol. + work_units += bve_project_batch_gpu(handle, cands, f_t(0)); + if (timer.check_time_limit()) break; + i_t committed = 0; + for (auto& cand : cands) { + if (timer.check_time_limit()) break; + // Valid for the block's rows regardless of the clause gates below, so harvest before them. + if (findings != nullptr) { + bve_extract_forcings(cand, *findings); + work_units += double(uint32_t(1) << cand.blk.nb) * double(cand.blk.nb); + } + work_units += + bve_commit_wall_ops(cand.blk.nb, cand.blk.n_rows + reducer.clause_growth_margin); + int64_t commit_ops = 0; + if (reducer.commit_projected(cand, &commit_ops)) ++committed; + work_units += double(commit_ops); + } + if (committed == 0) break; + } + return reducer.finalize(); +} + +// ---- implication adjacency from the probing cache (original-id -> current column) ---- +template +std::vector> bve_build_impl_adj( + const probing_cache_t& cache, + const std::vector& reverse_original_ids, + i_t n_vars, + const probe_findings_t* prior_original_id_findings) +{ + // original-id -> current column index (or -1 if the column no longer exists) + auto to_current = [&](i_t original_id) -> i_t { + if (original_id < 0 || original_id >= (i_t)reverse_original_ids.size()) return -1; + return reverse_original_ids[original_id]; + }; + std::vector> adj(n_vars); + auto add_edge = [&](i_t original_x, i_t original_y) { + const i_t x = to_current(original_x); + if (x < 0 || x >= n_vars) return; + const i_t y = to_current(original_y); + if (y < 0 || y >= n_vars || y == x) return; + adj[x].insert(y); + adj[y].insert(x); + }; + for (const auto& kv : cache.probing_cache) { + for (int p = 0; p < 2; ++p) { + for (const auto& yb : kv.second[p].var_to_cached_bound_map) + add_edge(kv.first, yb.first); + } + } + // Forcings mined from earlier projections. Pairs the cache never held become seed/absorb + // candidates, so a later round can grow blocks the first round could not see. + if (prior_original_id_findings != nullptr) { + for (const auto& forcing : prior_original_id_findings->forcings) + add_edge(forcing.var, forcing.forced_var); + } + std::vector> out(n_vars); + for (i_t v = 0; v < n_vars; ++v) + out[v].assign(adj[v].begin(), adj[v].end()); + return out; +} + +// Records every committed block on the unified append-only reconstruction log, translating +// detection-space column ids into the post-Papilo frame that postsolve replays in reverse. Commit +// order is preserved, which is what makes the reverse replay well-defined. +template +static void append_bve_reconstructions(const bve_plan_t& plan, + const std::vector& current_to_post_papilo, + presolve_data_t& presolve_data, + double& work_units) +{ + auto to_post_papilo = [&](i_t column) { + cuopt_assert(column >= 0 && column < (i_t)current_to_post_papilo.size(), + "block column out of variable_mapping range"); + return current_to_post_papilo[column]; + }; + + auto& recs = presolve_data.var_postsolve; + recs.reserve(recs.size() + plan.reductions.size()); + for (const auto& red : plan.reductions) { + work_units += double(red.interior.size() + red.boundary.size() + red.witness.size()); + var_postsolve_t rec; + rec.kind = reconstruction_kind_t::BlockBve; + rec.bve.interior.reserve(red.interior.size()); + for (i_t c : red.interior) + rec.bve.interior.push_back(to_post_papilo(c)); + rec.bve.boundary.reserve(red.boundary.size()); + for (i_t c : red.boundary) + rec.bve.boundary.push_back(to_post_papilo(c)); + rec.bve.witness = red.witness; + recs.push_back(std::move(rec)); + } +} + +// ---- the pass: detect (GPU-projected) -> install reduced model -> record reconstructions ---- +template +bool block_bve_presolve(problem_t& problem, + const std::vector>& impl_adj, + timer_t& timer, + double& work_units, + probe_findings_t* out_findings, + i_t boundary_cap, + i_t scope_cap, + i_t clause_growth_margin) +{ + work_units = 0.0; + // Local wall clock for the DEBUG total; `timer` is the caller's stage deadline. + timer_t wall(std::numeric_limits::infinity()); + [[maybe_unused]] double t_setup = 0.0, t_detect = 0.0, t_install = 0.0, t_compact = 0.0; + auto timer_raii_guard = cuopt::scope_guard([&]() { + CUOPT_LOG_DEBUG( + "Block-BVE phases: setup=%.2fs detect=%.2fs install=%.2fs compact=%.2fs total=%.2fs " + "work units: %.6g", + t_setup, + t_detect, + t_install, + t_compact, + wall.elapsed_time(), + work_units); + }); + + const raft::handle_t* handle = problem.handle_ptr; + auto stream = handle->get_stream(); + const i_t n_vars = problem.n_variables; + const i_t n_rows = problem.n_constraints; + const f_t tol = problem.tolerances.presolve_absolute_tolerance; + if (problem.empty || n_vars == 0 || n_rows == 0) return false; + + // ---- 1. host copy of the current (post-Papilo, post-initial-trivial-presolve) model ---- + auto h_off = cuopt::host_copy(problem.offsets, stream); + auto h_var = cuopt::host_copy(problem.variables, stream); + auto h_coef = cuopt::host_copy(problem.coefficients, stream); + auto h_clb = cuopt::host_copy(problem.constraint_lower_bounds, stream); + auto h_cub = cuopt::host_copy(problem.constraint_upper_bounds, stream); + auto h_vb = cuopt::host_copy(problem.variable_bounds, stream); + auto h_vtype = cuopt::host_copy(problem.variable_types, stream); + auto h_obj = cuopt::host_copy(problem.objective_coefficients, stream); + // variable_mapping maps current-space column -> post-Papilo index (the frame postsolve uses) + auto h_vmap = cuopt::host_copy(problem.presolve_data.variable_mapping, stream); + handle->sync_stream(); + + // Host mirror + reducer construction (each walks the CSR once). + const i_t nnz0 = (i_t)h_off.back(); + work_units = double(2 * nnz0) + double(2 * n_vars) + double(n_rows); + + if (timer.check_time_limit()) return false; + + // ---- 2. detector inputs (i_t CSR, f_t bounds/coeffs) ---- + std::vector offsets(h_off.begin(), h_off.end()); + std::vector variables(h_var.begin(), h_var.end()); + std::vector coefficients(h_coef.begin(), h_coef.end()); + std::vector row_lower(h_clb.begin(), h_clb.end()); + std::vector row_upper(h_cub.begin(), h_cub.end()); + std::vector col_lower(n_vars), col_upper(n_vars); + std::vector is_integer(n_vars); + for (i_t c = 0; c < n_vars; ++c) { + col_lower[c] = get_lower(h_vb[c]); + col_upper[c] = get_upper(h_vb[c]); + is_integer[c] = (h_vtype[c] == var_t::INTEGER) ? 1 : 0; + } + std::vector obj(h_obj.begin(), h_obj.end()); + + // ---- 3. detect + sanity check (probing-cache implication closure). Projection of each candidate + // block runs on the GPU: the batched detector stages scope-disjoint candidates per round and + // hands the whole batch to bve_project_batch_gpu (one enumeration-kernel launch per shape-bin), + // which fills feas/witness; commit (prime-implicate CNF + inline sanity check) then runs on the + // host. ---- + bve_reducer_t reducer(n_vars, + n_rows, + offsets, + variables, + coefficients, + row_lower, + row_upper, + col_lower, + col_upper, + is_integer, + obj, + tol, + boundary_cap, + scope_cap, + clause_growth_margin); + t_setup = wall.elapsed_time(); + probe_findings_t current_id_findings; + bve_plan_t plan = + bve_detect_closure_batched(*handle, + reducer, + impl_adj, + timer, + work_units, + out_findings != nullptr ? ¤t_id_findings : nullptr); + t_detect = wall.elapsed_time() - t_setup; + + // Projection findings hold for the block's rows whether or not the block was eliminated, so they + // are exported before the no-reduction exit; the rejected blocks are often the interesting ones. + if (out_findings != nullptr) { + auto to_original = [&](i_t column) { + cuopt_assert(column >= 0 && column < (i_t)h_vmap.size(), "column outside variable_mapping"); + return (i_t)h_vmap[column]; + }; + out_findings->forcings.reserve(out_findings->forcings.size() + + current_id_findings.forcings.size()); + for (const auto& forcing : current_id_findings.forcings) { + out_findings->forcings.push_back({to_original(forcing.var), + to_original(forcing.forced_var), + forcing.value, + forcing.forced_value}); + } + for (const auto& [column, value] : current_id_findings.fixings) + out_findings->fixings.emplace_back(to_original(column), value); + } + + if (plan.reductions.empty()) return false; + + // ---- 4. build the reduced forward CSR: keep original rows not removed, append clause rows ---- + const double t_install_begin = wall.elapsed_time(); + std::vector removed(n_rows, 0); + for (i_t r : plan.removed_rows) + removed[r] = 1; + std::vector new_off, new_var; + std::vector new_coef, new_clb, new_cub; + new_off.reserve(n_rows + plan.added_rows.size() + 1); + new_off.push_back(0); + for (i_t r = 0; r < n_rows; ++r) { + if (removed[r]) continue; + for (i_t k = offsets[r]; k < offsets[r + 1]; ++k) { + new_var.push_back(variables[k]); + new_coef.push_back(coefficients[k]); + } + new_off.push_back(new_var.size()); + new_clb.push_back(row_lower[r]); + new_cub.push_back(row_upper[r]); + } + for (const auto& ar : plan.added_rows) { + for (const auto& [var, coef] : ar.terms) { + new_var.push_back(var); + new_coef.push_back(coef); + } + new_off.push_back(new_var.size()); + new_clb.push_back(ar.lower); // eliminated interior cols become empty (only in removed rows) + // clause rows are >= no-goods; upper is +inf (problem_t convention) + new_cub.push_back(std::numeric_limits::infinity()); + } + // ---- 5. install the rewritten rows into problem_t (matrix + derived state) ---- + work_units += double(new_var.size()) + double(new_clb.size()); + problem.set_constraints_from_host_csr(new_off, new_var, new_coef, new_clb, new_cub, {}); + + // ---- 6. record reconstructions ---- + append_bve_reconstructions(plan, h_vmap, problem.presolve_data, work_units); + t_install = wall.elapsed_time() - t_install_begin; + + // ---- 7. compact the now-empty interior columns and update variable_mapping ---- + const double t_compact_begin = wall.elapsed_time(); + work_units += double(n_vars) + double(new_var.size()); + trivial_presolve(problem, /*remap_cache_ids=*/true); + handle->sync_stream(); + t_compact = wall.elapsed_time() - t_compact_begin; + const i_t reduced_cols = n_vars - problem.n_variables; + const i_t reduced_rows = n_rows - problem.n_constraints; + if (reduced_cols > 0 || reduced_rows > 0) { + CUOPT_LOG_DEBUG("Block-BVE reduced %d columns, %d rows", reduced_cols, reduced_rows); + } +#if (CUOPT_LOG_ACTIVE_LEVEL <= RAPIDS_LOGGER_LOG_LEVEL_DEBUG) + const i_t fractional_coefs = + thrust::count_if(handle->get_thrust_policy(), + problem.coefficients.begin(), + problem.coefficients.end(), + [] __device__(f_t v) -> bool { return floor(v) != v; }); + CUOPT_LOG_DEBUG("Block-BVE: %d fractional coefficients in A", fractional_coefs); +#endif + return true; +} + +#define INSTANTIATE(F_TYPE) \ + template double bve_project_batch_gpu( \ + const raft::handle_t&, std::vector>&, F_TYPE); \ + template std::vector> bve_build_impl_adj( \ + const probing_cache_t&, \ + const std::vector&, \ + int, \ + const probe_findings_t*); \ + template bool block_bve_presolve(problem_t&, \ + const std::vector>&, \ + timer_t&, \ + double&, \ + probe_findings_t*, \ + int, \ + int, \ + int) + +INSTANTIATE(double); +#ifdef MIP_INSTANTIATE_FLOAT +INSTANTIATE(float); +#endif +#undef INSTANTIATE + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/block_bve.cuh b/cpp/src/mip_heuristics/presolve/block_bve.cuh new file mode 100644 index 0000000000..355ff00a8d --- /dev/null +++ b/cpp/src/mip_heuristics/presolve/block_bve.cuh @@ -0,0 +1,145 @@ +/* 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 "probing_cache.cuh" + +#include + +#include +#include + +// Eliminates small blocks of zero-objective binary variables. A block is a set of columns to remove +// (the interior, na columns) together with every row they appear in; the other columns of those +// rows are the boundary (nb columns), which stays in the model and must also be binary. +// +// For each of the 2^nb boundary assignments the projection decides whether some interior +// assignment satisfies the block's rows. The ruled-out assignments are everything the block still +// forces on the rest of the model, so emitting them as prime-implicate no-goods over the boundary +// carries that force without the interior. Committing therefore deletes the interior columns and +// every block row, installing the no-goods in their place: interior variables disappear and the row +// count drops whenever the no-goods are fewer than the rows they replace, which the growth gate +// below requires. One feasible interior witness per surviving assignment is stored so postsolve +// can rebuild the deleted columns; since the interior carries no objective coefficients, any +// witness preserves the objective as well as feasibility. +// +// Candidate interiors are grown from the probing implication graph and committed only when the +// projected CNF satisfies the bounded-elimination growth limit of Eén and Biere, "Effective +// Preprocessing in SAT through Variable and Clause Elimination" (SAT 2005). Before commit, the +// emitted clauses are checked against the GPU-computed boundary feasibility table. + +namespace cuopt::mathematical_optimization::mip { + +// Caps for a single enumerated block. +static constexpr int BVE_MAX_BOUNDARY = 12; // nb <= 12 => 2^nb <= 4096 feasibility patterns +static constexpr int BVE_MAX_SCOPE = 16; // na + nb <= 16 +static constexpr int BVE_MAX_ROWS = 64; // rows spanned by the block; #clauses <= #rows +static constexpr int BVE_MAX_ROW_LEN = 24; // nnz within one block row (interior+boundary entries) +static constexpr int BVE_MAX_NNZ = BVE_MAX_ROWS * BVE_MAX_ROW_LEN; +static constexpr int BVE_MAX_CLAUSES = 64; // <= |rows| for any committed block +static constexpr int BVE_MAX_PATTERNS = 1 << BVE_MAX_BOUNDARY; + +// Packed projection block. Local ids [0, na) are interior and [na, na+nb) are boundary; rows use +// CSR layout and missing bounds are +/- infinity. +template +struct bve_block_t { + int na; // number of interior variables + int nb; // number of boundary variables + int n_rows; // rows spanned by the block + int row_off[BVE_MAX_ROWS + 1]; + int row_var[BVE_MAX_NNZ]; // local var id in [0, na+nb) + f_t row_coef[BVE_MAX_NNZ]; + f_t row_lo[BVE_MAX_ROWS]; // -inf if no lower bound + f_t row_up[BVE_MAX_ROWS]; // +inf if no upper bound +}; + +// Boundary clause forbidding patterns that match `bit_mask` at every position in `lit_mask`. +// It is emitted as sum_j (bit_j == 0 ? x_j : -x_j) >= 1 - popcount(bit_mask & lit_mask). +struct bve_clause_t { + uint32_t lit_mask; + uint32_t bit_mask; +}; + +// One bit per boundary pattern. The width tracks the block's own 2^nb, not BVE_MAX_PATTERNS, so +// raising BVE_MAX_BOUNDARY costs nothing on narrower blocks. +using bve_mask_t = std::vector; + +// Buffers the CNF construction reuses across blocks: `valid` alone is 4^nb bytes, so per-block +// allocation would dominate at wide boundaries. +struct bve_cover_scratch_t { + std::vector valid; // prime-cube validity table, grow-only + std::vector primes; + std::vector cover; // patterns matched by each prime + bve_mask_t uncovered; +}; + +// Derive a prime-implicate CNF from the boundary feasibility table by covering the infeasible +// patterns with a max-gain greedy over every prime forbidden cube; return -1 on cap overflow. +// Untemplated: the CNF is a Boolean computation over the feasibility table, and every dimension it +// touches is capped by the BVE_MAX_* constants above. +int bve_greedy_prime_cover(const uint8_t* feas, + int nb, + bve_clause_t* out, + int cap, + bve_cover_scratch_t& scratch, + int64_t* ops_out = nullptr); + +// Verify that the emitted clauses reproduce the boundary feasibility table exactly. +bool bve_sanity_check(const uint8_t* feas, int nb, const bve_clause_t* clauses, int n_clauses); + +// Exact existential projection of one block onto its boundary, filled by the projection backend. +// Both tables are sized to the block's own 2^nb rather than BVE_MAX_PATTERNS, so a narrow block +// does not carry the cost of raising BVE_MAX_BOUNDARY. +struct bve_projection_t { + std::vector feasible; // [2^nb] 1 iff the boundary pattern admits some interior + std::vector witness; // [2^nb] smallest feasible interior, 0 where infeasible +}; + +// Staged candidate. Vector fields use sorted current-problem ids; `blk` uses local ids. +template +struct bve_candidate_t { + std::vector interior; // sorted global column ids (to be eliminated) + std::vector boundary; // sorted global column ids (kept) + std::vector rows; // sorted global row ids spanned by the block + bve_block_t blk; // gathered block, local ids, for the projection + bve_projection_t projection; // sized and zeroed by stage(), filled by the projection backend +}; + +// Project shape-binned candidate batches on the GPU and return a deterministic work estimate. +template +double bve_project_batch_gpu(const raft::handle_t& handle, + std::vector>& cands, + f_t tol); + +// Build symmetric current-problem implication adjacency from the original-id keyed probing cache, +// optionally unioned with forcings harvested from earlier block projections (also original-id). +template +std::vector> bve_build_impl_adj( + const probing_cache_t& cache, + const std::vector& reverse_original_ids, + i_t n_vars, + const probe_findings_t* prior_original_id_findings = nullptr); + +// Run block BVE using caller-provided implication adjacency and deadline. Returns true iff at least +// one validated reduction was installed; `work_units` receives a deterministic unscaled estimate. +// `out_findings`, when given, is appended with the implications read off every projected block +// (original-id frame) -- including blocks that were not eliminated. +template +bool block_bve_presolve(problem_t& problem, + const std::vector>& impl_adj, + timer_t& timer, + double& work_units, + probe_findings_t* out_findings = nullptr, + i_t boundary_cap = BVE_MAX_BOUNDARY, + i_t scope_cap = BVE_MAX_SCOPE, + i_t clause_growth_margin = 0); + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/presolve/probing_cache.cu b/cpp/src/mip_heuristics/presolve/probing_cache.cu index fd4790479b..0d86fe66d4 100644 --- a/cpp/src/mip_heuristics/presolve/probing_cache.cu +++ b/cpp/src/mip_heuristics/presolve/probing_cache.cu @@ -161,10 +161,15 @@ void inline insert_current_probing_to_cache(i_t var_idx, const std::vector& modified_lb, const std::vector& modified_ub, const std::vector& h_integer_indices, + const std::vector& original_ids, std::atomic& n_implied_singletons) { f_t int_tol = bound_presolve.context.settings.tolerances.integrality_tolerance; + cuopt_assert(var_idx >= 0 && var_idx < (i_t)original_ids.size(), + "probe var out of original_ids range"); + const i_t var_original = original_ids[var_idx]; + cache_entry_t cache_item; cache_item.val_interval = probe_val; for (auto impacted_var_idx : h_integer_indices) { @@ -179,18 +184,21 @@ void inline insert_current_probing_to_cache(i_t var_idx, "Lower bound must be greater than or equal to original lower bound"); cuopt_assert(modified_ub[impacted_var_idx] <= get_upper(original_var_bounds), "Upper bound must be less than or equal to original upper bound"); + cuopt_assert(impacted_var_idx >= 0 && impacted_var_idx < (i_t)original_ids.size(), + "impacted var out of original_ids range"); cached_bound_t new_bound{modified_lb[impacted_var_idx], modified_ub[impacted_var_idx]}; - cache_item.var_to_cached_bound_map.insert({impacted_var_idx, new_bound}); + // Map keys are original-frame ids (same frame as reverse_original_ids / bve_build_impl_adj). + cache_item.var_to_cached_bound_map.insert({original_ids[impacted_var_idx], new_bound}); } } { std::lock_guard lock(bound_presolve.probing_cache.probing_cache_mutex); - if (!bound_presolve.probing_cache.probing_cache.count(var_idx) > 0) { + if (!bound_presolve.probing_cache.probing_cache.count(var_original) > 0) { std::array, 2> entries_per_var; entries_per_var[0] = cache_item; - bound_presolve.probing_cache.probing_cache.insert({var_idx, entries_per_var}); + bound_presolve.probing_cache.probing_cache.insert({var_original, entries_per_var}); } else { - bound_presolve.probing_cache.probing_cache[var_idx][1] = cache_item; + bound_presolve.probing_cache.probing_cache[var_original][1] = cache_item; } } } @@ -496,6 +504,7 @@ void compute_cache_for_var(i_t var_idx, h_improved_lower_bounds, h_improved_upper_bounds, h_integer_indices, + problem.original_ids, n_of_implied_singletons); } } @@ -703,36 +712,52 @@ void apply_substitution_queue_to_problem( std::vector offset_values; std::vector coefficient_values; - // Get variable_mapping to convert current indices to original indices + // Get variable_mapping to convert current indices to post-Papilo frame auto h_variable_mapping = host_copy(problem.presolve_data.variable_mapping, problem.handle_ptr->get_stream()); problem.handle_ptr->sync_stream(); + // Collect AffineSub reconstructions, then append in deterministic order (by substituted_var). + std::vector> batch_recs; + batch_recs.reserve(all_substitutions.size()); for (const auto& [substituting_var, substitutions] : all_substitutions) { for (const auto& [substituted_var, substitution] : substitutions) { CUOPT_LOG_TRACE("Applying substitution: %d -> %d", substitution.substituting_var, substitution.substituted_var); + cuopt_assert(substitution.substituted_var >= 0 && + substitution.substituted_var < (i_t)h_variable_mapping.size(), + "substituted_var out of variable_mapping range"); + cuopt_assert(substitution.substituting_var >= 0 && + substitution.substituting_var < (i_t)h_variable_mapping.size(), + "substituting_var out of variable_mapping range"); var_indices.push_back(substitution.substituted_var); substituting_var_indices.push_back(substitution.substituting_var); offset_values.push_back(substitution.offset); coefficient_values.push_back(substitution.coefficient); - // Store substitution for post-processing (convert to original variable IDs) - substitution_t sub; - sub.timestamp = substitution.timestamp; - sub.substituted_var = h_variable_mapping[substitution.substituted_var]; - sub.substituting_var = h_variable_mapping[substitution.substituting_var]; - sub.offset = substitution.offset; - sub.coefficient = substitution.coefficient; - problem.presolve_data.variable_substitutions.push_back(sub); - CUOPT_LOG_TRACE("Stored substitution for post-processing: x[%d] = %f + %f * x[%d]", - sub.substituted_var, - sub.offset, - sub.coefficient, - sub.substituting_var); + var_postsolve_t rec; + rec.kind = reconstruction_kind_t::AffineSub; + rec.sub = substitution; + rec.sub.substituted_var = h_variable_mapping[substitution.substituted_var]; + rec.sub.substituting_var = h_variable_mapping[substitution.substituting_var]; + batch_recs.push_back(std::move(rec)); + CUOPT_LOG_TRACE("Stored AffineSub for post-processing: x[%d] = %f + %f * x[%d]", + batch_recs.back().sub.substituted_var, + batch_recs.back().sub.offset, + batch_recs.back().sub.coefficient, + batch_recs.back().sub.substituting_var); } } + std::sort(batch_recs.begin(), + batch_recs.end(), + [](const var_postsolve_t& a, const var_postsolve_t& b) { + return a.sub.substituted_var < b.sub.substituted_var; + }); + auto& recs = problem.presolve_data.var_postsolve; + recs.insert(recs.end(), + std::make_move_iterator(batch_recs.begin()), + std::make_move_iterator(batch_recs.end())); if (!var_indices.empty()) { problem.substitute_variables( @@ -850,6 +875,15 @@ bool compute_probing_cache(bound_presolve_t& bound_presolve, timer_t timer) { raft::common::nvtx::range fun_scope("compute_probing_cache"); + + // Probing runs once per solve, ahead of the block-BVE rounds that consume the cache. A second + // call would drop everything those rounds folded back in, so refuse to start on a populated one. + cuopt_assert(bound_presolve.probing_cache.probing_cache.empty(), + "probing cache is built once per solve"); + // Entries are keyed by original id, so every caller must have compacted the problem with + // remap_cache_ids set. + cuopt_assert(problem.original_ids.size() == (size_t)problem.n_variables, + "probing cache needs id maps that match the current column set"); // we dont want to compute the probing cache for all variables for time and computation resources auto priority_indices = compute_priority_indices_by_implied_integers(problem); CUOPT_LOG_DEBUG("Computing probing cache"); @@ -951,6 +985,54 @@ bool compute_probing_cache(bound_presolve_t& bound_presolve, return problem_is_infeasible.load(); } +// incorporate implications discovered by block-BVE +template +void probing_cache_t::merge_forcings(const std::vector>& forcings, + std::vector>& fixings) +{ + i_t n_added = 0; + i_t n_tightened = 0; + i_t n_contradicted = 0; + for (const auto& forcing : forcings) { + cuopt_assert(forcing.var != forcing.forced_var, "self-forcing is not a projection finding"); + auto entry_it = probing_cache.find(forcing.var); + if (entry_it == probing_cache.end()) { continue; } + const f_t probed_val = forcing.value ? f_t(1) : f_t(0); + const f_t forced_val = forcing.forced_value ? f_t(1) : f_t(0); + for (cache_entry_t& entry : entry_it->second) { + if (entry.var_to_cached_bound_map.empty()) { continue; } + if (entry.val_interval.interval_type != interval_type_t::EQUALS) { continue; } + if (entry.val_interval.val != probed_val) { continue; } + auto [bound_it, inserted] = entry.var_to_cached_bound_map.insert( + {forcing.forced_var, cached_bound_t{forced_val, forced_val}}); + if (inserted) { + ++n_added; + continue; + } + cached_bound_t& bound = bound_it->second; + const f_t lb = std::max(bound.lb, forced_val); + const f_t ub = std::min(bound.ub, forced_val); + // Both the cached bound and the projection are valid and share the antecedent var == probed + // value, so an empty intersection proves only that the antecedent cannot hold. The slot is + // dead from here on, hence no tightening; the opposite value is the sound conclusion. + if (lb > ub) { + fixings.emplace_back(forcing.var, !forcing.value); + ++n_contradicted; + continue; + } + n_tightened += (lb != bound.lb || ub != bound.ub); + bound.lb = lb; + bound.ub = ub; + } + } + CUOPT_LOG_DEBUG( + "BVE forcings %zu: added %d and tightened %d probing cache bounds, %d contradicted a probe", + forcings.size(), + n_added, + n_tightened, + n_contradicted); +} + #define INSTANTIATE(F_TYPE) \ template bool compute_probing_cache(bound_presolve_t & bound_presolve, \ problem_t & problem, \ diff --git a/cpp/src/mip_heuristics/presolve/probing_cache.cuh b/cpp/src/mip_heuristics/presolve/probing_cache.cuh index ec532febb9..345b3b1f6c 100644 --- a/cpp/src/mip_heuristics/presolve/probing_cache.cuh +++ b/cpp/src/mip_heuristics/presolve/probing_cache.cuh @@ -66,6 +66,21 @@ struct cache_entry_t { std::unordered_map> var_to_cached_bound_map; }; +// A forcing read off an exactly projected block: var == value implies forced_var == forced_value. +template +struct probe_forcing_t { + i_t var; + i_t forced_var; + bool value; + bool forced_value; +}; + +template +struct probe_findings_t { + std::vector> forcings; + std::vector> fixings; // var forced to value by its block alone +}; + template class probing_cache_t { public: @@ -87,6 +102,12 @@ class probing_cache_t { f_t first_probe, f_t second_probe, f_t integrality_tolerance); + // Intersect block-BVE-derived forcings into the entries that already cover the same variable. + // Both sides are conditioned on the same antecedent, so an empty intersection disproves that + // antecedent, not the model: the variable is appended to fixings with the opposite value. Global + // infeasibility is the case where both polarities get fixed, which apply_bve_fixings detects. + void merge_forcings(const std::vector>& forcings, + std::vector>& fixings); // add the results of probing cache to secondary CG structure if not already in a gub constraint. // use the same activity computation that we will use in BP rounding. // use GUB constraints to find fixings in bulk rounding diff --git a/cpp/src/mip_heuristics/problem/presolve_data.cu b/cpp/src/mip_heuristics/problem/presolve_data.cu index e834ce8c21..e3976d022b 100644 --- a/cpp/src/mip_heuristics/problem/presolve_data.cu +++ b/cpp/src/mip_heuristics/problem/presolve_data.cu @@ -135,19 +135,45 @@ void presolve_data_t::post_process_assignment( } } - // Apply variable substitutions from probing: x_substituted = offset + coefficient * - // x_substituting - for (const auto& sub : variable_substitutions) { - cuopt_assert(sub.substituted_var < (i_t)h_assignment.size(), "substituted_var out of bounds"); - cuopt_assert(sub.substituting_var < (i_t)h_assignment.size(), "substituting_var out of bounds"); - h_assignment[sub.substituted_var] = - sub.offset + sub.coefficient * h_assignment[sub.substituting_var]; - CUOPT_LOG_DEBUG("Post-process substitution: x[%d] = %f + %f * x[%d] = %f", - sub.substituted_var, - sub.offset, - sub.coefficient, - sub.substituting_var, - h_assignment[sub.substituted_var]); + // Reverse-append undo of the unified GPU-presolve reconstruction log (probe AffineSub and BVE + // BlockBve interleaved in commit order across outer rounds). + for (auto it = var_postsolve.rbegin(); it != var_postsolve.rend(); ++it) { + const auto& rec = *it; + switch (rec.kind) { + case reconstruction_kind_t::BlockBve: { + cuopt_assert(rec.bve.witness.size() == (size_t{1} << rec.bve.boundary.size()), + "block witness size mismatch"); + uint32_t pattern = 0; + for (size_t j = 0; j < rec.bve.boundary.size(); ++j) { + cuopt_assert(rec.bve.boundary[j] < (i_t)h_assignment.size(), + "block boundary out of bounds"); + const int bit = (h_assignment[rec.bve.boundary[j]] > static_cast(0.5)) ? 1 : 0; + pattern |= (static_cast(bit) << j); + } + const uint32_t w = rec.bve.witness[pattern]; + for (size_t k = 0; k < rec.bve.interior.size(); ++k) { + cuopt_assert(rec.bve.interior[k] < (i_t)h_assignment.size(), + "block interior out of bounds"); + h_assignment[rec.bve.interior[k]] = static_cast((w >> k) & 1u); + } + break; + } + case reconstruction_kind_t::AffineSub: { + cuopt_assert(rec.sub.substituted_var < (i_t)h_assignment.size(), + "substituted_var out of bounds"); + cuopt_assert(rec.sub.substituting_var < (i_t)h_assignment.size(), + "substituting_var out of bounds"); + h_assignment[rec.sub.substituted_var] = + rec.sub.offset + rec.sub.coefficient * h_assignment[rec.sub.substituting_var]; + CUOPT_LOG_DEBUG("Post-process substitution: x[%d] = %f + %f * x[%d] = %f", + rec.sub.substituted_var, + rec.sub.offset, + rec.sub.coefficient, + rec.sub.substituting_var, + h_assignment[rec.sub.substituted_var]); + break; + } + } } // this separate resizing is needed because of the callback diff --git a/cpp/src/mip_heuristics/problem/presolve_data.cuh b/cpp/src/mip_heuristics/problem/presolve_data.cuh index 5f0b7f53c3..8ee4c7b6b1 100644 --- a/cpp/src/mip_heuristics/problem/presolve_data.cuh +++ b/cpp/src/mip_heuristics/problem/presolve_data.cuh @@ -34,6 +34,23 @@ struct substitution_t { f_t coefficient; }; +template +struct bve_postsolve_t { + std::vector interior; + std::vector boundary; + std::vector witness; // size 2^boundary.size() +}; + +enum class reconstruction_kind_t : uint8_t { AffineSub = 0, BlockBve = 1 }; + +// could be a tagged union, but alas non-trivial members +template +struct var_postsolve_t { + reconstruction_kind_t kind{}; + substitution_t sub{}; + bve_postsolve_t bve{}; +}; + template class presolve_data_t { public: @@ -62,7 +79,7 @@ class presolve_data_t { papilo_reduced_to_original_map(other.papilo_reduced_to_original_map), papilo_original_to_reduced_map(other.papilo_original_to_reduced_map), papilo_original_num_variables(other.papilo_original_num_variables), - variable_substitutions(other.variable_substitutions) + var_postsolve(other.var_postsolve) { } @@ -76,7 +93,7 @@ class presolve_data_t { fixed_var_assignment.begin(), fixed_var_assignment.end(), 0.); - variable_substitutions.clear(); + var_postsolve.clear(); } void reset_additional_vars(const problem_t& problem, const raft::handle_t* handle_ptr) @@ -128,9 +145,9 @@ class presolve_data_t { std::vector papilo_reduced_to_original_map{}; std::vector papilo_original_to_reduced_map{}; i_t papilo_original_num_variables{0}; - // Variable substitutions from probing: x_substituted = offset + coefficient * x_substituting - // Applied in post_process_assignment to recover substituted variable values - std::vector> variable_substitutions; + // Append-only GPU-presolve reconstruction log (AffineSub from probing, BlockBve from block-BVE). + // post_process_assignment replays in reverse append order. + std::vector> var_postsolve; }; } // namespace mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/problem/problem.cu b/cpp/src/mip_heuristics/problem/problem.cu index ccba2d5f2b..063b5c288d 100644 --- a/cpp/src/mip_heuristics/problem/problem.cu +++ b/cpp/src/mip_heuristics/problem/problem.cu @@ -12,6 +12,7 @@ #include #include +#include #include #include @@ -204,6 +205,7 @@ problem_t::problem_t(const problem_t& problem_) var_names(problem_.var_names), row_names(problem_.row_names), objective_name(problem_.objective_name), + objective_offset(problem_.presolve_data.objective_offset), is_scaled_(problem_.is_scaled_), preprocess_called(problem_.preprocess_called), objective_is_integral(problem_.objective_is_integral), @@ -263,6 +265,7 @@ problem_t::problem_t(const problem_t& problem_, var_names(problem_.var_names), row_names(problem_.row_names), objective_name(problem_.objective_name), + objective_offset(problem_.presolve_data.objective_offset), is_scaled_(problem_.is_scaled_), preprocess_called(problem_.preprocess_called), objective_is_integral(problem_.objective_is_integral), @@ -365,6 +368,9 @@ problem_t::problem_t(const problem_t& problem_, bool no_deep var_names(problem_.var_names), row_names(problem_.row_names), objective_name(problem_.objective_name), + // presolve_data above picks its source from no_deep_copy and presolve moves that offset, so + // read the member just built (declared ahead of this one) rather than problem_ again. + objective_offset(presolve_data.objective_offset), is_scaled_(problem_.is_scaled_), preprocess_called(problem_.preprocess_called), objective_is_integral(problem_.objective_is_integral), @@ -1229,124 +1235,6 @@ void problem_t::insert_constraints(constraints_delta_t& h_co pdlp::combine_constraint_bounds(*this, combined_bounds); } -// Best rational approximation p/q to x with q <= max_denom, via continued fractions. -// Returns the last valid convergent if the denominator limit is reached. -std::pair rational_approximation(double x, int64_t max_denom, double epsilon) -{ - double ax = std::abs(x); - if (ax < epsilon) { return {0, 1}; } - - if (x < 0) { - auto [p, q] = rational_approximation(-x, max_denom, epsilon); - return {-p, q}; - } - - int64_t p_prev2 = 1, q_prev2 = 0; - int64_t p_prev1 = (int64_t)std::floor(x), q_prev1 = 1; - - double remainder = x - std::floor(x); - - for (int iter = 0; iter < 100; ++iter) { - if (std::abs(remainder) < 1e-15) break; - - remainder = 1.0 / remainder; - int64_t a = (int64_t)std::floor(remainder); - remainder -= a; - - int64_t p_curr = a * p_prev1 + p_prev2; - int64_t q_curr = a * q_prev1 + q_prev2; - - if (q_curr > max_denom) break; - // overflow guard - if (std::abs(p_curr) < std::abs(p_prev1)) break; - - p_prev2 = p_prev1; - q_prev2 = q_prev1; - p_prev1 = p_curr; - q_prev1 = q_curr; - - double approx_err = x - (double)p_curr / (double)q_curr; - if (std::abs(approx_err) < epsilon) break; - } - - return {p_prev1, q_prev1}; -} - -// Brute-force: try scalars 1..max_brute and return the smallest that makes all coefficients -// integral. -double find_scaling_brute_force(const std::vector& coefficients, - int max_brute = 100, - double tol = 1e-6) -{ - for (int s = 1; s <= max_brute; ++s) { - bool ok = true; - for (double c : coefficients) { - double scaled = s * c; - if (std::abs(scaled - std::round(scaled)) > tol) { - ok = false; - break; - } - } - if (ok) return (double)s; - } - return std::numeric_limits::quiet_NaN(); -} - -// Continued-fractions approach: rationalize each coefficient, compute scm/gcd incrementally. -double find_scaling_rational(const std::vector& coefficients, - double maxscale = 1e6, - int64_t maxdnom = 10000000, - double maxfinal = 10000, - double intcheck_tol = 1e-6) -{ - constexpr double no_scaling = std::numeric_limits::quiet_NaN(); - double epsilon = 1.0 / maxscale; - - int64_t gcd = 0; - int64_t scm = 1; - - for (double c : coefficients) { - auto [num, den] = rational_approximation(c, maxdnom, epsilon); - if (den == 0 || num == 0) continue; - - int64_t abs_num = std::abs(num); - if (gcd == 0) { - gcd = abs_num; - scm = den; - } else { - gcd = std::gcd(gcd, abs_num); - int64_t factor = den / std::gcd(scm, den); - int64_t new_scm; - if (__builtin_mul_overflow(scm, factor, &new_scm)) return no_scaling; - scm = new_scm; - } - - if ((double)scm / (double)gcd > maxscale) return no_scaling; - } - - if (gcd == 0) return 1.0; - - double intscalar = (double)scm / (double)gcd; - if (intscalar > maxfinal) return no_scaling; - - for (double c : coefficients) { - double scaled = intscalar * c; - if (std::abs(scaled - std::round(scaled)) > intcheck_tol) return no_scaling; - } - - return intscalar; -} - -// Finds the smallest integer scaling factor s such that s * c_i is integral for all i. -// Tries a brute-force sweep first (cheap, numerically robust), then falls back to -// continued fractions for larger scalars. -double find_objective_scaling_factor(const std::vector& coefficients) -{ - double s = find_scaling_brute_force(coefficients); - if (!std::isnan(s)) return s; - return find_scaling_rational(coefficients); -} - template void problem_t::set_implied_integers(const std::vector& implied_integer_indices) { @@ -2168,36 +2056,29 @@ void problem_t::set_constraints_from_host_user_problem( raft::common::nvtx::range fun_scope("set_constraints_from_host_user_problem"); cuopt_assert(user_problem.handle_ptr == handle_ptr, "handle mismatch"); cuopt_assert(user_problem.num_cols == n_variables, "num cols mismatch"); - n_constraints = user_problem.num_rows; - cuopt_assert(user_problem.rhs.size() == static_cast(n_constraints), "rhs size mismatch"); - cuopt_assert(user_problem.row_sense.size() == static_cast(n_constraints), + const i_t num_rows = user_problem.num_rows; + cuopt_assert(user_problem.rhs.size() == static_cast(num_rows), "rhs size mismatch"); + cuopt_assert(user_problem.row_sense.size() == static_cast(num_rows), "row sense size mismatch"); cuopt_assert(user_problem.range_rows.size() == user_problem.range_value.size(), "range rows/value size mismatch"); - csr_matrix_t csr_A(n_constraints, n_variables, user_problem.A.nnz()); + csr_matrix_t csr_A(num_rows, n_variables, user_problem.A.nnz()); user_problem.A.to_compressed_row(csr_A); - nnz = csr_A.row_start[n_constraints]; - empty = (nnz == 0 && n_constraints == 0 && n_variables == 0); - auto stream = handle_ptr->get_stream(); - cuopt::device_copy(coefficients, csr_A.x, stream); - cuopt::device_copy(variables, csr_A.j, stream); - cuopt::device_copy(offsets, csr_A.row_start, stream); - - std::vector h_constraint_lower_bounds(n_constraints); - std::vector h_constraint_upper_bounds(n_constraints); - std::vector range_value_per_row(n_constraints, f_t{0}); - std::vector is_range_row(n_constraints, 0); + std::vector h_constraint_lower_bounds(num_rows); + std::vector h_constraint_upper_bounds(num_rows); + std::vector range_value_per_row(num_rows, f_t{0}); + std::vector is_range_row(num_rows, 0); for (size_t idx = 0; idx < user_problem.range_rows.size(); ++idx) { auto row = user_problem.range_rows[idx]; - cuopt_assert(row >= 0 && row < n_constraints, "range row out of bounds"); + cuopt_assert(row >= 0 && row < num_rows, "range row out of bounds"); is_range_row[row] = 1; range_value_per_row[row] = user_problem.range_value[idx]; } const auto inf = std::numeric_limits::infinity(); - for (i_t i = 0; i < n_constraints; ++i) { + for (i_t i = 0; i < num_rows; ++i) { const f_t rhs = user_problem.rhs[i]; const char sense = user_problem.row_sense[i]; if (sense == 'E') { @@ -2214,32 +2095,66 @@ void problem_t::set_constraints_from_host_user_problem( cuopt_assert(false, "Unsupported row sense"); } } + set_constraints_from_host_csr(csr_A.row_start, + csr_A.j, + csr_A.x, + h_constraint_lower_bounds, + h_constraint_upper_bounds, + user_problem.row_names); +} - cuopt::device_copy(constraint_lower_bounds, h_constraint_lower_bounds, stream); - cuopt::device_copy(constraint_upper_bounds, h_constraint_upper_bounds, stream); - - if (!user_problem.row_names.empty()) { - row_names = user_problem.row_names; - } else if (row_names.size() != static_cast(n_constraints)) { - row_names.clear(); +template +void problem_t::set_constraints_from_host_csr(const std::vector& offsets_in, + const std::vector& variables_in, + const std::vector& coefficients_in, + const std::vector& row_lower, + const std::vector& row_upper, + const std::vector& names) +{ + raft::common::nvtx::range fun_scope("set_constraints_from_host_csr"); + n_constraints = static_cast(row_lower.size()); + cuopt_assert(row_upper.size() == static_cast(n_constraints), "row bound size mismatch"); + cuopt_assert(offsets_in.size() == static_cast(n_constraints) + 1, + "offsets size mismatch"); + cuopt_assert(!offsets_in.empty() && offsets_in.front() == 0, "invalid CSR offsets"); + cuopt_assert(std::is_sorted(offsets_in.begin(), offsets_in.end()), "unsorted CSR offsets"); + cuopt_assert(variables_in.size() == coefficients_in.size(), "csr index/value size mismatch"); + cuopt_assert(static_cast(offsets_in.back()) == variables_in.size(), + "CSR offsets/entries size mismatch"); + cuopt_assert(names.empty() || names.size() == static_cast(n_constraints), + "row names size mismatch"); + for (i_t variable : variables_in) { + cuopt_assert(variable >= 0 && variable < n_variables, "CSR variable out of bounds"); } + nnz = static_cast(variables_in.size()); + empty = (nnz == 0 && n_constraints == 0 && n_variables == 0); + auto stream = handle_ptr->get_stream(); + cuopt::device_copy(coefficients, coefficients_in, stream); + cuopt::device_copy(variables, variables_in, stream); + cuopt::device_copy(offsets, offsets_in, stream); + cuopt::device_copy(constraint_lower_bounds, row_lower, stream); + cuopt::device_copy(constraint_upper_bounds, row_upper, stream); + + // the previous row set is gone: drop stale row names and any fixed-problem cache + row_names = names; integer_fixed_problem = nullptr; + fixing_helpers.reduction_in_rhs.resize(n_constraints, stream); - auto prev_dual_size = lp_state.prev_dual.size(); + thrust::fill(handle_ptr->get_thrust_policy(), + fixing_helpers.reduction_in_rhs.begin(), + fixing_helpers.reduction_in_rhs.end(), + f_t{0}); lp_state.prev_dual.resize(n_constraints, stream); - if (n_constraints > (i_t)prev_dual_size) { - thrust::fill(handle_ptr->get_thrust_policy(), - lp_state.prev_dual.begin() + prev_dual_size, - lp_state.prev_dual.end(), - f_t{0}); - } + thrust::fill( + handle_ptr->get_thrust_policy(), lp_state.prev_dual.begin(), lp_state.prev_dual.end(), f_t{0}); handle_ptr->sync_stream(); RAFT_CHECK_CUDA(stream); compute_transpose_of_problem(); combined_bounds.resize(n_constraints, stream); pdlp::combine_constraint_bounds(*this, combined_bounds); + recompute_auxilliary_data(false); } template diff --git a/cpp/src/mip_heuristics/problem/problem.cuh b/cpp/src/mip_heuristics/problem/problem.cuh index bcc3f06fc2..c3eedd4afb 100644 --- a/cpp/src/mip_heuristics/problem/problem.cuh +++ b/cpp/src/mip_heuristics/problem/problem.cuh @@ -143,6 +143,14 @@ class problem_t { cuopt::mathematical_optimization::simplex::user_problem_t& user_problem) const; void set_constraints_from_host_user_problem( const cuopt::mathematical_optimization::simplex::user_problem_t& user_problem); + // Replace the constraint matrix + row bounds in place from host CSR + // Used by presolve passes that rewrite rows in place (e.g. block-BVE) + void set_constraints_from_host_csr(const std::vector& offsets, + const std::vector& variables, + const std::vector& coefficients, + const std::vector& row_lower, + const std::vector& row_upper, + const std::vector& names); uint32_t get_fingerprint() const; @@ -326,7 +334,7 @@ class problem_t { std::vector row_names{}; /** name of the objective (only a single objective is currently allowed) */ std::string objective_name; - f_t objective_offset; + f_t objective_offset{0}; bool is_scaled_{false}; bool preprocess_called{false}; bool objective_is_integral{false}; diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 64d78efbc0..d60e3db81f 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -196,7 +196,7 @@ mip_solution_t run_mip_solver( scaled_problem.preprocess_problem(); scaled_problem.related_vars_time_limit = settings.heuristic_params.related_vars_time_limit; const i_t n_vars_before = scaled_problem.n_variables; - mip::trivial_presolve(scaled_problem); + mip::trivial_presolve(scaled_problem, /*remap_cache_ids=*/true); #ifdef DETECT_SYMMETRY_BEFORE_PRESOLVE // Trivial presolve may remove unused variables and renumber the remaining ones. diff --git a/cpp/src/utilities/integer_scaling.hpp b/cpp/src/utilities/integer_scaling.hpp new file mode 100644 index 0000000000..4456977698 --- /dev/null +++ b/cpp/src/utilities/integer_scaling.hpp @@ -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 */ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace cuopt { + +namespace detail { + +// Best rational approximation p/q to x with q <= max_denom, via continued fractions. Returns the +// last valid convergent if the denominator limit is reached. +inline std::pair rational_approximation(double x, + int64_t max_denom, + double epsilon) +{ + cuopt_assert(std::isfinite(x), "non-finite coefficient"); + if (!std::isfinite(x)) return {0, 0}; + + double ax = std::abs(x); + if (ax < epsilon) { return {0, 1}; } + + if (x < 0) { + auto [p, q] = rational_approximation(-x, max_denom, epsilon); + return {-p, q}; + } + + const double integer_part = std::floor(x); + if (integer_part >= (double)std::numeric_limits::max()) return {0, 0}; + + int64_t p_prev2 = 1, q_prev2 = 0; + int64_t p_prev1 = (int64_t)integer_part, q_prev1 = 1; + + double remainder = x - integer_part; + + for (int iter = 0; iter < 100; ++iter) { + if (std::abs(remainder) < 1e-15) break; + + remainder = 1.0 / remainder; + const double quotient = std::floor(remainder); + if (!std::isfinite(quotient) || quotient >= (double)std::numeric_limits::max()) { + return {0, 0}; + } + int64_t a = (int64_t)quotient; + remainder -= a; + + int64_t p_product; + int64_t q_product; + int64_t p_curr; + int64_t q_curr; + if (__builtin_mul_overflow(a, p_prev1, &p_product) || + __builtin_add_overflow(p_product, p_prev2, &p_curr) || + __builtin_mul_overflow(a, q_prev1, &q_product) || + __builtin_add_overflow(q_product, q_prev2, &q_curr)) { + return {0, 0}; + } + + if (q_curr > max_denom) break; + + p_prev2 = p_prev1; + q_prev2 = q_prev1; + p_prev1 = p_curr; + q_prev1 = q_curr; + + double approx_err = x - (double)p_curr / (double)q_curr; + if (std::abs(approx_err) < epsilon) break; + } + + return {p_prev1, q_prev1}; +} + +// Brute-force: try scalars 1..max_brute and return the smallest that makes all coefficients +// integral. +inline double find_scaling_brute_force(const std::vector& coefficients, + int max_brute = 100, + double tol = 1e-6) +{ + for (int s = 1; s <= max_brute; ++s) { + bool ok = true; + for (double c : coefficients) { + cuopt_assert(std::isfinite(c), "non-finite coefficient"); + if (!std::isfinite(c)) return std::numeric_limits::quiet_NaN(); + double scaled = s * c; + if (!std::isfinite(scaled) || std::abs(scaled - std::round(scaled)) > tol) { + ok = false; + break; + } + } + if (ok) return (double)s; + } + return std::numeric_limits::quiet_NaN(); +} + +} // namespace detail + +// Continued-fractions approach: rationalize each coefficient, compute scm/gcd incrementally. +// Returns the smallest positive multiplier s such that s * c is (near-)integer for every c, or NaN +// if no such multiplier exists within the caps. +inline double find_scaling_rational(const std::vector& coefficients, + double maxscale = 1e6, + int64_t maxdnom = 10000000, + double maxfinal = 10000, + double intcheck_tol = 1e-6) +{ + constexpr double no_scaling = std::numeric_limits::quiet_NaN(); + double epsilon = 1.0 / maxscale; + + int64_t gcd = 0; + int64_t scm = 1; + + for (double c : coefficients) { + auto [num, den] = detail::rational_approximation(c, maxdnom, epsilon); + if (den == 0) return no_scaling; + if (num == 0) continue; + + if (num == std::numeric_limits::min()) return no_scaling; + int64_t abs_num = std::abs(num); + if (gcd == 0) { + gcd = abs_num; + scm = den; + } else { + gcd = std::gcd(gcd, abs_num); + int64_t factor = den / std::gcd(scm, den); + int64_t new_scm; + if (__builtin_mul_overflow(scm, factor, &new_scm)) return no_scaling; + scm = new_scm; + } + + if ((double)scm / (double)gcd > maxscale) return no_scaling; + } + + if (gcd == 0) return 1.0; + + double intscalar = (double)scm / (double)gcd; + if (intscalar > maxfinal) return no_scaling; + + for (double c : coefficients) { + double scaled = intscalar * c; + if (!std::isfinite(scaled) || std::abs(scaled - std::round(scaled)) > intcheck_tol) + return no_scaling; + } + + return intscalar; +} + +// Finds the smallest integer scaling factor s such that s * c_i is integral for all i. Tries a +// brute-force sweep first (cheap, numerically robust), then falls back to continued fractions for +// larger scalars. +inline double find_objective_scaling_factor(const std::vector& coefficients) +{ + double s = detail::find_scaling_brute_force(coefficients); + if (!std::isnan(s)) return s; + return find_scaling_rational(coefficients); +} + +// A bound counts as "infinite" if non-finite or at/above the solver's large-bound sentinel. +template +inline bool scaling_bound_finite(f_t x) +{ + return std::isfinite(x) && std::abs(x) < f_t(1e30); +} + +// An exact subset sum of at most max_len integer terms, plus the bound compare, must stay inside +// the mantissa of the type that holds the sum for it to never round: 2^24 for fp32, 2^53 for fp64. +// Callers store the scaled row back as f_t and sum it as f_t, so the budget follows f_t rather than +// the double used internally to search for the multiplier. +template +inline constexpr double exact_subset_sum_budget = + (double)(uint64_t{1} << std::numeric_limits::digits); + +template +inline double row_int_scale(const f_t* coef, int n, f_t lo, f_t up, int max_len, int64_t scale_cap) +{ + static_assert(std::is_floating_point_v, "row scaling is defined for floating point rows"); + static_assert(std::numeric_limits::digits < 64, "mantissa wider than the budget shift"); + cuopt_assert(n >= 0, "negative row length"); + cuopt_assert(n <= max_len, "row length exceeds the exactness budget length"); + cuopt_assert(scale_cap > 0, "non-positive scale cap"); + + std::vector vals; + vals.reserve(n + 2); + for (int k = 0; k < n; ++k) + vals.push_back((double)coef[k]); + if (scaling_bound_finite(lo)) vals.push_back((double)lo); + if (scaling_bound_finite(up)) vals.push_back((double)up); + + const double scale = find_scaling_rational(vals, + /*maxscale=*/1e12, + /*maxdnom=*/scale_cap, + /*maxfinal=*/(double)scale_cap, + /*intcheck_tol=*/1e-9); + if (!std::isfinite(scale) || scale <= 0.0) return 0.0; + + // guard so the subset sum (<= max_len integer terms) stays within f_t's mantissa + double maxabs = 0.0; + for (double v : vals) + maxabs = std::max(maxabs, std::abs(v * scale)); + if (maxabs * (double)max_len >= exact_subset_sum_budget) return 0.0; + + return scale; +} + +} // namespace cuopt diff --git a/cpp/tests/internal/CMakeLists.txt b/cpp/tests/internal/CMakeLists.txt index c580e0117a..f13f57721d 100644 --- a/cpp/tests/internal/CMakeLists.txt +++ b/cpp/tests/internal/CMakeLists.txt @@ -27,6 +27,7 @@ ConfigureTest(NUMOPT_INTERNAL_TEST ${CUOPT_TEST_DIR}/mip/integer_with_real_bounds.cu ${CUOPT_TEST_DIR}/mip/empty_fixed_problems_test.cu ${CUOPT_TEST_DIR}/mip/presolve_test.cu + ${CUOPT_TEST_DIR}/mip/block_bve_test.cu ${CUOPT_TEST_DIR}/mip/termination_test.cu ${CUOPT_TEST_DIR}/mip/determinism_test.cu # socp diff --git a/cpp/tests/mip/block_bve_test.cu b/cpp/tests/mip/block_bve_test.cu new file mode 100644 index 0000000000..bc392bb7a0 --- /dev/null +++ b/cpp/tests/mip/block_bve_test.cu @@ -0,0 +1,879 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "../linear_programming/utilities/pdlp_test_utilities.cuh" // gtest + make_path_absolute (mip_utils.cuh deps) +#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 +#include +#include +#include + +namespace cuopt::mathematical_optimization::mip { + +// ---- host enumeration projection (the differential oracle) ---- + +template +inline bool bve_is_finite(f_t x) +{ + // finite iff it equals itself (rules out NaN) and is strictly within +/- inf + return (x == x) && (x < static_cast(INFINITY)) && (x > static_cast(-INFINITY)); +} + +// Feasibility of one packed row under a full local assignment `val` (length na+nb), with tolerance. +template +inline bool bve_row_sat(const bve_block_t& blk, int r, const int* val, f_t tol) +{ + f_t s = 0; + for (int k = blk.row_off[r]; k < blk.row_off[r + 1]; ++k) { + s += blk.row_coef[k] * static_cast(val[blk.row_var[k]]); + } + if (bve_is_finite(blk.row_up[r]) && s > blk.row_up[r] + tol) return false; + if (bve_is_finite(blk.row_lo[r]) && s < blk.row_lo[r] - tol) return false; + return true; +} + +// Project the block onto its boundary. `feas[m]` (length 2^nb) is set to 1 iff boundary pattern m +// (nb bits) admits SOME interior assignment satisfying every block row, and `witness[m]` receives +// the packed interior assignment (na bits) of the FIRST feasible completion. Both are left 0 for +// infeasible patterns. Mirrors the double loop in bve_blocks.cpp; the GPU kernel must match this. +template +inline void bve_project(const bve_block_t& blk, f_t tol, uint8_t* feas, uint32_t* witness) +{ + const int na = blk.na, nb = blk.nb; + int val[BVE_MAX_SCOPE]; + for (uint32_t m = 0; m < (1u << nb); ++m) { + for (int j = 0; j < nb; ++j) + val[na + j] = (m >> j) & 1u; + feas[m] = 0; + witness[m] = 0u; + for (uint32_t am = 0; am < (1u << na); ++am) { + for (int j = 0; j < na; ++j) + val[j] = (am >> j) & 1u; + bool ok = true; + for (int r = 0; r < blk.n_rows && ok; ++r) + ok = bve_row_sat(blk, r, val, tol); + if (ok) { + feas[m] = 1; + witness[m] = am; + break; + } + } + } +} + +enum class bve_status_t : int { + kReduced = 0, // sanity check passed; `clauses` is a sound replacement for the block rows + kSkipCaps = 1, // block violates a bound cap (defensive; detector should pre-filter) + kSkipGrowth = 2, // |clauses| > |rows| + margin (would grow the row count) + kSkipCheckFailed = + 3 // clauses did not reproduce feas (sanity check failed) => keep block verbatim +}; + +// Full per-block core on the host: project -> prime-implicate CNF -> growth gate -> inline sanity +// check. The production commit_projected does the same, but reads feas/witness from the GPU instead +// of the host bve_project above. +template +inline bve_status_t bve_project_and_check(const bve_block_t& blk, + f_t tol, + i_t margin, + bve_clause_t* clauses, + i_t* n_clauses, + uint32_t* witness) +{ + *n_clauses = 0; + if (blk.nb <= 0 || blk.nb > BVE_MAX_BOUNDARY) return bve_status_t::kSkipCaps; + if (blk.na < 0 || blk.na + blk.nb > BVE_MAX_SCOPE) return bve_status_t::kSkipCaps; + if (blk.n_rows < 0 || blk.n_rows > BVE_MAX_ROWS) return bve_status_t::kSkipCaps; + + uint8_t feas[BVE_MAX_PATTERNS]; + bve_project(blk, tol, feas, witness); + bve_cover_scratch_t scratch; + const int nc = bve_greedy_prime_cover(feas, blk.nb, clauses, BVE_MAX_CLAUSES, scratch); + if (nc < 0) return bve_status_t::kSkipGrowth; // clause explosion past cap + if (nc > blk.n_rows + margin) return bve_status_t::kSkipGrowth; + if (!bve_sanity_check(feas, blk.nb, clauses, nc)) return bve_status_t::kSkipCheckFailed; + *n_clauses = nc; + return bve_status_t::kReduced; +} + +} // namespace cuopt::mathematical_optimization::mip + +namespace cuopt::mathematical_optimization::test { + +namespace mip = cuopt::mathematical_optimization::mip; + +// A minimal "a = b OR c, with b+c <= 1 forced" block. `a` is the only zero-objective binary aux +// (b and c carry objective, so they stay on the boundary and are never absorbed into the interior). +// Eliminating `a` by exact projection leaves exactly ONE prime-implicate clause: b + c <= 1 (the +// boundary pattern b=c=1 is infeasible because it would force a=1 and violate a+b+c<=2). +static constexpr const char* kBlockLp = R"LP( +Minimize + obj: b + c +Subject To + r0: a - b >= 0 + r1: a - c >= 0 + r2: a + b + c <= 2 +Binaries + a + b + c +End +)LP"; + +// Same gadget with every row scaled by 1/2, so the block coefficients and bounds are FRACTIONAL. +// The feasible region (hence the reduction: b + c <= 1, `a` eliminated) is identical — positive +// row scaling preserves feasibility. This forces block-BVE's per-row integerization +// (row_int_scale) to recover integer coefficients before the exact tol-0 projection; if that +// path were wrong (N1), the reduction or its reconstruction would break. +static constexpr const char* kFractionalBlockLp = R"LP( +Minimize + obj: b + c +Subject To + r0: 0.5 a - 0.5 b >= 0 + r1: 0.5 a - 0.5 c >= 0 + r2: 0.5 a + 0.5 b + 0.5 c <= 1 +Binaries + a + b + c +End +)LP"; + +// solve_mip opens an OMP team before MIP internals that use taskloops; probing_cache sizes its +// pool from omp_get_num_threads()-1 (0 outside a parallel region → silent no-op). +template +static void with_mip_omp_team(F&& f) +{ + const int num_threads = std::max(2, omp_get_max_threads()); + const int saved_max_active_levels = omp_get_max_active_levels(); + if (saved_max_active_levels < 2) { omp_set_max_active_levels(2); } +#pragma omp parallel num_threads(num_threads) + { +#pragma omp masked + { + f(); + } + } + if (saved_max_active_levels < 2) { omp_set_max_active_levels(saved_max_active_levels); } +} + +// Production implication adjacency: bounds → probing cache → trivial compact → bve_build_impl_adj. +// If `out_infeasible` is non-null, probing infeasibility is reported there (empty adj returned); +// otherwise the caller is assumed to expect a feasible instance and we ASSERT that. +static std::vector> probing_impl_adj(mip::problem_t& problem, + bool* out_infeasible = nullptr) +{ + mip_solver_settings_t settings{}; + cuopt::timer_t timer(30.0); + mip::mip_solver_t solver(problem, settings, timer); + problem.tolerances = settings.get_tolerances(); + mip::bound_presolve_t bound_presolve(solver.context); + + bool infeasible = false; + with_mip_omp_team([&]() { + auto term_crit = bound_presolve.solve(problem); + if (term_crit != mip::termination_criterion_t::NO_UPDATE) { + bound_presolve.set_updated_bounds(problem); + } + cuopt::timer_t probing_timer(30.0); + infeasible = mip::compute_probing_cache(bound_presolve, problem, probing_timer); + if (!infeasible) { + constexpr bool remap_cache_ids = true; + mip::trivial_presolve(problem, remap_cache_ids); + } + }); + if (out_infeasible != nullptr) { + *out_infeasible = infeasible; + } else { + EXPECT_FALSE(infeasible); + } + if (infeasible) { return {}; } + return mip::bve_build_impl_adj( + bound_presolve.probing_cache, problem.reverse_original_ids, problem.n_variables); +} + +// Build one block by hand for the projection-core tests. Local ids: a=0 (interior), b=1, c=2. +static mip::bve_block_t make_block() +{ + const double INF = std::numeric_limits::infinity(); + mip::bve_block_t blk{}; + blk.na = 1; + blk.nb = 2; + blk.n_rows = 3; + int nz = 0; + auto row = [&](int r, std::initializer_list> terms, double lo, double up) { + blk.row_off[r] = nz; + for (const auto& t : terms) { + blk.row_var[nz] = t.first; + blk.row_coef[nz] = t.second; + ++nz; + } + blk.row_lo[r] = lo; + blk.row_up[r] = up; + }; + row(0, {{0, 1.0}, {1, -1.0}}, 0.0, INF); // a - b >= 0 + row(1, {{0, 1.0}, {2, -1.0}}, 0.0, INF); // a - c >= 0 + row(2, {{0, 1.0}, {1, 1.0}, {2, 1.0}}, -INF, 2.0); // a + b + c <= 2 + blk.row_off[blk.n_rows] = nz; + return blk; +} + +// --- 1. projection core: the block sanity checks, yields one clause and the right witness --- +TEST(block_bve_core, reduces_block_and_sanity_checks) +{ + auto blk = make_block(); + mip::bve_clause_t clauses[mip::BVE_MAX_CLAUSES]; + uint32_t witness[mip::BVE_MAX_PATTERNS]; + int n_clauses = 0; + auto st = mip::bve_project_and_check(blk, 1e-6, /*margin=*/0, clauses, &n_clauses, witness); + + EXPECT_EQ(st, mip::bve_status_t::kReduced); + ASSERT_EQ(n_clauses, 1); + // clause forbids boundary pattern b=1,c=1 (bits 0 and 1 both set): b + c <= 1 + EXPECT_EQ(clauses[0].lit_mask, 3u); + EXPECT_EQ(clauses[0].bit_mask, 3u); + // witness: (b=0,c=0)->a=0, (b=1,c=0)->a=1, (b=0,c=1)->a=1 + EXPECT_EQ(witness[0], 0u); + EXPECT_EQ(witness[1], 1u); + EXPECT_EQ(witness[2], 1u); +} + +// --- 2. sanity check safety: the INDEPENDENT clause evaluator rejects any clause set that +// misrepresents +// feas (the certifying-algorithm result check; not a machine-checkable certificate) --- +TEST(block_bve_core, sanity_check_rejects_corrupted_clauses) +{ + // feasible-pattern array for the block above (b=c=1 is the only infeasible pattern) + const uint8_t feas[4] = {1, 1, 1, 0}; + const mip::bve_clause_t correct[1] = {{3u, 3u}}; // b + c <= 1 + EXPECT_TRUE(mip::bve_sanity_check(feas, 2, correct, 1)); + + // dropping the clause entirely: the CNF would accept b=c=1, but feas forbids it -> rejected + EXPECT_FALSE(mip::bve_sanity_check(feas, 2, correct, 0)); + // a wrong clause (forbid b=1 only) makes a genuinely feasible pattern look infeasible -> rejected + const mip::bve_clause_t wrong[1] = {{1u, 1u}}; + EXPECT_FALSE(mip::bve_sanity_check(feas, 2, wrong, 1)); +} + +// --- N1 (numerical): the row integerization GATE. block-BVE scales each block row to integers via +// find_scaling_rational (strict caps mirroring row_int_scale) so the projection is exact at +// tolerance 0; a row that will not integerize within the caps must be REJECTED (NaN), never rounded +// into a different model. This pins the accept/reject decision that keeps large / non-rational +// coefficients off the exact-projection path. --- +TEST(block_bve_core, integer_scaling_accepts_rational_rejects_pathological) +{ + // Strict caps matching row_int_scale (maxdnom/maxfinal = BVE_INT_SCALE_MAX = 1e6). + const double kMaxScale = 1e12; + const int64_t kMaxDenom = 1000000; + const double kMaxFinal = 1e6; + const double kIntTol = 1e-9; + auto all_integer = [](double s, const std::vector& v) { + for (double c : v) + if (std::abs(s * c - std::round(s * c)) >= 1e-9) return false; + return true; + }; + + // Fractional-but-rational: {1/2, 1/4, -3/4, 1} integerize (expected multiplier 4). + { + std::vector v{0.5, 0.25, -0.75, 1.0}; + double s = cuopt::find_scaling_rational(v, kMaxScale, kMaxDenom, kMaxFinal, kIntTol); + ASSERT_TRUE(std::isfinite(s)) << "rational coefficients must integerize"; + EXPECT_GT(s, 0.0); + EXPECT_TRUE(all_integer(s, v)); + } + + // Large integer coefficients stay exact (already integer -> multiplier 1, no rounding). + { + std::vector v{1e9, -1e9, 3.0}; + double s = cuopt::find_scaling_rational(v, kMaxScale, kMaxDenom, kMaxFinal, kIntTol); + ASSERT_TRUE(std::isfinite(s)); + EXPECT_TRUE(all_integer(s, v)); + } + + // Pathological: distinct prime reciprocals need lcm(11,13,17,19,23) = 1062347 > maxfinal (1e6), + // so no bounded integer multiplier exists -> rejected (NaN), NOT silently rounded. + { + std::vector v{1.0 / 11, 1.0 / 13, 1.0 / 17, 1.0 / 19, 1.0 / 23}; + double s = cuopt::find_scaling_rational(v, kMaxScale, kMaxDenom, kMaxFinal, kIntTol); + EXPECT_TRUE(std::isnan(s)) << "un-integerizable coefficients must be rejected, got " << s; + } +} + +// A cached probe and a block projection are both valid, so they can only disagree when the +// antecedent they share is unsatisfiable. That fixes the variable to the opposite value; the model +// is infeasible only once both polarities are contradicted, which apply_bve_fixings derives from +// two fixings that disagree. Regression: the empty intersection used to be reported as global +// infeasibility outright, turning a feasible model into an INFEASIBLE answer. +TEST(block_bve_core, cache_contradiction_fixes_the_variable_instead_of_failing) +{ + constexpr int var = 7; + constexpr int forced = 9; + + // Probing has x7 = 0 => x9 = 0; the exact projection has x7 = 0 => x9 = 1. Slot 1 is left + // unpopulated, which also exercises the empty-bound-map guard. + { + mip::probing_cache_t cache; + std::array, 2> entries{}; + entries[0].val_interval = {0.0, mip::interval_type_t::EQUALS}; + entries[0].var_to_cached_bound_map[forced] = {0.0, 0.0}; + cache.probing_cache.insert({var, entries}); + + std::vector> fixings; + cache.merge_forcings({{var, forced, false, true}}, fixings); + + ASSERT_EQ(fixings.size(), 1u) << "a contradicted probe yields one fixing, not infeasibility"; + EXPECT_EQ(fixings[0].first, var); + EXPECT_TRUE(fixings[0].second) << "x7 = 0 is disproved, so x7 = 1"; + } + + // Both polarities contradicted: the two disagreeing fixings are what proves infeasibility. + { + mip::probing_cache_t cache; + std::array, 2> entries{}; + entries[0].val_interval = {0.0, mip::interval_type_t::EQUALS}; + entries[0].var_to_cached_bound_map[forced] = {0.0, 0.0}; + entries[1].val_interval = {1.0, mip::interval_type_t::EQUALS}; + entries[1].var_to_cached_bound_map[forced] = {0.0, 0.0}; + cache.probing_cache.insert({var, entries}); + + std::vector> fixings; + cache.merge_forcings({{var, forced, false, true}, {var, forced, true, true}}, fixings); + + ASSERT_EQ(fixings.size(), 2u); + std::sort(fixings.begin(), fixings.end()); + EXPECT_EQ(fixings[0], std::make_pair(var, false)); + EXPECT_EQ(fixings[1], std::make_pair(var, true)); + } + + // A forcing consistent with the cached interval tightens it and fixes nothing. + { + mip::probing_cache_t cache; + std::array, 2> entries{}; + entries[0].val_interval = {0.0, mip::interval_type_t::EQUALS}; + entries[0].var_to_cached_bound_map[forced] = {0.0, 1.0}; + cache.probing_cache.insert({var, entries}); + + std::vector> fixings; + cache.merge_forcings({{var, forced, false, true}}, fixings); + + EXPECT_TRUE(fixings.empty()) << "a consistent forcing must not fix anything"; + const auto& bound = cache.probing_cache.at(var)[0].var_to_cached_bound_map.at(forced); + EXPECT_EQ(bound.lb, 1.0); + EXPECT_EQ(bound.ub, 1.0); + } +} + +// Build a random block LAYOUT (na/nb/n_rows + sparsity pattern), coefficients/bounds left unset. +// Reps of one shape reuse the SAME layout so they land in one GPU shape-bin (exercising the num>1 +// path). +static mip::bve_block_t make_block_layout(std::mt19937& rng, int na, int nb, int n_rows) +{ + const int scope = na + nb; + mip::bve_block_t blk{}; + blk.na = na; + blk.nb = nb; + blk.n_rows = n_rows; + std::uniform_int_distribution present(0, 1); // is a var in this row + int nz = 0; + for (int r = 0; r < n_rows; ++r) { + blk.row_off[r] = nz; + for (int v = 0; v < scope; ++v) + if (present(rng)) blk.row_var[nz++] = v; + if (nz == blk.row_off[r]) blk.row_var[nz++] = r % scope; // never leave an empty row + } + blk.row_off[n_rows] = nz; + return blk; +} + +// Fill a layout's coefficients (small integers) and bounds (randomly ±inf), leaving the pattern +// fixed. +static void randomize_block_data(std::mt19937& rng, mip::bve_block_t& blk) +{ + const double INF = std::numeric_limits::infinity(); + const double coefs[4] = {-2.0, -1.0, 1.0, 2.0}; + std::uniform_int_distribution coef_pick(0, 3); + std::uniform_int_distribution bnd_pick(0, 2); // 0:[lo,inf] 1:[-inf,up] 2:[lo,up] + for (int k = 0; k < blk.row_off[blk.n_rows]; ++k) + blk.row_coef[k] = coefs[coef_pick(rng)]; + for (int r = 0; r < blk.n_rows; ++r) { + const int terms = blk.row_off[r + 1] - blk.row_off[r]; + // Activity under 0/1 vars and coefs in {-2,-1,1,2} lies in [-2*terms, 2*terms]. Pick finite + // uppers in [0, 2*terms] so they can bind (not always equal to the loose max activity). + const double lo = -static_cast(terms); + std::uniform_int_distribution up_pick(0, 2 * terms); + const double up = static_cast(up_pick(rng)); + const int kind = bnd_pick(rng); + blk.row_lo[r] = (kind == 1) ? -INF : lo; + blk.row_up[r] = (kind == 0) ? INF : up; + } +} + +// --- projection correctness: the GPU batch projection must equal the host enumeration oracle on a +// diverse batch (varied na/nb/rows, ±inf bounds, multiple distinct shapes, and >1-block bins). +// This is what pins projection correctness; the inline sanity check cannot (it trusts feas). +// Runs the same function two independent ways and asserts feas + witness agree everywhere. +TEST(block_bve_projection, gpu_batch_matches_host_oracle) +{ + const raft::handle_t handle_{}; + std::mt19937 rng(12345u); + + // several shapes, several blocks each; reps share a layout -> one shape-bin with num>1 + const int shapes[][3] = {{1, 2, 3}, {2, 2, 2}, {1, 3, 4}, {3, 3, 5}, {2, 4, 3}, {4, 2, 4}}; + std::vector> blocks; + for (const auto& s : shapes) { + const mip::bve_block_t layout = make_block_layout(rng, s[0], s[1], s[2]); + for (int rep = 0; rep < 6; ++rep) { + mip::bve_block_t blk = layout; + randomize_block_data(rng, blk); + blocks.push_back(blk); + } + } + + std::vector> cands(blocks.size()); + for (size_t i = 0; i < blocks.size(); ++i) + cands[i].blk = + blocks[i]; // the service reads only .blk; interior/boundary/rows are unused here + + mip::bve_project_batch_gpu(handle_, cands, 1e-6); + + for (size_t i = 0; i < blocks.size(); ++i) { + uint8_t exp_feas[mip::BVE_MAX_PATTERNS]; + uint32_t exp_wit[mip::BVE_MAX_PATTERNS]; + mip::bve_project(blocks[i], 1e-6, exp_feas, exp_wit); + const int patterns = 1 << blocks[i].nb; + for (int m = 0; m < patterns; ++m) { + EXPECT_EQ(cands[i].projection.feasible[m], exp_feas[m]) << "block " << i << " pattern " << m; + if (exp_feas[m]) // witness only defined for feasible patterns + EXPECT_EQ(cands[i].projection.witness[m], exp_wit[m]) << "block " << i << " pattern " << m; + } + } +} + +// Fill a layout with LARGE integer coefficients and integer bounds (all exact fp64 integers, well +// under 2^53), leaving the pattern fixed. This is the shape block-BVE feeds the projection after +// integerization, and the magnitude range where a 1e-6-tolerance fp test would be marginal but +// exact integer arithmetic is not. +static void randomize_block_data_integer(std::mt19937& rng, mip::bve_block_t& blk) +{ + const double INF = std::numeric_limits::infinity(); + const double coefs[] = {-2e6, -1e6, 1e6, 2e6, 5e6}; + std::uniform_int_distribution coef_pick(0, 4); + std::uniform_int_distribution bnd_pick(0, 2); // 0:[lo,inf] 1:[-inf,up] 2:[lo,up] + for (int k = 0; k < blk.row_off[blk.n_rows]; ++k) + blk.row_coef[k] = coefs[coef_pick(rng)]; + for (int r = 0; r < blk.n_rows; ++r) { + const int terms = blk.row_off[r + 1] - blk.row_off[r]; + // Activity lies in [-5e6*terms, 5e6*terms]; pick finite integer bounds (multiples of 1e6) that + // can bind. + const double lo = -5e6 * static_cast(terms); + std::uniform_int_distribution up_pick(0, 2 * terms); + const double up = 1e6 * static_cast(up_pick(rng)); + const int kind = bnd_pick(rng); + blk.row_lo[r] = (kind == 1) ? -INF : lo; + blk.row_up[r] = (kind == 0) ? INF : up; + } +} + +// --- N1: the EXACT projection path. Production integerizes each block and projects at tolerance 0; +// the 1e-6 differential test above never exercises that. On large-integer-coefficient blocks +// the GPU projection at tol 0 must still equal the host enumeration oracle at tol 0 everywhere. +// --- +TEST(block_bve_projection, exact_projection_matches_host_at_tol0) +{ + const raft::handle_t handle_{}; + std::mt19937 rng(2024u); + + const int shapes[][3] = {{1, 2, 3}, {2, 2, 2}, {1, 3, 4}, {3, 3, 5}, {2, 4, 3}}; + std::vector> blocks; + for (const auto& s : shapes) { + const mip::bve_block_t layout = make_block_layout(rng, s[0], s[1], s[2]); + for (int rep = 0; rep < 6; ++rep) { + mip::bve_block_t blk = layout; + randomize_block_data_integer(rng, blk); + blocks.push_back(blk); + } + } + + std::vector> cands(blocks.size()); + for (size_t i = 0; i < blocks.size(); ++i) + cands[i].blk = blocks[i]; + + mip::bve_project_batch_gpu(handle_, cands, 0.0); // exact: tol 0 + + for (size_t i = 0; i < blocks.size(); ++i) { + uint8_t exp_feas[mip::BVE_MAX_PATTERNS]; + uint32_t exp_wit[mip::BVE_MAX_PATTERNS]; + mip::bve_project(blocks[i], 0.0, exp_feas, exp_wit); + const int patterns = 1 << blocks[i].nb; + for (int m = 0; m < patterns; ++m) { + EXPECT_EQ(cands[i].projection.feasible[m], exp_feas[m]) << "block " << i << " pattern " << m; + if (exp_feas[m]) + EXPECT_EQ(cands[i].projection.witness[m], exp_wit[m]) << "block " << i << " pattern " << m; + } + } +} + +// --- 3. end-to-end: run the pass on a problem_t, then reconstruct through postsolve --- +TEST(block_bve_presolve, end_to_end_reduction_and_reconstruction) +{ + const raft::handle_t handle_{}; + auto model = io::read_lp_from_string(kBlockLp); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, model); + mip::problem_t problem(op_problem); + problem.preprocess_problem(); + problem.presolve_data.initialize_var_mapping(problem, problem.handle_ptr); + const int n_before = problem.n_variables; + + auto impl_adj = probing_impl_adj(problem); + + cuopt::timer_t bve_timer(10.0); + double bve_work_units = 0.0; + const bool applied = mip::block_bve_presolve(problem, impl_adj, bve_timer, bve_work_units); + EXPECT_TRUE(applied); + EXPECT_EQ(problem.n_variables, n_before - 1); // exactly `a` eliminated + + // Set a reduced solution with the first surviving (boundary) variable = 1; whichever of b/c it + // is, the block forces a = 1, so a correct reconstruction must satisfy the ORIGINAL constraints. + std::vector reduced(problem.n_variables, 0.0); + if (!reduced.empty()) reduced[0] = 1.0; + rmm::device_uvector assignment(problem.n_variables, handle_.get_stream()); + raft::copy(assignment.data(), reduced.data(), reduced.size(), handle_.get_stream()); + problem.presolve_data.post_process_assignment(problem, assignment, /*resize_to_original=*/true); + auto full = cuopt::host_copy(assignment, handle_.get_stream()); + handle_.sync_stream(); + + ASSERT_EQ(full.size(), static_cast(n_before)); // expanded back to all three variables + // The reconstructed full assignment must satisfy EVERY original constraint. This is order- + // independent (no assumption about which index is a/b/c): if the eliminated aux is reconstructed + // wrongly, a - b >= 0 or a - c >= 0 is violated. Since one boundary variable is set to 1, a + // correct reconstruction forces the aux to 1 — the feasibility check below is exactly that + // correctness test. + auto m_off = model.get_constraint_matrix_offsets(); + auto m_var = model.get_constraint_matrix_indices(); + auto m_val = model.get_constraint_matrix_values(); + auto m_rl = model.get_constraint_lower_bounds(); + auto m_ru = model.get_constraint_upper_bounds(); + for (size_t r = 0; r + 1 < m_off.size(); ++r) { + double s = 0.0; + for (int k = m_off[r]; k < m_off[r + 1]; ++k) + s += m_val[k] * full[m_var[k]]; + EXPECT_GE(s, m_rl[r] - 1e-6); + EXPECT_LE(s, m_ru[r] + 1e-6); + } +} + +// --- N1 end-to-end: the SAME gadget with fractional block coefficients (rows scaled by 1/2). The +// reduction and its reconstruction must be identical to the integer gadget — block-BVE has to +// integerize the 0.5 coefficients before the exact projection and undo it correctly at +// postsolve. +TEST(block_bve_presolve, fractional_gadget_reduces_and_reconstructs) +{ + const raft::handle_t handle_{}; + auto model = io::read_lp_from_string(kFractionalBlockLp); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, model); + mip::problem_t problem(op_problem); + problem.preprocess_problem(); + problem.presolve_data.initialize_var_mapping(problem, problem.handle_ptr); + const int n_before = problem.n_variables; + + auto impl_adj = probing_impl_adj(problem); + + cuopt::timer_t bve_timer(10.0); + double bve_work_units = 0.0; + const bool applied = mip::block_bve_presolve(problem, impl_adj, bve_timer, bve_work_units); + // Probing/trivial may already have eliminated the aux; either way exactly one variable is gone. + EXPECT_TRUE(applied || problem.n_variables < n_before); + ASSERT_EQ(problem.n_variables, n_before - 1) << "fractional gadget did not eliminate the aux"; + + // Set the first surviving (boundary) variable to 1; a correct reconstruction forces the aux so + // the full assignment satisfies every ORIGINAL (fractional) constraint. + std::vector reduced(problem.n_variables, 0.0); + if (!reduced.empty()) reduced[0] = 1.0; + rmm::device_uvector assignment(problem.n_variables, handle_.get_stream()); + raft::copy(assignment.data(), reduced.data(), reduced.size(), handle_.get_stream()); + problem.presolve_data.post_process_assignment(problem, assignment, /*resize_to_original=*/true); + auto full = cuopt::host_copy(assignment, handle_.get_stream()); + handle_.sync_stream(); + + ASSERT_EQ(full.size(), static_cast(n_before)); + auto m_off = model.get_constraint_matrix_offsets(); + auto m_var = model.get_constraint_matrix_indices(); + auto m_val = model.get_constraint_matrix_values(); + auto m_rl = model.get_constraint_lower_bounds(); + auto m_ru = model.get_constraint_upper_bounds(); + for (size_t r = 0; r + 1 < m_off.size(); ++r) { + double s = 0.0; + for (int k = m_off[r]; k < m_off[r + 1]; ++k) + s += m_val[k] * full[m_var[k]]; + EXPECT_GE(s, m_rl[r] - 1e-6); + EXPECT_LE(s, m_ru[r] + 1e-6); + } +} + +// Brute-force the (small, binary) reduced problem_t: enumerate all 2^n assignments, return whether +// any is feasible, the min solver-space objective, and its argmin. +struct bve_bf_t { + bool found; + double solver_obj; + std::vector x; +}; +static bve_bf_t brute_force_binary(mip::problem_t& problem) +{ + auto stream = problem.handle_ptr->get_stream(); + auto h_off = cuopt::host_copy(problem.offsets, stream); + auto h_var = cuopt::host_copy(problem.variables, stream); + auto h_coef = cuopt::host_copy(problem.coefficients, stream); + auto h_clb = cuopt::host_copy(problem.constraint_lower_bounds, stream); + auto h_cub = cuopt::host_copy(problem.constraint_upper_bounds, stream); + auto h_obj = cuopt::host_copy(problem.objective_coefficients, stream); + auto h_vb = cuopt::host_copy(problem.variable_bounds, stream); + problem.handle_ptr->sync_stream(); + + const int nv = problem.n_variables; + const int nr = problem.n_constraints; + for (int v = 0; v < nv; ++v) { // corpus is pure 0-1 + EXPECT_NEAR(get_lower(h_vb[v]), 0.0, 1e-9); + EXPECT_NEAR(get_upper(h_vb[v]), 1.0, 1e-9); + } + + bve_bf_t r{false, 0.0, {}}; + const double eps = 1e-6; + const uint64_t total = (nv >= 63) ? 0 : (uint64_t{1} << nv); + std::vector x(nv); + for (uint64_t mask = 0; mask < total; ++mask) { + for (int v = 0; v < nv; ++v) + x[v] = static_cast((mask >> v) & 1u); + bool ok = true; + for (int rr = 0; rr < nr && ok; ++rr) { + double s = 0.0; + for (int k = h_off[rr]; k < h_off[rr + 1]; ++k) + s += h_coef[k] * x[h_var[k]]; + if (s < h_clb[rr] - eps || s > h_cub[rr] + eps) ok = false; + } + if (!ok) continue; + double obj = 0.0; + for (int v = 0; v < nv; ++v) + obj += h_obj[v] * x[v]; + if (!r.found || obj < r.solver_obj - eps) { + r.found = true; + r.solver_obj = obj; + r.x = x; + } + } + return r; +} + +// Corpus of small 0-1 instances whose optima were cross-checked OFFLINE by brute force AND HiGHS. +// MPS live in datasets/mip/block_bve/ (generated by cpufj_sc22/bve_gen_fixtures.py); optima inlined +// here. Mix: gadget-rich (block-BVE fires), no-op/soundness (aux-with-objective, random feasible +// ILPs), and infeasible. +struct bve_case_t { + const char* file; + bool feasible; + double optimum; + bool expect_reduce; // gadget should shrink via probing and/or block-BVE +}; +static const bve_case_t kBveCases[] = { + {"mip/block_bve/or_used.mps", true, 1.0, true}, + {"mip/block_bve/and_used.mps", true, -2.0, true}, + {"mip/block_bve/neq_used.mps", true, -3.0, true}, + {"mip/block_bve/chain_or.mps", true, 1.0, true}, + {"mip/block_bve/two_gadgets.mps", true, 2.0, true}, + {"mip/block_bve/heavy_reduce.mps", true, 2.0, true}, + {"mip/block_bve/aux_with_obj.mps", true, 4.0, false}, + {"mip/block_bve/mixed.mps", true, -1.0, false}, + {"mip/block_bve/infeasible.mps", false, 0.0, false}, + {"mip/block_bve/random_a.mps", true, -3.0, false}, + {"mip/block_bve/random_b.mps", true, -5.0, false}, + {"mip/block_bve/random_c.mps", true, -1.0, false}, +}; + +// End-to-end equivalence: for each corpus instance, run the pass, brute-force the reduced model, +// and assert block-BVE preserved the answer. block-BVE is a PRIMAL, optimum-preserving reduction, +// so the bar is: reduced optimum == known optimum, the reduced optimum reconstructs to an +// ORIGINAL-feasible point with that objective, and infeasibility is preserved. This stresses the +// full detect -> project +// -> commit -> install -> reconstruct chain (incl. variable_mapping + witness replay), which the +// component tests above don't. +TEST(block_bve_equivalence, preserves_optimum_and_reconstruction_on_corpus) +{ + const raft::handle_t handle_{}; + bool any_reduced = false; + for (const auto& c : kBveCases) { + SCOPED_TRACE(c.file); + auto model = io::read_mps(make_path_absolute(c.file), /*fixed_format=*/false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, model); + mip::problem_t problem(op_problem); + problem.preprocess_problem(); + problem.presolve_data.initialize_var_mapping(problem, problem.handle_ptr); + const int n_before = problem.n_variables; + + bool probing_infeas = false; + auto impl_adj = probing_impl_adj(problem, &probing_infeas); + if (probing_infeas) { + EXPECT_FALSE(c.feasible) << "probing proved infeasible on a feasible instance"; + continue; + } + + cuopt::timer_t bve_timer(10.0); + double bve_work_units = 0.0; + const bool applied = mip::block_bve_presolve(problem, impl_adj, bve_timer, bve_work_units); + // Probing/trivial may already have eliminated the aux; BVE then correctly no-ops. + if (applied || problem.n_variables < n_before) { any_reduced = true; } + if (applied) { + EXPECT_LT(problem.n_variables, n_before) << "applied but variable count unchanged"; + } + if (c.expect_reduce) { + EXPECT_LT(problem.n_variables, n_before) + << "gadget fixture expected a reduction via probing and/or block-BVE"; + } + + ASSERT_LE(problem.n_variables, 24) << "brute force enumerates 2^n; keep the corpus small"; + auto bf = brute_force_binary(problem); + if (!c.feasible) { + // NOTE: if preprocess detects the infeasibility upstream and collapses the model, this may + // need to become a problem-status check instead of a no-feasible-point check. + EXPECT_FALSE(bf.found) << "reduced model is feasible but the instance is infeasible"; + continue; + } + ASSERT_TRUE(bf.found) << "reduced model is infeasible but the instance is feasible"; + + // The reduced optimum must reconstruct to an ORIGINAL-feasible point whose ORIGINAL objective + // equals the known optimum. This is offset/scaling-independent (evaluated directly on the + // original model) and catches both directions: a cut optimum -> recon_obj > optimum; a spurious + // better solution -> either the reconstruction is original-infeasible or recon_obj < optimum. + rmm::device_uvector assignment(problem.n_variables, handle_.get_stream()); + raft::copy(assignment.data(), bf.x.data(), bf.x.size(), handle_.get_stream()); + problem.presolve_data.post_process_assignment(problem, assignment, /*resize_to_original=*/true); + auto full = cuopt::host_copy(assignment, handle_.get_stream()); + handle_.sync_stream(); + + auto m_off = model.get_constraint_matrix_offsets(); + auto m_var = model.get_constraint_matrix_indices(); + auto m_val = model.get_constraint_matrix_values(); + auto m_rl = model.get_constraint_lower_bounds(); + auto m_ru = model.get_constraint_upper_bounds(); + for (size_t r = 0; r + 1 < m_off.size(); ++r) { + double s = 0.0; + for (int k = m_off[r]; k < m_off[r + 1]; ++k) + s += m_val[k] * full[m_var[k]]; + EXPECT_GE(s, m_rl[r] - 1e-6); + EXPECT_LE(s, m_ru[r] + 1e-6); + } + auto m_obj = model.get_objective_coefficients(); + ASSERT_EQ(full.size(), m_obj.size()) << "reconstruction is not in the original column frame"; + double recon_obj = 0.0; + for (size_t j = 0; j < m_obj.size(); ++j) + recon_obj += m_obj[j] * full[j]; + EXPECT_NEAR(recon_obj, c.optimum, 1e-6); + } + EXPECT_TRUE(any_reduced) << "corpus exercised no probing/block-BVE reduction path"; +} + +// Drive production MIP presolve (Papilo → cuOpt run_presolve) and optionally assert +// upper bounds on the reduced size. Pass std::numeric_limits::max() for a +// dimension to skip that check. +static void run_presolve_size_check(const char* relative_mps_path, + int max_vars = std::numeric_limits::max(), + int max_rows = std::numeric_limits::max()) +{ + const raft::handle_t handle_{}; + auto model = io::read_mps(make_path_absolute(relative_mps_path), + /*fixed_format=*/false); + auto op_problem = mps_data_model_to_optimization_problem(&handle_, model); + sort_csr(op_problem); + + mip_solver_settings_t settings{}; + settings.presolver = presolver_t::Papilo; + settings.probing = true; + settings.block_bve = true; + + auto papilo = std::make_unique>(); + auto result = papilo->apply_presolve_from_op_problem(op_problem, + problem_category_t::MIP, + settings.presolver, + /*dual_postsolve=*/false, + settings.tolerances.absolute_tolerance, + settings.tolerances.relative_tolerance, + /*time_limit=*/60.0, + /*num_cpu_threads=*/0); + ASSERT_NE(result.status, mip::third_party_presolve_status_t::INFEASIBLE) + << relative_mps_path << " infeasible after Papilo"; + ASSERT_NE(result.status, mip::third_party_presolve_status_t::UNBNDORINFEAS) + << relative_mps_path << " unbounded-or-infeasible after Papilo"; + ASSERT_NE(result.status, mip::third_party_presolve_status_t::UNBOUNDED) + << relative_mps_path << " unbounded after Papilo"; + + mip::problem_t problem(result.reduced_problem); + problem.set_papilo_presolve_data(papilo.get(), + result.reduced_to_original_map, + result.original_to_reduced_map, + op_problem.get_n_variables()); + problem.set_implied_integers(result.implied_integer_indices); + problem.preprocess_problem(); + mip::trivial_presolve(problem, /*remap_cache_ids=*/true); // mirrors solve.cu's setup + + cuopt::timer_t timer(120.0); + mip::mip_solver_t solver(problem, settings, timer); + problem.tolerances = settings.get_tolerances(); + mip::diversity_manager_t dm(solver.context); + + bool presolve_ok = false; + with_mip_omp_team([&]() { presolve_ok = dm.run_presolve(/*time_limit=*/60.0, timer); }); + + ASSERT_TRUE(presolve_ok) << relative_mps_path << " cuOpt run_presolve failed"; + if (max_vars != std::numeric_limits::max()) { + EXPECT_LT(problem.n_variables, max_vars) + << relative_mps_path << " reduced n_variables=" << problem.n_variables; + } + if (max_rows != std::numeric_limits::max()) { + EXPECT_LT(problem.n_constraints, max_rows) + << relative_mps_path << " reduced n_constraints=" << problem.n_constraints; + } +} + +TEST(block_bve_presolve, bnatt400_reduces_below_500_vars) +{ + run_presolve_size_check("mip/bnatt400.mps", /*max_vars=*/500); +} + +TEST(block_bve_presolve, bnatt500_reduces_below_500_vars) +{ + run_presolve_size_check("mip/bnatt500.mps", /*max_vars=*/500); +} + +} // namespace cuopt::mathematical_optimization::test diff --git a/datasets/mip/block_bve/and_used.mps b/datasets/mip/block_bve/and_used.mps new file mode 100644 index 0000000000..3708601858 --- /dev/null +++ b/datasets/mip/block_bve/and_used.mps @@ -0,0 +1,30 @@ +NAME +ROWS + N Obj + L r0 + L r1 + G r2 + G r3 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj -1 + c0 r0 -1 + c0 r2 -1 + c1 Obj -1 + c1 r1 -1 + c1 r2 -1 + c2 Obj 1 + c2 r3 -1 + c3 r0 1 + c3 r1 1 + c3 r2 1 + c3 r3 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r2 -1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 +ENDATA diff --git a/datasets/mip/block_bve/aux_with_obj.mps b/datasets/mip/block_bve/aux_with_obj.mps new file mode 100644 index 0000000000..ad097e0b9c --- /dev/null +++ b/datasets/mip/block_bve/aux_with_obj.mps @@ -0,0 +1,28 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + G r3 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c1 Obj 1 + c1 r1 -1 + c1 r2 -1 + c2 Obj 3 + c2 r0 1 + c2 r1 1 + c2 r2 1 + c2 r3 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r3 1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 +ENDATA diff --git a/datasets/mip/block_bve/chain_or.mps b/datasets/mip/block_bve/chain_or.mps new file mode 100644 index 0000000000..c1e1e772f1 --- /dev/null +++ b/datasets/mip/block_bve/chain_or.mps @@ -0,0 +1,40 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + G r3 + G r4 + L r5 + G r6 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c1 Obj 1 + c1 r1 -1 + c1 r2 -1 + c2 Obj 1 + c2 r4 -1 + c2 r5 -1 + c3 r0 1 + c3 r1 1 + c3 r2 1 + c3 r3 -1 + c3 r5 -1 + c4 r3 1 + c4 r4 1 + c4 r5 1 + c4 r6 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r6 1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 +ENDATA diff --git a/datasets/mip/block_bve/heavy_reduce.mps b/datasets/mip/block_bve/heavy_reduce.mps new file mode 100644 index 0000000000..4fd5d4ed45 --- /dev/null +++ b/datasets/mip/block_bve/heavy_reduce.mps @@ -0,0 +1,54 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + G r3 + G r4 + L r5 + L r6 + L r7 + G r8 + G r9 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c1 Obj 1 + c1 r1 -1 + c1 r2 -1 + c2 Obj 1 + c2 r3 -1 + c2 r5 -1 + c3 Obj 1 + c3 r4 -1 + c3 r5 -1 + c4 r0 1 + c4 r1 1 + c4 r2 1 + c4 r6 -1 + c4 r8 -1 + c5 r3 1 + c5 r4 1 + c5 r5 1 + c5 r7 -1 + c5 r8 -1 + c6 r6 1 + c6 r7 1 + c6 r8 1 + c6 r9 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r8 -1 + RHS_V r9 1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 + BV BOUND c5 + BV BOUND c6 +ENDATA diff --git a/datasets/mip/block_bve/infeasible.mps b/datasets/mip/block_bve/infeasible.mps new file mode 100644 index 0000000000..2ecd497664 --- /dev/null +++ b/datasets/mip/block_bve/infeasible.mps @@ -0,0 +1,31 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + G r3 + L r4 + L r5 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c0 r4 1 + c1 Obj 1 + c1 r1 -1 + c1 r2 -1 + c1 r5 1 + c2 r0 1 + c2 r1 1 + c2 r2 1 + c2 r3 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r3 1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 +ENDATA diff --git a/datasets/mip/block_bve/mixed.mps b/datasets/mip/block_bve/mixed.mps new file mode 100644 index 0000000000..424aeb95a5 --- /dev/null +++ b/datasets/mip/block_bve/mixed.mps @@ -0,0 +1,55 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + L r3 + L r4 + G r5 + L r6 + L r7 + L r8 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c0 r7 -1 + c0 r8 1 + c1 Obj -1 + c1 r1 -1 + c1 r2 -1 + c1 r3 -1 + c1 r5 -1 + c1 r8 1 + c2 Obj 2 + c2 r4 -1 + c2 r5 -1 + c2 r8 1 + c3 Obj -1 + c3 r6 1 + c3 r8 1 + c4 r0 1 + c4 r1 1 + c4 r2 1 + c4 r6 1 + c5 r3 1 + c5 r4 1 + c5 r5 1 + c5 r7 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r5 -1 + RHS_V r6 1 + RHS_V r8 3 +RANGES + RANGE r8 2 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 + BV BOUND c5 +ENDATA diff --git a/datasets/mip/block_bve/neq_used.mps b/datasets/mip/block_bve/neq_used.mps new file mode 100644 index 0000000000..3eff52e2c7 --- /dev/null +++ b/datasets/mip/block_bve/neq_used.mps @@ -0,0 +1,37 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + L r3 + L r4 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r1 1 + c0 r2 -1 + c0 r3 1 + c1 Obj 1 + c1 r0 1 + c1 r1 -1 + c1 r2 -1 + c1 r3 1 + c2 Obj -3 + c2 r4 1 + c3 r0 1 + c3 r1 1 + c3 r2 1 + c3 r3 1 + c3 r4 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r3 2 + RHS_V r4 1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 +ENDATA diff --git a/datasets/mip/block_bve/or_used.mps b/datasets/mip/block_bve/or_used.mps new file mode 100644 index 0000000000..18b3789a91 --- /dev/null +++ b/datasets/mip/block_bve/or_used.mps @@ -0,0 +1,34 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + L r3 + G r4 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c0 r4 1 + c1 Obj 1 + c1 r1 -1 + c1 r2 -1 + c1 r4 1 + c2 Obj -2 + c2 r3 1 + c3 r0 1 + c3 r1 1 + c3 r2 1 + c3 r3 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r3 1 + RHS_V r4 1 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 +ENDATA diff --git a/datasets/mip/block_bve/random_a.mps b/datasets/mip/block_bve/random_a.mps new file mode 100644 index 0000000000..b951a629e1 --- /dev/null +++ b/datasets/mip/block_bve/random_a.mps @@ -0,0 +1,37 @@ +NAME +ROWS + N Obj + G r0 + L r1 + G r2 + G r3 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj -1 + c0 r1 -2 + c0 r2 2 + c1 Obj 2 + c1 r3 -1 + c2 Obj -2 + c2 r2 -1 + c2 r3 2 + c3 r0 -2 + c3 r2 -2 + c4 Obj -2 + c4 r1 1 + c4 r2 -1 + c4 r3 -2 + c5 Obj 1 + c5 r0 2 + c5 r3 1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r1 -2 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 + BV BOUND c5 +ENDATA diff --git a/datasets/mip/block_bve/random_b.mps b/datasets/mip/block_bve/random_b.mps new file mode 100644 index 0000000000..f81788b761 --- /dev/null +++ b/datasets/mip/block_bve/random_b.mps @@ -0,0 +1,54 @@ +NAME +ROWS + N Obj + G r0 + L r1 + L r2 + G r3 + G r4 + G r5 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj -2 + c0 r0 2 + c1 Obj -2 + c1 r1 -1 + c2 Obj -2 + c2 r0 1 + c2 r1 -1 + c2 r3 2 + c3 r3 2 + c3 r4 2 + c4 Obj -1 + c4 r1 -2 + c4 r2 2 + c4 r5 2 + c5 r1 1 + c5 r2 -1 + c5 r4 1 + c6 r2 1 + c6 r3 -1 + c6 r4 -1 + c6 r5 1 + c7 Obj 2 + c7 r2 2 + c7 r4 2 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r0 -1 + RHS_V r1 -1 + RHS_V r2 5 + RHS_V r4 2 + RHS_V r5 1 +RANGES + RANGE r2 2 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 + BV BOUND c5 + BV BOUND c6 + BV BOUND c7 +ENDATA diff --git a/datasets/mip/block_bve/random_c.mps b/datasets/mip/block_bve/random_c.mps new file mode 100644 index 0000000000..83aaee35e0 --- /dev/null +++ b/datasets/mip/block_bve/random_c.mps @@ -0,0 +1,48 @@ +NAME +ROWS + N Obj + L r0 + L r1 + L r2 + G r3 + L r4 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj -1 + c0 r3 2 + c1 Obj 2 + c1 r1 -2 + c1 r3 -1 + c2 Obj 2 + c2 r4 1 + c3 Obj -1 + c3 r0 -1 + c3 r2 2 + c4 r0 -1 + c4 r1 1 + c4 r2 2 + c5 Obj 2 + c5 r0 2 + c5 r2 2 + c5 r4 2 + c6 Obj 1 + c6 r0 -1 + c6 r4 2 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r1 1 + RHS_V r2 3 + RHS_V r3 -2 + RHS_V r4 5 +RANGES + RANGE r0 2 + RANGE r4 3 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 + BV BOUND c5 + BV BOUND c6 +ENDATA diff --git a/datasets/mip/block_bve/two_gadgets.mps b/datasets/mip/block_bve/two_gadgets.mps new file mode 100644 index 0000000000..7bb2e51137 --- /dev/null +++ b/datasets/mip/block_bve/two_gadgets.mps @@ -0,0 +1,49 @@ +NAME +ROWS + N Obj + G r0 + G r1 + L r2 + L r3 + L r4 + G r5 + G r6 + G r7 +COLUMNS + MARK0000 'MARKER' 'INTORG' + c0 Obj 1 + c0 r0 -1 + c0 r2 -1 + c0 r7 1 + c1 Obj 1 + c1 r1 -1 + c1 r2 -1 + c1 r7 1 + c2 Obj 1 + c2 r3 -1 + c2 r5 -1 + c2 r7 1 + c3 Obj 1 + c3 r4 -1 + c3 r5 -1 + c3 r7 1 + c4 r0 1 + c4 r1 1 + c4 r2 1 + c4 r6 1 + c5 r3 1 + c5 r4 1 + c5 r5 1 + c5 r6 -1 + MARK0001 'MARKER' 'INTEND' +RHS + RHS_V r5 -1 + RHS_V r7 2 +BOUNDS + BV BOUND c0 + BV BOUND c1 + BV BOUND c2 + BV BOUND c3 + BV BOUND c4 + BV BOUND c5 +ENDATA diff --git a/skills/cuopt-developer/SKILL.md b/skills/cuopt-developer/SKILL.md index c856538856..5fe1b28464 100644 --- a/skills/cuopt-developer/SKILL.md +++ b/skills/cuopt-developer/SKILL.md @@ -171,6 +171,7 @@ cuopt/ - Keep operations stream-ordered - Follow existing RAFT/RMM patterns - 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) ## Build & Test @@ -224,7 +225,7 @@ For pre-commit setup, DCO sign-off (`git commit -s`), the fork-based PR workflow ## Coding Conventions -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`), test-impact rules, and 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), see [references/conventions.md](references/conventions.md). +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, and 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), see [references/conventions.md](references/conventions.md). ## OpenMP task/runtime compatibility diff --git a/skills/cuopt-developer/references/conventions.md b/skills/cuopt-developer/references/conventions.md index e9963d4824..5c5bcc133d 100644 --- a/skills/cuopt-developer/references/conventions.md +++ b/skills/cuopt-developer/references/conventions.md @@ -126,12 +126,51 @@ signed subtraction (`std::vector v(static_cast(hi - lo) + 2, 0)`), the narrowing `size_t`→`i_t` in `static_cast(x.size())` (established style; keep it) +### Integer widths — prefer fixed-width types + +Prefer `` fixed-width types (`int32_t`, `int64_t`, `uint32_t`, …) over +plain `int` / `long` / `long long` when the value range or ABI width matters +(counts that can exceed 32 bits, device grid math, work estimates, file offsets). + +Avoid multi-word functional casts such as `long long(x)` in `.cu`/`.cuh` — they +confuse CUDA-aware tooling (`type name is not allowed`). Use a C-style cast to a +fixed-width type instead: `(int64_t)x`. + +```cpp +const long long total = long long(num) * long long(patterns); // ❌ +const int64_t total = (int64_t)num * (int64_t)patterns; // ✅ +``` + +Keep `i_t` / `f_t` for problem-index and numeric template parameters; use +`int64_t` (etc.) for host-side wide counters outside that abstraction. + ### CUDA Error Checking ```cpp RAFT_CUDA_TRY(cudaMemcpy(...)); ``` +### Prefer modern CCCL utilities in device code + +When writing or editing CUDA kernels, prefer CCCL / libcu++ helpers over hand-rolled +bit math, reductions, or integer tricks. They encode the PTX-friendly form and avoid +boilerplate that compilers often fail to recover from runtime values. + +Examples (CUDA 13 / CCCL 3.x era — headers already used elsewhere in `cpp/src`): + +| Need | Prefer | Instead of | +|------|--------|------------| +| Extract a bitfield / decode packed indices | `cuda::bitfield_extract` (``) | `%` / `/` by a runtime `1 << k` (nvcc usually will not strength-reduce those to mask/shift) | +| Build a contiguous bit mask | `cuda::bitmask` | Hand-written `((1u << w) - 1u) << start` | +| Test / round to power of two | `cuda::is_power_of_two`, `next_power_of_two`, `prev_power_of_two` (``), or `cuda::std::has_single_bit` / `bit_ceil` / `bit_floor` (``) | Ad-hoc `(x & (x - 1)) == 0` / manual ceil loops | +| Divide/mod by a value that is constant for a launch (or across many ops) but not a compile-time constant | `cuda::fast_mod_div` (``) — construct on the host (or once), pass into the kernel, use `/` `%` / `cuda::div` | Hot-path `idiv` / handwritten libdivide magic | +| Warp/block algorithms | CUB / CCCL / RAFT primitives already used in-tree | Homegrown shared-memory reductions when an existing primitive fits | + +Docs: [CCCL bit extensions](https://nvidia.github.io/cccl/unstable/libcudacxx/extended_api/bit.html), +[pow2 helpers](https://nvidia.github.io/cccl/unstable/libcudacxx/extended_api/math/pow2.html), +[`cuda::fast_mod_div`](https://nvidia.github.io/cccl/unstable/libcudacxx/extended_api/math/fast_mod_div.html). +Check signatures in the installed headers rather than guessing — APIs evolve with CCCL. + ## Memory Management ```cpp @@ -158,3 +197,15 @@ Read existing code in `cpp/src/` for real examples of RMM allocation, stream-ord - Python pytest: `python/.../tests/` **Add at least one regression test for new behavior.** + +When a new MIP test loads a MIPLIB instance (e.g. via `make_path_absolute("mip/.mps")`), +that instance must appear in `datasets/mip/download_miplib_test_dataset.sh`'s `INSTANCES` +list. CI and local setups only fetch that allowlist — an unlisted name fails at parse time +with a missing-file error even though the test itself is correct. Add the basename there as part of the same change that introduces the test. + +Calling cuOpt MIP internals that use OpenMP taskloops (notably `diversity_manager::run_presolve` +→ `compute_probing_cache`) from a plain gtest must open an OMP team first, the same way +`solve_mip` does (`#pragma omp parallel num_threads(...)` + `#pragma omp masked`, with +`omp_set_max_active_levels(2)` if needed). Probing sizes its pool as +`omp_get_num_threads() - 1`; outside a parallel region that is 0 and probing becomes a silent +no-op (Papilo size unchanged, test finishes in a few hundred ms).