From 62481ca61b9c6eb9e826b83a023aa623df86bfb9 Mon Sep 17 00:00:00 2001 From: yboucher Date: Fri, 14 Aug 2026 13:49:24 -0700 Subject: [PATCH 1/9] initial impl --- .../diversity/diversity_manager.cu | 57 ++-- .../feasibility_jump/feasibility_jump.cuh | 6 + .../mip_heuristics/feasibility_jump/fj_cpu.cu | 304 +++++++++++++++++- .../feasibility_jump/fj_cpu.cuh | 14 + .../local_search/local_search.cu | 5 +- .../presolve/conflict_graph/clique_table.cu | 32 +- .../presolve/conflict_graph/clique_table.cuh | 44 +++ cpp/src/mip_heuristics/problem/problem.cu | 3 + 8 files changed, 397 insertions(+), 68 deletions(-) diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index 9d70ae17ee..d6c76df071 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -288,6 +288,9 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ raft::common::nvtx::range fun_scope("run_presolve"); CUOPT_LOG_INFO("\nRunning cuOpt presolve"); timer_t presolve_timer(time_limit); + // Presolve remaps variable ids, so a conflict graph from an earlier presolve of this problem no + // longer describes it. It is rebuilt at the end of this function. + problem_ptr->clique_table.reset(); auto term_crit = ls.constraint_prop.bounds_update.solve(*problem_ptr); if (ls.constraint_prop.bounds_update.infeas_constraints_count > 0) { @@ -328,32 +331,24 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ 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); } if (!problem_ptr->empty && !check_bounds_sanity(*problem_ptr)) { return false; } - // if (!presolve_timer.check_time_limit() && !context.settings.heuristics_only && - // !problem_ptr->empty) { - // f_t time_limit_for_clique_table = std::min(3., presolve_timer.remaining_time() / 5); - // timer_t clique_timer(time_limit_for_clique_table); - // simplex::user_problem_t host_problem(problem_ptr->handle_ptr); - // problem_ptr->get_host_user_problem(host_problem); - // std::shared_ptr> clique_table; - // constexpr bool modify_problem_with_cliques = false; - // find_initial_cliques(host_problem, - // context.settings.tolerances, - // &clique_table, - // clique_timer, - // modify_problem_with_cliques, - // nullptr); - // if (modify_problem_with_cliques) { - // problem_ptr->set_constraints_from_host_user_problem(host_problem); - // cuopt_assert(host_problem.lower.size() == static_cast(problem_ptr->n_variables), - // "host lower bound size mismatch"); - // cuopt_assert(host_problem.upper.size() == static_cast(problem_ptr->n_variables), - // "host upper bound size mismatch"); - // std::vector all_var_indices(problem_ptr->n_variables); - // std::iota(all_var_indices.begin(), all_var_indices.end(), 0); - // problem_ptr->update_variable_bounds(all_var_indices, host_problem.lower, - // host_problem.upper); trivial_presolve(*problem_ptr, remap_cache_ids); - // } - // } + // Build the conflict graph once the problem has its final shape: it is indexed by post-presolve + // variable ids. Consumed by the CPU FJ compound binary move and, through solver.cu, by the B&B + // clique and zero-half cut separators. The problem itself is left untouched. + if (!problem_ptr->empty && !presolve_timer.check_time_limit() && + !global_timer.check_time_limit()) { + const f_t clique_time_limit = std::min({(f_t)3., + (f_t)presolve_timer.remaining_time(), + (f_t)global_timer.remaining_time()}); + timer_t clique_timer(clique_time_limit); + simplex::user_problem_t host_problem(problem_ptr->handle_ptr); + problem_ptr->get_host_user_problem(host_problem); + find_initial_cliques( + host_problem, context.settings.tolerances, problem_ptr->clique_table, clique_timer); + // Publishing an empty table would make B&B skip its own build and lose clique cuts entirely. + if (problem_ptr->clique_table != nullptr && problem_ptr->clique_table->empty()) { + problem_ptr->clique_table.reset(); + } + } // May overconstrain if Papilo presolve has been run before if (context.settings.presolver == presolver_t::None) { if (!problem_ptr->empty) { @@ -392,19 +387,17 @@ void diversity_manager_t::generate_quick_feasible_solution() ls.generate_fast_solution(solution, sol_timer); if (solution.get_feasible()) { population.run_solution_callbacks(solution); - initial_sol_vector.emplace_back(std::move(solution)); problem_ptr->handle_ptr->sync_stream(); - solution_t searched_sol(initial_sol_vector.back()); + solution_t searched_sol(solution); ls_config_t ls_config; run_local_search(searched_sol, population.weights, sol_timer, ls_config); population.run_solution_callbacks(searched_sol); - initial_sol_vector.emplace_back(std::move(searched_sol)); - auto& feas_sol = initial_sol_vector.back().get_feasible() - ? initial_sol_vector.back() - : initial_sol_vector[initial_sol_vector.size() - 2]; + auto& feas_sol = searched_sol.get_feasible() ? searched_sol : solution; CUOPT_LOG_INFO("Generated fast solution in %f seconds with objective %f", timer.elapsed_time(), feas_sol.get_user_objective()); + population.add_solution(std::move(searched_sol)); + population.add_solution(std::move(solution)); } problem_ptr->handle_ptr->sync_stream(); } diff --git a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh index 8d1f39ce22..e9444c423a 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh @@ -75,6 +75,12 @@ struct fj_hyper_parameters_t { double small_move_tabu_threshold = 1e-6; int small_move_tabu_tenure = 4; + // binary 2-opt caps: violated rows sampled per local minimum, candidate variables taken from each + // sampled row, and total candidate pairs scored + int two_opt_max_rows = 4; + int two_opt_max_row_vars = 256; + int two_opt_max_pairs = 256; + // load-balancing related settings int old_codepath_total_var_to_relvar_ratio_threshold = 200; int load_balancing_codepath_min_varcount = 3200; diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 57a6a89479..8d92936786 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -15,6 +15,8 @@ #include "fj_cpu.cuh" #include "fj_cpu_worker.cuh" +#include + #include #include @@ -132,7 +134,6 @@ thrust::tuple get_mtm_for_constraint( template std::pair feas_score_constraint(const typename fj_t::climber_data_t::view_t& fj, - i_t var_idx, f_t delta, i_t cstr_idx, f_t cstr_coeff, @@ -613,7 +614,6 @@ static inline std::pair compute_score(fj_cpu_climber_t(fj_cpu.view, - var_idx, delta, cstr_idx, cstr_coeff, @@ -649,6 +649,273 @@ static inline std::pair compute_score(fj_cpu_climber_t::max()}; +}; + +// Deterministic total order over 2-opt candidates: higher score, then older last touch, then +// smaller variable indices. +static bool two_opt_cand_better(const two_opt_move_t& a, const two_opt_move_t& b) +{ + if (a.score != b.score) return a.score > b.score; + if (a.age != b.age) return a.age < b.age; + if (a.first.var_idx != b.first.var_idx) return a.first.var_idx < b.first.var_idx; + return a.second.var_idx < b.second.var_idx; +} + +/** + * @brief Score the combined effect of flipping two binaries at once. + * + * Scoring each flip on its own and summing would double count the rows both variables appear in, + * and would miss the case the compound move exists for: the two contributions cancelling out in a + * shared row. So the per-row LHS deltas of both flips are gathered, entries of the same row merged, + * and every touched row scored once from its aggregated delta. + */ +template +static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& fj_cpu, + i_t first, + f_t first_delta, + i_t second, + f_t second_delta) +{ + auto& row_deltas = fj_cpu.two_opt_row_deltas; + row_deltas.clear(); + auto collect = [&](i_t var_idx, f_t delta) { + const auto [offset_begin, offset_end] = reverse_range_for_var(fj_cpu, var_idx); + fj_cpu.nnz_processed_window += offset_end - offset_begin; + for (i_t i = offset_begin; i < offset_end; ++i) { + const i_t cstr_idx = fj_cpu.h_reverse_constraints[i]; + const f_t coeff = fj_cpu.h_reverse_coefficients[i]; + row_deltas.emplace_back(cstr_idx, coeff * delta); + } + }; + collect(first, first_delta); + collect(second, second_delta); + // Brings the entries of a shared row next to each other + std::sort(row_deltas.begin(), row_deltas.end()); + + f_t base_feas_sum = 0; + f_t bonus_robust_sum = 0; + for (size_t pos = 0; pos < row_deltas.size();) { + const i_t cstr_idx = row_deltas[pos].first; + f_t lhs_delta = 0; + do { + lhs_delta += row_deltas[pos++].second; + } while (pos < row_deltas.size() && row_deltas[pos].first == cstr_idx); + + // The coefficients are already folded into lhs_delta, hence the unit coefficient + auto [cstr_base_feas, cstr_bonus_robust] = + feas_score_constraint(fj_cpu.view, + lhs_delta, + cstr_idx, + 1, + fj_cpu.h_cstr_lb[cstr_idx], + fj_cpu.h_cstr_ub[cstr_idx], + fj_cpu.h_lhs[cstr_idx], + fj_cpu.h_cstr_left_weights[cstr_idx], + fj_cpu.h_cstr_right_weights[cstr_idx]); + base_feas_sum += cstr_base_feas; + bonus_robust_sum += cstr_bonus_robust; + } + + // Same staged score convention as the single variable path + const f_t obj_diff = + fj_cpu.h_obj_coeffs[first] * first_delta + fj_cpu.h_obj_coeffs[second] * second_delta; + f_t base_obj = 0; + if (obj_diff < 0) + base_obj = fj_cpu.h_objective_weight; + else if (obj_diff > 0) + base_obj = -fj_cpu.h_objective_weight; + + f_t bonus_breakthrough = 0; + bool old_obj_better = fj_cpu.h_incumbent_objective < fj_cpu.h_best_objective; + bool new_obj_better = fj_cpu.h_incumbent_objective + obj_diff < fj_cpu.h_best_objective; + if (!old_obj_better && new_obj_better) + bonus_breakthrough += fj_cpu.h_objective_weight; + else if (old_obj_better && !new_obj_better) + bonus_breakthrough -= fj_cpu.h_objective_weight; + + fj_staged_score_t score; + score.base = round(base_obj + base_feas_sum); + score.bonus = round(bonus_breakthrough + bonus_robust_sum); + return score; +} + +/** + * @brief Fill fj_cpu.two_opt_partners with candidates to flip together with `first`. + * + * Preferred source is the conflict graph: the literals conflicting with the one `first` is about to + * make true are exactly the ones that must become false, which gives both the partner and the value + * to move it to. A conflict with a negated literal asks for a same-direction pair, which is how + * covering rows contribute. When the graph has nothing for `first`, the variables sharing a row + * with it are the fallback, and there the only defensible direction is the swap. + */ +template +static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, + i_t first, + f_t first_delta, + size_t max_partners) +{ + auto& partners = fj_cpu.two_opt_partners; + const i_t n_variables = fj_cpu.view.pb.n_variables; + partners.clear(); + cuopt_assert(fj_cpu.h_is_binary_variable[first], "2-opt is only defined for binaries"); + + auto add_partner = [&](i_t var_idx, f_t target) { + if (var_idx == first) return; + const f_t val = fj_cpu.h_assignment[var_idx].get(); + // A partner between two integers has no opposite value to swap to + if (!fj_cpu.view.pb.is_integer(val)) return; + const f_t delta = target - val; + // Already at the value we would move it to, so there is no compound move to make + if (fabs(delta) < 0.5) return; + if (!check_variable_within_bounds(fj_cpu, var_idx, target)) return; + if (tabu_check(fj_cpu, var_idx, delta, true)) return; + partners.emplace_back(var_idx, delta); + }; + + if (fj_cpu.clique_table != nullptr) { + cuopt_assert(fj_cpu.clique_table->n_variables == n_variables, + "conflict graph belongs to another problem"); + const f_t new_val = fj_cpu.h_assignment[first].get() + first_delta; + // Literal `v` stands for "variable v is 1", literal `v + n_variables` for "variable v is 0" + const i_t literal = new_val > 0.5 ? first : first + n_variables; + fj_cpu.clique_table->for_each_conflict_partner(literal, [&](i_t partner_literal) { + // The conflicting literal has to become false + add_partner(partner_literal % n_variables, partner_literal < n_variables ? 0 : 1); + return partners.size() < max_partners; + }); + // A partner is yielded once per shared clique + std::sort(partners.begin(), partners.end()); + partners.erase(std::unique(partners.begin(), partners.end()), partners.end()); + } + if (!partners.empty()) return; + + const auto& related = fj_cpu.h_related_variables; + const auto& related_offsets = fj_cpu.h_related_variables_offsets; + if (related_offsets.size() != static_cast(n_variables) + 1) return; + // Row sharing carries no polarity, so take the value `first` is vacating: only a partner holding + // the opposite value then moves in the compensating direction, and one holding the same value is + // filtered out as a no-op. Same opposite-value rule as the GPU candidate build. + const f_t swap_target = fj_cpu.h_assignment[first].get(); + const i_t related_begin = related_offsets[first]; + const i_t related_end = related_offsets[first + 1]; + for (i_t i = related_begin; i < related_end && partners.size() < max_partners; ++i) { + const i_t var_idx = related[i]; + if (fj_cpu.h_is_binary_variable[var_idx]) { add_partner(var_idx, swap_target); } + } +} + +/** + * @brief Look for an improving simultaneous flip of two binaries (binary 2-opt). + * + * At a local minimum no single flip improves the score, but a pair often does: in set partitioning + * style rows the only way out is to turn one variable off and another on in the same step. The + * neighbourhood is sampled rather than enumerated, which would be quadratic. First variables come + * from a few violated rows, or from objective variables whose flip improves the objective when + * nothing is violated; partners come from two_opt_collect_partners. The search stops at + * two_opt_max_pairs candidates or once it has touched nnz_samples nonzeros, whichever comes first: + * the pair count bounds the neighbourhood, the nonzero count bounds the cost, which here is + * proportional to the degrees of both endpoints. Returns an invalid-scored move when no pair was + * worth applying. + */ +template +static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) +{ + CPUFJ_NVTX_RANGE("CPUFJ::find_two_opt_move"); + constexpr size_t max_obj_starts = 64; + // The GPU candidate table is indexed by pair, so it has no per-first-variable level to cap + constexpr size_t max_partners_per_var = 16; + + const auto& params = fj_cpu.settings.parameters; + const size_t max_target_rows = params.two_opt_max_rows; + const size_t max_first_vars = params.two_opt_max_row_vars; + const size_t max_pairs = params.two_opt_max_pairs; + + two_opt_move_t best; + // Climbers built from a host LP have neither partner source, so there is nothing to pair with + const bool partner_source_exists = + fj_cpu.clique_table != nullptr || + fj_cpu.h_related_variables_offsets.size() == + static_cast(fj_cpu.view.pb.n_variables) + 1; + if (fj_cpu.n_binary_vars == 0 || !partner_source_exists) return best; + + std::mt19937 rng(fj_cpu.settings.seed + fj_cpu.iterations); + auto& first_vars = fj_cpu.two_opt_first_vars; + first_vars.clear(); + + if (!fj_cpu.violated_constraints.empty()) { + auto& target_cstrs = fj_cpu.two_opt_target_cstrs; + target_cstrs.clear(); + std::sample(fj_cpu.violated_constraints.begin(), + fj_cpu.violated_constraints.end(), + std::back_inserter(target_cstrs), + max_target_rows, + rng); + for (i_t cstr_idx : target_cstrs) { + const auto [offset_begin, offset_end] = range_for_constraint(fj_cpu, cstr_idx); + for (i_t i = offset_begin; i < offset_end && first_vars.size() < max_first_vars; ++i) { + const i_t var_idx = fj_cpu.h_variables[i]; + if (fj_cpu.h_is_binary_variable[var_idx]) first_vars.push_back(var_idx); + } + } + } else { + std::sample(fj_cpu.h_objective_vars.underlying().begin(), + fj_cpu.h_objective_vars.underlying().end(), + std::back_inserter(first_vars), + max_obj_starts, + rng); + // Nothing is violated, so a pair can only help by improving the objective: keep the flips that + // move it down and let the pair scoring pay for the feasibility damage. + first_vars.erase(std::remove_if(first_vars.begin(), + first_vars.end(), + [&](i_t var_idx) { + if (!fj_cpu.h_is_binary_variable[var_idx]) return true; + const f_t delta = + round(1 - 2 * fj_cpu.h_assignment[var_idx].get()); + return fj_cpu.h_obj_coeffs[var_idx] * delta >= 0; + }), + first_vars.end()); + } + std::shuffle(first_vars.begin(), first_vars.end(), rng); + + const i_t nnz_at_entry = fj_cpu.nnz_processed_window; + size_t pairs_scored = 0; + for (i_t first : first_vars) { + if (pairs_scored >= max_pairs) break; + if (fj_cpu.nnz_processed_window - nnz_at_entry > fj_cpu.nnz_samples) break; + const f_t first_val = fj_cpu.h_assignment[first].get(); + if (!fj_cpu.view.pb.is_integer(first_val)) continue; + const f_t first_delta = round(1 - 2 * first_val); + if (tabu_check(fj_cpu, first, first_delta, true)) continue; + if (!check_variable_within_bounds(fj_cpu, first, first_val + first_delta)) continue; + const i_t first_inc = fj_cpu.h_tabu_lastinc[first]; + const i_t first_dec = fj_cpu.h_tabu_lastdec[first]; + const i_t first_touch = std::max(first_inc, first_dec); + + two_opt_collect_partners(fj_cpu, first, first_delta, max_partners_per_var); + for (const auto& [second, second_delta] : fj_cpu.two_opt_partners) { + const i_t second_inc = fj_cpu.h_tabu_lastinc[second]; + const i_t second_dec = fj_cpu.h_tabu_lastdec[second]; + two_opt_move_t cand; + cand.first = {first, first_delta}; + cand.second = {second, second_delta}; + cand.score = two_opt_compute_pair_score(fj_cpu, first, first_delta, second, second_delta); + cand.age = std::max(first_touch, std::max(second_inc, second_dec)); + if (two_opt_cand_better(cand, best)) { best = cand; } + ++pairs_scored; + if (pairs_scored >= max_pairs) return best; + if (fj_cpu.nnz_processed_window - nnz_at_entry > fj_cpu.nnz_samples) return best; + } + } + return best; +} + template static void smooth_weights(fj_cpu_climber_t& fj_cpu) { @@ -1272,6 +1539,17 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, fj_cpu.h_is_binary_variable = cuopt::host_copy(problem.is_binary_variable, handle_ptr->get_stream()); fj_cpu.h_binary_indices = cuopt::host_copy(problem.binary_indices, handle_ptr->get_stream()); + fj_cpu.h_related_variables = + cuopt::host_copy(problem.related_variables, handle_ptr->get_stream()); + fj_cpu.h_related_variables_offsets = + cuopt::host_copy(problem.related_variables_offsets, handle_ptr->get_stream()); + fj_cpu.clique_table = problem.clique_table; + // A problem built by fixing variables drops the graph, so a size mismatch means it outlived the + // index space it describes. Refuse to use it rather than index out of range. + if (fj_cpu.clique_table != nullptr && fj_cpu.clique_table->n_variables != problem.n_variables) { + cuopt_assert(false, "conflict graph does not match the problem"); + fj_cpu.clique_table.reset(); + } fj_cpu.h_cstr_left_weights = left_weights; fj_cpu.h_cstr_right_weights = right_weights; @@ -1691,11 +1969,23 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit, double w for (size_t i = 0; i < fj_cpu->cached_mtm_moves.size(); i++) fj_cpu->cached_mtm_moves[i].first = 0; } - thrust::tie(move, score) = - find_mtm_move_viol(*fj_cpu, 1, true); // pick a single random violated constraint - i_t var_idx = move.var_idx >= 0 ? move.var_idx : 0; - f_t delta = move.var_idx >= 0 ? move.value : 0; - apply_move(*fj_cpu, var_idx, delta, true); + // A simultaneous flip of two binaries escapes minima that no single flip can. Not attempted + // right after a perturbation, whose whole point is to leave the current region. + two_opt_move_t two_opt_move; + if (!should_perturb) two_opt_move = find_two_opt_move(*fj_cpu); + if (two_opt_move.score > fj_staged_score_t::zero()) { + // Applied as two moves; they were scored jointly, so the intermediate state after the first + // one may well be worse than the local minimum we came from + apply_move(*fj_cpu, two_opt_move.first.var_idx, two_opt_move.first.value, true); + apply_move(*fj_cpu, two_opt_move.second.var_idx, two_opt_move.second.value, true); + fj_cpu->n_mtm_viol_moves_window += 2; + } else { + thrust::tie(move, score) = + find_mtm_move_viol(*fj_cpu, 1, true); // pick a single random violated constraint + i_t var_idx = move.var_idx >= 0 ? move.var_idx : 0; + f_t delta = move.var_idx >= 0 ? move.value : 0; + apply_move(*fj_cpu, var_idx, delta, true); + } ++local_mins; ++fj_cpu->n_local_minima_window; } diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index 718c89615d..58d10b5c0d 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -44,6 +44,8 @@ struct fj_cpu_climber_t { ADD_INSTRUMENTED(h_is_binary_variable), ADD_INSTRUMENTED(h_objective_vars), ADD_INSTRUMENTED(h_binary_indices), + ADD_INSTRUMENTED(h_related_variables), + ADD_INSTRUMENTED(h_related_variables_offsets), ADD_INSTRUMENTED(h_tabu_nodec_until), ADD_INSTRUMENTED(h_tabu_noinc_until), ADD_INSTRUMENTED(h_tabu_lastdec), @@ -83,6 +85,12 @@ struct fj_cpu_climber_t { ins_vector h_is_binary_variable; ins_vector h_objective_vars; ins_vector h_binary_indices; + ins_vector h_related_variables; + ins_vector h_related_variables_offsets; + // Conflict graph of the problem this climber was built from, null when presolve built none. Held + // by shared_ptr and snapshotted alongside the host copies above because it is indexed by that + // problem's variable ids. + std::shared_ptr> clique_table; ins_vector h_tabu_nodec_until; ins_vector h_tabu_noinc_until; @@ -134,6 +142,12 @@ struct fj_cpu_climber_t { std::vector var_bitmap; ins_vector iter_mtm_vars; + // Scratch reused by the binary 2-opt search, which runs at every local minimum + std::vector two_opt_target_cstrs; + std::vector two_opt_first_vars; + std::vector> two_opt_partners; + std::vector> two_opt_row_deltas; + i_t mtm_viol_samples{25}; i_t mtm_sat_samples{15}; i_t nnz_samples{50000}; diff --git a/cpp/src/mip_heuristics/local_search/local_search.cu b/cpp/src/mip_heuristics/local_search/local_search.cu index 75c4185949..92be0822d3 100644 --- a/cpp/src/mip_heuristics/local_search/local_search.cu +++ b/cpp/src/mip_heuristics/local_search/local_search.cu @@ -145,8 +145,11 @@ void local_search_t::stop_cpufj_scratch_threads() { if (omp_get_num_threads() < CUOPT_MIP_FJ_REQUIRED_THREAD_COUNT) return; + for (auto& cpu_fj : scratch_cpu_fj) { + cuopt_assert(cpu_fj != nullptr, "scratch climbers must have been created"); + cpu_fj->halted = true; + } for (size_t i = 0; i < scratch_cpu_fj.size(); ++i) { - scratch_cpu_fj[i]->halted = true; #pragma omp taskwait depend(in : *scratch_cpu_fj[i]) // Wait for each scratch CPU FJ task to finish } diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu index a8e6997572..fb02c5bbc0 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu @@ -314,34 +314,10 @@ template std::unordered_set clique_table_t::get_adj_set_of_var(i_t var_idx) const { std::unordered_set adj_set; - - // First-clique edges: every member of each first-clique containing var_idx. - for (i_t clique_idx : var_clique_first.slice(var_idx)) { - const auto& c = first[clique_idx]; - adj_set.insert(c.begin(), c.end()); - } - - // Addtl-clique edges. - for (i_t addtl_idx : var_clique_addtl.slice(var_idx)) { - const auto& a = addtl_cliques[addtl_idx]; - if (a.vertex_idx == var_idx) { - // var_idx is the extension vertex; new neighbors are the base suffix. - const auto& base = first[a.clique_idx]; - adj_set.insert(base.begin() + a.start_pos_on_clique, base.end()); - } else { - // var_idx is a base member; only new edge is to the extension vertex. - adj_set.insert(a.vertex_idx); - } - } - - for (i_t adj_var : small_clique_adj.slice(var_idx)) { - adj_set.insert(adj_var); - } - - // Add the complement of var_idx to the adjacency set - i_t complement_idx = (var_idx >= n_variables) ? (var_idx - n_variables) : (var_idx + n_variables); - adj_set.insert(complement_idx); - adj_set.erase(var_idx); + for_each_conflict_partner(var_idx, [&](i_t partner) { + adj_set.insert(partner); + return true; + }); return adj_set; } diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh index fe6db31c26..a12cfb2c1c 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh @@ -9,6 +9,7 @@ #include #include +#include #include #include @@ -174,6 +175,12 @@ struct clique_table_t { clique_table_t& operator=(clique_table_t&&) = delete; std::unordered_set get_adj_set_of_var(i_t var_idx) const; + // Yields the conflict partners of `literal`, also as literals, until `fn` returns false. A + // partner is yielded once per clique it shares with `literal`, so a caller that cannot tolerate + // repeats has to filter them. Nothing is materialized: a caller that wants only a few partners + // returns false early and pays for those only. + template + void for_each_conflict_partner(i_t literal, callable_t&& fn) const; i_t get_degree_of_var(i_t var_idx); bool check_adjacency(i_t var_idx1, i_t var_idx2) const; bool empty() const @@ -201,6 +208,43 @@ struct clique_table_t { typename mip_solver_settings_t::tolerances_t tolerances; }; +template +template +void clique_table_t::for_each_conflict_partner(i_t literal, callable_t&& fn) const +{ + cuopt_assert(literal >= 0 && literal < 2 * n_variables, "literal out of range"); + const i_t complement = + (literal >= n_variables) ? (literal - n_variables) : (literal + n_variables); + if (!fn(complement)) { return; } + + for (i_t clique_idx : var_clique_first.slice(literal)) { + for (i_t member : first[clique_idx]) { + if (member != literal && !fn(member)) { return; } + } + } + + for (i_t addtl_idx : var_clique_addtl.slice(literal)) { + const auto& addtl = addtl_cliques[addtl_idx]; + if (addtl.vertex_idx == literal) { + // literal is the extension vertex; its new neighbors are the base suffix + const auto& base = first[addtl.clique_idx]; + cuopt_assert(addtl.start_pos_on_clique >= 0 && + static_cast(addtl.start_pos_on_clique) <= base.size(), + "extension start position out of range"); + for (auto it = base.begin() + addtl.start_pos_on_clique; it != base.end(); ++it) { + if (*it != literal && !fn(*it)) { return; } + } + } else if (!fn(addtl.vertex_idx)) { + // literal is a base member; the only new edge is to the extension vertex + return; + } + } + + for (i_t member : small_clique_adj.slice(literal)) { + if (member != literal && !fn(member)) { return; } + } +} + // Builds the conflict-graph clique table for `problem`. The base cliques are // published to `clique_table_out` before the (optional, signal-gated) extension // phase begins, so cut generation can pick up the table while extension keeps diff --git a/cpp/src/mip_heuristics/problem/problem.cu b/cpp/src/mip_heuristics/problem/problem.cu index ccba2d5f2b..b8214831de 100644 --- a/cpp/src/mip_heuristics/problem/problem.cu +++ b/cpp/src/mip_heuristics/problem/problem.cu @@ -1862,6 +1862,9 @@ void problem_t::remove_given_variables(problem_t& original_p thrust::fill(handle_ptr->get_thrust_policy(), offsets.begin(), offsets.end(), 0); cuopt_assert(assignment.size() == n_variables, "Variable size mismatch"); cuopt_assert(variable_map.size() < n_variables, "Too many variables to fix"); + + // the clique table could be propagated after fixings. but to err on the safe side, reset it here + clique_table.reset(); rmm::device_uvector tmp_assignment(assignment, handle_ptr->get_stream()); // first remove the assignment and variable related vectors From 6cea550521781febf76f4c1ddda9b9a4c8ce0b5f Mon Sep 17 00:00:00 2001 From: yboucher Date: Tue, 18 Aug 2026 00:37:04 -0700 Subject: [PATCH 2/9] cleanup --- cpp/src/mip_heuristics/diversity/diversity_manager.cu | 2 -- cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu | 2 +- .../mip_heuristics/presolve/conflict_graph/clique_table.cuh | 5 +---- 3 files changed, 2 insertions(+), 7 deletions(-) diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index d6c76df071..5472d97f37 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -288,8 +288,6 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ raft::common::nvtx::range fun_scope("run_presolve"); CUOPT_LOG_INFO("\nRunning cuOpt presolve"); timer_t presolve_timer(time_limit); - // Presolve remaps variable ids, so a conflict graph from an earlier presolve of this problem no - // longer describes it. It is rebuilt at the end of this function. problem_ptr->clique_table.reset(); auto term_crit = ls.constraint_prop.bounds_update.solve(*problem_ptr); diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 8d92936786..5964fb346e 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -838,7 +838,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) const size_t max_pairs = params.two_opt_max_pairs; two_opt_move_t best; - // Climbers built from a host LP have neither partner source, so there is nothing to pair with + const bool partner_source_exists = fj_cpu.clique_table != nullptr || fj_cpu.h_related_variables_offsets.size() == diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh index a12cfb2c1c..d58c99e304 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh @@ -175,10 +175,7 @@ struct clique_table_t { clique_table_t& operator=(clique_table_t&&) = delete; std::unordered_set get_adj_set_of_var(i_t var_idx) const; - // Yields the conflict partners of `literal`, also as literals, until `fn` returns false. A - // partner is yielded once per clique it shares with `literal`, so a caller that cannot tolerate - // repeats has to filter them. Nothing is materialized: a caller that wants only a few partners - // returns false early and pays for those only. + // Yields the conflict partners of literal, also as literals, until fn returns false template void for_each_conflict_partner(i_t literal, callable_t&& fn) const; i_t get_degree_of_var(i_t var_idx); From 1960c1f91fe60623292806a349bcea123ba461c5 Mon Sep 17 00:00:00 2001 From: yboucher Date: Tue, 18 Aug 2026 04:06:52 -0700 Subject: [PATCH 3/9] cleanup --- .../diversity/diversity_manager.cu | 55 ++++---- .../feasibility_jump/feasibility_jump.cuh | 4 + .../mip_heuristics/feasibility_jump/fj_cpu.cu | 120 +++++++++++++----- .../feasibility_jump/fj_cpu.cuh | 22 +++- .../local_search/local_search.cu | 11 +- .../presolve/conflict_graph/clique_table.cu | 32 ++++- .../presolve/conflict_graph/clique_table.cuh | 41 ------ .../mip_heuristics/presolve/probing_cache.cuh | 2 +- cpp/src/mip_heuristics/problem/problem.cu | 3 - 9 files changed, 177 insertions(+), 113 deletions(-) diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index 5472d97f37..9d70ae17ee 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -288,7 +288,6 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ raft::common::nvtx::range fun_scope("run_presolve"); CUOPT_LOG_INFO("\nRunning cuOpt presolve"); timer_t presolve_timer(time_limit); - problem_ptr->clique_table.reset(); auto term_crit = ls.constraint_prop.bounds_update.solve(*problem_ptr); if (ls.constraint_prop.bounds_update.infeas_constraints_count > 0) { @@ -329,24 +328,32 @@ bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_ 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); } if (!problem_ptr->empty && !check_bounds_sanity(*problem_ptr)) { return false; } - // Build the conflict graph once the problem has its final shape: it is indexed by post-presolve - // variable ids. Consumed by the CPU FJ compound binary move and, through solver.cu, by the B&B - // clique and zero-half cut separators. The problem itself is left untouched. - if (!problem_ptr->empty && !presolve_timer.check_time_limit() && - !global_timer.check_time_limit()) { - const f_t clique_time_limit = std::min({(f_t)3., - (f_t)presolve_timer.remaining_time(), - (f_t)global_timer.remaining_time()}); - timer_t clique_timer(clique_time_limit); - simplex::user_problem_t host_problem(problem_ptr->handle_ptr); - problem_ptr->get_host_user_problem(host_problem); - find_initial_cliques( - host_problem, context.settings.tolerances, problem_ptr->clique_table, clique_timer); - // Publishing an empty table would make B&B skip its own build and lose clique cuts entirely. - if (problem_ptr->clique_table != nullptr && problem_ptr->clique_table->empty()) { - problem_ptr->clique_table.reset(); - } - } + // if (!presolve_timer.check_time_limit() && !context.settings.heuristics_only && + // !problem_ptr->empty) { + // f_t time_limit_for_clique_table = std::min(3., presolve_timer.remaining_time() / 5); + // timer_t clique_timer(time_limit_for_clique_table); + // simplex::user_problem_t host_problem(problem_ptr->handle_ptr); + // problem_ptr->get_host_user_problem(host_problem); + // std::shared_ptr> clique_table; + // constexpr bool modify_problem_with_cliques = false; + // find_initial_cliques(host_problem, + // context.settings.tolerances, + // &clique_table, + // clique_timer, + // modify_problem_with_cliques, + // nullptr); + // if (modify_problem_with_cliques) { + // problem_ptr->set_constraints_from_host_user_problem(host_problem); + // cuopt_assert(host_problem.lower.size() == static_cast(problem_ptr->n_variables), + // "host lower bound size mismatch"); + // cuopt_assert(host_problem.upper.size() == static_cast(problem_ptr->n_variables), + // "host upper bound size mismatch"); + // std::vector all_var_indices(problem_ptr->n_variables); + // std::iota(all_var_indices.begin(), all_var_indices.end(), 0); + // problem_ptr->update_variable_bounds(all_var_indices, host_problem.lower, + // host_problem.upper); trivial_presolve(*problem_ptr, remap_cache_ids); + // } + // } // May overconstrain if Papilo presolve has been run before if (context.settings.presolver == presolver_t::None) { if (!problem_ptr->empty) { @@ -385,17 +392,19 @@ void diversity_manager_t::generate_quick_feasible_solution() ls.generate_fast_solution(solution, sol_timer); if (solution.get_feasible()) { population.run_solution_callbacks(solution); + initial_sol_vector.emplace_back(std::move(solution)); problem_ptr->handle_ptr->sync_stream(); - solution_t searched_sol(solution); + solution_t searched_sol(initial_sol_vector.back()); ls_config_t ls_config; run_local_search(searched_sol, population.weights, sol_timer, ls_config); population.run_solution_callbacks(searched_sol); - auto& feas_sol = searched_sol.get_feasible() ? searched_sol : solution; + initial_sol_vector.emplace_back(std::move(searched_sol)); + auto& feas_sol = initial_sol_vector.back().get_feasible() + ? initial_sol_vector.back() + : initial_sol_vector[initial_sol_vector.size() - 2]; CUOPT_LOG_INFO("Generated fast solution in %f seconds with objective %f", timer.elapsed_time(), feas_sol.get_user_objective()); - population.add_solution(std::move(searched_sol)); - population.add_solution(std::move(solution)); } problem_ptr->handle_ptr->sync_stream(); } diff --git a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh index e9444c423a..f1af5744d6 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh @@ -204,6 +204,9 @@ struct fj_move_candidate_t { template struct fj_cpu_climber_t; +template +class probing_cache_t; + template class fj_t { public: @@ -221,6 +224,7 @@ class fj_t { const std::vector& right_weights, f_t objective_weight, std::atomic& preemption_flag, + const probing_cache_t* probing_cache, fj_settings_t settings = fj_settings_t{}, bool randomize_params = false); i_t alloc_max_climbers(i_t desired_climbers); diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 5964fb346e..24aeeb2b1e 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -15,7 +15,7 @@ #include "fj_cpu.cuh" #include "fj_cpu_worker.cuh" -#include +#include #include @@ -746,14 +746,33 @@ static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& return score; } +// Translate a variable of this problem into the pre-trivial-presolve id the probing cache is keyed +// by. An empty map means presolve removed nothing, so the ids coincide. +template +static inline i_t two_opt_probed_id(const fj_cpu_climber_t& fj_cpu, i_t var_idx) +{ + if (fj_cpu.h_original_ids.size() == 0) { return var_idx; } + cuopt_assert(var_idx < (i_t)fj_cpu.h_original_ids.size(), "variable has no original id"); + return fj_cpu.h_original_ids[var_idx]; +} + +// Reverse of two_opt_probed_id. Returns -1 for a variable presolve has since removed. +template +static inline i_t two_opt_problem_id(const fj_cpu_climber_t& fj_cpu, i_t probed_id) +{ + if (fj_cpu.h_reverse_original_ids.size() == 0) { return probed_id; } + if (probed_id < 0 || probed_id >= (i_t)fj_cpu.h_reverse_original_ids.size()) { return -1; } + return fj_cpu.h_reverse_original_ids[probed_id]; +} + /** * @brief Fill fj_cpu.two_opt_partners with candidates to flip together with `first`. * - * Preferred source is the conflict graph: the literals conflicting with the one `first` is about to - * make true are exactly the ones that must become false, which gives both the partner and the value - * to move it to. A conflict with a negated literal asks for a same-direction pair, which is how - * covering rows contribute. When the graph has nothing for `first`, the variables sharing a row - * with it are the fallback, and there the only defensible direction is the swap. + * Preferred source is the probing cache: it recorded, for each probed variable and value, the + * bounds propagation implies on every other variable. An implied bound pinning a binary to a value + * names both the partner and the value it has to take once `first` moves, so a pair moving in the + * same direction is reached as naturally as a swap. When probing has nothing for `first`, the + * variables sharing a row with it are the fallback, and the only defensible direction is the swap. */ template static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, @@ -779,20 +798,36 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, partners.emplace_back(var_idx, delta); }; - if (fj_cpu.clique_table != nullptr) { - cuopt_assert(fj_cpu.clique_table->n_variables == n_variables, - "conflict graph belongs to another problem"); - const f_t new_val = fj_cpu.h_assignment[first].get() + first_delta; - // Literal `v` stands for "variable v is 1", literal `v + n_variables` for "variable v is 0" - const i_t literal = new_val > 0.5 ? first : first + n_variables; - fj_cpu.clique_table->for_each_conflict_partner(literal, [&](i_t partner_literal) { - // The conflicting literal has to become false - add_partner(partner_literal % n_variables, partner_literal < n_variables ? 0 : 1); - return partners.size() < max_partners; - }); - // A partner is yielded once per shared clique - std::sort(partners.begin(), partners.end()); - partners.erase(std::unique(partners.begin(), partners.end()), partners.end()); + if (fj_cpu.probing_cache != nullptr) { + const auto& cache = fj_cpu.probing_cache->probing_cache; + const auto cached_probe = cache.find(two_opt_probed_id(fj_cpu, first)); + if (cached_probe != cache.end()) { + const f_t new_val = fj_cpu.h_assignment[first].get() + first_delta; + // Pick the probed interval that covers the value `first` is moving to; the two entries per + // variable are the two values or intervals it was probed at. + i_t hit_interval = -1; + i_t unused_hit = -1; + for (i_t interval = 0; interval < 2; ++interval) { + const auto& entry = cached_probe->second[interval]; + if (entry.var_to_cached_bound_map.empty()) { continue; } + entry.val_interval.fill_cache_hits(interval, new_val, new_val, hit_interval, unused_hit); + } + if (hit_interval != -1) { + const auto& implications = cached_probe->second[hit_interval].var_to_cached_bound_map; + for (const auto& [probed_id, implied] : implications) { + if (partners.size() >= max_partners) break; + const i_t var_idx = two_opt_problem_id(fj_cpu, probed_id); + // -1 means presolve removed the variable after the probe recorded it + if (var_idx < 0) { continue; } + cuopt_assert(var_idx < n_variables, "implied variable out of range"); + if (!fj_cpu.h_is_binary_variable[var_idx]) { continue; } + // Only an implication that pins the partner names a value to move it to; a bound that + // still admits both says nothing about what the partner should do. + if (!fj_cpu.view.pb.integer_equal(implied.lb, implied.ub)) { continue; } + add_partner(var_idx, round(implied.lb)); + } + } + } } if (!partners.empty()) return; @@ -840,7 +875,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) two_opt_move_t best; const bool partner_source_exists = - fj_cpu.clique_table != nullptr || + (fj_cpu.probing_cache != nullptr && !fj_cpu.probing_cache->probing_cache.empty()) || fj_cpu.h_related_variables_offsets.size() == static_cast(fj_cpu.view.pb.n_variables) + 1; if (fj_cpu.n_binary_vars == 0 || !partner_source_exists) return best; @@ -850,6 +885,9 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) first_vars.clear(); if (!fj_cpu.violated_constraints.empty()) { + cuopt_assert(fj_cpu.h_binrow_offsets.size() == + static_cast(fj_cpu.view.pb.n_constraints) + 1, + "binary row table missing"); auto& target_cstrs = fj_cpu.two_opt_target_cstrs; target_cstrs.clear(); std::sample(fj_cpu.violated_constraints.begin(), @@ -858,10 +896,10 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) max_target_rows, rng); for (i_t cstr_idx : target_cstrs) { - const auto [offset_begin, offset_end] = range_for_constraint(fj_cpu, cstr_idx); - for (i_t i = offset_begin; i < offset_end && first_vars.size() < max_first_vars; ++i) { - const i_t var_idx = fj_cpu.h_variables[i]; - if (fj_cpu.h_is_binary_variable[var_idx]) first_vars.push_back(var_idx); + const i_t bin_begin = fj_cpu.h_binrow_offsets[cstr_idx]; + const i_t bin_end = fj_cpu.h_binrow_offsets[cstr_idx + 1]; + for (i_t i = bin_begin; i < bin_end && first_vars.size() < max_first_vars; ++i) { + first_vars.push_back(fj_cpu.h_binrow_vars[i].get()); } } } else { @@ -1510,7 +1548,8 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, solution_t& solution, const std::vector& left_weights, const std::vector& right_weights, - f_t objective_weight) + f_t objective_weight, + const probing_cache_t* probing_cache) { auto& problem = *solution.problem_ptr; auto handle_ptr = solution.handle_ptr; @@ -1543,13 +1582,9 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, cuopt::host_copy(problem.related_variables, handle_ptr->get_stream()); fj_cpu.h_related_variables_offsets = cuopt::host_copy(problem.related_variables_offsets, handle_ptr->get_stream()); - fj_cpu.clique_table = problem.clique_table; - // A problem built by fixing variables drops the graph, so a size mismatch means it outlived the - // index space it describes. Refuse to use it rather than index out of range. - if (fj_cpu.clique_table != nullptr && fj_cpu.clique_table->n_variables != problem.n_variables) { - cuopt_assert(false, "conflict graph does not match the problem"); - fj_cpu.clique_table.reset(); - } + fj_cpu.probing_cache = probing_cache; + fj_cpu.h_original_ids = problem.original_ids; + fj_cpu.h_reverse_original_ids = problem.reverse_original_ids; fj_cpu.h_cstr_left_weights = left_weights; fj_cpu.h_cstr_right_weights = right_weights; @@ -1686,6 +1721,18 @@ void finalize_fj_cpu_host_initialization( } } + fj_cpu.h_binrow_offsets.resize(n_constraints + 1); + fj_cpu.h_binrow_vars.clear(); + for (i_t cstr_idx = 0; cstr_idx < n_constraints; ++cstr_idx) { + fj_cpu.h_binrow_offsets[cstr_idx] = fj_cpu.h_binrow_vars.size(); + auto [offset_begin, offset_end] = range_for_constraint(fj_cpu, cstr_idx); + for (i_t i = offset_begin; i < offset_end; ++i) { + const i_t var_idx = fj_cpu.h_variables[i]; + if (fj_cpu.h_is_binary_variable[var_idx]) { fj_cpu.h_binrow_vars.push_back(var_idx); } + } + } + fj_cpu.h_binrow_offsets[n_constraints] = fj_cpu.h_binrow_vars.size(); + fj_cpu.flip_move_computed.resize(n_variables, false); fj_cpu.var_bitmap.resize(n_variables, false); fj_cpu.iter_mtm_vars.reserve(n_variables); @@ -1862,6 +1909,7 @@ std::unique_ptr> fj_t::create_cpu_climber( const std::vector& right_weights, f_t objective_weight, std::atomic& preemption_flag, + const probing_cache_t* probing_cache, fj_settings_t settings, bool randomize_params) { @@ -1870,7 +1918,7 @@ std::unique_ptr> fj_t::create_cpu_climber( auto fj_cpu = std::make_unique>(preemption_flag); // Initialize fj_cpu with all the data - init_fj_cpu(*fj_cpu, solution, left_weights, right_weights, objective_weight); + init_fj_cpu(*fj_cpu, solution, left_weights, right_weights, objective_weight, probing_cache); fj_cpu->settings = settings; if (randomize_params) { auto rng = std::mt19937(cuopt::seed_generator::get_seed()); @@ -2074,7 +2122,9 @@ std::unique_ptr> init_fj_cpu_standalone( auto fj_cpu = std::make_unique>(preemption_flag); std::vector default_weights(problem.n_constraints, 1.0); - init_fj_cpu(*fj_cpu, solution, default_weights, default_weights, 0.0); + // Early CPUFJ runs while presolve is still probing, so there are no implications to hand it + const probing_cache_t* no_implications = nullptr; + init_fj_cpu(*fj_cpu, solution, default_weights, default_weights, 0.0, no_implications); fj_cpu->settings = settings; fj_cpu->settings.seed = cuopt::seed_generator::get_seed(); diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index 58d10b5c0d..609b08618f 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -20,6 +20,9 @@ namespace cuopt::mathematical_optimization::mip { +template +class probing_cache_t; + // NOTE: this seems an easy pick for reflection/xmacros once this is available (C++26?) // Maintaining a single source of truth for all members would be nice template @@ -46,6 +49,10 @@ struct fj_cpu_climber_t { ADD_INSTRUMENTED(h_binary_indices), ADD_INSTRUMENTED(h_related_variables), ADD_INSTRUMENTED(h_related_variables_offsets), + ADD_INSTRUMENTED(h_binrow_offsets), + ADD_INSTRUMENTED(h_binrow_vars), + ADD_INSTRUMENTED(h_original_ids), + ADD_INSTRUMENTED(h_reverse_original_ids), ADD_INSTRUMENTED(h_tabu_nodec_until), ADD_INSTRUMENTED(h_tabu_noinc_until), ADD_INSTRUMENTED(h_tabu_lastdec), @@ -87,10 +94,17 @@ struct fj_cpu_climber_t { ins_vector h_binary_indices; ins_vector h_related_variables; ins_vector h_related_variables_offsets; - // Conflict graph of the problem this climber was built from, null when presolve built none. Held - // by shared_ptr and snapshotted alongside the host copies above because it is indexed by that - // problem's variable ids. - std::shared_ptr> clique_table; + + // precompute the binary variables per row for bin 2opt + ins_vector h_binrow_offsets; + ins_vector h_binrow_vars; + // Implications recorded by probing: for a probed variable set to a value, the bounds propagation + // implies on everything else. The 2-opt reads it to learn which value a partner has to take. Null + // when probing was disabled or had not run yet when this climber was created. + const probing_cache_t* probing_cache{nullptr}; + // Probing cache keys are pre-trivial-presolve variable ids; these translate to and from them + ins_vector h_original_ids; + ins_vector h_reverse_original_ids; ins_vector h_tabu_nodec_until; ins_vector h_tabu_noinc_until; diff --git a/cpp/src/mip_heuristics/local_search/local_search.cu b/cpp/src/mip_heuristics/local_search/local_search.cu index 92be0822d3..bdeabee0a6 100644 --- a/cpp/src/mip_heuristics/local_search/local_search.cu +++ b/cpp/src/mip_heuristics/local_search/local_search.cu @@ -71,6 +71,7 @@ void local_search_t::start_cpufj_scratch_threads(population_t 0); @@ -117,8 +118,12 @@ void local_search_t::start_cpufj_lptopt_scratch_threads( solution_lp.copy_new_assignment( host_copy(lp_optimal_solution, context.problem_ptr->handle_ptr->get_stream())); solution_lp.round_random_nearest(500); - scratch_cpu_fj_on_lp_opt = fj.create_cpu_climber( - solution_lp, default_weights, default_weights, 0., context.preempt_heuristic_solver_); + scratch_cpu_fj_on_lp_opt = fj.create_cpu_climber(solution_lp, + default_weights, + default_weights, + 0., + context.preempt_heuristic_solver_, + &constraint_prop.bounds_update.probing_cache); scratch_cpu_fj_on_lp_opt->log_prefix = "******* scratch on LP optimal: "; scratch_cpu_fj_on_lp_opt->improvement_callback = [this, &population](f_t obj, const std::vector& h_vec, double /*work_units*/) { @@ -186,6 +191,7 @@ void local_search_t::start_cpufj_deterministic(mip::branch_and_bound_t default_weights, 0., context.preempt_heuristic_solver_, + &constraint_prop.bounds_update.probing_cache, fj_settings_t{}, /*randomize=*/true); @@ -261,6 +267,7 @@ bool local_search_t::do_fj_solve(solution_t& solution, h_weights, h_objective_weight, context.preempt_heuristic_solver_, + &constraint_prop.bounds_update.probing_cache, fj_settings_t{}, true); } diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu index fb02c5bbc0..a8e6997572 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu @@ -314,10 +314,34 @@ template std::unordered_set clique_table_t::get_adj_set_of_var(i_t var_idx) const { std::unordered_set adj_set; - for_each_conflict_partner(var_idx, [&](i_t partner) { - adj_set.insert(partner); - return true; - }); + + // First-clique edges: every member of each first-clique containing var_idx. + for (i_t clique_idx : var_clique_first.slice(var_idx)) { + const auto& c = first[clique_idx]; + adj_set.insert(c.begin(), c.end()); + } + + // Addtl-clique edges. + for (i_t addtl_idx : var_clique_addtl.slice(var_idx)) { + const auto& a = addtl_cliques[addtl_idx]; + if (a.vertex_idx == var_idx) { + // var_idx is the extension vertex; new neighbors are the base suffix. + const auto& base = first[a.clique_idx]; + adj_set.insert(base.begin() + a.start_pos_on_clique, base.end()); + } else { + // var_idx is a base member; only new edge is to the extension vertex. + adj_set.insert(a.vertex_idx); + } + } + + for (i_t adj_var : small_clique_adj.slice(var_idx)) { + adj_set.insert(adj_var); + } + + // Add the complement of var_idx to the adjacency set + i_t complement_idx = (var_idx >= n_variables) ? (var_idx - n_variables) : (var_idx + n_variables); + adj_set.insert(complement_idx); + adj_set.erase(var_idx); return adj_set; } diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh index d58c99e304..fe6db31c26 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh @@ -9,7 +9,6 @@ #include #include -#include #include #include @@ -175,9 +174,6 @@ struct clique_table_t { clique_table_t& operator=(clique_table_t&&) = delete; std::unordered_set get_adj_set_of_var(i_t var_idx) const; - // Yields the conflict partners of literal, also as literals, until fn returns false - template - void for_each_conflict_partner(i_t literal, callable_t&& fn) const; i_t get_degree_of_var(i_t var_idx); bool check_adjacency(i_t var_idx1, i_t var_idx2) const; bool empty() const @@ -205,43 +201,6 @@ struct clique_table_t { typename mip_solver_settings_t::tolerances_t tolerances; }; -template -template -void clique_table_t::for_each_conflict_partner(i_t literal, callable_t&& fn) const -{ - cuopt_assert(literal >= 0 && literal < 2 * n_variables, "literal out of range"); - const i_t complement = - (literal >= n_variables) ? (literal - n_variables) : (literal + n_variables); - if (!fn(complement)) { return; } - - for (i_t clique_idx : var_clique_first.slice(literal)) { - for (i_t member : first[clique_idx]) { - if (member != literal && !fn(member)) { return; } - } - } - - for (i_t addtl_idx : var_clique_addtl.slice(literal)) { - const auto& addtl = addtl_cliques[addtl_idx]; - if (addtl.vertex_idx == literal) { - // literal is the extension vertex; its new neighbors are the base suffix - const auto& base = first[addtl.clique_idx]; - cuopt_assert(addtl.start_pos_on_clique >= 0 && - static_cast(addtl.start_pos_on_clique) <= base.size(), - "extension start position out of range"); - for (auto it = base.begin() + addtl.start_pos_on_clique; it != base.end(); ++it) { - if (*it != literal && !fn(*it)) { return; } - } - } else if (!fn(addtl.vertex_idx)) { - // literal is a base member; the only new edge is to the extension vertex - return; - } - } - - for (i_t member : small_clique_adj.slice(literal)) { - if (member != literal && !fn(member)) { return; } - } -} - // Builds the conflict-graph clique table for `problem`. The base cliques are // published to `clique_table_out` before the (optional, signal-gated) extension // phase begins, so cut generation can pick up the table while extension keeps diff --git a/cpp/src/mip_heuristics/presolve/probing_cache.cuh b/cpp/src/mip_heuristics/presolve/probing_cache.cuh index 24d9a9cfc1..079c99edbd 100644 --- a/cpp/src/mip_heuristics/presolve/probing_cache.cuh +++ b/cpp/src/mip_heuristics/presolve/probing_cache.cuh @@ -44,7 +44,7 @@ struct val_interval_t { f_t first_probe, f_t second_probe, i_t& hit_interval_for_first_probe, - i_t& hit_interval_for_second_probe) + i_t& hit_interval_for_second_probe) const { if (interval_type == interval_type_t::EQUALS) { if (val == first_probe) { hit_interval_for_first_probe = interval; } diff --git a/cpp/src/mip_heuristics/problem/problem.cu b/cpp/src/mip_heuristics/problem/problem.cu index b8214831de..ccba2d5f2b 100644 --- a/cpp/src/mip_heuristics/problem/problem.cu +++ b/cpp/src/mip_heuristics/problem/problem.cu @@ -1862,9 +1862,6 @@ void problem_t::remove_given_variables(problem_t& original_p thrust::fill(handle_ptr->get_thrust_policy(), offsets.begin(), offsets.end(), 0); cuopt_assert(assignment.size() == n_variables, "Variable size mismatch"); cuopt_assert(variable_map.size() < n_variables, "Too many variables to fix"); - - // the clique table could be propagated after fixings. but to err on the safe side, reset it here - clique_table.reset(); rmm::device_uvector tmp_assignment(assignment, handle_ptr->get_stream()); // first remove the assignment and variable related vectors From 39f8d84d4dfbaf1b649b99e16429c1c24e412818 Mon Sep 17 00:00:00 2001 From: yboucher Date: Tue, 18 Aug 2026 04:30:33 -0700 Subject: [PATCH 4/9] cleanup --- .../feasibility_jump/feasibility_jump.cuh | 2 - .../mip_heuristics/feasibility_jump/fj_cpu.cu | 45 ++++++++----------- .../feasibility_jump/fj_cpu.cuh | 5 +-- 3 files changed, 21 insertions(+), 31 deletions(-) diff --git a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh index f1af5744d6..ac1da031e3 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cuh @@ -75,8 +75,6 @@ struct fj_hyper_parameters_t { double small_move_tabu_threshold = 1e-6; int small_move_tabu_tenure = 4; - // binary 2-opt caps: violated rows sampled per local minimum, candidate variables taken from each - // sampled row, and total candidate pairs scored int two_opt_max_rows = 4; int two_opt_max_row_vars = 256; int two_opt_max_pairs = 256; diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 24aeeb2b1e..7c3be243e1 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -649,8 +649,6 @@ static inline std::pair compute_score(fj_cpu_climber_t static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& fj_cpu, @@ -877,16 +870,16 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) const bool partner_source_exists = (fj_cpu.probing_cache != nullptr && !fj_cpu.probing_cache->probing_cache.empty()) || fj_cpu.h_related_variables_offsets.size() == - static_cast(fj_cpu.view.pb.n_variables) + 1; + fj_cpu.view.pb.n_variables + 1; + if (fj_cpu.n_binary_vars == 0 || !partner_source_exists) return best; - std::mt19937 rng(fj_cpu.settings.seed + fj_cpu.iterations); auto& first_vars = fj_cpu.two_opt_first_vars; first_vars.clear(); if (!fj_cpu.violated_constraints.empty()) { cuopt_assert(fj_cpu.h_binrow_offsets.size() == - static_cast(fj_cpu.view.pb.n_constraints) + 1, + fj_cpu.view.pb.n_constraints + 1, "binary row table missing"); auto& target_cstrs = fj_cpu.two_opt_target_cstrs; target_cstrs.clear(); @@ -894,7 +887,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) fj_cpu.violated_constraints.end(), std::back_inserter(target_cstrs), max_target_rows, - rng); + fj_cpu.rng); for (i_t cstr_idx : target_cstrs) { const i_t bin_begin = fj_cpu.h_binrow_offsets[cstr_idx]; const i_t bin_end = fj_cpu.h_binrow_offsets[cstr_idx + 1]; @@ -907,7 +900,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) fj_cpu.h_objective_vars.underlying().end(), std::back_inserter(first_vars), max_obj_starts, - rng); + fj_cpu.rng); // Nothing is violated, so a pair can only help by improving the objective: keep the flips that // move it down and let the pair scoring pay for the feasibility damage. first_vars.erase(std::remove_if(first_vars.begin(), @@ -920,10 +913,11 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) }), first_vars.end()); } - std::shuffle(first_vars.begin(), first_vars.end(), rng); + std::shuffle(first_vars.begin(), first_vars.end(), fj_cpu.rng); const i_t nnz_at_entry = fj_cpu.nnz_processed_window; size_t pairs_scored = 0; + // find a (first, second) pair for the 2opt for (i_t first : first_vars) { if (pairs_scored >= max_pairs) break; if (fj_cpu.nnz_processed_window - nnz_at_entry > fj_cpu.nnz_samples) break; @@ -932,21 +926,20 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) const f_t first_delta = round(1 - 2 * first_val); if (tabu_check(fj_cpu, first, first_delta, true)) continue; if (!check_variable_within_bounds(fj_cpu, first, first_val + first_delta)) continue; - const i_t first_inc = fj_cpu.h_tabu_lastinc[first]; - const i_t first_dec = fj_cpu.h_tabu_lastdec[first]; - const i_t first_touch = std::max(first_inc, first_dec); + const i_t first_touch = std::max(fj_cpu.h_tabu_lastinc[first], fj_cpu.h_tabu_lastdec[first]); + // look for potential other binary vars to flip alongside the first var two_opt_collect_partners(fj_cpu, first, first_delta, max_partners_per_var); for (const auto& [second, second_delta] : fj_cpu.two_opt_partners) { - const i_t second_inc = fj_cpu.h_tabu_lastinc[second]; - const i_t second_dec = fj_cpu.h_tabu_lastdec[second]; + const i_t second_touch = std::max(fj_cpu.h_tabu_lastinc[second], fj_cpu.h_tabu_lastdec[second]); two_opt_move_t cand; cand.first = {first, first_delta}; cand.second = {second, second_delta}; cand.score = two_opt_compute_pair_score(fj_cpu, first, first_delta, second, second_delta); - cand.age = std::max(first_touch, std::max(second_inc, second_dec)); + cand.age = std::max(first_touch, second_touch); if (two_opt_cand_better(cand, best)) { best = cand; } ++pairs_scored; + if (pairs_scored >= max_pairs) return best; if (fj_cpu.nnz_processed_window - nnz_at_entry > fj_cpu.nnz_samples) return best; } @@ -1335,7 +1328,7 @@ static thrust::tuple find_mtm_move_viol( fj_cpu.violated_constraints.end(), std::back_inserter(sampled_cstrs), sample_size, - std::mt19937(fj_cpu.settings.seed + fj_cpu.iterations)); + fj_cpu.rng); return find_mtm_move(fj_cpu, sampled_cstrs, localmin); } @@ -1353,7 +1346,7 @@ static thrust::tuple find_mtm_move_sat( fj_cpu.satisfied_constraints.end(), std::back_inserter(sampled_cstrs), sample_size, - std::mt19937(fj_cpu.settings.seed + fj_cpu.iterations)); + fj_cpu.rng); return find_mtm_move(fj_cpu, sampled_cstrs); } @@ -1521,7 +1514,7 @@ static void perturb(fj_cpu_climber_t& fj_cpu) fj_cpu.h_objective_vars.end(), std::back_inserter(sampled_vars), 2, - std::mt19937(fj_cpu.settings.seed + fj_cpu.iterations)); + fj_cpu.rng); raft::random::PCGenerator rng(fj_cpu.settings.seed + fj_cpu.iterations, 0, 0); for (auto var_idx : sampled_vars) { @@ -1721,6 +1714,7 @@ void finalize_fj_cpu_host_initialization( } } + // precompute the binvars-pre-row tables for 2opt fj_cpu.h_binrow_offsets.resize(n_constraints + 1); fj_cpu.h_binrow_vars.clear(); for (i_t cstr_idx = 0; cstr_idx < n_constraints; ++cstr_idx) { @@ -1939,6 +1933,8 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit, double w auto time_limit = std::chrono::milliseconds(static_cast(std::floor(in_time_limit * 1000.0))); auto loop_time_start = std::chrono::high_resolution_clock::now(); + fj_cpu->rng.seed(fj_cpu->settings.seed); + // Initialize feature tracking fj_cpu->last_feature_log_time = loop_start; fj_cpu->prev_best_objective = fj_cpu->h_best_objective; @@ -2017,13 +2013,10 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, f_t in_time_limit, double w for (size_t i = 0; i < fj_cpu->cached_mtm_moves.size(); i++) fj_cpu->cached_mtm_moves[i].first = 0; } - // A simultaneous flip of two binaries escapes minima that no single flip can. Not attempted - // right after a perturbation, whose whole point is to leave the current region. + two_opt_move_t two_opt_move; if (!should_perturb) two_opt_move = find_two_opt_move(*fj_cpu); if (two_opt_move.score > fj_staged_score_t::zero()) { - // Applied as two moves; they were scored jointly, so the intermediate state after the first - // one may well be worse than the local minimum we came from apply_move(*fj_cpu, two_opt_move.first.var_idx, two_opt_move.first.value, true); apply_move(*fj_cpu, two_opt_move.second.var_idx, two_opt_move.second.value, true); fj_cpu->n_mtm_viol_moves_window += 2; diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index 609b08618f..411b4083f7 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -76,6 +77,7 @@ struct fj_cpu_climber_t { problem_t* pb_ptr; fj_settings_t settings; + std::mt19937 rng; typename fj_t::climber_data_t::view_t view; // Host copies of device data as struct members ins_vector h_reverse_coefficients; @@ -98,9 +100,6 @@ struct fj_cpu_climber_t { // precompute the binary variables per row for bin 2opt ins_vector h_binrow_offsets; ins_vector h_binrow_vars; - // Implications recorded by probing: for a probed variable set to a value, the bounds propagation - // implies on everything else. The 2-opt reads it to learn which value a partner has to take. Null - // when probing was disabled or had not run yet when this climber was created. const probing_cache_t* probing_cache{nullptr}; // Probing cache keys are pre-trivial-presolve variable ids; these translate to and from them ins_vector h_original_ids; From 14cc888c79ec1ed0728960ec51ce99ada24878ee Mon Sep 17 00:00:00 2001 From: yboucher Date: Tue, 18 Aug 2026 04:48:15 -0700 Subject: [PATCH 5/9] more cleanup --- .../mip_heuristics/feasibility_jump/fj_cpu.cu | 74 +++++++------------ 1 file changed, 25 insertions(+), 49 deletions(-) diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 7c3be243e1..4dfc9993a2 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -654,17 +654,15 @@ struct two_opt_move_t { fj_move_t second{-1, 0}; fj_staged_score_t score{fj_staged_score_t::invalid()}; int age{std::numeric_limits::max()}; -}; -// Deterministic total order over 2-opt candidates: higher score, then older last touch, then -// smaller variable indices. -static bool two_opt_cand_better(const two_opt_move_t& a, const two_opt_move_t& b) -{ - if (a.score != b.score) return a.score > b.score; - if (a.age != b.age) return a.age < b.age; - if (a.first.var_idx != b.first.var_idx) return a.first.var_idx < b.first.var_idx; - return a.second.var_idx < b.second.var_idx; -} + bool operator>(const two_opt_move_t& other) const + { + if (score != other.score) return score > other.score; + if (age != other.age) return age < other.age; + if (first.var_idx != other.first.var_idx) return first.var_idx < other.first.var_idx; + return second.var_idx < other.second.var_idx; + } +}; /** * @brief Score the combined effect of flipping two binaries at once. @@ -716,7 +714,6 @@ static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& bonus_robust_sum += cstr_bonus_robust; } - // Same staged score convention as the single variable path const f_t obj_diff = fj_cpu.h_obj_coeffs[first] * first_delta + fj_cpu.h_obj_coeffs[second] * second_delta; f_t base_obj = 0; @@ -739,25 +736,6 @@ static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& return score; } -// Translate a variable of this problem into the pre-trivial-presolve id the probing cache is keyed -// by. An empty map means presolve removed nothing, so the ids coincide. -template -static inline i_t two_opt_probed_id(const fj_cpu_climber_t& fj_cpu, i_t var_idx) -{ - if (fj_cpu.h_original_ids.size() == 0) { return var_idx; } - cuopt_assert(var_idx < (i_t)fj_cpu.h_original_ids.size(), "variable has no original id"); - return fj_cpu.h_original_ids[var_idx]; -} - -// Reverse of two_opt_probed_id. Returns -1 for a variable presolve has since removed. -template -static inline i_t two_opt_problem_id(const fj_cpu_climber_t& fj_cpu, i_t probed_id) -{ - if (fj_cpu.h_reverse_original_ids.size() == 0) { return probed_id; } - if (probed_id < 0 || probed_id >= (i_t)fj_cpu.h_reverse_original_ids.size()) { return -1; } - return fj_cpu.h_reverse_original_ids[probed_id]; -} - /** * @brief Fill fj_cpu.two_opt_partners with candidates to flip together with `first`. * @@ -777,6 +755,12 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, const i_t n_variables = fj_cpu.view.pb.n_variables; partners.clear(); cuopt_assert(fj_cpu.h_is_binary_variable[first], "2-opt is only defined for binaries"); + cuopt_assert(fj_cpu.probing_cache == nullptr || + fj_cpu.h_original_ids.size() == (size_t)n_variables, + "original id map does not cover every variable"); + cuopt_assert(fj_cpu.probing_cache == nullptr || + fj_cpu.h_reverse_original_ids.size() >= fj_cpu.h_original_ids.size(), + "reverse original id map smaller than the problem"); auto add_partner = [&](i_t var_idx, f_t target) { if (var_idx == first) return; @@ -792,14 +776,12 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, }; if (fj_cpu.probing_cache != nullptr) { - const auto& cache = fj_cpu.probing_cache->probing_cache; - const auto cached_probe = cache.find(two_opt_probed_id(fj_cpu, first)); + const auto& cache = fj_cpu.probing_cache->probing_cache; + const auto cached_probe = cache.find(fj_cpu.h_original_ids[first]); if (cached_probe != cache.end()) { const f_t new_val = fj_cpu.h_assignment[first].get() + first_delta; - // Pick the probed interval that covers the value `first` is moving to; the two entries per - // variable are the two values or intervals it was probed at. - i_t hit_interval = -1; - i_t unused_hit = -1; + i_t hit_interval = -1; + i_t unused_hit = -1; for (i_t interval = 0; interval < 2; ++interval) { const auto& entry = cached_probe->second[interval]; if (entry.var_to_cached_bound_map.empty()) { continue; } @@ -809,13 +791,11 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, const auto& implications = cached_probe->second[hit_interval].var_to_cached_bound_map; for (const auto& [probed_id, implied] : implications) { if (partners.size() >= max_partners) break; - const i_t var_idx = two_opt_problem_id(fj_cpu, probed_id); + const i_t var_idx = fj_cpu.h_reverse_original_ids[probed_id]; // -1 means presolve removed the variable after the probe recorded it if (var_idx < 0) { continue; } cuopt_assert(var_idx < n_variables, "implied variable out of range"); if (!fj_cpu.h_is_binary_variable[var_idx]) { continue; } - // Only an implication that pins the partner names a value to move it to; a bound that - // still admits both says nothing about what the partner should do. if (!fj_cpu.view.pb.integer_equal(implied.lb, implied.ub)) { continue; } add_partner(var_idx, round(implied.lb)); } @@ -826,10 +806,7 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, const auto& related = fj_cpu.h_related_variables; const auto& related_offsets = fj_cpu.h_related_variables_offsets; - if (related_offsets.size() != static_cast(n_variables) + 1) return; - // Row sharing carries no polarity, so take the value `first` is vacating: only a partner holding - // the opposite value then moves in the compensating direction, and one holding the same value is - // filtered out as a no-op. Same opposite-value rule as the GPU candidate build. + if (related_offsets.size() != (size_t)n_variables + 1) return; const f_t swap_target = fj_cpu.h_assignment[first].get(); const i_t related_begin = related_offsets[first]; const i_t related_end = related_offsets[first + 1]; @@ -856,8 +833,7 @@ template static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) { CPUFJ_NVTX_RANGE("CPUFJ::find_two_opt_move"); - constexpr size_t max_obj_starts = 64; - // The GPU candidate table is indexed by pair, so it has no per-first-variable level to cap + constexpr size_t max_obj_starts = 64; constexpr size_t max_partners_per_var = 16; const auto& params = fj_cpu.settings.parameters; @@ -869,7 +845,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) const bool partner_source_exists = (fj_cpu.probing_cache != nullptr && !fj_cpu.probing_cache->probing_cache.empty()) || - fj_cpu.h_related_variables_offsets.size() == + (int64_t)fj_cpu.h_related_variables_offsets.size() == fj_cpu.view.pb.n_variables + 1; if (fj_cpu.n_binary_vars == 0 || !partner_source_exists) return best; @@ -877,6 +853,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) auto& first_vars = fj_cpu.two_opt_first_vars; first_vars.clear(); + // target binvars in violated constraints for flips if (!fj_cpu.violated_constraints.empty()) { cuopt_assert(fj_cpu.h_binrow_offsets.size() == fj_cpu.view.pb.n_constraints + 1, @@ -896,13 +873,12 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) } } } else { + // target objective-bearing binary vars in satisfied constraints std::sample(fj_cpu.h_objective_vars.underlying().begin(), fj_cpu.h_objective_vars.underlying().end(), std::back_inserter(first_vars), max_obj_starts, fj_cpu.rng); - // Nothing is violated, so a pair can only help by improving the objective: keep the flips that - // move it down and let the pair scoring pay for the feasibility damage. first_vars.erase(std::remove_if(first_vars.begin(), first_vars.end(), [&](i_t var_idx) { @@ -937,7 +913,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) cand.second = {second, second_delta}; cand.score = two_opt_compute_pair_score(fj_cpu, first, first_delta, second, second_delta); cand.age = std::max(first_touch, second_touch); - if (two_opt_cand_better(cand, best)) { best = cand; } + if (cand > best) { best = cand; } ++pairs_scored; if (pairs_scored >= max_pairs) return best; From c296e6d98c2b6286d472ed16a51fa4607fa485ce Mon Sep 17 00:00:00 2001 From: yboucher Date: Tue, 18 Aug 2026 05:10:00 -0700 Subject: [PATCH 6/9] would you look at that, even more cleannup! --- .../feasibility_jump_impl_common.cuh | 2 +- .../mip_heuristics/feasibility_jump/fj_cpu.cu | 67 +++++++++---------- 2 files changed, 32 insertions(+), 37 deletions(-) diff --git a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump_impl_common.cuh b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump_impl_common.cuh index 98267f117c..046e138c5b 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump_impl_common.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/feasibility_jump_impl_common.cuh @@ -171,7 +171,7 @@ HDI std::pair feas_score_constraint( base_feas += (i_t)(cstr_weight * fj.settings->parameters.excess_improvement_weight); } // simple worsening - else if (!old_sat && !new_sat && old_lhs <= new_lhs) { + else if (!old_sat && !new_sat && old_lhs < new_lhs) { cuopt_assert(old_viol && new_viol, ""); base_feas -= (i_t)(cstr_weight * fj.settings->parameters.excess_improvement_weight); } diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 4dfc9993a2..d0e3609dcb 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -206,7 +206,7 @@ std::pair feas_score_constraint(const typename fj_t::climber base_feas += (i_t)(cstr_weight * fj.settings->parameters.excess_improvement_weight); } // simple worsening - else if (!old_sat && !new_sat && old_lhs <= new_lhs) { + else if (!old_sat && !new_sat && old_lhs < new_lhs) { cuopt_assert(old_viol && new_viol, ""); base_feas -= (i_t)(cstr_weight * fj.settings->parameters.excess_improvement_weight); } @@ -664,9 +664,8 @@ struct two_opt_move_t { } }; -/** - * @brief Score the combined effect of flipping two binaries at once. - */ + +// returns the combined score of a joint 2opt move template static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& fj_cpu, i_t first, @@ -676,7 +675,8 @@ static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& { auto& row_deltas = fj_cpu.two_opt_row_deltas; row_deltas.clear(); - auto collect = [&](i_t var_idx, f_t delta) { + const fj_move_t endpoints[2] = {{first, first_delta}, {second, second_delta}}; + for (const auto& [var_idx, delta] : endpoints) { const auto [offset_begin, offset_end] = reverse_range_for_var(fj_cpu, var_idx); fj_cpu.nnz_processed_window += offset_end - offset_begin; for (i_t i = offset_begin; i < offset_end; ++i) { @@ -684,9 +684,7 @@ static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& const f_t coeff = fj_cpu.h_reverse_coefficients[i]; row_deltas.emplace_back(cstr_idx, coeff * delta); } - }; - collect(first, first_delta); - collect(second, second_delta); + } // Brings the entries of a shared row next to each other std::sort(row_deltas.begin(), row_deltas.end()); @@ -736,6 +734,24 @@ static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& return score; } +template +static void two_opt_add_partner(fj_cpu_climber_t& fj_cpu, + i_t first, + i_t var_idx, + f_t target) +{ + if (var_idx == first) return; + const f_t val = fj_cpu.h_assignment[var_idx].get(); + // A partner between two integers has no opposite value to swap to + if (!fj_cpu.view.pb.is_integer(val)) return; + const f_t delta = target - val; + // Already at the value we would move it to, so there is no compound move to make + if (fabs(delta) < 0.5) return; + if (!check_variable_within_bounds(fj_cpu, var_idx, target)) return; + if (tabu_check(fj_cpu, var_idx, delta, true)) return; + fj_cpu.two_opt_partners.emplace_back(var_idx, delta); +} + /** * @brief Fill fj_cpu.two_opt_partners with candidates to flip together with `first`. * @@ -762,19 +778,6 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, fj_cpu.h_reverse_original_ids.size() >= fj_cpu.h_original_ids.size(), "reverse original id map smaller than the problem"); - auto add_partner = [&](i_t var_idx, f_t target) { - if (var_idx == first) return; - const f_t val = fj_cpu.h_assignment[var_idx].get(); - // A partner between two integers has no opposite value to swap to - if (!fj_cpu.view.pb.is_integer(val)) return; - const f_t delta = target - val; - // Already at the value we would move it to, so there is no compound move to make - if (fabs(delta) < 0.5) return; - if (!check_variable_within_bounds(fj_cpu, var_idx, target)) return; - if (tabu_check(fj_cpu, var_idx, delta, true)) return; - partners.emplace_back(var_idx, delta); - }; - if (fj_cpu.probing_cache != nullptr) { const auto& cache = fj_cpu.probing_cache->probing_cache; const auto cached_probe = cache.find(fj_cpu.h_original_ids[first]); @@ -797,7 +800,7 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, cuopt_assert(var_idx < n_variables, "implied variable out of range"); if (!fj_cpu.h_is_binary_variable[var_idx]) { continue; } if (!fj_cpu.view.pb.integer_equal(implied.lb, implied.ub)) { continue; } - add_partner(var_idx, round(implied.lb)); + two_opt_add_partner(fj_cpu, first, var_idx, round(implied.lb)); } } } @@ -812,23 +815,15 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, const i_t related_end = related_offsets[first + 1]; for (i_t i = related_begin; i < related_end && partners.size() < max_partners; ++i) { const i_t var_idx = related[i]; - if (fj_cpu.h_is_binary_variable[var_idx]) { add_partner(var_idx, swap_target); } + if (fj_cpu.h_is_binary_variable[var_idx]) { + two_opt_add_partner(fj_cpu, first, var_idx, swap_target); + } } } -/** - * @brief Look for an improving simultaneous flip of two binaries (binary 2-opt). - * - * At a local minimum no single flip improves the score, but a pair often does: in set partitioning - * style rows the only way out is to turn one variable off and another on in the same step. The - * neighbourhood is sampled rather than enumerated, which would be quadratic. First variables come - * from a few violated rows, or from objective variables whose flip improves the objective when - * nothing is violated; partners come from two_opt_collect_partners. The search stops at - * two_opt_max_pairs candidates or once it has touched nnz_samples nonzeros, whichever comes first: - * the pair count bounds the neighbourhood, the nonzero count bounds the cost, which here is - * proportional to the degrees of both endpoints. Returns an invalid-scored move when no pair was - * worth applying. - */ +// Look for binary 2opt moves at a local minimum. by definition no 1opt move can improve, but combined moves may +// especially in the case of set partitioning constraints / cliques. Use information from the probing cache +// to find potential good 2opt moves. template static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) { From 5093f2dfe99b4f4e0ccc7be454050c8593f374fc Mon Sep 17 00:00:00 2001 From: yboucher Date: Wed, 19 Aug 2026 06:59:18 -0700 Subject: [PATCH 7/9] style --- .../mip_heuristics/feasibility_jump/fj_cpu.cu | 37 ++++++++----------- .../local_search/local_search.cu | 2 +- 2 files changed, 17 insertions(+), 22 deletions(-) diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index d0e3609dcb..9915073a96 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -664,14 +664,10 @@ struct two_opt_move_t { } }; - // returns the combined score of a joint 2opt move template -static fj_staged_score_t two_opt_compute_pair_score(fj_cpu_climber_t& fj_cpu, - i_t first, - f_t first_delta, - i_t second, - f_t second_delta) +static fj_staged_score_t two_opt_compute_pair_score( + fj_cpu_climber_t& fj_cpu, i_t first, f_t first_delta, i_t second, f_t second_delta) { auto& row_deltas = fj_cpu.two_opt_row_deltas; row_deltas.clear(); @@ -771,9 +767,9 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, const i_t n_variables = fj_cpu.view.pb.n_variables; partners.clear(); cuopt_assert(fj_cpu.h_is_binary_variable[first], "2-opt is only defined for binaries"); - cuopt_assert(fj_cpu.probing_cache == nullptr || - fj_cpu.h_original_ids.size() == (size_t)n_variables, - "original id map does not cover every variable"); + cuopt_assert( + fj_cpu.probing_cache == nullptr || fj_cpu.h_original_ids.size() == (size_t)n_variables, + "original id map does not cover every variable"); cuopt_assert(fj_cpu.probing_cache == nullptr || fj_cpu.h_reverse_original_ids.size() >= fj_cpu.h_original_ids.size(), "reverse original id map smaller than the problem"); @@ -821,9 +817,9 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, } } -// Look for binary 2opt moves at a local minimum. by definition no 1opt move can improve, but combined moves may -// especially in the case of set partitioning constraints / cliques. Use information from the probing cache -// to find potential good 2opt moves. +// Look for binary 2opt moves at a local minimum. by definition no 1opt move can improve, but +// combined moves may especially in the case of set partitioning constraints / cliques. Use +// information from the probing cache to find potential good 2opt moves. template static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) { @@ -840,8 +836,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) const bool partner_source_exists = (fj_cpu.probing_cache != nullptr && !fj_cpu.probing_cache->probing_cache.empty()) || - (int64_t)fj_cpu.h_related_variables_offsets.size() == - fj_cpu.view.pb.n_variables + 1; + (int64_t)fj_cpu.h_related_variables_offsets.size() == fj_cpu.view.pb.n_variables + 1; if (fj_cpu.n_binary_vars == 0 || !partner_source_exists) return best; @@ -850,8 +845,7 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) // target binvars in violated constraints for flips if (!fj_cpu.violated_constraints.empty()) { - cuopt_assert(fj_cpu.h_binrow_offsets.size() == - fj_cpu.view.pb.n_constraints + 1, + cuopt_assert(fj_cpu.h_binrow_offsets.size() == fj_cpu.view.pb.n_constraints + 1, "binary row table missing"); auto& target_cstrs = fj_cpu.two_opt_target_cstrs; target_cstrs.clear(); @@ -902,12 +896,13 @@ static two_opt_move_t find_two_opt_move(fj_cpu_climber_t& fj_cpu) // look for potential other binary vars to flip alongside the first var two_opt_collect_partners(fj_cpu, first, first_delta, max_partners_per_var); for (const auto& [second, second_delta] : fj_cpu.two_opt_partners) { - const i_t second_touch = std::max(fj_cpu.h_tabu_lastinc[second], fj_cpu.h_tabu_lastdec[second]); + const i_t second_touch = + std::max(fj_cpu.h_tabu_lastinc[second], fj_cpu.h_tabu_lastdec[second]); two_opt_move_t cand; cand.first = {first, first_delta}; cand.second = {second, second_delta}; - cand.score = two_opt_compute_pair_score(fj_cpu, first, first_delta, second, second_delta); - cand.age = std::max(first_touch, second_touch); + cand.score = two_opt_compute_pair_score(fj_cpu, first, first_delta, second, second_delta); + cand.age = std::max(first_touch, second_touch); if (cand > best) { best = cand; } ++pairs_scored; @@ -1546,8 +1541,8 @@ static void init_fj_cpu(fj_cpu_climber_t& fj_cpu, cuopt::host_copy(problem.related_variables, handle_ptr->get_stream()); fj_cpu.h_related_variables_offsets = cuopt::host_copy(problem.related_variables_offsets, handle_ptr->get_stream()); - fj_cpu.probing_cache = probing_cache; - fj_cpu.h_original_ids = problem.original_ids; + fj_cpu.probing_cache = probing_cache; + fj_cpu.h_original_ids = problem.original_ids; fj_cpu.h_reverse_original_ids = problem.reverse_original_ids; fj_cpu.h_cstr_left_weights = left_weights; diff --git a/cpp/src/mip_heuristics/local_search/local_search.cu b/cpp/src/mip_heuristics/local_search/local_search.cu index bdeabee0a6..23edf555cd 100644 --- a/cpp/src/mip_heuristics/local_search/local_search.cu +++ b/cpp/src/mip_heuristics/local_search/local_search.cu @@ -118,7 +118,7 @@ void local_search_t::start_cpufj_lptopt_scratch_threads( solution_lp.copy_new_assignment( host_copy(lp_optimal_solution, context.problem_ptr->handle_ptr->get_stream())); solution_lp.round_random_nearest(500); - scratch_cpu_fj_on_lp_opt = fj.create_cpu_climber(solution_lp, + scratch_cpu_fj_on_lp_opt = fj.create_cpu_climber(solution_lp, default_weights, default_weights, 0., From e044bed82b41c84540b2edc1e07a2d763294ff9c Mon Sep 17 00:00:00 2001 From: yboucher Date: Thu, 20 Aug 2026 00:42:37 -0700 Subject: [PATCH 8/9] some cleanup --- cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 9915073a96..5bd0a288de 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -754,8 +754,8 @@ static void two_opt_add_partner(fj_cpu_climber_t& fj_cpu, * Preferred source is the probing cache: it recorded, for each probed variable and value, the * bounds propagation implies on every other variable. An implied bound pinning a binary to a value * names both the partner and the value it has to take once `first` moves, so a pair moving in the - * same direction is reached as naturally as a swap. When probing has nothing for `first`, the - * variables sharing a row with it are the fallback, and the only defensible direction is the swap. + * same direction is reached as naturally as a swap. The + * variables sharing a row with it are used as fallback. */ template static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, From 087c3aadeac51fa691a139abedea90940a8aa947 Mon Sep 17 00:00:00 2001 From: yboucher Date: Thu, 20 Aug 2026 00:48:36 -0700 Subject: [PATCH 9/9] consider both 2opt sources --- cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu index 5bd0a288de..5d9b0267b1 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu @@ -801,7 +801,6 @@ static void two_opt_collect_partners(fj_cpu_climber_t& fj_cpu, } } } - if (!partners.empty()) return; const auto& related = fj_cpu.h_related_variables; const auto& related_offsets = fj_cpu.h_related_variables_offsets;