diff --git a/benchmarks/linear_programming/cuopt/run_cpufj.cu b/benchmarks/linear_programming/cuopt/run_cpufj.cu new file mode 100644 index 0000000000..e259c2ffc0 --- /dev/null +++ b/benchmarks/linear_programming/cuopt/run_cpufj.cu @@ -0,0 +1,811 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "miplib2017_bks.hpp" + +#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 +#include +#include +#include +#include +#include +#include +#include + +#include "initial_problem_check.hpp" + +namespace { + +using i_t = int; +using f_t = double; +namespace mip = cuopt::mathematical_optimization::mip; + +using clk = std::chrono::high_resolution_clock; +double since(clk::time_point t0) +{ + return std::chrono::duration_cast>(clk::now() - t0).count(); +} + +struct climber_result_t { + bool crossed{false}; + double t_first{-1.0}; + f_t best_objective{std::numeric_limits::infinity()}; + i_t iterations{0}; + double seconds{0.0}; +}; + +void pin_to_core(int core) +{ + cpu_set_t set; + CPU_ZERO(&set); + CPU_SET(core, &set); + pthread_setaffinity_np(pthread_self(), sizeof(set), &set); +} + +// The CPUs this process is actually permitted to run on. A cgroup mask can be non-contiguous, so +// indexing hardware_concurrency() directly would collide several climbers onto one core. +std::vector allowed_cpus() +{ + std::vector allowed; + cpu_set_t set; + CPU_ZERO(&set); + if (sched_getaffinity(0, sizeof(set), &set) == 0) { + for (int cpu = 0; cpu < CPU_SETSIZE; ++cpu) { + if (CPU_ISSET(cpu, &set)) allowed.push_back(cpu); + } + } + if (allowed.empty()) allowed.push_back(0); + return allowed; +} + +void run_climber(mip::fj_cpu_climber_t* climber, + f_t time_limit, + int core, + climber_result_t& result) +{ + pin_to_core(core); + const auto t0 = clk::now(); + + climber->improvement_callback = [&result, t0](f_t objective, const std::vector&, double) { + if (!result.crossed) { + result.crossed = true; + result.t_first = since(t0); + } + result.best_objective = objective; + }; + + mip::cpufj_solve(climber, time_limit); + + result.seconds = since(t0); + result.iterations = climber->iterations; +} + +std::vector uncrush_assignment(mip::problem_t& problem, + const std::vector& assignment, + rmm::cuda_stream_view stream) +{ + rmm::device_uvector d_assignment(assignment.size(), stream); + raft::copy(d_assignment.data(), assignment.data(), assignment.size(), stream); + problem.post_process_assignment(d_assignment, true, stream); + if (problem.has_papilo_presolve_data()) { + problem.papilo_uncrush_assignment(d_assignment, stream); + } + return cuopt::host_copy(d_assignment, stream); +} + +bool write_lane_solutions( + const std::string& dir, + const std::string& instance, + const std::vector& var_names, + const std::vector>>& climbers, + const std::vector& results, + mip::problem_t& problem, + rmm::cuda_stream_view stream, + int& written) +{ + std::error_code ec; + std::filesystem::create_directories(dir, ec); + if (ec) { return false; } + + const std::string stem = std::filesystem::path(instance).stem().string(); + for (size_t k = 0; k < climbers.size(); ++k) { + const auto& c = *climbers[k]; + if (!c.feasible_found) { continue; } + + const std::vector solver(c.h_best_assignment.data(), + c.h_best_assignment.data() + c.h_best_assignment.size()); + const std::vector user = uncrush_assignment(problem, solver, stream); + if (user.size() != var_names.size()) { return false; } + + const std::string file = + (std::filesystem::path(dir) / (stem + ".lane" + std::to_string(k) + ".sol")).string(); + cuopt::mathematical_optimization::solution_writer_t::write_solution_to_sol_file( + file, + "Feasible", + problem.get_user_obj_from_solver_obj(results[k].best_objective), + var_names, + user); + ++written; + } + return true; +} + +} // namespace + +int main(int argc, char** argv) +{ + argparse::ArgumentParser program("solve_CPUFJ"); + + program.add_argument("instance").help("input .mps path"); + + program.add_argument("time_limit") + .help("time limit in seconds") + .scan<'g', double>() + .default_value(60.0) + .nargs(argparse::nargs_pattern::optional); + + program.add_argument("climbers") + .help("number of climbers") + .scan<'i', int>() + .default_value(16) + .nargs(argparse::nargs_pattern::optional); + + program.add_argument("seed") + .help("base seed") + .scan<'i', unsigned>() + .default_value(12345u) + .nargs(argparse::nargs_pattern::optional); + + program.add_argument("--sol-dir") + .help("directory to write one .sol per feasible lane into") + .default_value(std::string("")); + + program.add_argument("--no-related-vars") + .help("skip the related-variable structure build in trivial presolve") + .flag(); + + program.add_argument("--low-latency") + .help("skip bound propagation and the binary fast path scan") + .flag(); + + program.add_argument("--probing") + .help("compute the probing cache before the solve window and hand it to the portfolio") + .flag(); + + program.add_argument("--probing-time-limit") + .help("wall clock cap for the probing cache, in seconds; not charged to the solve window") + .scan<'g', double>() + .default_value(60.0); + + try { + program.parse_args(argc, argv); + } catch (const std::exception& err) { + std::cerr << err.what() << std::endl; + std::cerr << program; + return 2; + } + + const std::string path = program.get("instance"); + const f_t time_limit = program.get("time_limit"); + const int n_climbers = program.get("climbers"); + const unsigned base_seed = program.get("seed"); + const bool no_related_vars = program.get("--no-related-vars"); + const std::string sol_dir = program.get("--sol-dir"); + const bool low_latency = program.get("--low-latency"); + const bool run_probing = program.get("--probing"); + const double probing_time = program.get("--probing-time-limit"); + + // Console sink so the engine's end-of-solve incumbent audit is visible, as solve_MIP does it. + cuopt::init_logger_t log_guard("", true); + + raft::handle_t handle; + + const auto mps_data_model = cuopt::mathematical_optimization::io::read_mps(path, false); + const auto op_problem = + cuopt::mathematical_optimization::mps_data_model_to_optimization_problem( + &handle, mps_data_model); + mip::problem_t problem(op_problem); + + // Anonymise the instance before anything under evolution can see it. + // + // problem_t exposes var_names, row_names and objective_name as public members, and + // the FJ code receives problem_t&. For a fixed benchmark set those strings are an + // exact fingerprint -- row_names[0] alone identifies most MIPLIB instances -- so a + // candidate could branch on identity and return a memorised objective. Reading the + // MODEL is intended and useful: coefficients, bounds, variable types, sparsity and + // row structure are all untouched here, so recognising set-packing rows, knapsack + // substructure or GUB constraints still works exactly as before. Only the labels go. + // + // Each string is cleared in place rather than the vectors being emptied, so size() + // and indexing stay valid and any code that walks names by variable index still + // works -- it just gets empty strings. + // + // This file is outside target_code and is sha256-gated by evaluate.py's FROZEN_FILES, + // so a candidate cannot restore the names. Do not move this below the solve. + for (auto& name : problem.var_names) + name.clear(); + for (auto& name : problem.row_names) + name.clear(); + problem.objective_name.clear(); + + // The same pair solve.cu runs before handing the model to the heuristics. Without it the harness + // solves a model the production path never produces: trivial presolve is what drops explicit zero + // coefficients, and every consumer downstream is entitled to assume a nonzero is nonzero. + // trivial_presolve requires preprocess_problem first and says so. + problem.preprocess_problem(); + mip::trivial_presolve( + problem, /*remap_cache_ids=*/true, /*compute_related_vars=*/!no_related_vars); + + std::unique_ptr> probing_context; + std::unique_ptr> probing_presolve; + if (run_probing) { + cuopt::mathematical_optimization::mip_solver_settings_t probing_settings; + probing_settings.seed = (i_t)base_seed; + probing_context = + std::make_unique>(&handle, &problem, probing_settings); + probing_presolve = std::make_unique>(*probing_context); + + const auto probing_features = mip::probing_presolve_features(problem); + const auto probing_budget = + mip::evaluate_presolve_budget(probing_settings.heuristic_params, probing_features); + const int probing_threads = std::max(1, (int)allowed_cpus().size()); + probing_presolve->settings.num_tasks = std::max(1, probing_threads - 1); + + const auto probing_t0 = clk::now(); + bool infeasible = false; + 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(probing_threads) + { +#pragma omp masked + { + infeasible = mip::compute_probing_cache(*probing_presolve, + problem, + cuopt::timer_t{probing_time}, + probing_budget.probing_work_limit, + (size_t)probing_budget.probing_step_size); + } + } + if (saved_max_active_levels < 2) { omp_set_max_active_levels(saved_max_active_levels); } + handle.sync_stream(); + + if (infeasible) { + std::printf("probing: problem proved infeasible\n"); + return 1; + } + mip::trivial_presolve( + problem, /*remap_cache_ids=*/true, /*compute_related_vars=*/!no_related_vars); + std::printf("probing: %zu cached vars %.3fs threads=%d work_limit=%.3f step=%d\n", + probing_presolve->probing_cache.probing_cache.size(), + since(probing_t0), + probing_threads, + probing_budget.probing_work_limit, + probing_budget.probing_step_size); + } + + std::printf("instance: %s n_vars=%d n_cstrs=%d nnz=%d\n", + path.c_str(), + problem.n_variables, + problem.n_constraints, + problem.nnz); + + // Taken from the host-side parse, so it is independent of everything under target_code. + { + const auto& col_indices = mps_data_model.get_constraint_matrix_indices(); + const auto& row_lb = mps_data_model.get_constraint_lower_bounds(); + const auto& row_ub = mps_data_model.get_constraint_upper_bounds(); + const int64_t nnz = (int64_t)col_indices.size(); + + const i_t n_cols = mps_data_model.get_n_variables(); + std::vector degree(n_cols, 0); + for (i_t index : col_indices) { + if (index >= 0 && index < n_cols) ++degree[index]; + } + std::sort(degree.begin(), degree.end()); + + const i_t max_degree = degree.empty() ? 0 : degree.back(); + auto quantile = [&](double q) { + return degree.empty() + ? 0 + : degree[std::min(degree.size() - 1, (size_t)(q * degree.size()))]; + }; + int64_t top10 = 0; + for (size_t k = 0; k < 10 && k < degree.size(); ++k) + top10 += degree[degree.size() - 1 - k]; + const double mean_degree = n_cols > 0 ? (double)nnz / n_cols : 0.0; + std::printf( + "census cols: n=%d degree max=%d p99=%d p90=%d median=%d mean=%.1f" + " widest=%.1f%% top10=%.1f%% of nnz hub=%.0fx mean\n", + n_cols, + max_degree, + quantile(0.99), + quantile(0.90), + quantile(0.50), + mean_degree, + nnz > 0 ? 100.0 * max_degree / nnz : 0.0, + nnz > 0 ? 100.0 * top10 / nnz : 0.0, + mean_degree > 0 ? max_degree / mean_degree : 0.0); + + const i_t n_rows = (i_t)std::min(row_lb.size(), row_ub.size()); + i_t lb_only = 0, ub_only = 0, equality = 0, ranged = 0, free_rows = 0; + for (i_t r = 0; r < n_rows; ++r) { + const bool has_lb = std::isfinite((double)row_lb[r]); + const bool has_ub = std::isfinite((double)row_ub[r]); + if (has_lb && has_ub) { + ++(row_lb[r] == row_ub[r] ? equality : ranged); + } else if (has_lb) { + ++lb_only; + } else if (has_ub) { + ++ub_only; + } else { + ++free_rows; + } + } + std::printf( + "census rows: n=%d lb_only=%d ub_only=%d equality=%d ranged=%d free=%d" + " one_sided=%.1f%%\n", + n_rows, + lb_only, + ub_only, + equality, + ranged, + free_rows, + n_rows > 0 ? 100.0 * (lb_only + ub_only) / n_rows : 0.0); + } + + // FROZEN -- defines t=0 for the benchmark. Everything above it (the MPS parse, + // problem construction under problem/, preprocess and trivial presolve, and the + // name anonymisation) is outside target_code; everything below it is editable. A + // marker any later would leave editable code ahead of the clock, which is somewhere + // to do unmeasured work; any earlier would charge the budget for a parse and a CUDA + // context no candidate can influence. + CUOPT_LOG_INFO("CPUFJ solve window start"); + + std::vector> preemption_flags(n_climbers); + std::vector>> climbers(n_climbers); + mip::build_climber_portfolio( + problem, preemption_flags, climbers, base_seed, low_latency); + if (probing_presolve != nullptr) { + for (int k = 0; k < n_climbers; ++k) + const_cast*>(climbers[k]->problem.get())->probing_cache = + &probing_presolve->probing_cache; + } + for (int k = 0; k < n_climbers; ++k) { + climbers[k]->log_prefix = "[climber " + std::to_string(k) + "] "; + climbers[k]->log_interval = 1000; + } + + const std::vector cpus = allowed_cpus(); + std::printf("running %d climbers x %.0fs, base seed %u, %zu allowed CPUs (%d..%d)\n", + n_climbers, + (double)time_limit, + base_seed, + cpus.size(), + cpus.front(), + cpus.back()); + + std::vector results(n_climbers); + std::vector threads; + threads.reserve(n_climbers); + const auto wall0 = clk::now(); + for (int k = 0; k < n_climbers; ++k) { + threads.emplace_back( + run_climber, climbers[k].get(), time_limit, cpus[k % cpus.size()], std::ref(results[k])); + } + for (auto& t : threads) { + t.join(); + } + const double wall = since(wall0); + + int crossed = 0; + double sum_iters = 0; + f_t best_overall = std::numeric_limits::infinity(); + std::printf("\n climber | crossed | t_first(s) | obj | iters | iters/s\n"); + std::printf("---------+---------+------------+--------------+----------+---------\n"); + for (int k = 0; k < n_climbers; ++k) { + const auto& r = results[k]; + sum_iters += r.iterations; + if (r.crossed) { + ++crossed; + best_overall = std::min(best_overall, r.best_objective); + } + std::printf(" %7d | %7s | %10s | %12.6g | %8d | %8.0f\n", + k, + r.crossed ? "YES" : "no", + r.crossed ? std::to_string(r.t_first).c_str() : "-", + r.crossed ? (double)problem.get_user_obj_from_solver_obj(r.best_objective) : 0.0, + r.iterations, + r.seconds > 0 ? r.iterations / r.seconds : 0.0); + } + // Runs after the measured window closes, so its cost is off the clock. + // Solver space is always a minimisation, so beating the best known is always a smaller value. + const auto bks_user = cuopt_bench::lookup_miplib_bks(path); + const double bks = bks_user ? (double)problem.get_solver_obj_from_user_obj((f_t)*bks_user) : 0.0; + const double bks_slack = std::max(1e-6, std::fabs(bks) * 1e-9); + + int audited = 0, invalid = 0; + int64_t escalated_rows = 0; + std::printf( + "\n climber | viol rows worst/tol | bnd viol worst/tol | int viol worst/tol |" + " obj drift rel | vs bks\n"); + std::printf( + "---------+----------------------+---------------------+---------------------+" + "----------------------+----------\n"); + for (int k = 0; k < n_climbers; ++k) { + auto& c = *climbers[k]; + if (c.feasible_found != results[k].crossed) { + std::printf(" %7d | feasible_found=%d disagrees with a reported incumbent=%d\n", + k, + (int)c.feasible_found, + (int)results[k].crossed); + ++invalid; + continue; + } + if (!c.feasible_found) continue; + ++audited; + + const auto& cpu_problem = *c.problem; + const double int_tol = cpu_problem.tolerances.integrality_tolerance; + + i_t rows_over = 0; + i_t rows_exact = 0; + double worst_row_ratio = 0.0; + // Both counters and the max are order-independent, so any schedule leaves the verdict + // identical. Per-row cost tracks the row's nnz, which is skewed, hence guided. +#pragma omp parallel for num_threads(n_climbers) schedule(guided) \ + reduction(+ : rows_over, rows_exact) reduction(max : worst_row_ratio) + for (i_t r = 0; r < cpu_problem.n_constraints; ++r) { + const i_t begin = cpu_problem.offsets[r]; + const i_t end = cpu_problem.offsets[r + 1]; + const f_t lb = cpu_problem.cstr_lb[r]; + const f_t ub = cpu_problem.cstr_ub[r]; + + const double row_tol = + mip::get_cstr_tolerance(lb, + ub, + cpu_problem.tolerances.absolute_tolerance, + cpu_problem.tolerances.relative_tolerance); + const double tol = std::max(row_tol, 1e-12); + + // Naive summation over w products: each product carries eps/2 and each of the w-1 additions + // carries eps, both against the running magnitude, so the row's error is within + // (w + 1) * eps * abs_sum. A verdict farther from the tolerance than that cannot flip. + // Only rows the double pass cannot place on one side of the tolerance pay for _Float128, + // which is soft-float on x86-64. Includes rows whose double excess is zero but whose error + // bound alone exceeds the tolerance: a real violation can hide there. + const auto verdict = check_row( + cpu_problem.coefficients.data(), + cpu_problem.variables.data(), + (int64_t)begin, + (int64_t)end, + c.h_best_assignment.data(), + (double)lb, + (double)ub, + [&](double) { return std::pair{(double)lb - tol, (double)ub + tol}; }); + if (verdict.escalated) ++rows_exact; + if (verdict.raw_excess <= 0.0) continue; + + const double ratio = + tol > 0 ? verdict.raw_excess / tol : std::numeric_limits::infinity(); + if (ratio > 1.0) ++rows_over; + worst_row_ratio = std::max(worst_row_ratio, ratio); + } + escalated_rows += rows_exact; + + i_t bounds_over = 0; + i_t integers_over = 0; + double worst_bound_ratio = 0.0; + double worst_integer_ratio = 0.0; + _Float128 objective = 0; + for (i_t v = 0; v < cpu_problem.n_variables; ++v) { + const auto bounds = c.h_var_bounds[v].get(); + const double x = (double)c.h_best_assignment[v]; + const double out = std::max( + std::max((double)cuopt::get_lower(bounds) - x, x - (double)cuopt::get_upper(bounds)), 0.0); + if (out > int_tol) ++bounds_over; + worst_bound_ratio = std::max(worst_bound_ratio, int_tol > 0 ? out / int_tol : 0.0); + + if (cpu_problem.h_var_types[v] == cuopt::mathematical_optimization::var_t::INTEGER) { + const double residual = std::fabs(x - std::round(x)); + if (residual > int_tol) ++integers_over; + worst_integer_ratio = std::max(worst_integer_ratio, int_tol > 0 ? residual / int_tol : 0.0); + } + const double coefficient = cpu_problem.h_obj_coeffs[v]; + objective += (_Float128)coefficient * (_Float128)x; + } + + // Differenced before narrowing; the drift is smaller than a double ulp of the sum. + const _Float128 difference = objective - (_Float128)results[k].best_objective; + const double drift = (double)(difference < 0 ? -difference : difference); + const double exact = (double)objective; + const double scale = std::max(std::fabs(exact), 1.0); + const bool below_bks = bks_user && exact < bks - bks_slack; + const bool bad = rows_over > 0 || bounds_over > 0 || integers_over > 0 || below_bks; + if (bad) ++invalid; + std::printf(" %7d | %9d %10.3g | %8d %10.3g | %8d %10.3g | %12.3g %6.1e | %9.3g%s%s\n", + k, + rows_over, + worst_row_ratio, + bounds_over, + worst_bound_ratio, + integers_over, + worst_integer_ratio, + drift, + drift / scale, + bks_user ? exact - bks : 0.0, + below_bks ? " BELOW BKS" : "", + bad ? " INVALID" : ""); + } + std::printf( + "AUDIT: %d/%d reporting climbers checked, %d invalid, %lld rows re-summed exactly," + " bks %s\n", + audited, + crossed, + invalid, + (long long)escalated_rows, + bks_user ? std::to_string(*bks_user).c_str() + : (cuopt_bench::is_known_infeasible(path) ? "known infeasible" : "unknown")); + + // checked the uncrushed solution against the original model + { + const auto& A_val = mps_data_model.get_constraint_matrix_values(); + const auto& A_idx = mps_data_model.get_constraint_matrix_indices(); + const auto& A_off = mps_data_model.get_constraint_matrix_offsets(); + const auto& row_lb = mps_data_model.get_constraint_lower_bounds(); + const auto& row_ub = mps_data_model.get_constraint_upper_bounds(); + const auto& col_lb = mps_data_model.get_variable_lower_bounds(); + const auto& col_ub = mps_data_model.get_variable_upper_bounds(); + const auto& v_type = mps_data_model.get_variable_types(); + const i_t n_orig_rows = (i_t)A_off.size() - 1; + const double abs_tol = problem.tolerances.absolute_tolerance; + const double int_tol = problem.tolerances.integrality_tolerance; + + int lifted_checked = 0, lifted_bad = 0; + for (int k = 0; k < n_climbers; ++k) { + const auto& c = *climbers[k]; + if (!c.feasible_found) continue; + const std::vector solver(c.h_best_assignment.data(), + c.h_best_assignment.data() + c.h_best_assignment.size()); + const std::vector user = uncrush_assignment(problem, solver, handle.get_stream()); + if ((i_t)user.size() != (i_t)col_lb.size()) { + std::printf("LIFTED AUDIT: climber %d produced %d values for %d original columns\n", + k, + (int)user.size(), + (int)col_lb.size()); + ++lifted_bad; + continue; + } + ++lifted_checked; + + i_t bad_rows = 0, bad_bnd = 0, bad_int = 0; + double worst_row = 0.0; + i_t worst_row_id = -1; + for (i_t r = 0; r < n_orig_rows; ++r) { + const double lb = (double)row_lb[r]; + const double ub = (double)row_ub[r]; + const auto verdict = + check_row(A_val.data(), + A_idx.data(), + (int64_t)A_off[r], + (int64_t)A_off[r + 1], + user.data(), + lb, + ub, + [&](double positive) { return scaled_row_limits(abs_tol, positive, lb, ub); }); + if (verdict.excess > 0.0) { + ++bad_rows; + if (verdict.excess > worst_row) { + worst_row = verdict.excess; + worst_row_id = r; + } + } + } + for (i_t v = 0; v < (i_t)user.size(); ++v) { + const double x = (double)user[v]; + if (x < (double)col_lb[v] - int_tol || x > (double)col_ub[v] + int_tol) ++bad_bnd; + if ((v_type[v] == 'I' || v_type[v] == 'B') && std::fabs(x - std::round(x)) > int_tol) + ++bad_int; + } + if (bad_rows || bad_bnd || bad_int) { + ++lifted_bad; + std::printf( + "LIFTED AUDIT: climber %d INVALID on the original model -- %d rows, %d bounds," + " %d integrality; worst row %d by %.6g\n", + k, + (int)bad_rows, + (int)bad_bnd, + (int)bad_int, + (int)worst_row_id, + worst_row); + } + } + std::printf( + "LIFTED AUDIT: %d/%d crossing climbers verified against the original model," + " %d INVALID, %d original rows\n", + lifted_checked, + crossed, + lifted_bad, + (int)n_orig_rows); + } + + std::printf( + "\n climber | moves | apply nnz | nnz/move | bitmap elems | ratio |" + " bump/apply | bump/weight | mtm inval | cache hit%%\n"); + std::printf( + "---------+-----------+------------+----------+--------------+-------+" + "------------+-------------+-----------+-----------\n"); + for (int k = 0; k < n_climbers; ++k) { + const auto& c = *climbers[k]; + const int64_t bitmap = 2 * c.n_moves_applied * (int64_t)c.problem->n_variables; + const int64_t probes = c.hit_count + c.miss_count; + std::printf( + " %7d | %9lld | %10lld | %8.1f | %12lld | %5.0f | %10lld | %11lld | %9lld |" + " %9.2f\n", + k, + (long long)c.n_moves_applied, + (long long)c.apply_move_nnz, + c.n_moves_applied > 0 ? (double)c.apply_move_nnz / c.n_moves_applied : 0.0, + (long long)bitmap, + c.apply_move_nnz > 0 ? (double)bitmap / c.apply_move_nnz : 0.0, + (long long)c.n_version_bumps_apply, + (long long)c.n_version_bumps_weights, + (long long)c.n_mtm_cache_invalidations, + probes > 0 ? 100.0 * c.hit_count / probes : 0.0); + } + + std::printf( + "\n climber | mtm calls | row entries | ent/call | capped ent | capped/call |" + " score calls | score nnz | nnz/score | nnz budget\n"); + std::printf( + "---------+-----------+-------------+----------+-------------+-------------+" + "-------------+-----------+-----------+-----------\n"); + for (int k = 0; k < n_climbers; ++k) { + const auto& c = *climbers[k]; + std::printf( + " %7d | %9lld | %11lld | %8.0f | %11lld | %11.0f | %11lld | %9lld | %9.1f |" + " %10d\n", + k, + (long long)c.n_mtm_calls, + (long long)c.mtm_row_entries, + c.n_mtm_calls > 0 ? (double)c.mtm_row_entries / c.n_mtm_calls : 0.0, + (long long)c.mtm_entries_capped, + c.n_mtm_calls > 0 ? (double)c.mtm_entries_capped / c.n_mtm_calls : 0.0, + (long long)c.n_compute_score_calls, + (long long)c.compute_score_nnz, + c.n_compute_score_calls > 0 ? (double)c.compute_score_nnz / c.n_compute_score_calls : 0.0, + c.nnz_samples); + } + + std::printf( + "\n climber | refresh period | lhs total | periodic | bigval | perturb | restart |" + " epi vars | epi projections\n"); + std::printf( + "---------+----------------+-----------+----------+--------+---------+---------+" + "----------+----------------\n"); + for (int k = 0; k < n_climbers; ++k) { + const auto& c = *climbers[k]; + std::printf(" %7d | %14d | %9lld | %8lld | %6lld | %7lld | %7lld | %8zu | %15lld\n", + k, + c.lhs_refresh_period_used, + (long long)c.n_lhs_recompute_total, + (long long)c.n_lhs_recompute_periodic, + (long long)c.n_lhs_recompute_bigval, + (long long)c.n_lhs_recompute_perturb, + (long long)c.n_lhs_recompute_restart, + c.epigraph_vars.size(), + (long long)c.n_epigraph_projections); + } + + // Everything a climber spends outside the search loop. lp solve is the simplex share of the LP + // start, so it is shown for attribution and left out of the total. A phase a lane does not run + // reads 0. + std::printf( + "\n climber | start | bnd prop | lp start | (lp solve) | colouring | features |" + " init lhs | bin setup | total\n"); + std::printf( + "---------+----------+----------+----------+------------+-----------+----------+" + "----------+-----------+---------\n"); + for (int k = 0; k < n_climbers; ++k) { + const auto& c = *climbers[k]; + const double total = c.t_start + c.t_bound_prop + c.t_lp_start + c.t_coloring + c.t_features + + c.t_init_lhs + c.bin_setup.total(); + std::printf(" %7d | %8.4f | %8.4f | %8.4f | %10.4f | %9.4f | %8.4f | %8.4f | %9.4f | %8.4f\n", + k, + c.t_start, + c.t_bound_prop, + c.t_lp_start, + c.t_lp_relaxation, + c.t_coloring, + c.t_features, + c.t_init_lhs, + c.bin_setup.total(), + total); + } + + // The bin setup column above, by phase. Charged even when the fast path declines, so a scan that + // only produces a rejection still shows. narrow and transpose are the all-binary path; encode is + // the general-integer one and runs twice when int8 is enough. + std::printf( + "\n climber | bin scan | bin narrow | transpose | bin encode | engine init | bin total\n"); + std::printf( + "---------+----------+------------+-----------+------------+-------------+----------\n"); + for (int k = 0; k < n_climbers; ++k) { + const auto& b = climbers[k]->bin_setup; + std::printf(" %7d | %8.4f | %10.4f | %9.4f | %10.4f | %11.4f | %9.4f\n", + k, + b.scan, + b.narrow, + b.transpose, + b.encode, + b.engine_init, + b.total()); + } + + std::printf("\nSUMMARY: %d/%d crossed (%.0f%%) wall=%.1fs total_iters=%.0f agg_iters/s=%.0f\n", + crossed, + n_climbers, + 100.0 * crossed / n_climbers, + wall, + sum_iters, + wall > 0 ? sum_iters / wall : 0.0); + if (crossed > 0) { + std::printf("BEST OBJECTIVE: %.10g\n", + (double)problem.get_user_obj_from_solver_obj(best_overall)); + } + + if (!sol_dir.empty()) { + int written = 0; + const bool ok = write_lane_solutions(sol_dir, + path, + mps_data_model.get_variable_names(), + climbers, + results, + problem, + handle.get_stream(), + written); + std::printf("SOLUTIONS: %d/%d lanes -> %s (%s)\n", + written, + n_climbers, + sol_dir.c_str(), + ok ? "written" : "WRITE FAILED"); + } + + return 0; +} diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 555b55273d..5b1dfe7d5e 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -1543,6 +1543,30 @@ if (BUILD_MIP_BENCHMARKS AND NOT BUILD_LP_ONLY) "${CMAKE_CURRENT_SOURCE_DIR}/src" ) + add_executable(solve_CPUFJ ../benchmarks/linear_programming/cuopt/run_cpufj.cu) + set_target_properties(solve_CPUFJ PROPERTIES CXX_SCAN_FOR_MODULES OFF) + target_compile_options(solve_CPUFJ + PRIVATE "$<$:${CUOPT_CXX_FLAGS}>" + "$<$:${CUOPT_CUDA_FLAGS}>" + "$<$:-fopenmp>" + ) + target_link_libraries(solve_CPUFJ + PUBLIC + cuopt_static + OpenMP::OpenMP_CXX + OpenMP::OpenMP_CUDA + ) + target_include_directories(solve_CPUFJ + PRIVATE + "${CMAKE_CURRENT_SOURCE_DIR}/src" + "${papilo_SOURCE_DIR}/src" + "${papilo_BINARY_DIR}" + ) + target_include_directories(solve_CPUFJ SYSTEM PRIVATE + "${pslp_SOURCE_DIR}/include" + "${dejavu_SOURCE_DIR}" + ) + endif () option(BUILD_LP_BENCHMARKS "Build LP benchmarks" OFF) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 7d82f12134..c185fc580e 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -3064,8 +3064,14 @@ void branch_and_bound_t::launch_root_heuristics( [this](f_t obj, const std::vector& assignment, double work_units) { set_solution_from_cpu_fj(obj, assignment, work_units); }; - current_heuristic->fj_cpu_worker_.create_worker( - lp, var_types_, original_problem_.num_cols, lp_solution.x, settings_, "[RootCut CPUFJ] "); + current_heuristic->fj_cpu_worker_.create_worker(lp, + var_types_, + original_problem_.num_cols, + lp_solution.x, + settings_, + "[RootCut CPUFJ] ", + -1, + cut_pass); ++(*worker_count); ++current_heuristic->active_workers_; diff --git a/cpp/src/mip_heuristics/CMakeLists.txt b/cpp/src/mip_heuristics/CMakeLists.txt index 5771b69fe6..336b4486be 100644 --- a/cpp/src/mip_heuristics/CMakeLists.txt +++ b/cpp/src/mip_heuristics/CMakeLists.txt @@ -50,6 +50,11 @@ set(MIP_NON_LP_FILES ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/loop.cpp ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/portfolio.cpp ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/search/escape.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/starts/affine.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/starts/cardinality.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/starts/chain.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/starts/covering.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/setup/bounds.cpp ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/setup/lp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/cpu/setup/structure.cpp ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/fj_cpu_bridge.cu diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/climber.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/climber.cpp index e132eb8327..89c84c02ad 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/climber.cpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/climber.cpp @@ -9,6 +9,7 @@ #include "internal.hpp" #include "problem.hpp" #include "search/api.hpp" +#include "setup/bounds.hpp" #include "setup/lp.hpp" #include "setup/structure.hpp" @@ -93,6 +94,7 @@ void wire_fj_cpu_host_views( set_host_data_view(fj_cpu, n_variables, n_constraints, n_integer_vars, nnz, tolerances); + cap_integer_domains(fj_cpu, n_variables); fj_cpu.h_best_objective = +std::numeric_limits::infinity(); // cached_mtm_moves, cached_mtm_moves_version and h_cstr_version are indexed by search row and @@ -119,6 +121,8 @@ void finalize_fj_cpu_host_initialization( detect_implied_integers(fj_cpu, problem); wire_fj_cpu_host_views(fj_cpu, n_variables, n_constraints, n_integer_vars, nnz, tolerances); + build_cardinality_index(fj_cpu, problem); + detect_free_equality_singletons(fj_cpu); problem.h_objective_vars.resize(n_variables); auto end = std::copy_if( @@ -151,6 +155,8 @@ void finalize_fj_cpu_host_initialization( phase_timer_t timer(fj_cpu.t_init_lhs); recompute_lhs(fj_cpu); } + + precompute_problem_features(fj_cpu, problem); } template diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/loop.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/loop.cpp index bafb3bbf72..7cdb40446b 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/loop.cpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/loop.cpp @@ -13,16 +13,165 @@ #include "search/moves.hpp" #include "search/score.hpp" #include "search/update.hpp" +#include "setup/bounds.hpp" +#include "setup/lp.hpp" #include "setup/structure.hpp" +#include "starts/starts.hpp" #include namespace cuopt::mathematical_optimization::mip { +template +static std::vector lift_equality_substituted_assignment( + fj_cpu_climber_t& c, + const std::vector& assignment, + const std::vector& retained, + const std::vector>& substitutions) +{ + std::vector lifted(c.problem->n_variables, 0); + for (size_t j = 0; j < retained.size(); ++j) + lifted[retained[j]] = assignment[j]; + for (auto it = substitutions.rbegin(); it != substitutions.rend(); ++it) { + const auto coefficients = thrust::make_transform_iterator( + it->terms.begin(), [](const auto& term) { return term.second; }); + const auto values = thrust::make_transform_iterator( + it->terms.begin(), [&lifted](const auto& term) { return lifted[term.first]; }); + lifted[it->variable] = it->constant + compensated_dot2(coefficients, values, it->terms.size()); + } + for (const auto& sub : substitutions) { + const auto bounds = c.h_var_bounds[sub.variable].get(); + f_t value = std::clamp(lifted[sub.variable], get_lower(bounds), get_upper(bounds)); + if (is_integer_var(c, sub.variable)) value = std::round(value); + lifted[sub.variable] = value; + } + // Substitution and the final bound/integrality repair must both survive a check in the unchanged + // parent model + const auto& p = *c.problem; + for (i_t v = 0; v < p.n_variables; ++v) { + if (!std::isfinite(lifted[v]) || !check_variable_within_bounds(c, v, lifted[v]) || + (is_integer_var(c, v) && !p.is_integer(lifted[v]))) + return {}; + } + for (i_t r = 0; r < p.n_constraints; ++r) { + const f_t activity = compensated_dot2_csr(p, lifted, r); + const f_t tol = p.tolerances.absolute_tolerance; + if (!std::isfinite(activity) || activity < p.cstr_lb[r] - tol || activity > p.cstr_ub[r] + tol) + return {}; + } + return lifted; +} + +// Solve a lane in equality-reduced coordinates, then lift every candidate back into the unchanged +// parent model before reporting it. +template +bool try_equality_substituted_solve(fj_cpu_climber_t& c, + double time_limit, + double work_unit_limit) +{ + if (!c.use_equality_substitution || c.feasible_found || c.producer_sync || time_limit <= 0) + return false; + const auto started = std::chrono::steady_clock::now(); + std::vector> substitutions; + std::vector retained; + std::unique_ptr> child; + { + phase_timer_t timer(c.t_start); + child = + make_equality_reduced_climber(c, std::min(0.75, 0.15 * time_limit), substitutions, retained); + } + if (!child) return false; + const double remaining = + time_limit - std::chrono::duration(std::chrono::steady_clock::now() - started).count(); + if (remaining <= 0) return false; + + bool rejected_lift = false; + child->work_unit_bias = c.work_unit_bias; + child->improvement_callback = + [&](f_t child_objective, const std::vector& assignment, double work) { + if (assignment.size() != retained.size()) { + rejected_lift = true; + child->halted = true; + return; + } + const std::vector lifted = + lift_equality_substituted_assignment(c, assignment, retained, substitutions); + if (lifted.empty()) { + rejected_lift = true; + child->halted = true; + return; + } + const f_t objective = child_objective + child->problem->objective_offset; + if (!c.feasible_found || objective < c.h_best_objective) { + c.h_best_assignment = lifted; + c.h_best_objective = objective; + c.feasible_found = true; + } + report_cpu_incumbent(c, objective, lifted, work); + }; + + const auto setup_stats = static_cast&>(c); + cpufj_solve(child.get(), remaining, work_unit_limit); + static_cast&>(c) = static_cast&>(*child); + c.t_start += setup_stats.t_start; + c.t_bound_prop += setup_stats.t_bound_prop; + c.t_lp_start += setup_stats.t_lp_start; + c.t_features += setup_stats.t_features; + c.t_init_lhs += setup_stats.t_init_lhs; + c.iterations = child->iterations; + c.work_units_elapsed.store(child->work_units_elapsed.load(std::memory_order_relaxed), + std::memory_order_relaxed); + + if (child->feasible_found && child->h_best_assignment.size() == retained.size()) { + std::vector lifted = lift_equality_substituted_assignment( + c, child->h_best_assignment.underlying(), retained, substitutions); + if (lifted.empty()) return c.feasible_found; + const f_t objective = compensated_dot2( + thrust::make_permutation_iterator(c.problem->h_obj_coeffs.data(), + c.problem->h_objective_vars.data()), + thrust::make_permutation_iterator(lifted.data(), c.problem->h_objective_vars.data()), + c.problem->h_objective_vars.size()); + if (!c.feasible_found || objective < c.h_best_objective) { + c.h_best_assignment = std::move(lifted); + c.h_best_objective = objective; + c.feasible_found = true; + } + } + return !rejected_lift || c.feasible_found; +} + template void cpufj_solve(fj_cpu_climber_t* fj_cpu, double time_limit, double work_unit_limit) { - if (try_cpufj_binary_solve(*fj_cpu, time_limit, work_unit_limit)) return; + const auto solve_start = std::chrono::steady_clock::now(); + if (fj_cpu->use_precedence_start) apply_precedence_completion_start(*fj_cpu); + apply_bound_propagation(*fj_cpu); + if (fj_cpu->use_equality_substitution) { + const double elapsed = + std::chrono::duration(std::chrono::steady_clock::now() - solve_start).count(); + if (try_equality_substituted_solve(*fj_cpu, time_limit - elapsed, work_unit_limit)) return; + } + const double setup_time_left = + fj_cpu->use_equality_substitution + ? std::max( + 0.0, + time_limit - + std::chrono::duration(std::chrono::steady_clock::now() - solve_start).count()) + : time_limit; + if (!fj_cpu->feasible_found) { apply_lp_rounded_start(*fj_cpu, setup_time_left); } + + const bool paid_setup = fj_cpu->use_bound_prop || fj_cpu->use_lp_start || + fj_cpu->use_precedence_start || fj_cpu->use_equality_substitution; + const double setup_seconds = + paid_setup + ? std::chrono::duration(std::chrono::steady_clock::now() - solve_start).count() + : 0.0; + const double remaining = std::max(0.0, time_limit - setup_seconds); + if (remaining <= 0.0) return; + + if (try_cpufj_binary_solve(*fj_cpu, remaining, work_unit_limit)) return; + + clamp_start_magnitude(*fj_cpu, fj_cpu->problem->n_variables); // Past this point the search runs on one-sided rows. Everything above reasons about the original // model, which is why the rows are built here and not at construction. @@ -38,7 +187,8 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, double time_limit, double w } [[maybe_unused]] i_t local_mins = 0; - const auto loop_start = std::chrono::steady_clock::now(); + const auto loop_start = paid_setup ? solve_start : std::chrono::steady_clock::now(); + bool first_cross_needs_polish = fj_cpu->use_lp_polish; fj_cpu->rng.set_seed(fj_cpu->settings.seed); @@ -86,6 +236,13 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, double time_limit, double w break; } + if (first_cross_needs_polish && fj_cpu->feasible_found) { + first_cross_needs_polish = false; + const double elapsed = + std::chrono::duration(std::chrono::steady_clock::now() - loop_start).count(); + apply_lp_polish(*fj_cpu, fj_cpu->hp.lp_polish_budget_share * (time_limit - elapsed)); + } + // periodically recompute the slacks and violation scores // to correct any accumulated numerical errors if (fj_cpu->trigger_early_lhs_recomputation) { @@ -152,6 +309,11 @@ void cpufj_solve(fj_cpu_climber_t* fj_cpu, double time_limit, double w should_perturb = true; // Without this the counter stays above the interval and every later iteration perturbs. fj_cpu->iterations_since_best = 0; + if (fj_cpu->use_lp_polish && fj_cpu->feasible_found) { + const double elapsed = + std::chrono::duration(std::chrono::steady_clock::now() - loop_start).count(); + apply_lp_polish(*fj_cpu, fj_cpu->hp.lp_polish_budget_share * (time_limit - elapsed)); + } } if (score > fj_staged_score_t::zero() && !should_perturb) { diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/portfolio.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/portfolio.cpp index b68d7460ec..10a153c490 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/portfolio.cpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/portfolio.cpp @@ -8,6 +8,8 @@ #include "climber.hpp" #include "internal.hpp" #include "problem.hpp" +#include "setup/bounds.hpp" +#include "starts/starts.hpp" namespace cuopt::mathematical_optimization::mip { @@ -25,7 +27,8 @@ void fj_cpu_worker_t::create_worker( const std::vector& start_assignment, const simplex_solver_settings_t& settings, std::string log_prefix, - int64_t seed) + int64_t seed, + int lane) { auto new_climber = init_fj_cpu_from_host_lp( problem, variable_types, n_structural, start_assignment, settings, preemption_flag, seed); @@ -35,6 +38,7 @@ void fj_cpu_worker_t::create_worker( fj_cpu->halted = false; preemption_flag = false; is_initialized = true; + if (lane >= 0) apply_lane_diversification(*fj_cpu, lane, fj_cpu->settings.seed); } template @@ -78,12 +82,184 @@ void fj_cpu_worker_t::send_stop_signal() preemption_flag = true; } +template +void apply_lane_diversification(fj_cpu_climber_t& c, int lane, int64_t base_seed) +{ + cuopt_assert(lane >= 0, "CPUFJ lane must be nonnegative"); + + const bool extreme_hub = + c.problem->max_var_degree > 512 && c.problem->max_var_degree > 64.0 * c.problem->avg_var_degree; + const bool low_rank_integer_equalities = + c.n_binary_vars == 0 && c.n_integer_vars == c.problem->n_variables && + c.problem->equality_fraction > 0.99 && + int64_t{32} * c.problem->n_constraints < c.problem->n_variables; + + c.use_lp_start = lane == 4 || lane == 10 || lane == 11 || lane == 14; + c.lp_start_feasibility_objective = lane == 4 || lane == 10 || lane == 11; + c.use_deep_lp_pump = lane == 4; + c.use_integer_bit_encoding = lane != 0 && lane != 7 && c.n_binary_vars > 0; + c.use_lp_polish = lane == 3 || lane == 5 || lane == 8 || lane == 9 || lane == 10 || lane == 11 || + lane == 12 || lane == 13 || lane == 14; + c.use_precedence_start = lane % 8 == 0 && !c.low_latency; + c.use_equality_substitution = lane % 4 == 0 && !c.low_latency; + c.use_bound_prop = lane % 2 == 0 && !c.low_latency; + + if (lane == 13 && low_rank_integer_equalities) { + c.use_lp_start = true; + c.lp_start_feasibility_objective = true; + c.use_deep_lp_pump = false; + } + + { + phase_timer_t timer(c.t_start); + switch (lane % 8) { + case 1: apply_structural_completion_start(c); break; + case 3: apply_greedy_covering_start(c); break; + case 4: + if (!c.use_lp_start) { + apply_ambiguous_lock_start(c); + apply_greedy_covering_start(c); + } + break; + case 0: + if (lane == 8) { + apply_exact_k_start(c); + apply_greedy_covering_start(c); + } + break; + default: break; + } + if (lane == 10) apply_greedy_covering_start(c); + if (lane == 12 || lane == 15) apply_structural_completion_start(c); + if (lane == 11) { + apply_structural_completion_start(c); + apply_greedy_covering_start(c); + } + } + + std::mt19937 rng(base_seed + 7919u * lane); + c.mtm_viol_samples = std::uniform_int_distribution(10, 80)(rng); + c.mtm_sat_samples = std::uniform_int_distribution(5, 50)(rng); + c.nnz_samples = std::uniform_int_distribution(1000, 20000)(rng); + c.perturb_interval = std::uniform_int_distribution(10, 2000)(rng); + + static constexpr double smoothing[8] = {0.0003, 0.0, 0.001, 0.003, 0.0001, 0.0006, 0.002, 0.0003}; + static constexpr int tabu_min[8] = {3, 1, 5, 3, 2, 6, 4, 3}; + static constexpr int tabu_max[8] = {13, 7, 21, 13, 10, 25, 17, 13}; + const int policy = lane % 8; + c.settings.parameters.weight_smoothing_probability = smoothing[policy]; + c.settings.parameters.tabu_tenure_min = tabu_min[policy]; + c.settings.parameters.tabu_tenure_max = tabu_max[policy]; + + if (lane == 7 || lane == 8) { + c.mtm_viol_samples = 192; + c.mtm_sat_samples = 64; + c.nnz_samples = 100000; + } + if (lane == 0) { + c.mtm_viol_samples = std::uniform_int_distribution(40, 100)(rng); + c.mtm_sat_samples = std::uniform_int_distribution(20, 60)(rng); + c.nnz_samples = std::uniform_int_distribution(10000, 30000)(rng); + } + if (lane == 12) { + c.mtm_viol_samples = std::uniform_int_distribution(50, 120)(rng); + c.mtm_sat_samples = std::uniform_int_distribution(25, 70)(rng); + } + if (lane == 11) { + c.mtm_viol_samples = std::uniform_int_distribution(30, 100)(rng); + c.mtm_sat_samples = std::uniform_int_distribution(15, 50)(rng); + } + if (lane == 9 && extreme_hub) { + c.mtm_viol_samples = 8; + c.mtm_sat_samples = 3; + c.nnz_samples = 2000; + } + if (lane == 10) { + c.use_bound_prop = false; + c.use_lp_start = true; + c.lp_start_feasibility_objective = true; + c.settings.seed += 224737; + } + if (extreme_hub && lane == 5) { + c.use_lp_start = c.use_deep_lp_pump = true; + c.lp_start_feasibility_objective = false; + } + + static constexpr f_t objective_weight[4] = {2, 8, 32, 1}; + i_t continuous_objective_vars = 0; + for (i_t var : c.problem->h_objective_vars) + continuous_objective_vars += !is_integer_var(c, var); + const int64_t objective_var_count = c.problem->h_objective_vars.size(); + const bool continuous_objective_model = + objective_var_count > 0 && + int64_t{10} * continuous_objective_vars >= int64_t{9} * objective_var_count && + int64_t{10} * objective_var_count >= int64_t{c.problem->n_variables}; + c.h_objective_weight = !continuous_objective_model ? f_t{0} + : lane == 3 ? f_t{4} + : lane == 9 ? f_t{8} + : lane == 15 ? f_t{16} + : f_t{0}; + c.seed_objective_weight = lane == 1 ? f_t{32} + : lane == 4 ? f_t{8} + : lane == 5 ? f_t{16} + : lane == 15 ? f_t{8} + : objective_weight[lane % 4]; +} + +template +void complete_climber_portfolio(std::unique_ptr> first_climber, + const std::vector& lane_seed, + std::vector>& preemption_flags, + std::vector>>& climbers, + int64_t base_seed, + bool low_latency) +{ + const int n_climbers = climbers.size(); + cuopt_assert(n_climbers > 0, "a CPUFJ portfolio needs at least one climber"); + cuopt_assert(preemption_flags.size() == climbers.size(), "preemption flag count mismatch"); + cuopt_assert(lane_seed.size() == climbers.size(), "lane seed count mismatch"); + + climbers[0] = std::move(first_climber); + cuopt_assert(climbers[0] != nullptr, "missing first CPUFJ climber"); + climbers[0]->low_latency = low_latency; + apply_exact_k_start(*climbers[0]); + repair_difficult_anchor(*climbers[0]); + apply_lane_diversification(*climbers[0], 0, base_seed); + +#ifdef _OPENMP +#pragma omp parallel for num_threads(std::max(1, n_climbers - 1)) schedule(static) +#endif + for (int k = 1; k < n_climbers; ++k) { + fj_settings_t settings; + settings.seed = lane_seed[k]; + climbers[k] = init_fj_cpu_clone(*climbers[0], preemption_flags[k], settings); + climbers[k]->low_latency = low_latency; + apply_lane_diversification(*climbers[k], k, base_seed); + } +} + #if MIP_INSTANTIATE_FLOAT template struct fj_cpu_worker_t; +template void apply_lane_diversification(fj_cpu_climber_t&, int, int64_t); +template void complete_climber_portfolio( + std::unique_ptr>, + const std::vector&, + std::vector>&, + std::vector>>&, + int64_t, + bool); #endif #if MIP_INSTANTIATE_DOUBLE template struct fj_cpu_worker_t; +template void apply_lane_diversification(fj_cpu_climber_t&, int, int64_t); +template void complete_climber_portfolio( + std::unique_ptr>, + const std::vector&, + std::vector>&, + std::vector>>&, + int64_t, + bool); #endif } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/bounds.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/bounds.cpp new file mode 100644 index 0000000000..a07eb9a7f5 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/bounds.cpp @@ -0,0 +1,332 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "bounds.hpp" +#include "../audit.hpp" +#include "../internal.hpp" +#include "../problem.hpp" +#include "../search/api.hpp" + +namespace cuopt::mathematical_optimization::mip { + +template +void cap_integer_domains(fj_cpu_climber_t& fj_cpu, i_t n_variables) +{ + cuopt_assert(fj_cpu.h_best_assignment.size() == static_cast(n_variables), + "best assignment size mismatch"); + + for (i_t var = 0; var < n_variables; ++var) { + if (var_t::INTEGER != fj_cpu.problem->h_var_types[var]) continue; + + const auto bounds = fj_cpu.h_var_bounds[var].get(); + const f_t lower = std::max(get_lower(bounds), (f_t)-fj_cpu.hp.integer_domain_limit); + const f_t upper = std::min(get_upper(bounds), (f_t)fj_cpu.hp.integer_domain_limit); + if (lower > upper) continue; + if (lower == get_lower(bounds) && upper == get_upper(bounds)) continue; + + // Both assignments, since the from-template path inherits h_lhs instead of recomputing it. + fj_cpu.h_var_bounds[var] = typename type_2::type{lower, upper}; + fj_cpu.h_assignment[var] = std::clamp((f_t)fj_cpu.h_assignment[var], lower, upper); + fj_cpu.h_best_assignment[var] = std::clamp((f_t)fj_cpu.h_best_assignment[var], lower, upper); + } +} + +template +void clamp_start_magnitude(fj_cpu_climber_t& fj_cpu, i_t n_variables) +{ + cuopt_assert(fj_cpu.h_assignment.size() == static_cast(n_variables), + "assignment size mismatch"); + cuopt_assert(fj_cpu.h_best_assignment.size() == static_cast(n_variables), + "best assignment size mismatch"); + + for (i_t var = 0; var < n_variables; ++var) { + const auto bounds = fj_cpu.h_var_bounds[var].get(); + const f_t lower = std::max(get_lower(bounds), (f_t)-fj_cpu.hp.start_magnitude_limit); + const f_t upper = std::min(get_upper(bounds), (f_t)fj_cpu.hp.start_magnitude_limit); + if (lower > upper) continue; + + fj_cpu.h_assignment[var] = std::clamp((f_t)fj_cpu.h_assignment[var], lower, upper); + fj_cpu.h_best_assignment[var] = std::clamp((f_t)fj_cpu.h_best_assignment[var], lower, upper); + } +} + +template +bool tighten_lower_bound(fj_cpu_climber_t& fj_cpu, + std::vector& lower, + const std::vector& upper, + i_t var, + f_t limit, + f_t commit_threshold) +{ + if (!std::isfinite(limit)) return false; + if (is_integer_var(fj_cpu, var)) + limit = std::ceil(limit - fj_cpu.problem->tolerances.integrality_tolerance); + if (limit > upper[var]) return false; + if (limit <= lower[var] + commit_threshold) return false; + lower[var] = limit; + return true; +} + +template +bool tighten_upper_bound(fj_cpu_climber_t& fj_cpu, + const std::vector& lower, + std::vector& upper, + i_t var, + f_t limit, + f_t commit_threshold) +{ + if (!std::isfinite(limit)) return false; + if (is_integer_var(fj_cpu, var)) + limit = std::floor(limit + fj_cpu.problem->tolerances.integrality_tolerance); + if (limit < lower[var]) return false; + if (limit >= upper[var] - commit_threshold) return false; + upper[var] = limit; + return true; +} + +// a light bounds propagation phase that runs much faster than the full scale presolve +// really helps on some instances. +template +void apply_bound_propagation(fj_cpu_climber_t& fj_cpu) +{ + if (!fj_cpu.use_bound_prop) return; + CPUFJ_NVTX_RANGE("CPUFJ::apply_bound_propagation"); + phase_timer_t timer(fj_cpu.t_bound_prop); + + const i_t n_variables = fj_cpu.problem->n_variables; + const i_t n_constraints = fj_cpu.problem->n_constraints; + const f_t commit = + (f_t)fj_cpu.hp.bound_prop_commit_scale * fj_cpu.problem->tolerances.absolute_tolerance; + + std::vector lower(n_variables); + std::vector upper(n_variables); + for (i_t var = 0; var < n_variables; ++var) { + auto bounds = fj_cpu.h_var_bounds[var].get(); + lower[var] = get_lower(bounds); + upper[var] = get_upper(bounds); + } + + bool changed = true; + int32_t pass = 0; + for (; changed && pass < fj_cpu.hp.bound_prop_rounds; ++pass) { + changed = false; + for (i_t row = 0; row < n_constraints; ++row) { + const f_t row_lb = fj_cpu.problem->cstr_lb[row]; + const f_t row_ub = fj_cpu.problem->cstr_ub[row]; + const bool has_lb = std::isfinite(row_lb); + const bool has_ub = std::isfinite(row_ub); + if (!has_lb && !has_ub) continue; + + const i_t begin = fj_cpu.problem->offsets[row]; + const i_t end = fj_cpu.problem->offsets[row + 1]; + + f_t min_activity = 0; + f_t max_activity = 0; + bool finite_min = true; + bool finite_max = true; + for (i_t p = begin; p < end; ++p) { + const f_t coeff = fj_cpu.problem->coefficients[p]; + if (coeff == f_t{0}) continue; + const i_t var = fj_cpu.problem->variables[p]; + const f_t min_x = coeff > 0 ? lower[var] : upper[var]; + const f_t max_x = coeff > 0 ? upper[var] : lower[var]; + finite_min &= std::isfinite(min_x); + finite_max &= std::isfinite(max_x); + } + const auto indices = thrust::make_counting_iterator(begin); + if (finite_min) { + const auto min_values = thrust::make_transform_iterator(indices, [&](i_t p) { + const f_t coeff = fj_cpu.problem->coefficients[p]; + if (coeff == f_t{0}) return f_t{0}; + const i_t var = fj_cpu.problem->variables[p]; + return coeff > 0 ? lower[var] : upper[var]; + }); + min_activity = + compensated_dot2(fj_cpu.problem->coefficients.data() + begin, min_values, end - begin); + } + if (finite_max) { + const auto max_values = thrust::make_transform_iterator(indices, [&](i_t p) { + const f_t coeff = fj_cpu.problem->coefficients[p]; + if (coeff == f_t{0}) return f_t{0}; + const i_t var = fj_cpu.problem->variables[p]; + return coeff > 0 ? upper[var] : lower[var]; + }); + max_activity = + compensated_dot2(fj_cpu.problem->coefficients.data() + begin, max_values, end - begin); + } + + const bool from_row_ub = finite_min && has_ub; + const bool from_row_lb = finite_max && has_lb; + if (!from_row_ub && !from_row_lb) continue; + + // The activities are not refreshed as the loop below narrows the row's own variables, and a + // stale bound is the looser one, so a deduction taken against it is the weaker one. + for (i_t p = begin; p < end; ++p) { + const f_t coeff = fj_cpu.problem->coefficients[p]; + if (coeff == f_t{0}) continue; + const i_t var = fj_cpu.problem->variables[p]; + + if (from_row_ub) { + const f_t rest = min_activity - coeff * (coeff > 0 ? lower[var] : upper[var]); + const f_t limit = (row_ub - rest) / coeff; + changed |= coeff > 0 ? tighten_upper_bound(fj_cpu, lower, upper, var, limit, commit) + : tighten_lower_bound(fj_cpu, lower, upper, var, limit, commit); + } + if (from_row_lb) { + const f_t rest = max_activity - coeff * (coeff > 0 ? upper[var] : lower[var]); + const f_t limit = (row_lb - rest) / coeff; + changed |= coeff > 0 ? tighten_lower_bound(fj_cpu, lower, upper, var, limit, commit) + : tighten_upper_bound(fj_cpu, lower, upper, var, limit, commit); + } + } + } + } + + fj_cpu.h_binary_indices.clear(); + fj_cpu.n_binary_vars = 0; + fj_cpu.n_integer_vars = 0; + [[maybe_unused]] i_t tightened = 0; + bool clamped = false; + for (i_t var = 0; var < n_variables; ++var) { + auto bounds = fj_cpu.h_var_bounds[var].get(); + cuopt_assert(!(lower[var] < get_lower(bounds)), "propagation widened a lower bound"); + cuopt_assert(!(upper[var] > get_upper(bounds)), "propagation widened an upper bound"); + cuopt_assert(!(lower[var] > upper[var]), "propagation emptied a domain"); + const bool moved = lower[var] != get_lower(bounds) || upper[var] != get_upper(bounds); + + // Same rule as problem_t::compute_binary_var_table, fixed binaries included: a domain narrowed + // to a point is no longer binary. + const bool integer = is_integer_var(fj_cpu, var); + const bool binary = integer && fj_cpu.problem->integer_equal(lower[var], (f_t)0) && + fj_cpu.problem->integer_equal(upper[var], (f_t)1); + fj_cpu.h_is_binary_variable[var] = binary; + if (binary) { + fj_cpu.h_binary_indices.push_back(var); + ++fj_cpu.n_binary_vars; + } else if (integer) { + ++fj_cpu.n_integer_vars; + } + if (!moved) continue; + + ++tightened; + fj_cpu.h_var_bounds[var] = typename type_2::type{lower[var], upper[var]}; + + const f_t value = fj_cpu.h_assignment[var]; + const f_t clamped_value = std::clamp(value, lower[var], upper[var]); + if (clamped_value != value) { + cuopt_assert(!integer || fj_cpu.problem->is_integer(clamped_value), + "bound clamp broke integrality"); + fj_cpu.h_assignment[var] = clamped_value; + clamped = true; + } + fj_cpu.h_best_assignment[var] = + std::clamp((f_t)fj_cpu.h_best_assignment[var], lower[var], upper[var]); + } + + if (clamped) recompute_lhs(fj_cpu); + cuopt_func_call(audit_assignment_bounds(fj_cpu, "bound prop")); + + CUOPT_LOG_DEBUG("%sCPUFJ bound prop: %d passes, %d domains tightened, %d binary of %d integer", + fj_cpu.log_prefix.c_str(), + pass, + tightened, + fj_cpu.n_binary_vars, + fj_cpu.n_binary_vars + fj_cpu.n_integer_vars); +} + +template +void apply_lock_weighted_start(fj_cpu_climber_t& fj_cpu) +{ + if (fj_cpu.problem->nnz > fj_cpu.hp.start_nnz_limit) return; + + const i_t n_variables = fj_cpu.problem->n_variables; + for (i_t var_idx = 0; var_idx < n_variables; ++var_idx) { + const f_t lb = get_lower(fj_cpu.h_var_bounds[var_idx].get()); + const f_t ub = get_upper(fj_cpu.h_var_bounds[var_idx].get()); + if (!std::isfinite(lb) || !std::isfinite(ub) || lb >= ub) continue; + + i_t lock_up = 0; + i_t lock_down = 0; + const auto range = model_range_for_var(fj_cpu, var_idx); + for (i_t i = range.first; i < range.second; ++i) { + const f_t coeff = fj_cpu.problem->reverse_coefficients[i]; + const i_t cstr_idx = fj_cpu.problem->reverse_constraints[i]; + const bool has_lb = std::isfinite((f_t)fj_cpu.problem->cstr_lb[cstr_idx]); + const bool has_ub = std::isfinite((f_t)fj_cpu.problem->cstr_ub[cstr_idx]); + if (coeff > 0) { + lock_up += has_ub; + lock_down += has_lb; + } else if (coeff < 0) { + lock_up += has_lb; + lock_down += has_ub; + } + } + + f_t new_val = lock_up <= lock_down ? ub : lb; + if (is_integer_var(fj_cpu, var_idx)) new_val = std::round(new_val); + fj_cpu.h_assignment[var_idx] = new_val; + } + + recompute_lhs(fj_cpu); + fj_cpu.h_best_assignment = fj_cpu.h_assignment; +} + +template +void apply_ambiguous_lock_start(fj_cpu_climber_t& fj_cpu) +{ + if (fj_cpu.problem->nnz > fj_cpu.hp.start_nnz_limit) return; + + std::mt19937 rng((uint32_t)fj_cpu.settings.seed ^ 0x9e3779b9u); + std::bernoulli_distribution flip(0.5); + for (i_t var = 0; var < fj_cpu.problem->n_variables; ++var) { + const f_t lower = get_lower(fj_cpu.h_var_bounds[var].get()); + const f_t upper = get_upper(fj_cpu.h_var_bounds[var].get()); + if (!std::isfinite(lower) || !std::isfinite(upper) || lower >= upper) continue; + + i_t up = 0, down = 0; + const auto [begin, end] = model_range_for_var(fj_cpu, var); + for (i_t p = begin; p < end; ++p) { + const f_t coeff = fj_cpu.problem->reverse_coefficients[p]; + const i_t row = fj_cpu.problem->reverse_constraints[p]; + const bool lb = std::isfinite((f_t)fj_cpu.problem->cstr_lb[row]); + const bool ub = std::isfinite((f_t)fj_cpu.problem->cstr_ub[row]); + if (coeff > 0) { + up += ub; + down += lb; + } else if (coeff < 0) { + up += lb; + down += ub; + } + } + + const bool ambiguous = std::abs(up - down) <= 1; + const bool choose_up = up < down || (ambiguous && flip(rng)); + f_t value = choose_up ? upper : lower; + if (is_integer_var(fj_cpu, var)) value = std::round(value); + fj_cpu.h_assignment[var] = value; + } + recompute_lhs(fj_cpu); + fj_cpu.h_best_assignment = fj_cpu.h_assignment; +} + +#if MIP_INSTANTIATE_FLOAT +template void cap_integer_domains(fj_cpu_climber_t&, int); +template void clamp_start_magnitude(fj_cpu_climber_t&, int); +template void apply_bound_propagation(fj_cpu_climber_t&); +template void apply_lock_weighted_start(fj_cpu_climber_t&); +template void apply_ambiguous_lock_start(fj_cpu_climber_t&); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void cap_integer_domains(fj_cpu_climber_t&, int); +template void clamp_start_magnitude(fj_cpu_climber_t&, int); +template void apply_bound_propagation(fj_cpu_climber_t&); +template void apply_lock_weighted_start(fj_cpu_climber_t&); +template void apply_ambiguous_lock_start(fj_cpu_climber_t&); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/bounds.hpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/bounds.hpp new file mode 100644 index 0000000000..55ea7f6cc9 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/bounds.hpp @@ -0,0 +1,18 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +#pragma once +#include "../state.hpp" +namespace cuopt::mathematical_optimization::mip { +template +void cap_integer_domains(fj_cpu_climber_t&, i_t); +template +void clamp_start_magnitude(fj_cpu_climber_t&, i_t); +template +void apply_bound_propagation(fj_cpu_climber_t&); +template +void apply_lock_weighted_start(fj_cpu_climber_t&); +template +void apply_ambiguous_lock_start(fj_cpu_climber_t&); +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.cpp index 37911cc9fc..ab67c2e699 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.cpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.cpp @@ -71,12 +71,366 @@ void eliminate_slacks(const lp_problem_t& problem, csr_A.n = n_structural; } +template +bool solve_lp_relaxation(const simplex::user_problem_t& relaxation, + double time_limit, + std::vector& x, + double& lp_seconds) +{ + simplex::lp_status_t status = simplex::lp_status_t::UNSET; + double seconds = 0; + + // solve_linear_program_advanced, whose status separates a limit -- which leaves a usable vertex + // behind -- from infeasibility. Guarded on f_t because dual simplex is only built for double. + if constexpr (std::is_same_v) { + simplex_solver_settings_t lp_settings; + lp_settings.relaxation = true; + lp_settings.time_limit = time_limit; + lp_settings.log.log = false; + // The portfolio already pins one CPU per lane, and the simplex default is + // omp_get_max_threads() - 1, which would open a second portfolio inside this lane's worker. + lp_settings.num_threads = 1; + + const f_t lp_start = tic(); + simplex::lp_solution_t lp_solution(relaxation.num_rows, relaxation.num_cols); + status = simplex::solve_linear_program(relaxation, lp_settings, lp_start, lp_solution); + x = std::move(lp_solution.x); + seconds = toc(lp_start); + } + lp_seconds += seconds; + + const bool usable = + status == simplex::lp_status_t::OPTIMAL || status == simplex::lp_status_t::TIME_LIMIT || + status == simplex::lp_status_t::ITERATION_LIMIT || + status == simplex::lp_status_t::CONCURRENT_LIMIT || status == simplex::lp_status_t::WORK_LIMIT; + CUOPT_LOG_DEBUG("CPUFJ LP relaxation: %s after %.3fs of %.3fs%s", + simplex::lp_status_to_string(status).c_str(), + seconds, + time_limit, + usable ? "" : ", discarded"); + return usable; +} + +template +static simplex::user_problem_t make_lp_distance_problem( + const simplex::user_problem_t& base, + fj_cpu_climber_t& fj_cpu, + const std::vector& rounded) +{ + std::vector integer_vars; + for (i_t var = 0; var < fj_cpu.problem->n_variables; ++var) + if (is_integer_var(fj_cpu, var)) integer_vars.push_back(var); + const i_t n_distance = (i_t)integer_vars.size(); + + simplex::user_problem_t result(base.handle_ptr); + result.num_rows = base.num_rows + 2 * n_distance; + result.num_cols = base.num_cols + n_distance; + + // The model's own objective is dropped: this LP measures distance alone. + result.objective.assign(result.num_cols, f_t{0}); + for (i_t k = 0; k < n_distance; ++k) + result.objective[base.num_cols + k] = f_t{1}; + + result.lower = base.lower; + result.upper = base.upper; + result.lower.resize(result.num_cols, f_t{0}); + result.upper.resize(result.num_cols, std::numeric_limits::infinity()); + + result.rhs = base.rhs; + result.row_sense = base.row_sense; + result.rhs.reserve(result.num_rows); + result.row_sense.reserve(result.num_rows); + for (i_t k = 0; k < n_distance; ++k) { + result.rhs.push_back(rounded[integer_vars[k]]); + result.row_sense.push_back('L'); + result.rhs.push_back(-rounded[integer_vars[k]]); + result.row_sense.push_back('L'); + } + result.range_rows = base.range_rows; + result.range_value = base.range_value; + result.num_range_rows = base.num_range_rows; + + const i_t base_nnz = base.A.col_start[base.A.n]; + csc_matrix_t matrix(result.num_rows, result.num_cols, base_nnz + 4 * n_distance); + i_t out = 0; + i_t next_integer = 0; + for (i_t j = 0; j < base.num_cols; ++j) { + matrix.col_start[j] = out; + for (i_t p = base.A.col_start[j]; p < base.A.col_start[j + 1]; ++p) { + matrix.i[out] = base.A.i[p]; + matrix.x[out++] = base.A.x[p]; + } + if (next_integer < n_distance && integer_vars[next_integer] == j) { + const i_t row = base.num_rows + 2 * next_integer++; + matrix.i[out] = row; + matrix.x[out++] = f_t{1}; + matrix.i[out] = row + 1; + matrix.x[out++] = f_t{-1}; + } + } + for (i_t k = 0; k < n_distance; ++k) { + matrix.col_start[base.num_cols + k] = out; + const i_t row = base.num_rows + 2 * k; + matrix.i[out] = row; + matrix.x[out++] = f_t{-1}; + matrix.i[out] = row + 1; + matrix.x[out++] = f_t{-1}; + } + matrix.col_start[result.num_cols] = out; + cuopt_assert(out == base_nnz + 4 * n_distance, "distance problem nonzero count mismatch"); + result.A = std::move(matrix); + return result; +} + +template +static simplex::user_problem_t make_fixed_integer_lp( + const simplex::user_problem_t& base, + fj_cpu_climber_t& fj_cpu, + const std::vector& rounded) +{ + simplex::user_problem_t fixed(base.handle_ptr); + fixed.num_rows = base.num_rows; + fixed.num_cols = base.num_cols; + fixed.objective.assign(base.num_cols, f_t{0}); + fixed.rhs = base.rhs; + fixed.row_sense = base.row_sense; + fixed.range_rows = base.range_rows; + fixed.range_value = base.range_value; + fixed.num_range_rows = base.num_range_rows; + fixed.A = base.A; + fixed.lower = base.lower; + fixed.upper = base.upper; + for (i_t var = 0; var < (i_t)rounded.size(); ++var) { + if (!is_integer_var(fj_cpu, var)) continue; + fixed.lower[var] = fixed.upper[var] = rounded[var]; + } + return fixed; +} + +template +void apply_lp_rounded_start(fj_cpu_climber_t& fj_cpu, f_t lane_time_limit) +{ + if (!fj_cpu.use_lp_start || !fj_cpu.problem->host_lp) return; + if (fj_cpu.problem->nnz > fj_cpu.hp.lp_start_nnz_limit) return; + + // In a nonnegative all-integer equality system, flooring a feasible relaxation is a particularly + // useful FJ start: it cannot overshoot any equality, so the residual search only has to fill + // deficits instead of simultaneously undoing stochastic round-ups. Give one feasibility-LP + // persona enough time to obtain a real basic solution on this broad structural class. The + // ordinary LP personas retain their tiny opportunistic budget, and deep pumps retain theirs. + bool monotone_integer_equalities = + fj_cpu.lp_start_feasibility_objective && fj_cpu.problem->equality_fraction == 1.0 && + fj_cpu.n_integer_vars + fj_cpu.n_binary_vars == fj_cpu.problem->n_variables; + if (monotone_integer_equalities) { + for (i_t var = 0; var < fj_cpu.problem->n_variables; ++var) { + const auto bounds = fj_cpu.h_var_bounds[var].get(); + if (get_lower(bounds) < f_t{0} || fj_cpu.problem->h_obj_coeffs[var] != f_t{0}) { + monotone_integer_equalities = false; + break; + } + } + } + if (monotone_integer_equalities) { + for (i_t p = 0; p < fj_cpu.problem->nnz; ++p) { + if (fj_cpu.problem->coefficients[p] < f_t{0}) { + monotone_integer_equalities = false; + break; + } + } + } + + const double budget = fj_cpu.use_deep_lp_pump ? std::min(3.5, 0.7 * (double)lane_time_limit) + : monotone_integer_equalities + ? std::min(2.0, 0.45 * (double)lane_time_limit) + : std::min(fj_cpu.hp.lp_pump_max_budget_s, + fj_cpu.hp.lp_pump_budget_share * (double)lane_time_limit); + if (budget <= 0) return; + + CPUFJ_NVTX_RANGE("CPUFJ::apply_lp_rounded_start"); + phase_timer_t timer(fj_cpu.t_lp_start); + + simplex::user_problem_t base = *fj_cpu.problem->host_lp; + if (fj_cpu.lp_start_feasibility_objective) + std::fill(base.objective.begin(), base.objective.end(), f_t{0}); + + // A zero-objective relaxation returns the same arbitrary basic solution in every feasibility-LP + // lane. On a nonnegative integer equality master, all of those lanes then floor the same vertex + // and spend most of the budget repairing the same residual. Positive random costs keep the LP + // bounded and feasible while selecting independent vertices for the portfolio. This is gated by + // the certificate above; ordinary objective-bearing and mixed-sign models are unchanged. + if (monotone_integer_equalities) { + cuopt::pcgenerator_t objective_rng(fj_cpu.settings.seed ^ 0xd1b54a32d192ed03ULL); + for (i_t var = 0; var < fj_cpu.problem->n_variables; ++var) + base.objective[var] = f_t{1} + (f_t)objective_rng.next_double(); + } + + const auto started = std::chrono::steady_clock::now(); + const i_t n_variables = fj_cpu.problem->n_variables; + + std::vector rounded; + std::vector selected; + // Keep the least-infeasible rounded LP projection as the FJ starting point. + // total_violations sums negative excesses, so the greatest value is the least infeasible. + f_t selected_violation = -std::numeric_limits::infinity(); + + const int32_t projections = fj_cpu.use_deep_lp_pump ? 100 : fj_cpu.hp.lp_pump_projections; + for (int32_t projection = 0; projection < projections; ++projection) { + const double remaining = + budget - std::chrono::duration(std::chrono::steady_clock::now() - started).count(); + if (remaining <= 0) break; + + // Projection 0 is the plain relaxation; the rest chase the previous rounding. + const auto distance = projection == 0 ? simplex::user_problem_t(base.handle_ptr) + : make_lp_distance_problem(base, fj_cpu, rounded); + const auto& relaxation = projection == 0 ? base : distance; + + std::vector x; + if (!solve_lp_relaxation(relaxation, remaining, x, fj_cpu.t_lp_relaxation)) break; + // convert_user_problem appends slacks, so the model's own variables are the leading columns. + if ((i_t)x.size() < n_variables) break; + + rounded.resize(n_variables); + cuopt::pcgenerator_t rng(fj_cpu.settings.seed + 0x9e3779b9ULL * (uint64_t)projection); + bool valid = true; + for (i_t var = 0; var < n_variables && valid; ++var) { + const auto bounds = fj_cpu.h_var_bounds[var].get(); + const f_t lower = get_lower(bounds); + const f_t upper = get_upper(bounds); + f_t value = std::clamp(x[var], lower, upper); + if (!std::isfinite(value)) { + valid = false; + break; + } + if (is_integer_var(fj_cpu, var)) { + if (monotone_integer_equalities) { + // Every coefficient is nonnegative, hence this preserves every equality's upper side. + // Clamp once more because an LP value can sit a few ulps below an integral lower bound. + value = std::clamp(std::floor(value), std::ceil(lower), std::floor(upper)); + } else { + // Rounded up with probability equal to the fractional part, so successive projections of + // the same point explore different corners. + const f_t fraction = value - std::floor(value); + value = rng.next_double() < fraction ? std::ceil(value) : std::floor(value); + } + // A variable with no integral value inside its bounds cannot form a valid start without + // breaking the engine's integrality invariant. + valid = value >= lower && value <= upper; + } + rounded[var] = value; + } + if (!valid) break; + + std::vector candidate = rounded; + if (fj_cpu.use_deep_lp_pump) { + const double repair_budget = + budget - std::chrono::duration(std::chrono::steady_clock::now() - started).count(); + if (repair_budget > 0.01) { + auto fixed = make_fixed_integer_lp(base, fj_cpu, rounded); + std::vector repaired; + if (solve_lp_relaxation(fixed, repair_budget, repaired, fj_cpu.t_lp_relaxation) && + (i_t)repaired.size() >= n_variables) { + for (i_t var = 0; var < n_variables; ++var) + if (!is_integer_var(fj_cpu, var)) candidate[var] = repaired[var]; + } + } + } + std::copy(candidate.begin(), candidate.end(), fj_cpu.h_assignment.begin()); + recompute_lhs(fj_cpu); + cuopt_assert(fj_cpu.total_violations <= f_t{0}, "total_violations should be nonpositive"); + if (fj_cpu.total_violations > selected_violation) { + selected_violation = fj_cpu.total_violations; + selected = candidate; + } + + // The rounded point can already be integral-feasible. It never passed through apply_move, so + // the incumbent is recorded here through the same contract that path uses. + if (fj_cpu.violated_constraints.empty() && check_variable_feasibility(fj_cpu)) { + std::copy(candidate.begin(), candidate.end(), fj_cpu.h_best_assignment.begin()); + fj_cpu.h_best_objective = + fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; + fj_cpu.feasible_found = true; + report_cpu_incumbent(fj_cpu); + return; + } + } + + if (selected.empty()) return; + std::copy(selected.begin(), selected.end(), fj_cpu.h_assignment.begin()); + std::copy(selected.begin(), selected.end(), fj_cpu.h_best_assignment.begin()); + recompute_lhs(fj_cpu); + cuopt_func_call(audit_assignment_bounds(fj_cpu, "lp pump")); +} + +template +bool apply_lp_polish(fj_cpu_climber_t& fj_cpu, double budget_s) +{ + if (!fj_cpu.problem->host_lp || fj_cpu.problem->nnz > fj_cpu.hp.lp_polish_nnz_limit || + budget_s < fj_cpu.hp.lp_polish_min_budget_s) + return false; + + const i_t n = fj_cpu.problem->n_variables; + std::vector incumbent(fj_cpu.h_best_assignment.begin(), + fj_cpu.h_best_assignment.begin() + n); + if (fj_cpu.shared_incumbent) + fj_cpu.shared_incumbent->adopt(fj_cpu.h_best_objective + f_t{1}, incumbent); + + bool has_continuous = false; + for (i_t var = 0; var < n; ++var) + has_continuous |= !is_integer_var(fj_cpu, var); + if (!has_continuous) return false; + + phase_timer_t timer(fj_cpu.t_lp_start); + simplex::user_problem_t relaxation = *fj_cpu.problem->host_lp; + if ((i_t)relaxation.lower.size() < n || (i_t)relaxation.upper.size() < n) return false; + + for (i_t var = 0; var < n; ++var) { + if (!is_integer_var(fj_cpu, var)) continue; + relaxation.lower[var] = relaxation.upper[var] = std::round(incumbent[var]); + } + + std::vector x; + if (!solve_lp_relaxation(relaxation, budget_s, x, fj_cpu.t_lp_relaxation) || (i_t)x.size() < n) + return false; + + std::vector completion(n); + for (i_t var = 0; var < n; ++var) { + const auto bounds = fj_cpu.h_var_bounds[var].get(); + const f_t value = is_integer_var(fj_cpu, var) ? std::round(incumbent[var]) : x[var]; + if (!std::isfinite(value)) return false; + completion[var] = std::clamp(value, get_lower(bounds), get_upper(bounds)); + } + + std::copy(completion.begin(), completion.end(), fj_cpu.h_assignment.begin()); + recompute_slack(fj_cpu); + const bool improved = fj_cpu.h_incumbent_objective < fj_cpu.h_best_objective && + fj_cpu.violated_constraints.empty() && + check_variable_feasibility(fj_cpu); + if (!improved) { + std::copy(fj_cpu.h_best_assignment.begin(), + fj_cpu.h_best_assignment.begin() + n, + fj_cpu.h_assignment.begin()); + recompute_slack(fj_cpu); + return false; + } + + std::copy(completion.begin(), completion.end(), fj_cpu.h_best_assignment.begin()); + fj_cpu.h_best_objective = + fj_cpu.h_incumbent_objective - fj_cpu.settings.parameters.breakthrough_move_epsilon; + fj_cpu.iterations_since_best = 0; + fj_cpu.perturb_streak = 0; + fj_cpu.feasible_found = true; + report_cpu_incumbent(fj_cpu); + return true; +} + #if MIP_INSTANTIATE_FLOAT template void eliminate_slacks(const lp_problem_t&, int, csr_matrix_t&, std::vector&, std::vector&); +template void apply_lp_rounded_start(fj_cpu_climber_t&, float); +template bool apply_lp_polish(fj_cpu_climber_t&, double); #endif #if MIP_INSTANTIATE_DOUBLE @@ -85,6 +439,8 @@ template void eliminate_slacks(const lp_problem_t&, csr_matrix_t&, std::vector&, std::vector&); +template void apply_lp_rounded_start(fj_cpu_climber_t&, double); +template bool apply_lp_polish(fj_cpu_climber_t&, double); #endif } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.hpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.hpp index d2dedc08c2..5fb43806fe 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.hpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/lp.hpp @@ -27,4 +27,8 @@ void eliminate_slacks(const simplex::lp_problem_t&, std::vector&, std::vector&); +template +void apply_lp_rounded_start(fj_cpu_climber_t&, f_t); +template +bool apply_lp_polish(fj_cpu_climber_t&, double); } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.cpp index 1ee1173118..71ba3c81c3 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.cpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.cpp @@ -151,6 +151,158 @@ void detect_implied_integers(fj_cpu_climber_t& fj_cpu, "CPUFJ implied integrality: %d forced, %d integral-completable", n_forced, n_completable); } +template +void precompute_problem_features(fj_cpu_climber_t& fj_cpu, + fj_cpu_problem_t& problem) +{ + phase_timer_t timer(fj_cpu.t_features); + fj_cpu.n_binary_vars = 0; + fj_cpu.n_integer_vars = 0; + for (i_t i = 0; i < (i_t)fj_cpu.h_is_binary_variable.size(); i++) { + if (fj_cpu.h_is_binary_variable[i]) { + fj_cpu.n_binary_vars++; + } else if (problem.h_var_types[i] == var_t::INTEGER) { + fj_cpu.n_integer_vars++; + } + } + + i_t total_nnz = problem.reverse_offsets.back(); + i_t n_vars = problem.reverse_offsets.size() - 1; + i_t n_cstrs = problem.offsets.size() - 1; + + problem.avg_var_degree = (double)total_nnz / n_vars; + + problem.max_var_degree = 0; + std::vector var_degrees(n_vars); + for (i_t i = 0; i < n_vars; i++) { + i_t degree = problem.reverse_offsets[i + 1] - problem.reverse_offsets[i]; + var_degrees[i] = degree; + problem.max_var_degree = std::max(problem.max_var_degree, degree); + } + + i_t equalities = 0; + for (i_t row = 0; row < n_cstrs; ++row) + equalities += problem.cstr_lb[row] == problem.cstr_ub[row]; + problem.equality_fraction = n_cstrs ? (double)equalities / n_cstrs : 0.0; +} + +template +void detect_free_equality_singletons(fj_cpu_climber_t& c) +{ + const auto& model = *c.problem; + const i_t n = model.reverse_offsets.size() - 1, m = model.offsets.size() - 1; + c.bin_ignore_row.assign(m, 0); + c.bin_ignore_var.assign(n, 0); + for (i_t row = 0; row < m; ++row) { + const f_t rhs = model.cstr_lb[row]; + if (!std::isfinite(rhs) || rhs != model.cstr_ub[row]) continue; + typename fj_cpu_climber_t::bin_eliminated_row_t rec{row, rhs}; + // A single nonnegative continuous slack can be substituted without fill-in. + // Keep the original model intact; the encoder applies its bounds and cost. + i_t singleton = -1; + for (i_t entry = model.offsets[row]; entry < model.offsets[row + 1]; ++entry) { + const i_t var = model.variables[entry]; + if (model.h_var_types[var] != var_t::CONTINUOUS) continue; + const auto bounds = c.h_var_bounds[var].get(); + if (singleton >= 0 || model.reverse_offsets[var + 1] - model.reverse_offsets[var] != 1 || + std::abs(model.coefficients[entry]) != f_t{1} || get_lower(bounds) != f_t{0} || + std::isfinite(get_upper(bounds)) || !std::isfinite(model.h_obj_coeffs[var])) { + singleton = -1; + break; + } + singleton = var; + } + if (singleton >= 0) { + c.bin_singletons.emplace_back(row, singleton); + c.bin_ignore_var[singleton] = 1; + rec.all.push_back(singleton); + c.bin_eliminated_rows.push_back(std::move(rec)); + continue; + } + bool valid = true, raises = false, lowers = false; + for (i_t p = model.offsets[row]; p < model.offsets[row + 1]; ++p) { + const i_t var = model.variables[p]; + const f_t a = model.coefficients[p]; + if (!a || c.h_is_binary_variable[var]) continue; + if (c.problem->h_var_types[var] != var_t::CONTINUOUS || + model.reverse_offsets[var + 1] - model.reverse_offsets[var] != 1) { + valid = false; + break; + } + const auto bounds = c.h_var_bounds[var].get(); + const bool up = (a > 0 && !std::isfinite(get_upper(bounds))) || + (a < 0 && !std::isfinite(get_lower(bounds))); + const bool down = (a > 0 && !std::isfinite(get_lower(bounds))) || + (a < 0 && !std::isfinite(get_upper(bounds))); + if (!up && !down) { + valid = false; + break; + } + rec.all.push_back(var); + if (up) { + raises = true; + rec.positive.push_back(var); + rec.positive_coeff.push_back(a); + } + if (down) { + lowers = true; + rec.negative.push_back(var); + rec.negative_coeff.push_back(a); + } + } + if (!valid || rec.all.empty() || !raises || !lowers) continue; + c.bin_ignore_row[row] = 1; + for (i_t var : rec.all) + c.bin_ignore_var[var] = 1; + c.bin_eliminated_rows.push_back(std::move(rec)); + } + c.has_bin_elimination = !c.bin_eliminated_rows.empty(); +} + +template +void build_cardinality_index(fj_cpu_climber_t& c, fj_cpu_problem_t& problem) +{ + auto& offsets = problem.card_row_offsets; + auto& members = problem.card_variables; + auto& cardinalities = problem.card_cardinalities; + offsets.assign(1, 0); + members.clear(); + cardinalities.clear(); + problem.card_group_of_variable.assign(problem.n_variables, -1); + + for (i_t row = 0; row < problem.n_constraints; ++row) { + const f_t lb = problem.cstr_lb[row], ub = problem.cstr_ub[row]; + const i_t begin = problem.offsets[row], end = problem.offsets[row + 1]; + if (!std::isfinite(lb) || !std::isfinite(ub) || std::abs(lb - ub) > 1e-6 || end - begin < 2 || + end - begin > 20000) + continue; + + const f_t common = problem.coefficients[begin]; + if (std::abs(common) <= 1e-6) continue; + const f_t cardinality = lb / common; + if (std::abs(cardinality - std::round(cardinality)) > 1e-6 || cardinality < 0 || + cardinality > end - begin) + continue; + + bool valid = true; + for (i_t p = begin; p < end && valid; ++p) { + const i_t var = problem.variables[p]; + valid = c.h_is_binary_variable[var] && std::abs(problem.coefficients[p] - common) <= + 1e-6 * std::max(1, std::abs(common)); + } + if (!valid) continue; + const i_t group = (i_t)offsets.size() - 1; + for (i_t p = begin; p < end; ++p) { + const i_t var = problem.variables[p]; + members.push_back(var); + i_t& owner = problem.card_group_of_variable[var]; + owner = owner == -1 ? group : (owner == group ? group : -2); + } + offsets.push_back((i_t)members.size()); + cardinalities.push_back((i_t)std::llround(cardinality)); + } +} + template void certify_epigraph_variables(fj_cpu_climber_t& fj_cpu, i_t n_variables) { @@ -299,18 +451,333 @@ void build_one_sided_rows(fj_cpu_climber_t& fj_cpu) fj_cpu.release_setup_structures(); } +// Eliminate coordinates through exact equalities while retaining each pivot's domain as a row. +// Integer pivots are accepted only when divisibility proves that every lifted value stays integral. +// FJ usually struggles with equality-heavy models since every move may result in equality rows +// being violated and repair having to be applied to many other variables to "compensate". Rewriting +// the problem may help in some cases. +template +std::unique_ptr> make_equality_reduced_climber( + fj_cpu_climber_t& c, + double budget, + std::vector>& substitutions, + std::vector& retained) +{ + using term_t = std::pair; + const auto started = std::chrono::steady_clock::now(); + auto expired = [&] { + return c.preemption_flag.load(std::memory_order_relaxed) || + std::chrono::duration(std::chrono::steady_clock::now() - started).count() >= + budget; + }; + const auto& p = *c.problem; + std::vector equalities; + for (i_t r = 0; r < p.n_constraints; ++r) + if (std::isfinite(p.cstr_lb[r]) && p.cstr_lb[r] == p.cstr_ub[r]) equalities.push_back(r); + if (equalities.empty() || budget <= 0) return nullptr; + + std::vector> rows(p.n_constraints); + std::vector> incidence(p.n_variables); + std::vector lower = p.cstr_lb, upper = p.cstr_ub, objective = p.h_obj_coeffs; + f_t objective_offset = 0; + std::vector eliminated(p.n_variables, 0); + int64_t nnz = 0; + for (i_t r = 0; r < p.n_constraints; ++r) { + auto& row = rows[r]; + for (i_t k = p.offsets[r]; k < p.offsets[r + 1]; ++k) + row.emplace_back(p.variables[k], p.coefficients[k]); + std::sort(row.begin(), row.end()); + size_t out = 0; + for (size_t k = 0; k < row.size();) { + const i_t v = row[k].first; + f_t a = 0; + do { + a += row[k++].second; + } while (k < row.size() && row[k].first == v); + if (a != 0) { + row[out++] = {v, a}; + incidence[v].push_back(r); + } + } + row.resize(out); + nnz += out; + } + std::stable_sort(equalities.begin(), equalities.end(), [&](i_t a, i_t b) { + return rows[a].size() < rows[b].size(); + }); + + struct staged_row_t { + i_t index; + std::vector terms; + f_t lower; + f_t upper; + }; + // Sparse binary scheduling models can encode resource usage as a chain of unit + // differences. Partial elimination leaves the same equality barrier in place; + // allow enough fill to expose the cumulative capacity rows on this class only. + i_t chain_rows = 0; + if (c.n_binary_vars > p.n_variables / 2 && p.max_var_degree <= 4 && + equalities.size() > 0.75 * p.n_constraints) { + for (i_t r : equalities) { + if (lower[r] != f_t{0}) continue; + i_t count = 0; + f_t sum = 0; + bool unit = true; + for (const auto& [v, a] : rows[r]) { + if (c.h_is_binary_variable[v]) continue; + ++count; + sum += a; + const auto bounds = c.h_var_bounds[v].get(); + unit &= std::fabs(a) == f_t{1} && objective[v] == f_t{0} && incidence[v].size() <= 2 && + std::isfinite(get_lower(bounds)) && std::isfinite(get_upper(bounds)); + } + if (unit && count == 2 && sum == f_t{0}) ++chain_rows; + } + } + const bool cumulative_chain = + chain_rows >= 16 && 2 * chain_rows >= p.n_variables - c.n_binary_vars; + const int64_t fill_limit = (cumulative_chain ? 8 : 2) * std::max(1, nnz); + for (i_t r : equalities) { + if (expired()) break; + const auto& equation = rows[r]; + bool integer_row = std::isfinite(lower[r]) && lower[r] == std::round(lower[r]) && + std::fabs(lower[r]) <= f_t{1e12}; + int64_t row_gcd = integer_row ? (int64_t)std::fabs(lower[r]) : 0; + if (integer_row) { + for (const auto& [v, a] : equation) { + if (p.h_var_types[v] != var_t::INTEGER || !std::isfinite(a) || a != std::round(a) || + std::fabs(a) > f_t{1e12}) { + integer_row = false; + break; + } + row_gcd = std::gcd(row_gcd, (int64_t)std::fabs(a)); + } + } + + i_t pivot = -1; + f_t divisor = 0; + uint64_t best_work = std::numeric_limits::max(); + for (const auto& [v, a] : equation) { + const auto bounds = c.h_var_bounds[v].get(); + if (eliminated[v] || c.h_is_binary_variable[v] || get_lower(bounds) == get_upper(bounds)) + continue; + bool admissible = std::fabs(a) == f_t{1}; + if (p.h_var_types[v] == var_t::INTEGER) + admissible = integer_row && row_gcd > 0 && std::fabs(a) == (f_t)row_gcd; + if (!admissible) continue; + const uint64_t work = (uint64_t)incidence[v].size() * (equation.size() - 1); + if (work < best_work) { + pivot = v; + divisor = a; + best_work = work; + } + } + if (pivot < 0) continue; + + fj_equality_substitution_t sub{pivot, lower[r] / divisor, {}}; + if (!std::isfinite(sub.constant)) continue; + for (const auto& [v, a] : equation) + if (v != pivot) sub.terms.emplace_back(v, -a / divisor); + const auto domain = c.h_var_bounds[pivot].get(); + const f_t bound_lower = get_lower(domain) - sub.constant; + const f_t bound_upper = get_upper(domain) - sub.constant; + if ((std::isfinite(get_lower(domain)) && !std::isfinite(bound_lower)) || + (std::isfinite(get_upper(domain)) && !std::isfinite(bound_upper))) + continue; + + auto affected = incidence[pivot]; + std::sort(affected.begin(), affected.end()); + affected.erase(std::unique(affected.begin(), affected.end()), affected.end()); + std::vector staged; + int64_t next_nnz = nnz - (int64_t)equation.size() + (int64_t)sub.terms.size(); + bool rejected = false; + for (i_t row_index : affected) { + if (row_index == r) continue; + if (expired()) { + rejected = true; + break; + } + const auto& old = rows[row_index]; + auto entry = std::lower_bound( + old.begin(), old.end(), pivot, [](const term_t& t, i_t v) { return t.first < v; }); + if (entry == old.end() || entry->first != pivot) continue; + const f_t factor = entry->second; + staged_row_t replacement{row_index, + {}, + std::fma(-factor, sub.constant, lower[row_index]), + std::fma(-factor, sub.constant, upper[row_index])}; + if ((std::isfinite(lower[row_index]) && !std::isfinite(replacement.lower)) || + (std::isfinite(upper[row_index]) && !std::isfinite(replacement.upper))) { + rejected = true; + break; + } + size_t i = 0, j = 0; + replacement.terms.reserve(old.size() + sub.terms.size()); + while (i < old.size() || j < sub.terms.size()) { + if (i < old.size() && old[i].first == pivot) { + ++i; + continue; + } + const i_t v = std::min(i < old.size() ? old[i].first : p.n_variables, + j < sub.terms.size() ? sub.terms[j].first : p.n_variables); + f_t a = 0; + if (i < old.size() && old[i].first == v) a = old[i++].second; + if (j < sub.terms.size() && sub.terms[j].first == v) + a = std::fma(factor, sub.terms[j++].second, a); + if (!std::isfinite(a) || std::fabs(a) > f_t{1e12}) { + rejected = true; + break; + } + if (a != 0) replacement.terms.emplace_back(v, a); + } + if (rejected) break; + next_nnz += (int64_t)replacement.terms.size() - (int64_t)old.size(); + if (next_nnz > fill_limit) { + rejected = true; + break; + } + staged.push_back(std::move(replacement)); + } + if (rejected) continue; + + for (auto& replacement : staged) { + auto& old = rows[replacement.index]; + size_t k = 0; + for (const auto& [v, a] : replacement.terms) { + while (k < old.size() && old[k].first < v) + ++k; + if (k == old.size() || old[k].first != v) incidence[v].push_back(replacement.index); + } + old = std::move(replacement.terms); + lower[replacement.index] = replacement.lower; + upper[replacement.index] = replacement.upper; + } + rows[r] = sub.terms; + lower[r] = bound_lower; + upper[r] = bound_upper; + for (const auto& [v, a] : sub.terms) { + objective[v] = std::fma(objective[pivot], a, objective[v]); + if (!std::isfinite(objective[v])) return nullptr; + } + objective_offset = std::fma(objective[pivot], sub.constant, objective_offset); + if (!std::isfinite(objective_offset)) return nullptr; + objective[pivot] = 0; + eliminated[pivot] = 1; + incidence[pivot].clear(); + substitutions.push_back(std::move(sub)); + nnz = next_nnz; + } + if (substitutions.empty()) return nullptr; + + std::vector mapping(p.n_variables, -1), offsets{0}, variables; + std::vector coefficients, row_lower, row_upper, var_lower, var_upper, costs; + std::vector types; + for (i_t v = 0; v < p.n_variables; ++v) { + if (eliminated[v]) continue; + mapping[v] = retained.size(); + retained.push_back(v); + const auto bounds = c.h_var_bounds[v].get(); + var_lower.push_back(get_lower(bounds)); + var_upper.push_back(get_upper(bounds)); + costs.push_back(objective[v]); + types.push_back(p.h_var_types[v]); + } + for (i_t r = 0; r < p.n_constraints; ++r) { + if (!std::isfinite(lower[r]) && !std::isfinite(upper[r])) continue; + if (rows[r].empty()) { + if (lower[r] > 0 || upper[r] < 0) return nullptr; + continue; + } + for (const auto& [v, a] : rows[r]) { + cuopt_assert(mapping[v] >= 0, "eliminated column survived substitution"); + variables.push_back(mapping[v]); + coefficients.push_back(a); + } + offsets.push_back(variables.size()); + row_lower.push_back(lower[r]); + row_upper.push_back(upper[r]); + } + if (retained.empty() || row_lower.empty() || coefficients.size() > (size_t)INT32_MAX) + return nullptr; + + const i_t reduced_nnz = coefficients.size(); + const i_t reduced_rows = row_lower.size(); + auto child = init_fj_cpu_from_host_model((i_t)retained.size(), + reduced_rows, + reduced_nnz, + false, + f_t{1}, + objective_offset, + std::move(coefficients), + std::move(variables), + std::move(offsets), + std::move(costs), + std::move(var_lower), + std::move(var_upper), + std::move(row_lower), + std::move(row_upper), + {}, + {}, + std::move(types), + p.tolerances, + c.preemption_flag, + c.settings); + static_cast&>(*child) = static_cast&>(c); + child->use_equality_substitution = false; + child->use_lp_start = child->use_lp_polish = false; + child->use_bound_prop = true; + child->use_move_batching &= child->n_colors > 0; + child->log_prefix = c.log_prefix; + child->suppress_incumbent_log = true; + for (i_t j = 0; j < (i_t)retained.size(); ++j) { + const auto bounds = child->h_var_bounds[j].get(); + f_t value = c.h_assignment[retained[j]]; + if (is_integer_var(*child, j)) value = std::round(value); + child->h_assignment[j] = std::clamp(value, get_lower(bounds), get_upper(bounds)); + } + child->h_best_assignment = child->h_assignment; + recompute_lhs(*child); + CUOPT_LOG_DEBUG("%sCPUFJ equality substitution: %zu pivots, %d columns, %d rows, %d nonzeros", + c.log_prefix.c_str(), + substitutions.size(), + (int)retained.size(), + (int)reduced_rows, + (int)reduced_nnz); + return child; +} + #if MIP_INSTANTIATE_FLOAT template void detect_implied_integers(fj_cpu_climber_t&, fj_cpu_problem_t&); +template void detect_free_equality_singletons(fj_cpu_climber_t&); +template void precompute_problem_features(fj_cpu_climber_t&, + fj_cpu_problem_t&); +template void build_cardinality_index(fj_cpu_climber_t&, + fj_cpu_problem_t&); template void certify_epigraph_variables(fj_cpu_climber_t&, int); template void build_one_sided_rows(fj_cpu_climber_t&); +template std::unique_ptr> make_equality_reduced_climber( + fj_cpu_climber_t&, + double, + std::vector>&, + std::vector&); #endif #if MIP_INSTANTIATE_DOUBLE template void detect_implied_integers(fj_cpu_climber_t&, fj_cpu_problem_t&); +template void detect_free_equality_singletons(fj_cpu_climber_t&); +template void precompute_problem_features(fj_cpu_climber_t&, + fj_cpu_problem_t&); +template void build_cardinality_index(fj_cpu_climber_t&, + fj_cpu_problem_t&); template void certify_epigraph_variables(fj_cpu_climber_t&, int); template void build_one_sided_rows(fj_cpu_climber_t&); +template std::unique_ptr> make_equality_reduced_climber( + fj_cpu_climber_t&, + double, + std::vector>&, + std::vector&); #endif } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.hpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.hpp index 6870468cae..4579d4bff9 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.hpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/setup/structure.hpp @@ -6,10 +6,30 @@ #include "../state.hpp" namespace cuopt::mathematical_optimization::mip { +// An eliminated coordinate, recorded in original variable indices. Records are lifted in reverse. +template +struct fj_equality_substitution_t { + i_t variable; + f_t constant; + std::vector> terms; +}; + template void detect_implied_integers(fj_cpu_climber_t&, fj_cpu_problem_t&); template +void detect_free_equality_singletons(fj_cpu_climber_t&); +template +void precompute_problem_features(fj_cpu_climber_t&, fj_cpu_problem_t&); +template +void build_cardinality_index(fj_cpu_climber_t&, fj_cpu_problem_t&); +template void certify_epigraph_variables(fj_cpu_climber_t&, i_t); template void build_one_sided_rows(fj_cpu_climber_t&); +template +std::unique_ptr> make_equality_reduced_climber( + fj_cpu_climber_t&, + double, + std::vector>&, + std::vector&); } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/affine.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/affine.cpp new file mode 100644 index 0000000000..7b4a762ee8 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/affine.cpp @@ -0,0 +1,30 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "../setup/bounds.hpp" +#include "starts.hpp" + +namespace cuopt::mathematical_optimization::mip { + +template +void apply_structural_completion_start(fj_cpu_climber_t& fj_cpu) +{ + apply_lock_weighted_start(fj_cpu); + apply_exact_k_start(fj_cpu); + apply_greedy_covering_start(fj_cpu); + repair_difficult_anchor(fj_cpu); +} + +#if MIP_INSTANTIATE_FLOAT +template void apply_structural_completion_start(fj_cpu_climber_t&); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void apply_structural_completion_start(fj_cpu_climber_t&); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/cardinality.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/cardinality.cpp new file mode 100644 index 0000000000..ffad900515 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/cardinality.cpp @@ -0,0 +1,175 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "../internal.hpp" +#include "../problem.hpp" +#include "../search/api.hpp" +#include "starts.hpp" + +namespace cuopt::mathematical_optimization::mip { + +template +void apply_exact_k_start(fj_cpu_climber_t& fj_cpu) +{ + if (fj_cpu.problem->nnz > fj_cpu.hp.start_nnz_limit) return; + + const auto started = std::chrono::steady_clock::now(); + auto timed_out = [&] { + return std::chrono::duration(std::chrono::steady_clock::now() - started).count() > + fj_cpu.hp.exact_k_budget_s; + }; + + struct exact_k_row_t { + i_t k, begin, end; + }; + std::vector rows; + for (i_t row = 0; row < fj_cpu.problem->n_constraints; ++row) { + if ((row & 0xFFF) == 0 && timed_out()) return; + + const f_t lb = fj_cpu.problem->cstr_lb[row]; + const f_t ub = fj_cpu.problem->cstr_ub[row]; + if (!std::isfinite(lb) || !std::isfinite(ub) || std::abs(lb - ub) > fj_cpu.hp.exact_k_tol) + continue; + + const i_t begin = fj_cpu.problem->offsets[row]; + const i_t end = fj_cpu.problem->offsets[row + 1]; + if (end - begin < 2 || end - begin > fj_cpu.hp.exact_k_max_width) continue; + + const f_t scale = fj_cpu.problem->coefficients[begin]; + if (scale <= 0) continue; + bool uniform_binary = true; + for (i_t p = begin; p < end && uniform_binary; ++p) { + const i_t var = fj_cpu.problem->variables[p]; + const f_t coeff = fj_cpu.problem->coefficients[p]; + const f_t agreement = fj_cpu.hp.exact_k_tol * std::max((f_t)1, std::abs(scale)); + uniform_binary = + fj_cpu.h_is_binary_variable[var] && coeff > 0 && std::abs(coeff - scale) <= agreement; + } + if (!uniform_binary) continue; + + const double cardinality = (double)lb / scale; + const i_t k = (i_t)std::lround(cardinality); + if (std::abs(cardinality - k) <= 1e-4 && k >= 0 && k <= end - begin) + rows.push_back({k, begin, end}); + } + if (rows.empty()) return; + + std::sort(rows.begin(), rows.end(), [](const exact_k_row_t& a, const exact_k_row_t& b) { + return a.end - a.begin < b.end - b.begin; + }); + + const i_t n_variables = fj_cpu.problem->n_variables; + std::vector degree(n_variables, 0); + for (const auto& row : rows) + for (i_t p = row.begin; p < row.end; ++p) + ++degree[fj_cpu.problem->variables[p]]; + + std::vector state(n_variables, -1); + std::vector free_vars; + for (size_t index = 0; index < rows.size(); ++index) { + if ((index & 0xFFF) == 0 && timed_out()) break; + const auto& row = rows[index]; + + i_t selected = 0; + free_vars.clear(); + for (i_t p = row.begin; p < row.end; ++p) { + const i_t var = fj_cpu.problem->variables[p]; + selected += state[var] == 1; + if (state[var] < 0) free_vars.push_back(var); + } + const i_t needed = row.k - selected; + if (needed < 0 || (i_t)free_vars.size() < needed) continue; + + std::sort(free_vars.begin(), free_vars.end(), [°ree](i_t a, i_t b) { + return degree[a] < degree[b]; + }); + for (i_t p = 0; p < (i_t)free_vars.size(); ++p) + state[free_vars[p]] = (int8_t)(p < needed); + } + + recompute_lhs(fj_cpu); + const i_t baseline = fj_cpu.violated_constraints.size(); + const auto anchor = fj_cpu.h_assignment; + for (i_t var = 0; var < n_variables; ++var) + if (state[var] >= 0) fj_cpu.h_assignment[var] = state[var]; + + recompute_lhs(fj_cpu); + if ((i_t)fj_cpu.violated_constraints.size() >= baseline) { + fj_cpu.h_assignment = anchor; + recompute_lhs(fj_cpu); + } + fj_cpu.h_best_assignment = fj_cpu.h_assignment; +} + +template +void repair_difficult_anchor(fj_cpu_climber_t& fj_cpu) +{ + recompute_lhs(fj_cpu); + const i_t baseline = fj_cpu.violated_constraints.size(); + if (baseline == 0 || + baseline <= fj_cpu.problem->n_constraints / fj_cpu.hp.anchor_repair_violated_share) + return; + + const auto started = std::chrono::steady_clock::now(); + const auto anchor = fj_cpu.h_assignment; + const std::vector violated(fj_cpu.violated_constraints.begin(), + fj_cpu.violated_constraints.end()); + std::vector> candidates; + + for (i_t row : violated) { + if (std::chrono::duration(std::chrono::steady_clock::now() - started).count() > + fj_cpu.hp.anchor_repair_budget_s) + break; + + const f_t lb = fj_cpu.problem->cstr_lb[row]; + const f_t ub = fj_cpu.problem->cstr_ub[row]; + f_t sum = fj_cpu.h_lhs[row]; + f_t target = 0; + f_t direction = 0; + if (sum < lb) { + direction = 1; + target = lb; + } else if (sum > ub) { + direction = -1; + target = ub; + } else { + continue; + } + + collect_row_repair_moves(fj_cpu, + fj_cpu.problem->offsets[row], + fj_cpu.problem->offsets[row + 1], + direction, + fj_cpu.hp.exact_k_tol, + candidates); + for (const auto& move : candidates) { + if (direction > 0 ? sum >= target : sum <= target) break; + const f_t delta = move.new_val - (f_t)fj_cpu.h_assignment[move.var]; + sum += move.coeff * delta; + fj_cpu.h_assignment[move.var] = move.new_val; + } + } + + recompute_lhs(fj_cpu); + if ((i_t)fj_cpu.violated_constraints.size() >= baseline) { + fj_cpu.h_assignment = anchor; + recompute_lhs(fj_cpu); + } + fj_cpu.h_best_assignment = fj_cpu.h_assignment; +} + +#if MIP_INSTANTIATE_FLOAT +template void apply_exact_k_start(fj_cpu_climber_t&); +template void repair_difficult_anchor(fj_cpu_climber_t&); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void apply_exact_k_start(fj_cpu_climber_t&); +template void repair_difficult_anchor(fj_cpu_climber_t&); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/chain.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/chain.cpp new file mode 100644 index 0000000000..a2c414f190 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/chain.cpp @@ -0,0 +1,125 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "../audit.hpp" +#include "../internal.hpp" +#include "../problem.hpp" +#include "../search/api.hpp" +#include "starts.hpp" + +namespace cuopt::mathematical_optimization::mip { + +template +void apply_precedence_completion_start(fj_cpu_climber_t& fj_cpu) +{ + phase_timer_t timer(fj_cpu.t_start); + const i_t n_constraints = fj_cpu.problem->n_constraints; + + i_t lower_only = 0; + for (i_t row = 0; row < n_constraints; ++row) + lower_only += + std::isfinite(fj_cpu.problem->cstr_lb[row]) && !std::isfinite(fj_cpu.problem->cstr_ub[row]); + if (lower_only * fj_cpu.hp.precedence_lower_den < n_constraints * fj_cpu.hp.precedence_lower_num) + return; + + const auto started = std::chrono::steady_clock::now(); + auto timed_out = [&] { + return std::chrono::duration(std::chrono::steady_clock::now() - started).count() > + fj_cpu.hp.precedence_budget_s; + }; + + recompute_lhs(fj_cpu); + const auto anchor = fj_cpu.h_assignment; + auto best = anchor; + const i_t anchor_count = fj_cpu.violated_constraints.size(); + const f_t anchor_severity = -fj_cpu.total_violations; + i_t best_count = anchor_count; + f_t best_severity = anchor_severity; + + for (i_t pass = 0; pass < fj_cpu.hp.precedence_passes; ++pass) { + bool changed = false; + for (i_t row = 0; row < n_constraints; ++row) { + if ((row & 0x1FF) == 0 && timed_out()) break; + const f_t lb = fj_cpu.problem->cstr_lb[row]; + if (!std::isfinite(lb) || std::isfinite(fj_cpu.problem->cstr_ub[row])) continue; + const f_t lhs = fj_cpu.h_lhs[row]; + const f_t deficit = lb - lhs; + if (deficit <= fj_cpu.row_tolerance) continue; + + // The head is the row's only positive continuous coefficient. A row with none, or with + // several, is not a precedence row and is left to the search. + i_t head = -1; + f_t head_coeff = 0; + const auto [begin, end] = model_range_for_row(fj_cpu, row); + for (i_t p = begin; p < end; ++p) { + const i_t var = fj_cpu.problem->variables[p]; + const f_t coeff = fj_cpu.problem->coefficients[p]; + if (coeff <= 0 || is_integer_var(fj_cpu, var)) continue; + if (head >= 0 && head != var) { + head = -2; + break; + } + head = var; + head_coeff = coeff; + } + if (head < 0 || head_coeff == 0) continue; + + const auto bounds = fj_cpu.h_var_bounds[head].get(); + const f_t old_val = fj_cpu.h_assignment[head]; + const f_t value = std::min(get_upper(bounds), old_val + deficit / head_coeff); + const f_t delta = value - old_val; + if (!(delta > 0) || !std::isfinite(value)) continue; + + fj_cpu.h_assignment[head] = value; + const auto [col_begin, col_end] = model_range_for_var(fj_cpu, head); + for (i_t q = col_begin; q < col_end; ++q) { + const i_t touched = fj_cpu.problem->reverse_constraints[q]; + const f_t patched = fj_cpu.h_lhs[touched]; + fj_cpu.h_lhs[touched] = patched + fj_cpu.problem->reverse_coefficients[q] * delta; + } + changed = true; + + cuopt_assert( + (f_t)fj_cpu.h_lhs[row] >= lb - fj_cpu.row_tolerance || value >= get_upper(bounds), + "precedence step neither repaired the row nor saturated its head"); + } + + recompute_lhs(fj_cpu); + const i_t count = fj_cpu.violated_constraints.size(); + const f_t severity = -fj_cpu.total_violations; + if (count < best_count || (count == best_count && severity < best_severity)) { + best_count = count; + best_severity = severity; + best = fj_cpu.h_assignment; + if (count == 0) break; + } + if (!changed || timed_out()) break; + } + + cuopt_assert(fj_cpu.h_assignment.size() == anchor.size(), + "incumbent_assignment span would be invalidated"); + const bool keep = + best_count < anchor_count || (best_count == anchor_count && best_severity < anchor_severity); + if (keep) { + fj_cpu.h_assignment = best; + } else { + fj_cpu.h_assignment = anchor; + } + recompute_lhs(fj_cpu); + cuopt_func_call(audit_assignment_bounds(fj_cpu, "precedence start")); + fj_cpu.h_best_assignment = fj_cpu.h_assignment; +} + +#if MIP_INSTANTIATE_FLOAT +template void apply_precedence_completion_start(fj_cpu_climber_t&); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void apply_precedence_completion_start(fj_cpu_climber_t&); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/covering.cpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/covering.cpp new file mode 100644 index 0000000000..7470254678 --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/covering.cpp @@ -0,0 +1,152 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include "../internal.hpp" +#include "../problem.hpp" +#include "../search/api.hpp" +#include "starts.hpp" + +namespace cuopt::mathematical_optimization::mip { + +template +void collect_row_repair_moves(fj_cpu_climber_t& fj_cpu, + i_t row_begin, + i_t row_end, + f_t direction, + f_t tol, + std::vector>& out) +{ + out.clear(); + for (i_t i = row_begin; i < row_end; ++i) { + const i_t var = fj_cpu.problem->variables[i]; + if (!is_integer_var(fj_cpu, var)) continue; + + const f_t coeff = fj_cpu.problem->coefficients[i]; + const f_t val = fj_cpu.h_assignment[var]; + const f_t lb = get_lower(fj_cpu.h_var_bounds[var].get()); + const f_t ub = get_upper(fj_cpu.h_var_bounds[var].get()); + const bool is_bin = fj_cpu.h_is_binary_variable[var] != 0; + + // Raising the variable shifts the sum by `direction * coeff`; lowering it by the negation. + const f_t raise = direction * coeff; + if (raise > 0 && val < ub - tol) { + const f_t new_val = is_bin ? (f_t)1 : std::floor(val) + 1; + if (new_val > val && new_val <= ub + tol) out.push_back({raise, var, coeff, new_val}); + } else if (raise < 0 && val > lb + tol) { + const f_t new_val = is_bin ? (f_t)0 : std::ceil(val) - 1; + if (new_val < val && new_val >= lb - tol) out.push_back({-raise, var, coeff, new_val}); + } + } + std::sort(out.begin(), + out.end(), + [](const row_repair_move_t& a, const row_repair_move_t& b) { + return a.effect > b.effect; + }); +} + +template +void apply_greedy_covering_start(fj_cpu_climber_t& fj_cpu) +{ + if (fj_cpu.problem->nnz > fj_cpu.hp.start_nnz_limit) return; + + recompute_lhs(fj_cpu); + const i_t baseline_violated = fj_cpu.violated_constraints.size(); + const auto anchor_assignment = fj_cpu.h_assignment; + + const i_t n_constraints = fj_cpu.problem->n_constraints; + std::vector row_order(n_constraints); + for (i_t i = 0; i < n_constraints; ++i) + row_order[i] = i; + std::sort(row_order.begin(), row_order.end(), [&](i_t a, i_t b) { + return (fj_cpu.problem->offsets[a + 1] - fj_cpu.problem->offsets[a]) < + (fj_cpu.problem->offsets[b + 1] - fj_cpu.problem->offsets[b]); + }); + + const auto started = std::chrono::steady_clock::now(); + const double time_budget_s = fj_cpu.hp.covering_budget_s; + const f_t tol = 1e-6; + const i_t max_passes = 2; + std::vector> candidates; + bool out_of_time = false; + + for (i_t pass = 0; pass < max_passes && !out_of_time; ++pass) { + for (i_t k = 0; k < n_constraints; ++k) { + if ((k & 0xFFF) == 0 && + std::chrono::duration(std::chrono::steady_clock::now() - started).count() > + time_budget_s) { + out_of_time = true; + break; + } + const i_t cstr_idx = row_order[k]; + const i_t row_begin = fj_cpu.problem->offsets[cstr_idx]; + const i_t row_end = fj_cpu.problem->offsets[cstr_idx + 1]; + if (row_begin == row_end) continue; + + const f_t lb = fj_cpu.problem->cstr_lb[cstr_idx]; + const f_t ub = fj_cpu.problem->cstr_ub[cstr_idx]; + const bool has_lb = std::isfinite(lb); + const bool has_ub = std::isfinite(ub); + if (!has_lb && !has_ub) continue; + + f_t sum = compensated_dot2_csr(*fj_cpu.problem, fj_cpu.h_assignment, cstr_idx); + + // Equality rows are driven to their bound; one-sided rows only to the side they violate. + const bool is_equality = has_lb && has_ub && std::abs(lb - ub) < tol; + f_t direction = 0; + f_t target = 0; + if (is_equality && std::abs(sum - lb) > tol) { + direction = sum < lb ? (f_t)1 : (f_t)-1; + target = lb; + } else if (has_lb && sum < lb - tol) { + direction = 1; + target = lb; + } else if (has_ub && sum > ub + tol) { + direction = -1; + target = ub; + } else { + continue; + } + + collect_row_repair_moves(fj_cpu, row_begin, row_end, direction, tol, candidates); + for (const auto& m : candidates) { + if (direction > 0 ? sum >= target - tol : sum <= target + tol) break; + const f_t delta = m.new_val - (f_t)fj_cpu.h_assignment[m.var]; + sum += m.coeff * delta; + fj_cpu.h_assignment[m.var] = m.new_val; + } + } + } + + recompute_lhs(fj_cpu); + if ((i_t)fj_cpu.violated_constraints.size() >= baseline_violated) { + fj_cpu.h_assignment = anchor_assignment; + recompute_lhs(fj_cpu); + } + fj_cpu.h_best_assignment = fj_cpu.h_assignment; +} + +#if MIP_INSTANTIATE_FLOAT +template void collect_row_repair_moves(fj_cpu_climber_t&, + int, + int, + float, + float, + std::vector>&); +template void apply_greedy_covering_start(fj_cpu_climber_t&); +#endif + +#if MIP_INSTANTIATE_DOUBLE +template void collect_row_repair_moves(fj_cpu_climber_t&, + int, + int, + double, + double, + std::vector>&); +template void apply_greedy_covering_start(fj_cpu_climber_t&); +#endif + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/starts.hpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/starts.hpp new file mode 100644 index 0000000000..3c0d907e4c --- /dev/null +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/starts/starts.hpp @@ -0,0 +1,43 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include "../internal.hpp" + +namespace cuopt::mathematical_optimization::mip { + +template +struct row_repair_move_t { + f_t effect; + i_t var; + f_t coeff; + f_t new_val; +}; + +template +void collect_row_repair_moves(fj_cpu_climber_t& c, + i_t row_begin, + i_t row_end, + f_t direction, + f_t tolerance, + std::vector>& out); + +template +void apply_greedy_covering_start(fj_cpu_climber_t& c); + +template +void apply_exact_k_start(fj_cpu_climber_t& c); + +template +void repair_difficult_anchor(fj_cpu_climber_t& c); + +template +void apply_precedence_completion_start(fj_cpu_climber_t& c); + +template +void apply_structural_completion_start(fj_cpu_climber_t& c); + +} // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/cpu/state.hpp b/cpp/src/mip_heuristics/feasibility_jump/cpu/state.hpp index 4ce9496331..c1fafea66e 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/cpu/state.hpp +++ b/cpp/src/mip_heuristics/feasibility_jump/cpu/state.hpp @@ -372,9 +372,7 @@ template struct fj_runtime_t { explicit fj_runtime_t(std::atomic& flag) : preemption_flag(flag) {} i_t log_interval{0}; - i_t diversity_callback_interval{3000}; std::function&, double)> improvement_callback{nullptr}; - std::function&)> diversity_callback{nullptr}; std::string log_prefix; std::shared_ptr> shared_incumbent; std::atomic work_units_elapsed{0.0}; @@ -512,4 +510,15 @@ std::unique_ptr> init_fj_cpu_clone( std::atomic& preemption_flag, fj_settings_t settings = fj_settings_t{}); +template +void apply_lane_diversification(fj_cpu_climber_t& climber, int lane, int64_t base_seed); + +template +void complete_climber_portfolio(std::unique_ptr> first_climber, + const std::vector& lane_seeds, + std::vector>& preemption_flags, + std::vector>>& climbers, + int64_t base_seed, + bool low_latency = false); + } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cu b/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cu index 0547a12354..b0d9f091ed 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cu @@ -8,9 +8,13 @@ #include "early_cpufj.cuh" #include +#include #include +#include +#include + namespace cuopt::mathematical_optimization::mip { template @@ -33,11 +37,11 @@ early_cpufj_t::~early_cpufj_t() } template -void early_cpufj_t::start(bool low_latency) +void early_cpufj_t::start(int n_lanes, bool low_latency) { const bool threaded = !omp_in_parallel(); // 1: presolve, 1: early GPU FJ, 1: early CPU FJ - if (climber_ || + if (!climbers_.empty() || (!threaded && omp_get_num_threads() < CUOPT_MIP_EARLY_CPUFJ_REQUIRED_THREAD_COUNT)) { return; } @@ -45,47 +49,84 @@ void early_cpufj_t::start(bool low_latency) this->preemption_flag_.store(false); this->start_time_ = std::chrono::steady_clock::now(); + // Tasks are not preempted, so a lane posted beyond the team size would sit in the queue for the + // whole of presolve without running an iteration. + n_lanes = threaded ? 1 : std::clamp(n_lanes, 1, omp_get_num_threads()); + const int64_t base_seed = cuopt::seed_generator::get_seed(); + climbers_.resize(n_lanes); + auto report_incumbent = [this](f_t solver_obj, const std::vector& assignment, double) { + std::lock_guard guard(incumbent_mutex_); this->try_update_best(solver_obj, assignment); }; - fj_settings_t settings; - settings.seed = (int)seed_; - climber_ = init_fj_cpu_from_optimization_problem( - *this->problem_ptr_, tolerances_, preemption_flag_, settings); - climber_->low_latency = low_latency; - climber_->log_prefix = "[Early CPUFJ] "; - climber_->improvement_callback = report_incumbent; + // Lane 0 builds the host problem representation and every other lane copies it. All of it + // finishes before the first task is posted, so no lane reads a template another lane is running + // on. seed_generator steps a non-atomic global, which is why the draws stay on this thread. + for (int k = 0; k < n_lanes; ++k) { + if (k == 0) { + climbers_[0] = + init_fj_cpu_from_optimization_problem(*this->problem_ptr_, tolerances_, preemption_flag_); + } else { + fj_settings_t settings; + settings.seed = (int)cuopt::seed_generator::get_seed(); + climbers_[k] = init_fj_cpu_clone(*climbers_[0], preemption_flag_, settings); + } + climbers_[k]->low_latency = low_latency; + apply_lane_diversification(*climbers_[k], k, base_seed); + climbers_[k]->log_prefix = "[Early CPUFJ " + std::to_string(k) + "] "; + climbers_[k]->improvement_callback = report_incumbent; + } - CUOPT_LOG_DEBUG("Launching early CPUFJ %s", threaded ? "thread" : "task"); - auto* climber = climber_.get(); + auto shared = std::make_shared>(); + for (int k = 0; k < n_lanes; ++k) + climbers_[k]->shared_incumbent = shared; + + CUOPT_LOG_DEBUG("Launching %d early CPUFJ %s", n_lanes, threaded ? "thread" : "tasks"); if (threaded) { - worker_ = std::thread([climber] { cpufj_solve(climber); }); + auto* climber = climbers_[0].get(); + worker_ = std::thread([climber] { cpufj_solve(climber); }); return; } + for (int k = 0; k < n_lanes; ++k) { + auto* climber = climbers_[k].get(); #pragma omp task firstprivate(climber) priority(CUOPT_DEFAULT_TASK_PRIORITY) \ depend(out : *climber) default(none) - cpufj_solve(climber); + cpufj_solve(climber); + } } template void early_cpufj_t::stop() { - if (!climber_) { return; } + if (climbers_.empty()) { return; } preemption_flag_.store(true); - climber_->halted = true; + + // Every lane is told to stop before any wait, otherwise the first wait blocks on a lane that has + // not been asked to exit yet. + for (auto& climber : climbers_) { + climber->halted = true; + } if (worker_.joinable()) { worker_.join(); } else { -#pragma omp taskwait depend(in : *climber_) // Wait for the early CPUFJ task to finish + for (size_t k = 0; k < climbers_.size(); ++k) { +#pragma omp taskwait depend(in : *climbers_[k]) // Wait for each early CPUFJ task to finish + } + } + + [[maybe_unused]] i_t total_iterations = 0; + for (const auto& climber : climbers_) { + total_iterations += climber->iterations; } - CUOPT_LOG_DEBUG("[Early CPUFJ] Stopped after %d iterations, solution_found=%d", - climber_->iterations, + CUOPT_LOG_DEBUG("[Early CPUFJ] Stopped after %d iterations over %d climbers, solution_found=%d", + total_iterations, + (int)climbers_.size(), this->solution_found_); - climber_.reset(); + climbers_.clear(); } template diff --git a/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cuh b/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cuh index 531e2ea704..2d6ceee04c 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/early_cpufj.cuh @@ -12,6 +12,7 @@ #include #include +#include #include #include @@ -29,9 +30,13 @@ class early_cpufj_t : public early_heuristic_t static constexpr const char* name() { return "CPUFJ"; } - void start(bool low_latency = false); + // Lanes are OMP tasks that never yield, so n_lanes threads are unavailable to anything else + // until stop(). Callers sharing the team with other work size it accordingly. + void start(int n_lanes, bool low_latency = false); void stop(); + int lane_count() const { return (int)climbers_.size(); } + private: friend class early_heuristic_t>; @@ -39,12 +44,15 @@ class early_cpufj_t : public early_heuristic_t const optimization_problem_t* problem_ptr_{nullptr}; typename mip_solver_settings_t::tolerances_t tolerances_; - std::unique_ptr> climber_; + std::vector>> climbers_; std::thread worker_; std::atomic preemption_flag_{false}; // Explicit seed for this climber's FJ RNG, resolved once from the solve's base seed (see // mip_solver_context_t::base_seed) since this heuristic runs before that context exists. uint64_t seed_; + // try_update_best and the incumbent callback behind it are not thread-safe, and every lane + // reports into them from its own task. + std::mutex incumbent_mutex_; }; } // namespace cuopt::mathematical_optimization::mip diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh index b803191642..316e15a38f 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cuh @@ -18,6 +18,20 @@ namespace cuopt::mathematical_optimization::mip { template class problem_t; +template +std::unique_ptr> init_fj_cpu_standalone( + problem_t& problem, + std::atomic& preemption_flag, + uint64_t seed, + fj_settings_t settings = fj_settings_t{}); + +template +void build_climber_portfolio(problem_t& problem, + std::vector>& preemption_flags, + std::vector>>& climbers, + int64_t base_seed, + bool low_latency = false); + template std::unique_ptr> init_fj_cpu_from_optimization_problem( const optimization_problem_t& problem, diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cu index 93581e55c8..99bbddbe45 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cu @@ -42,6 +42,49 @@ static const char* fj_binary_reject_name(fj_binary_reject_t reason) return "unknown"; } +template +static void audit_binary_incumbent(const fj_cpu_climber_t& climber, + const ins_vector& assignment, + f_t reported_objective, + f_t engine_objective) +{ + const auto& problem = *climber.problem; + cuopt_assert(assignment.size() == (size_t)problem.n_variables, "binary incumbent size mismatch"); + for (i_t variable = 0; variable < problem.n_variables; ++variable) { + const f_t value = assignment[variable]; + cuopt_assert(std::isfinite(value), "binary incumbent contains a non-finite value"); + cuopt_assert(climber.check_variable_within_bounds(variable, value), + "binary incumbent violates original variable bounds"); + if (problem.h_var_types[variable] == var_t::INTEGER) { + cuopt_assert(problem.is_integer(value), "binary incumbent violates original integrality"); + } + } + const f_t row_tolerance = problem.tolerances.absolute_tolerance; + for (i_t row = 0; row < problem.n_constraints; ++row) { + const f_t activity = compensated_dot2_csr(problem, assignment, row); + cuopt_assert(std::isfinite(activity), "binary incumbent has non-finite row activity"); + cuopt_assert( + !std::isfinite(problem.cstr_lb[row]) || activity >= problem.cstr_lb[row] - row_tolerance, + "binary incumbent violates original row lower bound"); + cuopt_assert( + !std::isfinite(problem.cstr_ub[row]) || activity <= problem.cstr_ub[row] + row_tolerance, + "binary incumbent violates original row upper bound"); + } + const f_t objective = + compensated_dot2(problem.h_obj_coeffs.data(), assignment.data(), problem.n_variables); + const f_t objective_scale = + std::max(f_t{1}, std::max(std::fabs(objective), std::fabs(reported_objective))); + const f_t objective_tolerance = + std::max(problem.tolerances.absolute_tolerance, + f_t{64} * std::numeric_limits::epsilon() * objective_scale); + cuopt_assert(std::fabs(objective - reported_objective) <= objective_tolerance, + "binary incumbent objective disagrees with the original assignment"); + if (!climber.bin_singletons.empty()) { + cuopt_assert(std::fabs(objective - engine_objective) <= objective_tolerance, + "binary singleton substitution changed the objective"); + } +} + // work unit proxy. will likely require a lot of tuning constexpr double fj_bin_bytes_per_nnz = 16.0; // restarts can help a lot on some smaller combinatorial instances @@ -237,9 +280,8 @@ struct fj_bin_engine_t { row_slack[r] = slack; if (slack < 0) set_violated(r); } - incumbent_objective = objective_offset; - for (int32_t v = 0; v < pb.n_variables; ++v) - incumbent_objective += pb.objective[v] * assign[v]; + incumbent_objective = + objective_offset + compensated_dot2(pb.objective.data(), assign.data(), pb.n_variables); nnz_touched += pb.nnz; rebuild_scores(); } @@ -421,11 +463,14 @@ struct fj_bin_engine_t { ? get_lower(bounds) : (isfinite(get_upper(bounds)) ? get_upper(bounds) : f_t{0}); } - f_t lhs = 0; - for (i_t p = climber.problem->offsets[rec.row]; p < climber.problem->offsets[rec.row + 1]; - ++p) - lhs += climber.problem->coefficients[p] * values[climber.problem->variables[p]]; + const f_t lhs = compensated_dot2_csr(*climber.problem, values, rec.row); const f_t residual = rec.rhs - lhs; + if (rec.all.size() == 1) { + const i_t var = rec.all[0]; + values[var] += + residual / climber.problem->reverse_coefficients[climber.problem->reverse_offsets[var]]; + continue; + } if (residual > 0 && !rec.positive.empty()) values[rec.positive[0]] += residual / rec.positive_coeff[0]; else if (residual < 0 && !rec.negative.empty()) @@ -433,6 +478,13 @@ struct fj_bin_engine_t { } } + static f_t original_objective(const fj_cpu_climber_t& climber, + const ins_vector& values) + { + return compensated_dot2( + climber.problem->h_obj_coeffs.data(), values.data(), climber.problem->h_obj_coeffs.size()); + } + void report_incumbent(fj_cpu_climber_t& climber) { auto& h_assign = climber.h_assignment; @@ -454,24 +506,21 @@ struct fj_bin_engine_t { uncrush(climber, h_assign); uncrush(climber, h_best); } - auto objective = [&](const auto& values) { - f_t result = 0; - for (i_t var = 0; var < (i_t)climber.problem->h_obj_coeffs.size(); ++var) - result += climber.problem->h_obj_coeffs[var] * values[var]; - return result; - }; - const f_t reported = climber.has_bin_elimination ? objective(h_best) : (f_t)best_objective; - climber.h_incumbent_objective = - climber.has_bin_elimination ? objective(h_assign) : (f_t)incumbent_objective; - climber.h_best_objective = reported; - climber.feasible_found = true; + const f_t reported = + climber.has_bin_elimination ? original_objective(climber, h_best) : (f_t)best_objective; + climber.h_incumbent_objective = climber.has_bin_elimination + ? original_objective(climber, h_assign) + : (f_t)incumbent_objective; + climber.h_best_objective = reported; + climber.feasible_found = true; + cuopt_func_call(audit_binary_incumbent(climber, h_best, reported, (f_t)best_objective)); if (shared_incumbent) shared_incumbent->publish(reported, climber.get_user_objective(reported), h_best); CUOPT_LOG_DEBUG("%sCPUFJ[bin%d] new incumbent: objective %.17g", climber.log_prefix.c_str(), coefficient_bits(), - reported); + climber.get_user_objective(reported)); if (climber.improvement_callback) { const double work_units = climber.work_units_elapsed; climber.improvement_callback(reported, h_best, work_units); @@ -925,9 +974,9 @@ struct fj_bin_engine_t { cuopt_assert(std::isfinite(obj_magnitude) && obj_magnitude > 0, "objective magnitude unit must be finite and positive"); - objective_offset = 0; - for (int32_t v = 0; v < pb.n_original; ++v) - objective_offset += pb.orig_objective[v] * pb.var_offset[v]; + objective_offset = + pb.substitution_offset + + compensated_dot2(pb.orig_objective.data(), pb.var_offset.data(), pb.n_original); argmax_tile = fj_bin_argmax_tile(); set_objective_weight(seeded_weight > 0 ? seeded_weight : 0); @@ -1012,18 +1061,9 @@ struct fj_bin_engine_t { coefficient_bits(), iters, violated_list.size(), - best_objective, + climber.get_user_objective((f_t)best_objective), max_weight); } - if (iters % climber.diversity_callback_interval == 0 && climber.diversity_callback) { - auto& h_assign = climber.h_assignment; - for (int32_t v = 0; v < pb.n_original; ++v) - h_assign[v] = (f_t)pb.var_offset[v]; - for (int32_t b = 0; b < pb.n_variables; ++b) - if (assign[b]) h_assign[pb.bit_owner[b]] += (f_t)pb.bit_weight[b]; - climber.diversity_callback((f_t)incumbent_objective, h_assign); - } - // Work-unit proxy. nnz_touched is cumulative, reproducing the accumulation shape the general // path gets from its cumulative byte counters. if (iters % 100 == 0 && iters > 0) { @@ -1046,7 +1086,7 @@ struct fj_bin_engine_t { climber.log_prefix.c_str(), coefficient_bits(), iters, - best_objective, + climber.get_user_objective((f_t)best_objective), max_weight, max_aggregate_base, fj_bin_base_limit, diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cuh index f90c064a90..4e68d6d942 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary.cuh @@ -95,6 +95,7 @@ struct fj_bin_problem_t { std::vector original_to_bin_mapping; std::vector bit_weight; std::vector orig_objective; + double substitution_offset{0}; }; // Result of the width-independent eligibility scan. diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary_preprocess.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary_preprocess.cu index c44a93d119..d333baef6e 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary_preprocess.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_binary_preprocess.cu @@ -98,6 +98,11 @@ fj_bin_scan_t fj_bin_scan(const fj_cpu_climber_t& c, fj_bin_setup_time { phase_timer_t timer(times.scan); fj_bin_scan_t out; + // Singleton row/objective substitution is performed by the encoded path. + if (!c.bin_singletons.empty()) { + out.reject = fj_binary_reject_t::non_binary_var; + return out; + } const int32_t n_cols = c.problem->n_variables; const int32_t n_rows = c.problem->n_constraints; if (n_cols <= 0 || n_rows <= 0) { @@ -399,6 +404,24 @@ bool fj_bin_encode(const fj_cpu_climber_t& c, const auto& left_w = c.h_initial_left_weights; const auto& right_w = c.h_initial_right_weights; const auto& obj = c.problem->h_obj_coeffs; + std::vector encoded_obj(obj.begin(), obj.end()); + std::vector singleton(n_rows, -1); + std::vector substituted(n_cols, 0); + pb.substitution_offset = 0; + for (const auto& [row, var] : c.bin_singletons) { + singleton[row] = var; + substituted[var] = 1; + const double pivot = c.problem->reverse_coefficients[c.problem->reverse_offsets[var]]; + const double cost = obj[var] / pivot; + pb.substitution_offset += cost * cstr_lb[row]; + encoded_obj[var] = 0; + for (int32_t entry = offsets[row]; entry < offsets[row + 1]; ++entry) { + if (variables[entry] != var) encoded_obj[variables[entry]] -= cost * coeffs[entry]; + } + } + if (!std::isfinite(pb.substitution_offset)) return false; + for (double cost : encoded_obj) + if (!std::isfinite(cost)) return false; std::vector lower(n_cols); std::vector upper(n_cols); @@ -407,6 +430,7 @@ bool fj_bin_encode(const fj_cpu_climber_t& c, int64_t total_bits = 0; // count the total bits that'd be required to encode this model as pure-binary for (int32_t v = 0; v < n_cols; ++v) { + if (substituted[v]) continue; if (var_types[v] != var_t::INTEGER) return false; auto bounds = var_bounds[v]; const double x = (double)cuopt::get_lower(bounds); @@ -462,14 +486,14 @@ bool fj_bin_encode(const fj_cpu_climber_t& c, // emit onesided a row auto emit = [&](int32_t r, double side_bound, long side, double weight) -> bool { - double fixed = 0; - for (int32_t k = offsets[r]; k < offsets[r + 1]; ++k) - fixed += coeffs[k] * lower[variables[k]]; + const double fixed = + compensated_dot2_csr(offsets.data(), variables.data(), coeffs.data(), lower.data(), r); const double folded_bound = side_bound - fixed; row_values.clear(); bool integral = true; for (int32_t k = offsets[r]; k < offsets[r + 1]; ++k) { + if (substituted[variables[k]]) continue; row_values.push_back(coeffs[k]); if (!is_integer(coeffs[k], tol)) integral = false; } @@ -485,7 +509,8 @@ bool fj_bin_encode(const fj_cpu_climber_t& c, double row_abs_sum = 0; for (int32_t k = offsets[r]; k < offsets[r + 1]; ++k) { const int32_t v = variables[k]; - const double a = s * coeffs[k]; + if (substituted[v]) continue; + const double a = s * coeffs[k]; if (!is_integer(a, tol)) return false; const long ai = std::lround(a); for (int32_t bk = 0; bk < nbits[v]; ++bk) { @@ -523,8 +548,15 @@ bool fj_bin_encode(const fj_cpu_climber_t& c, const uint8_t* ignored_row = c.has_bin_elimination ? c.bin_ignore_row.data() : nullptr; for (int32_t r = 0; r < n_rows; ++r) { if (ignored_row && ignored_row[r]) continue; - const double lb = cstr_lb[r]; - const double ub = cstr_ub[r]; + double lb = cstr_lb[r], ub = cstr_ub[r]; + if (singleton[r] >= 0) { + const int32_t var = singleton[r]; + const double pivot = c.problem->reverse_coefficients[c.problem->reverse_offsets[var]]; + const auto bounds = var_bounds[var]; + const double left = pivot * get_lower(bounds), right = pivot * get_upper(bounds); + lb = cstr_lb[r] - std::max(left, right); + ub = cstr_lb[r] - std::min(left, right); + } if (std::isfinite(lb) && !emit(r, lb, -1, left_w[r])) return false; if (std::isfinite(ub) && !emit(r, ub, 1, right_w[r])) return false; } @@ -549,11 +581,11 @@ bool fj_bin_encode(const fj_cpu_climber_t& c, pb.objective.assign(n_bits, 0.0); pb.objective_vars.clear(); for (int32_t v = 0; v < n_cols; ++v) { - pb.orig_objective[v] = obj[v]; - if (obj[v] == 0.0) continue; + pb.orig_objective[v] = encoded_obj[v]; + if (encoded_obj[v] == 0.0) continue; for (int32_t bk = 0; bk < nbits[v]; ++bk) { const int32_t bit = bit_start[v] + bk; - pb.objective[bit] = obj[v] * pb.bit_weight[bit]; + pb.objective[bit] = encoded_obj[v] * pb.bit_weight[bit]; if (pb.objective[bit] != 0.0) pb.objective_vars.push_back(bit); } } diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_bridge.cu b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_bridge.cu index 517d78876b..a6cd4a5e40 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_bridge.cu +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_bridge.cu @@ -155,6 +155,53 @@ std::unique_ptr> init_fj_cpu_from_optimization_proble settings); } +template +std::unique_ptr> init_fj_cpu_standalone( + problem_t& problem, + std::atomic& preemption_flag, + uint64_t seed, + fj_settings_t settings) +{ + raft::common::nvtx::range scope("init_fj_cpu_standalone"); + auto fj_cpu = std::make_unique>(preemption_flag); + std::vector default_weights(problem.n_constraints, f_t{1}); + const probing_cache_t* no_implications = nullptr; + init_fj_cpu_from_problem(*fj_cpu, + problem, + problem.handle_ptr, + std::vector{}, + default_weights, + default_weights, + f_t{0}, + no_implications); + fj_cpu->settings = settings; + fj_cpu->settings.seed = seed; + return fj_cpu; +} + +template +void build_climber_portfolio(problem_t& problem, + std::vector>& preemption_flags, + std::vector>>& climbers, + int64_t base_seed, + bool low_latency) +{ + cuopt_assert(!climbers.empty(), "a CPUFJ portfolio needs at least one climber"); + cuopt_assert(preemption_flags.size() == climbers.size(), "preemption flag count mismatch"); + + std::vector lane_seeds(climbers.size()); + for (size_t k = 0; k < climbers.size(); ++k) { + preemption_flags[k].store(false); + lane_seeds[k] = base_seed + k; + } + + fj_settings_t settings; + settings.seed = lane_seeds[0]; + auto first = init_fj_cpu_standalone(problem, preemption_flags[0], lane_seeds[0], settings); + complete_climber_portfolio( + std::move(first), lane_seeds, preemption_flags, climbers, base_seed, low_latency); +} + template std::unique_ptr> fj_t::create_cpu_climber( solution_t& solution, @@ -198,6 +245,14 @@ template std::unique_ptr> init_fj_cpu_from_optimiza const typename mip_solver_settings_t::tolerances_t&, std::atomic&, fj_settings_t); +template std::unique_ptr> init_fj_cpu_standalone( + problem_t&, std::atomic&, uint64_t, fj_settings_t); +template void build_climber_portfolio( + problem_t&, + std::vector>&, + std::vector>>&, + int64_t, + bool); template std::unique_ptr> fj_t::create_cpu_climber( solution_t&, const std::vector&, @@ -215,6 +270,14 @@ template std::unique_ptr> init_fj_cpu_from_optimiz const typename mip_solver_settings_t::tolerances_t&, std::atomic&, fj_settings_t); +template std::unique_ptr> init_fj_cpu_standalone( + problem_t&, std::atomic&, uint64_t, fj_settings_t); +template void build_climber_portfolio( + problem_t&, + std::vector>&, + std::vector>>&, + int64_t, + bool); template std::unique_ptr> fj_t::create_cpu_climber( solution_t&, const std::vector&, diff --git a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh index 521df10b8b..50a8f9f423 100644 --- a/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh +++ b/cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh @@ -50,7 +50,8 @@ struct fj_cpu_worker_t { const std::vector& start_assignment, const simplex::simplex_solver_settings_t& settings, std::string log_prefix, - int64_t seed = -1); + int64_t seed = -1, + int lane = -1); // Run the worker asynchronously (i.e., launch an openmp task and then continue the // execution). Call `stop()` for stopping the worker diff --git a/cpp/src/mip_heuristics/mip_constants.hpp b/cpp/src/mip_heuristics/mip_constants.hpp index 11a67aef3f..d9cacd9b4a 100644 --- a/cpp/src/mip_heuristics/mip_constants.hpp +++ b/cpp/src/mip_heuristics/mip_constants.hpp @@ -29,6 +29,10 @@ #define CUOPT_MIP_BATCH_PDLP_REQUIRED_THREAD_COUNT 3 #define CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT 3 +/* @brief Threads the early CPUFJ portfolio leaves to the rest of the team. Every lane holds its + * own host copy of the problem and occupies an OMP task for the whole of presolve. */ +#define CUOPT_MIP_EARLY_CPUFJ_RESERVED_THREADS 4 + /* @brief Priority classes for the omp tasks. Highest value = higher priority. * Note that this only gives a hint to the runtime, such that the high priority * is not guarantee to be executed before a low priority one (i.e., do not rely on diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 07e0080206..263122bd03 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -320,7 +320,7 @@ mip_solution_t run_mip_solver( if (std::isfinite(initial_upper_bound)) { early_cpufj->set_best_objective(problem.get_solver_obj_from_user_obj(initial_upper_bound)); } - early_cpufj->start(); + early_cpufj->start(omp_get_num_threads() - CUOPT_MIP_EARLY_CPUFJ_RESERVED_THREADS); solver.context.early_cpufj_ptr = early_cpufj.get(); CUOPT_LOG_DEBUG("Started early CPUFJ on papilo-presolved problem during cuOpt presolve"); @@ -567,8 +567,10 @@ mip_solution_t solve_mip_helper( if (pre_solve_heuristics && pre_solve_heuristics->solution_found()) { early_cpufj->set_best_objective(pre_solve_heuristics->get_best_objective()); } - early_cpufj->start(); - CUOPT_LOG_DEBUG("Started early CPUFJ on original problem"); + // Papilo runs on its own threads, so the team is otherwise idle here. + early_cpufj->start(omp_get_num_threads() - CUOPT_MIP_EARLY_CPUFJ_RESERVED_THREADS); + CUOPT_LOG_DEBUG("Started early CPUFJ on original problem with %d lanes", + early_cpufj->lane_count()); } auto early_cpufj_guard = cuopt::scope_guard([&]() { @@ -965,7 +967,7 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, std::chrono::duration(std::chrono::steady_clock::now() - probe_start).count()); }, mip::derive_seed(settings_const.seed, mip::rng_id_t::early_cpufj)); - pre_solve_heuristics->start(/*low_latency=*/true); + pre_solve_heuristics->start(1, /*low_latency=*/true); } cuopt::scope_guard release_probe([&pre_solve_heuristics] { if (pre_solve_heuristics) { pre_solve_heuristics->stop(); }