diff --git a/cpp/include/cuopt/mathematical_optimization/constants.h b/cpp/include/cuopt/mathematical_optimization/constants.h index 024fae7583..5e720a645e 100644 --- a/cpp/include/cuopt/mathematical_optimization/constants.h +++ b/cpp/include/cuopt/mathematical_optimization/constants.h @@ -72,6 +72,7 @@ #define CUOPT_MIP_FLOW_COVER_CUTS "mip_flow_cover_cuts" #define CUOPT_MIP_IMPLIED_BOUND_CUTS "mip_implied_bound_cuts" #define CUOPT_MIP_CLIQUE_CUTS "mip_clique_cuts" +#define CUOPT_MIP_ZERO_HALF_CUTS "mip_zero_half_cuts" #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS "mip_strong_chvatal_gomory_cuts" #define CUOPT_MIP_REDUCED_COST_STRENGTHENING "mip_reduced_cost_strengthening" #define CUOPT_MIP_OBJECTIVE_STEP "mip_objective_step" diff --git a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp index ba4dce20dc..6d2dc8e694 100644 --- a/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp +++ b/cpp/include/cuopt/mathematical_optimization/mip/solver_settings.hpp @@ -127,6 +127,7 @@ class mip_solver_settings_t { i_t knapsack_cuts = -1; i_t flow_cover_cuts = -1; i_t clique_cuts = -1; + i_t zero_half_cuts = -1; i_t implied_bound_cuts = -1; i_t strong_chvatal_gomory_cuts = -1; i_t reduced_cost_strengthening = -1; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a17bc6716c..3a57ba49ff 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2493,7 +2493,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut omp_atomic_t* clique_signal = &signal_extend_cliques_; - if (settings_.clique_cuts != 0 && clique_table_ == nullptr && + if ((settings_.clique_cuts != 0 || settings_.zero_half_cuts != 0) && clique_table_ == nullptr && omp_get_num_threads() >= CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT) { signal_extend_cliques_.store(false, std::memory_order_release); typename mip_solver_settings_t::tolerances_t tolerances_for_clique{}; @@ -2510,7 +2510,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut user_problem_t problem_copy = original_problem_; timer_t timer(std::numeric_limits::infinity()); mip::find_initial_cliques( - problem_copy, tolerances_for_clique, &clique_table_, timer, clique_signal); + problem_copy, tolerances_for_clique, clique_table_, timer, clique_signal); } } diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index 40efbe6185..6c4fa51144 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -49,22 +50,38 @@ using simplex::vector_norm_inf; namespace { -#define DEBUG_CLIQUE_CUTS 0 -#define CHECK_WORKSPACE 0 +#define DEBUG_CLIQUE_CUTS 0 +#define DEBUG_ZERO_HALF_CUTS 0 +#define CHECK_WORKSPACE 0 enum class clique_cut_build_status_t : int8_t { NO_CUT = 0, CUT_ADDED = 1, INFEASIBLE = 2 }; -#if DEBUG_CLIQUE_CUTS -#define CLIQUE_CUTS_DEBUG(...) \ - do { \ - std::fprintf(stderr, "[DEBUG_CLIQUE_CUTS] "); \ - std::fprintf(stderr, __VA_ARGS__); \ - std::fprintf(stderr, "\n"); \ +// Shared crash-tolerant debug logger: writes a prefixed line to stderr and +// flushes immediately so the last line is visible even if the process +// aborts/terminates right after. Each channel below enables it through its own +// DEBUG_* flag and supplies its own prefix; when the flag is 0 the call expands +// to a no-op that still consumes its arguments. +#define CUTS_DEBUG_LOG(prefix, ...) \ + do { \ + std::fprintf(stderr, prefix " "); \ + std::fprintf(stderr, __VA_ARGS__); \ + std::fprintf(stderr, "\n"); \ + std::fflush(stderr); \ } while (0) -#else -#define CLIQUE_CUTS_DEBUG(...) \ - do { \ +#define CUTS_DEBUG_NOOP(...) \ + do { \ } while (0) + +#if DEBUG_CLIQUE_CUTS +#define CLIQUE_CUTS_DEBUG(...) CUTS_DEBUG_LOG("[DEBUG_CLIQUE_CUTS]", __VA_ARGS__) +#else +#define CLIQUE_CUTS_DEBUG(...) CUTS_DEBUG_NOOP(__VA_ARGS__) +#endif + +#if DEBUG_ZERO_HALF_CUTS +#define ZERO_HALF_DEBUG(...) CUTS_DEBUG_LOG("[zero_half]", __VA_ARGS__) +#else +#define ZERO_HALF_DEBUG(...) CUTS_DEBUG_NOOP(__VA_ARGS__) #endif template @@ -113,12 +130,12 @@ clique_cut_build_status_t build_clique_cut(const std::vector& clique_vertic std::unordered_set complement_pairs; seen_original.reserve(clique_vertices.size()); seen_complement.reserve(clique_vertices.size()); - for (const auto vertex_idx : clique_vertices) { + for (const i_t vertex_idx : clique_vertices) { cuopt_assert(vertex_idx >= 0 && vertex_idx < 2 * num_vars, "Clique vertex out of range"); - const i_t var_idx = vertex_idx % num_vars; - const bool complement = vertex_idx >= num_vars; - const f_t lower_bound = lower_bounds[var_idx]; - const f_t upper_bound = upper_bounds[var_idx]; + const i_t var_idx = vertex_idx % num_vars; + const bool complement = vertex_idx >= num_vars; + [[maybe_unused]] const f_t lower_bound = lower_bounds[var_idx]; + [[maybe_unused]] const f_t upper_bound = upper_bounds[var_idx]; cuopt_assert(var_types[var_idx] != variable_type_t::CONTINUOUS, "Clique contains continuous variable"); cuopt_assert(lower_bound >= -bound_tol, "Clique variable lower bound below zero"); @@ -149,7 +166,7 @@ clique_cut_build_status_t build_clique_cut(const std::vector& clique_vertic // variable from the support and bump rhs by 1, producing a fixing cut. const bool has_pair = complement_pairs.size() == 1; i_t num_complements = 0; - for (const auto vertex_idx : clique_vertices) { + for (const i_t vertex_idx : clique_vertices) { const i_t var_idx = vertex_idx % num_vars; const bool complement = vertex_idx >= num_vars; if (has_pair && complement_pairs.count(var_idx) > 0) { continue; } @@ -229,7 +246,7 @@ inline size_t bitset_words(size_t n) { return (n + 63) / 64; } inline bool bitset_any(const std::vector& bs) { - for (auto word : bs) { + for (const uint64_t word : bs) { if (word != 0) { return true; } } return false; @@ -343,7 +360,7 @@ void bron_kerbosch(bk_bitset_context_t& ctx, ctx.add_work(static_cast(2 * ctx.words + num_candidates)); ctx.add_work(static_cast(num_candidates) * static_cast(7 * ctx.words + 6)); // note that candidates will include pivot if it is in P - for (auto v : candidates) { + for (const i_t v : candidates) { if (ctx.over_call_limit()) { ctx.call_limit_reached = true; return; @@ -370,6 +387,85 @@ void bron_kerbosch(bk_bitset_context_t& ctx, } } +// ---- Shared helpers for greedy CG-based set extension (clique & odd-wheel) ---- + +// Pick the seed vertex with the smallest conflict-graph degree. Returns -1 if +// the seed is empty or the time limit is hit while scanning. +template +i_t min_degree_anchor(const std::vector& seed, + mip::clique_table_t& graph, + f_t start_time, + f_t time_limit) +{ + i_t smallest_degree = std::numeric_limits::max(); + i_t smallest_degree_var = -1; + for (const i_t v : seed) { + if (toc(start_time) >= time_limit) { return -1; } + i_t degree = graph.get_degree_of_var(v); + if (degree < smallest_degree) { + smallest_degree = degree; + smallest_degree_var = v; + } + } + return smallest_degree_var; +} + +// Reduced-cost key for a CG vertex. A complement vertex (idx >= num_vars) maps +// to the original variable and flips the sign. Sorting candidates by this key +// keeps xstar minimally disturbed so the resulting cut stays binding and the +// dual simplex resolve stays cheap. +template +f_t cg_reduced_cost(i_t vertex_idx, const std::vector& reduced_costs, i_t num_vars) +{ + i_t var_idx = vertex_idx % num_vars; + cuopt_assert(var_idx >= 0 && var_idx < static_cast(reduced_costs.size()), + "Reduced cost index out of range"); + f_t rc = reduced_costs[var_idx]; + if (!std::isfinite(rc)) { rc = 0.0; } + return vertex_idx >= num_vars ? -rc : rc; +} + +template +void sort_candidates_by_reduced_cost(std::vector& candidates, + const std::vector& reduced_costs, + i_t num_vars) +{ + std::sort(candidates.begin(), candidates.end(), [&](i_t a, i_t b) { + return cg_reduced_cost(a, reduced_costs, num_vars) < + cg_reduced_cost(b, reduced_costs, num_vars); + }); +} + +// Greedily grow `selected` by appending candidates (assumed already ordered by +// reduced cost) that are adjacent to every current member of `selected`. The +// resulting `selected` is therefore a clique. Stops early when the time or work +// budget is exhausted. +template +void greedy_extend_clique(std::vector& selected, + const std::vector& candidates, + mip::clique_table_t& graph, + f_t adj_check_cost, + f_t start_time, + f_t time_limit, + f_t* work_estimate, + f_t max_work_estimate) +{ + for (const i_t candidate : candidates) { + if (toc(start_time) >= time_limit) { return; } + bool add = true; + i_t checks = 0; + for (const i_t v : selected) { + checks++; + if (!graph.check_adjacency(candidate, v)) { + add = false; + break; + } + } + if (add_work_estimate(adj_check_cost * checks, work_estimate, max_work_estimate)) { break; } + if (add) { selected.push_back(candidate); } + } +} + template void extend_clique_vertices(std::vector& clique_vertices, mip::clique_table_t& graph, @@ -391,35 +487,25 @@ void extend_clique_vertices(std::vector& clique_vertices, static_cast(clique_vertices.size())); const f_t initial_clique_size = static_cast(clique_vertices.size()); - i_t smallest_degree = std::numeric_limits::max(); - i_t smallest_degree_var = -1; - for (auto v : clique_vertices) { - if (toc(start_time) >= time_limit) { return; } - i_t degree = graph.get_degree_of_var(v); - if (degree < smallest_degree) { - smallest_degree = degree; - smallest_degree_var = v; - } - } + const i_t smallest_degree_var = min_degree_anchor(clique_vertices, graph, start_time, time_limit); + if (smallest_degree_var < 0) { return; } - auto adj_set = graph.get_adj_set_of_var(smallest_degree_var); + std::unordered_set adj_set = graph.get_adj_set_of_var(smallest_degree_var); std::unordered_set clique_members(clique_vertices.begin(), clique_vertices.end()); std::vector candidates; candidates.reserve(adj_set.size()); // the candidate list if only the integer valued vertices - for (const auto& candidate : adj_set) { + for (const i_t& candidate : adj_set) { if (toc(start_time) >= time_limit) { return; } if (clique_members.count(candidate) != 0) { continue; } i_t var_idx = candidate % num_vars; f_t value = candidate >= num_vars ? (1.0 - xstar[var_idx]) : xstar[var_idx]; if (std::abs(value - std::round(value)) <= integer_tol) { candidates.push_back(candidate); } } - CLIQUE_CUTS_DEBUG( - "extend_clique_vertices anchor=%lld degree=%lld adj_size=%lld integer_candidates=%lld", - static_cast(smallest_degree_var), - static_cast(smallest_degree), - static_cast(adj_set.size()), - static_cast(candidates.size())); + CLIQUE_CUTS_DEBUG("extend_clique_vertices anchor=%lld adj_size=%lld integer_candidates=%lld", + static_cast(smallest_degree_var), + static_cast(adj_set.size()), + static_cast(candidates.size())); const f_t candidate_size = static_cast(candidates.size()); const f_t sort_work = candidate_size > 0.0 ? 2.0 * candidate_size * std::log2(candidate_size + 1.0) : 0.0; @@ -443,44 +529,481 @@ void extend_clique_vertices(std::vector& clique_vertices, // less refactors and less iterations after resolve. // it also increases the cut's effectiveness by keeping xstar not disturbed much // if it is disturbed too much, the cut might become non-binding - auto reduced_cost = [&](i_t vertex_idx) -> f_t { - i_t var_idx = vertex_idx % num_vars; - cuopt_assert(var_idx >= 0 && var_idx < static_cast(reduced_costs.size()), - "Variable index out of range"); - f_t rc = reduced_costs[var_idx]; - if (!std::isfinite(rc)) { rc = 0.0; } - return vertex_idx >= num_vars ? -rc : rc; + sort_candidates_by_reduced_cost(candidates, reduced_costs, num_vars); + + // adj_check_cost folds in addtl_cliques_scan_cost so each check_adjacency + // charges its own addtl scan cost as the clique grows. + greedy_extend_clique(clique_vertices, + candidates, + graph, + adj_check_cost, + start_time, + time_limit, + work_estimate, + max_work_estimate); + CLIQUE_CUTS_DEBUG("extend_clique_vertices done start=%lld final=%lld added=%lld", + static_cast(initial_clique_vertices), + static_cast(clique_vertices.size()), + static_cast(clique_vertices.size() - initial_clique_vertices)); +} + +// Build a zero-half (odd-cycle / odd-wheel) cut from a cycle and optional wheel +// centers. cycle_vertices is a simple odd cycle in the conflict graph using the +// 2*num_vars vertex indexing (var j and complement j+num_vars). wheel_centers +// are extra vertices each adjacent to every vertex in cycle_vertices. The +// resulting cut is stored in the form a^T x >= rhs to match cut_pool_t. +template +clique_cut_build_status_t build_zero_half_cut(const std::vector& cycle_vertices, + const std::vector& wheel_centers, + i_t num_vars, + const std::vector& var_types, + const std::vector& lower_bounds, + const std::vector& upper_bounds, + const std::vector& xstar, + f_t bound_tol, + f_t min_violation, + sparse_vector_t& cut, + f_t& cut_rhs, + f_t* work_estimate, + f_t max_work_estimate) +{ + const size_t cycle_size = cycle_vertices.size(); + if (cycle_size < 5 || (cycle_size % 2) == 0) { + ZERO_HALF_DEBUG("build_zero_half_cut reject cycle_size=%zu", cycle_size); + return clique_cut_build_status_t::NO_CUT; + } + cuopt_assert(num_vars > 0, "Zero-half cut num_vars must be positive"); + cuopt_assert(static_cast(num_vars) <= lower_bounds.size(), + "Zero-half cut lower bounds size mismatch"); + cuopt_assert(static_cast(num_vars) <= xstar.size(), "Zero-half cut xstar size mismatch"); + + const i_t m = static_cast((cycle_size - 1) / 2); + const f_t f_m = static_cast(m); + // The guard above rejects even or <5 cycles, so the cycle decomposes as + // exactly 2m+1 literals with m >= 2. The whole zero-half lift (rhs = -m, + // unit cycle coefficients, m-weighted wheel centers) depends on this. + cuopt_assert(2 * m + 1 == static_cast(cycle_size), + "Zero-half cut: cycle_size must equal 2m+1 (odd cycle)"); + cuopt_assert(m >= 2, "Zero-half cut: odd cycle must have length >= 5 (m >= 2)"); + const f_t total_size = static_cast(cycle_size + wheel_centers.size()); + const f_t estimated_work = 8.0 * total_size + 2.0 * total_size * std::log2(total_size + 1.0); + if (add_work_estimate(estimated_work, work_estimate, max_work_estimate)) { + ZERO_HALF_DEBUG("build_zero_half_cut work_limit hit"); + return clique_cut_build_status_t::NO_CUT; + } + + cut.i.clear(); + cut.x.clear(); + + std::unordered_map coeff_by_var; + std::unordered_set seen_original; + std::unordered_set seen_complement; + coeff_by_var.reserve(cycle_size + wheel_centers.size()); + seen_original.reserve(cycle_size + wheel_centers.size()); + seen_complement.reserve(cycle_size + wheel_centers.size()); + + f_t rhs_acc = -f_m; + + auto accumulate = + [&](const std::vector& verts, f_t weight, bool is_cycle) -> clique_cut_build_status_t { + ZERO_HALF_DEBUG("build_zero_half_cut accumulate verts.size=%zu weight=%g is_cycle=%d", + verts.size(), + static_cast(weight), + static_cast(is_cycle)); + for (const i_t vertex_idx : verts) { + ZERO_HALF_DEBUG(" acc vertex_idx=%lld (range [0, %lld))", + static_cast(vertex_idx), + static_cast(2 * num_vars)); + cuopt_assert(vertex_idx >= 0 && vertex_idx < 2 * num_vars, "Zero-half vertex out of range"); + const i_t var_idx = vertex_idx % num_vars; + const bool complement = vertex_idx >= num_vars; + const f_t lower_bound = lower_bounds[var_idx]; + const f_t upper_bound = upper_bounds[var_idx]; + cuopt_assert(var_types[var_idx] != variable_type_t::CONTINUOUS, + "Zero-half cut contains continuous variable"); + cuopt_assert(lower_bound >= -bound_tol, "Zero-half variable lower bound below zero"); + cuopt_assert(upper_bound <= 1 + bound_tol, "Zero-half variable upper bound above one"); + + if (complement) { + if (seen_original.count(var_idx) > 0) { return clique_cut_build_status_t::NO_CUT; } + seen_complement.insert(var_idx); + coeff_by_var[var_idx] += weight; + rhs_acc += weight; + } else { + if (seen_complement.count(var_idx) > 0) { return clique_cut_build_status_t::NO_CUT; } + seen_original.insert(var_idx); + coeff_by_var[var_idx] -= weight; + } + } + return clique_cut_build_status_t::CUT_ADDED; }; - std::sort(candidates.begin(), candidates.end(), [&](i_t a, i_t b) { - return reduced_cost(a) < reduced_cost(b); - }); + if (accumulate(cycle_vertices, static_cast(1), true) != + clique_cut_build_status_t::CUT_ADDED) { + ZERO_HALF_DEBUG("build_zero_half_cut cycle accumulate failed"); + return clique_cut_build_status_t::NO_CUT; + } + if (m > 0 && !wheel_centers.empty()) { + if (accumulate(wheel_centers, f_m, false) != clique_cut_build_status_t::CUT_ADDED) { + ZERO_HALF_DEBUG("build_zero_half_cut wheel accumulate failed"); + return clique_cut_build_status_t::NO_CUT; + } + } - for (const auto candidate : candidates) { - bool add = true; - i_t checks = 0; - for (const auto v : clique_vertices) { - checks++; + const f_t coeff_zero_tol = static_cast(1e-12); + cut.i.reserve(coeff_by_var.size()); + cut.x.reserve(coeff_by_var.size()); + for (const std::pair& kv : coeff_by_var) { + if (std::abs(kv.second) <= coeff_zero_tol) { continue; } + // Each variable appears at most once on the cycle (contributing +/-1) and + // at most once among the wheel centers (contributing +/-m), so no final + // coefficient can exceed 1 + m in magnitude. A larger value means a vertex + // was double-counted in accumulation. + cuopt_assert(std::abs(kv.second) <= f_m + static_cast(1) + bound_tol, + "Zero-half coefficient exceeds 1 + m (vertex double-counted?)"); + cut.i.push_back(kv.first); + cut.x.push_back(kv.second); + } + // Support is bounded by the number of distinct accumulated vertices. + cuopt_assert(cut.i.size() <= cycle_size + wheel_centers.size(), + "Zero-half cut support exceeds accumulated vertex count"); + + if (cut.i.empty()) { + ZERO_HALF_DEBUG("build_zero_half_cut empty support after accumulation"); + return clique_cut_build_status_t::NO_CUT; + } + + cut_rhs = rhs_acc; + cut.sort(); + + const f_t dot = cut.dot(xstar); + const f_t violation = cut_rhs - dot; + ZERO_HALF_DEBUG( + "build_zero_half_cut nz=%lld rhs=%g dot=%g violation=%g threshold=%g cycle=%lld wheel=%lld", + static_cast(cut.i.size()), + static_cast(cut_rhs), + static_cast(dot), + static_cast(violation), + static_cast(min_violation), + static_cast(cycle_size), + static_cast(wheel_centers.size())); + cuopt_assert(violation > -bound_tol, "Zero-half cut violation flipped sign unexpectedly"); + if (violation > min_violation) { return clique_cut_build_status_t::CUT_ADDED; } + return clique_cut_build_status_t::NO_CUT; +} + +// Reusable scratch for dijkstra_odd_cycle. The separation loop runs Dijkstra +// once per source vertex; re-allocating and re-initializing dist/prev. Instead +// we allocate the buffers once and reset them in O(1) using a generation stamp: +// dist[v]/prev[v] are considered valid for the current call only when +// stamp[v] == gen. +template +struct dijkstra_scratch_t { + std::vector dist; + std::vector prev; + std::vector stamp; // stamp[v] == gen <=> dist[v]/prev[v] valid this call + std::uint64_t gen{0}; + + void ensure_size(std::size_t n) + { + if (stamp.size() < n) { + dist.resize(n); + prev.resize(n); + stamp.assign(n, 0); + gen = 0; + } + } +}; + +template +bool dijkstra_odd_cycle(i_t source_local, + const std::vector>& local_adj, + const std::vector& weights, + f_t cutoff, + std::vector& path, + f_t& total_weight, + f_t* work_estimate, + f_t max_work_estimate, + dijkstra_scratch_t& scratch) +{ + const i_t num_local = static_cast(local_adj.size()); + if (source_local < 0 || source_local >= num_local) { return false; } + if (weights.size() != static_cast(num_local)) { return false; } + cuopt_assert(source_local >= 0 && source_local < num_local, + "Zero-half Dijkstra source out of range"); + cuopt_assert(weights.size() == static_cast(num_local), + "Zero-half Dijkstra weights size mismatch"); + + const i_t source_idx = source_local; + const i_t target_idx = source_local + num_local; + const i_t total_idx = 2 * num_local; + const f_t f_inf = std::numeric_limits::infinity(); + + scratch.ensure_size(static_cast(total_idx)); + ++scratch.gen; + const std::uint64_t gen = scratch.gen; + std::vector& dist = scratch.dist; + std::vector& prev = scratch.prev; + std::vector& stamp = scratch.stamp; + // dist[v]/prev[v] are valid only if last written this call (stamp[v] == gen); + // otherwise the node is unreached, i.e. distance infinity. + auto cur_dist = [&](i_t v) -> f_t { return stamp[v] == gen ? dist[v] : f_inf; }; + + dist[source_idx] = 0; + prev[source_idx] = -1; + stamp[source_idx] = gen; + + using node_t = std::pair; + std::priority_queue, std::greater> pq; + pq.emplace(static_cast(0), source_idx); + + i_t pops = 0; + while (!pq.empty()) { + auto [d, u] = pq.top(); + pq.pop(); + ++pops; + if (d > cur_dist(u)) { continue; } + if (u == target_idx) { break; } + if (cutoff > 0 && d >= cutoff) { break; } + + const i_t u_local = u % num_local; + const i_t u_part = u / num_local; + const i_t v_part = 1 - u_part; + cuopt_assert(u_part == 0 || u_part == 1, "Bipartite part out of range"); + + const std::vector& neigh = local_adj[u_local]; + if (add_work_estimate(static_cast(neigh.size()) + 4.0, work_estimate, max_work_estimate)) { + ZERO_HALF_DEBUG("dijkstra_odd_cycle work_limit hit pops=%lld", static_cast(pops)); + return false; + } + for (const i_t v_local : neigh) { + cuopt_assert(v_local >= 0 && v_local < num_local, "Zero-half Dijkstra neighbor out of range"); + f_t edge_w = (static_cast(1) - weights[u_local] - weights[v_local]) / 2; + if (edge_w < 0) { edge_w = 0; } + const i_t v = v_local + v_part * num_local; + const f_t nd = d + edge_w; + if (nd < cur_dist(v)) { + dist[v] = nd; + prev[v] = u; + stamp[v] = gen; + pq.emplace(nd, v); + } + } + } + + const f_t target_dist = cur_dist(target_idx); + if (!std::isfinite(target_dist)) { + ZERO_HALF_DEBUG("dijkstra_odd_cycle no path pops=%lld", static_cast(pops)); + return false; + } + total_weight = target_dist; + // All G' edge weights are clamped to >= 0, so the shortest-path distance must + // be non-negative; a negative total means the clamp/relaxation invariant broke. + cuopt_assert(total_weight >= -static_cast(1e-9), + "Zero-half Dijkstra shortest-path distance must be non-negative"); + if (cutoff > 0 && total_weight >= cutoff) { + ZERO_HALF_DEBUG("dijkstra_odd_cycle path too long total=%g cutoff=%g", + static_cast(total_weight), + static_cast(cutoff)); + return false; + } + + path.clear(); + for (i_t cur = target_idx; cur != -1; cur = prev[cur]) { + path.push_back(cur); + if (cur == source_idx) { break; } + } + cuopt_assert(!path.empty(), "Zero-half Dijkstra path empty"); + cuopt_assert(path.back() == source_idx, "Zero-half Dijkstra path missing source"); + std::reverse(path.begin(), path.end()); + cuopt_assert(path.front() == source_idx, "Zero-half Dijkstra path must start at source"); + cuopt_assert(path.back() == target_idx, "Zero-half Dijkstra path must end at target"); + // bipartite path from j1 to j2 must have odd number of edges + cuopt_assert((path.size() % 2) == 0, "Zero-half bipartite path must have even node count"); +#ifdef ASSERT_MODE + // Every G' edge crosses between the two bipartite copies, so consecutive path + // nodes must live in opposite parts (part = bipartite_idx / num_local). + for (size_t k = 0; k + 1 < path.size(); ++k) { + cuopt_assert((path[k] / num_local) != (path[k + 1] / num_local), + "Zero-half Dijkstra path must alternate bipartite parts"); + } +#endif + ZERO_HALF_DEBUG("dijkstra_odd_cycle done path.size=%zu total_weight=%g pops=%lld", + path.size(), + static_cast(total_weight), + static_cast(pops)); + return true; +} + +template +bool path_to_odd_cycle(const std::vector& bipartite_path, + const std::vector& vertices, + i_t num_local, + i_t num_vars, + std::vector& cycle_vertices, + f_t* work_estimate, + f_t max_work_estimate) +{ + ZERO_HALF_DEBUG( + "path_to_odd_cycle enter bipartite_path.size=%zu vertices.size=%zu num_local=%lld " + "num_vars=%lld", + bipartite_path.size(), + vertices.size(), + static_cast(num_local), + static_cast(num_vars)); + cycle_vertices.clear(); + if (bipartite_path.size() < 4) { + ZERO_HALF_DEBUG("path_to_odd_cycle reject short path"); + return false; + } + if (add_work_estimate( + static_cast(bipartite_path.size()) * 2.0, work_estimate, max_work_estimate)) { + ZERO_HALF_DEBUG("path_to_odd_cycle work_limit hit"); + return false; + } + + std::vector local_seq; + local_seq.reserve(bipartite_path.size()); + for (const i_t bv : bipartite_path) { + local_seq.push_back(bv % num_local); + } + cuopt_assert(local_seq.front() == local_seq.back(), "Zero-half cycle path endpoints must match"); + + // Drop the duplicate end so we have a sequence covering each cycle vertex once + local_seq.pop_back(); + + std::unordered_set seen_local; + seen_local.reserve(local_seq.size()); + for (const i_t lv : local_seq) { + if (!seen_local.insert(lv).second) { + // Same CG vertex appears twice in the path; reject (degenerate cycle) + ZERO_HALF_DEBUG("path_to_odd_cycle duplicate local vertex lv=%lld", + static_cast(lv)); + return false; + } + } + + cycle_vertices.reserve(local_seq.size()); + std::unordered_set seen_var; + seen_var.reserve(local_seq.size()); + for (const i_t lv : local_seq) { + const i_t global = vertices[lv]; + cuopt_assert(global >= 0 && global < 2 * num_vars, "Zero-half global vertex out of range"); + const i_t var_idx = global % num_vars; + if (!seen_var.insert(var_idx).second) { + // Variable appears as both x and ¯x in the cycle; reject (degenerate) + ZERO_HALF_DEBUG("path_to_odd_cycle duplicate var_idx=%lld", static_cast(var_idx)); + return false; + } + cycle_vertices.push_back(global); + } + cuopt_assert(cycle_vertices.size() == local_seq.size(), + "Zero-half cycle dropped vertices during global mapping"); + cuopt_assert((cycle_vertices.size() % 2) == 1, "Zero-half extracted cycle must have odd length"); + ZERO_HALF_DEBUG("path_to_odd_cycle done cycle_vertices.size=%zu", cycle_vertices.size()); + return cycle_vertices.size() >= 5; +} + +// Greedy lifting: extend an odd cycle by attaching a clique of "wheel center" +// vertices that are adjacent (in CG) to every vertex of the cycle. +template +void extend_to_odd_wheel(const std::vector& cycle_vertices, + std::vector& wheel_centers, + mip::clique_table_t& graph, + const std::vector& reduced_costs, + i_t num_vars, + f_t start_time, + f_t time_limit, + f_t* work_estimate, + f_t max_work_estimate) +{ + ZERO_HALF_DEBUG( + "extend_to_odd_wheel enter cycle.size=%zu num_vars=%lld reduced_costs.size=%zu " + "graph.n_variables=%lld", + cycle_vertices.size(), + static_cast(num_vars), + reduced_costs.size(), + static_cast(graph.n_variables)); + wheel_centers.clear(); + if (cycle_vertices.empty()) { return; } + if (toc(start_time) >= time_limit) { return; } + + const i_t smallest_degree_var = min_degree_anchor(cycle_vertices, graph, start_time, time_limit); + ZERO_HALF_DEBUG("extend_to_odd_wheel smallest_degree_var=%lld", + static_cast(smallest_degree_var)); + if (smallest_degree_var < 0) { return; } + + std::unordered_set adj_set = graph.get_adj_set_of_var(smallest_degree_var); + ZERO_HALF_DEBUG("extend_to_odd_wheel adj_set.size=%zu", adj_set.size()); + std::vector cycle_members(2 * num_vars, 0); + for (const i_t v : cycle_vertices) { + cuopt_assert(v >= 0 && v < 2 * num_vars, "Zero-half cycle vertex out of range"); + cycle_members[v] = 1; + } + std::vector candidates; + candidates.reserve(adj_set.size()); + for (const i_t candidate : adj_set) { + if (toc(start_time) >= time_limit) { return; } + if (cycle_members[candidate] != 0) { continue; } + bool adj_to_all = true; + for (const i_t v : cycle_vertices) { + if (candidate == v) { + adj_to_all = false; + break; + } if (!graph.check_adjacency(candidate, v)) { - add = false; + adj_to_all = false; break; } } - // Each check_adjacency now charges its own addtl_cliques_scan_cost - // term so the per-iteration budget reflects the addtl scan cost. - if (add_work_estimate( - adj_check_cost * static_cast(checks), work_estimate, max_work_estimate)) { - break; + if (adj_to_all) { candidates.push_back(candidate); } + } + ZERO_HALF_DEBUG("extend_to_odd_wheel candidates.size=%zu", candidates.size()); + if (candidates.empty()) { return; } + + const f_t candidate_size = static_cast(candidates.size()); + const f_t cycle_size_f = static_cast(cycle_vertices.size()); + const f_t adj_set_cost = 2.0 * static_cast(adj_set.size()); + const f_t sort_cost = + candidate_size > 0.0 ? 2.0 * candidate_size * std::log2(candidate_size + 1.0) : 0.0; + if (add_work_estimate(adj_set_cost + cycle_size_f * candidate_size + sort_cost, + work_estimate, + max_work_estimate)) { + ZERO_HALF_DEBUG("extend_to_odd_wheel work_limit hit pre-sort"); + return; + } + + sort_candidates_by_reduced_cost(candidates, reduced_costs, num_vars); + + // Candidates are already adjacent to every cycle vertex (filtered above), so + // growing a clique among them yields centers adjacent to the whole cycle and + // to each other. + const f_t adj_check_cost = 5.0; + greedy_extend_clique(wheel_centers, + candidates, + graph, + adj_check_cost, + start_time, + time_limit, + work_estimate, + max_work_estimate); +#ifdef ASSERT_MODE + // Post-condition: the selected centers must form a clique that is fully + // adjacent to the cycle — each center adjacent to every cycle vertex and to + // every other center. This is exactly what makes the m-weighted wheel lift a + // valid zero-half inequality. + for (size_t a = 0; a < wheel_centers.size(); ++a) { + for (const i_t cv : cycle_vertices) { + cuopt_assert(graph.check_adjacency(wheel_centers[a], cv), + "Zero-half wheel center not adjacent to every cycle vertex"); } - if (add) { - clique_vertices.push_back(candidate); - clique_members.insert(candidate); + for (size_t b = a + 1; b < wheel_centers.size(); ++b) { + cuopt_assert(graph.check_adjacency(wheel_centers[a], wheel_centers[b]), + "Zero-half wheel centers must be mutually adjacent (clique)"); } } - CLIQUE_CUTS_DEBUG("extend_clique_vertices done start=%lld final=%lld added=%lld", - static_cast(initial_clique_vertices), - static_cast(clique_vertices.size()), - static_cast(clique_vertices.size() - initial_clique_vertices)); +#endif + ZERO_HALF_DEBUG("extend_to_odd_wheel done wheel_centers.size=%zu", wheel_centers.size()); } } // namespace @@ -510,7 +1033,7 @@ std::vector> find_maximal_cliques_for_test( const size_t words = bitset_words(n_vertices); std::vector> adj_bitset(n_vertices, std::vector(words, 0)); for (size_t v = 0; v < n_vertices; ++v) { - for (const auto& nbr : adjacency_list[v]) { + for (const int& nbr : adjacency_list[v]) { cuopt_assert(nbr >= 0 && static_cast(nbr) < n_vertices, "Neighbor index out of range in clique test helper"); bitset_set(adj_bitset[v], static_cast(nbr)); @@ -541,6 +1064,95 @@ std::vector> find_maximal_cliques_for_test( return ctx.cliques; } +// This function is only used in tests +std::vector> find_violated_odd_cycles_for_test( + const std::vector>& adjacency_list, + const std::vector& x_values, + double min_violation, + double time_limit) +{ + const size_t n_vertices = adjacency_list.size(); + if (n_vertices == 0) { return {}; } + cuopt_assert(x_values.size() == n_vertices, "x_values size mismatch in odd-cycle test helper"); + + const int num_local = static_cast(n_vertices); + std::vector> adj_local(n_vertices); + for (size_t v = 0; v < n_vertices; ++v) { + adj_local[v].reserve(adjacency_list[v].size()); + for (const int nbr : adjacency_list[v]) { + cuopt_assert(nbr >= 0 && static_cast(nbr) < n_vertices, + "Neighbor index out of range in odd-cycle test helper"); + adj_local[v].push_back(nbr); + } + } + + double work_estimate = 0.0; + const double max_work_estimate = std::numeric_limits::infinity(); + const double start_time = tic(); + const double cutoff = 0.5 - min_violation; + + std::vector> result; + std::vector bipartite_path; + dijkstra_scratch_t dijkstra_scratch; + + for (int s = 0; s < num_local; ++s) { + if (toc(start_time) >= time_limit) { break; } + + double total_weight = 0; + if (!dijkstra_odd_cycle(s, + adj_local, + x_values, + cutoff, + bipartite_path, + total_weight, + &work_estimate, + max_work_estimate, + dijkstra_scratch)) { + continue; + } + if (bipartite_path.size() < 4) { continue; } + std::vector seq; + seq.reserve(bipartite_path.size()); + for (const int bv : bipartite_path) { + seq.push_back(bv % num_local); + } + cuopt_assert(seq.front() == seq.back(), "Odd-cycle test helper path endpoints must match"); + seq.pop_back(); + if ((seq.size() % 2) == 0 || seq.size() < 5) { continue; } + bool simple = true; + std::unordered_set seen; + seen.reserve(seq.size()); + for (const int v : seq) { + if (!seen.insert(v).second) { + simple = false; + break; + } + } + if (!simple) { continue; } + result.push_back(seq); + } + return result; +} + +namespace { + +// 64-bit integer mixer (SplitMix64). Used as the building block for the +// cousin filter's per-slot independent hash family. +inline uint64_t splitmix64_mix(uint64_t x) +{ + x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; + x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; + x = x ^ (x >> 31); + return x; +} + +inline uint64_t hash64_with_seed(uint64_t value, uint64_t seed) +{ + return splitmix64_mix(value ^ (seed * 0xbf58476d1ce4e5b9ULL + 0x9e3779b97f4a7c15ULL)); +} + +} // namespace + template void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { @@ -2577,6 +3189,148 @@ void cut_generation_t::generate_implied_bound_cuts( } } +template +void cut_generation_t::prepare_fractional_sub_conflict_graph( + const simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time) +{ + sub_cg_.clear(); + + if (settings.clique_cuts == 0 && settings.zero_half_cuts == 0) { return; } + if (toc(start_time) >= settings.time_limit) { return; } + + // The clique table is produced by a background OpenMP task (spawned in + // branch_and_bound). Its base cliques are published early, but the extension + // phase keeps mutating the same table. Now that cut generation needs the + // table, signal the task to stop extending and block until it finishes, so the + // table is no longer mutated while we read it below. The wait must happen + // unconditionally (even when the pointer is already non-null), otherwise we + // would race the still-running extension. + if (signal_extend_) { + signal_extend_->store(true, std::memory_order_release); +#pragma omp taskwait depend(in : *signal_extend_) + } + + if (clique_table_ == nullptr) { return; } + // small_clique_adj may carry pairwise CG edges from cliques demoted by + // remove_small_cliques; clique_table_t::empty() accounts for that. + if (clique_table_->empty()) { return; } + + const i_t num_vars = user_problem_.num_cols; + cuopt_assert(clique_table_->n_variables == num_vars, + "prepare_fractional_sub_conflict_graph clique table variable count mismatch"); + cuopt_assert(static_cast(num_vars) <= xstar.size(), + "prepare_fractional_sub_conflict_graph xstar size mismatch"); + cuopt_assert(user_problem_.var_types.size() == static_cast(num_vars), + "prepare_fractional_sub_conflict_graph user problem var_types size mismatch"); + + const f_t bound_tol = settings.primal_tol; + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e7; + + sub_cg_.num_vars = num_vars; + sub_cg_.vertices.reserve(static_cast(num_vars) * 2); + sub_cg_.weights.reserve(static_cast(num_vars) * 2); + + for (i_t j = 0; j < num_vars; ++j) { + if (user_problem_.var_types[j] == variable_type_t::CONTINUOUS) { continue; } + const f_t lower_bound = user_problem_.lower[j]; + const f_t upper_bound = user_problem_.upper[j]; + if (lower_bound < -bound_tol || upper_bound > 1 + bound_tol) { continue; } + const f_t xj = xstar[j]; + if (std::abs(xj - std::round(xj)) <= settings.integer_tol) { continue; } + sub_cg_.vertices.push_back(j); + sub_cg_.weights.push_back(xj); + sub_cg_.vertices.push_back(j + num_vars); + sub_cg_.weights.push_back(static_cast(1.0) - xj); + } + work_estimate += + 4.0 * static_cast(num_vars) + 2.0 * static_cast(sub_cg_.vertices.size()); + if (work_estimate > max_work_estimate) { + sub_cg_.clear(); + return; + } + + if (sub_cg_.vertices.empty()) { + // No fractional binaries — both separators have nothing to do, but the + // build itself succeeded. Mark ready so callers don't keep retrying. + sub_cg_.ready = true; + return; + } + + sub_cg_.vertex_to_local.assign(static_cast(2 * num_vars), -1); + sub_cg_.in_subgraph.assign(static_cast(2 * num_vars), 0); + for (size_t idx = 0; idx < sub_cg_.vertices.size(); ++idx) { + if (toc(start_time) >= settings.time_limit) { + sub_cg_.clear(); + return; + } + const i_t v_idx = sub_cg_.vertices[idx]; + sub_cg_.vertex_to_local[v_idx] = static_cast(idx); + sub_cg_.in_subgraph[v_idx] = 1; + } + work_estimate += 3.0 * static_cast(sub_cg_.vertices.size()); + if (work_estimate > max_work_estimate) { + sub_cg_.clear(); + return; + } + + sub_cg_.adj_local.assign(sub_cg_.vertices.size(), {}); + size_t total_adj_entries = 0; + size_t kept_adj_entries = 0; + for (size_t idx = 0; idx < sub_cg_.vertices.size(); ++idx) { + if (toc(start_time) >= settings.time_limit) { + sub_cg_.clear(); + return; + } + const i_t v_idx = sub_cg_.vertices[idx]; + std::unordered_set adj_set = clique_table_->get_adj_set_of_var(v_idx); + total_adj_entries += adj_set.size(); + std::vector& adj = sub_cg_.adj_local[idx]; + adj.reserve(adj_set.size()); + for (const i_t neighbor : adj_set) { + cuopt_assert(neighbor >= 0 && neighbor < 2 * num_vars, + "prepare_fractional_sub_conflict_graph neighbor out of range"); + if (!sub_cg_.in_subgraph[neighbor]) { continue; } + const i_t local_neighbor = sub_cg_.vertex_to_local[neighbor]; + cuopt_assert(local_neighbor >= 0, + "prepare_fractional_sub_conflict_graph local_neighbor out of range"); + adj.push_back(local_neighbor); + } + kept_adj_entries += adj.size(); +#ifdef ASSERT_MODE + { + std::unordered_set adj_global; + adj_global.reserve(adj.size()); + for (const i_t neighbor : adj) { + const i_t v = sub_cg_.vertices[neighbor]; + cuopt_assert(adj_global.insert(v).second, + "Duplicate neighbor in fractional sub-CG adjacency list"); + } + } +#endif + } + work_estimate += static_cast(sub_cg_.vertices.size()) + static_cast(total_adj_entries) + + 2.0 * static_cast(kept_adj_entries); + if (work_estimate > max_work_estimate) { + sub_cg_.clear(); + return; + } + + sub_cg_.ready = true; + CLIQUE_CUTS_DEBUG( + "prepare_fractional_sub_conflict_graph ready vertices=%lld raw_adj=%lld kept_adj=%lld", + static_cast(sub_cg_.vertices.size()), + static_cast(total_adj_entries), + static_cast(kept_adj_entries)); + ZERO_HALF_DEBUG( + "prepare_fractional_sub_conflict_graph ready vertices=%lld raw_adj=%lld kept_adj=%lld", + static_cast(sub_cg_.vertices.size()), + static_cast(total_adj_entries), + static_cast(kept_adj_entries)); +} + template bool cut_generation_t::generate_cuts(const lp_problem_t& lp, const simplex_solver_settings_t& settings, @@ -2633,6 +3387,24 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } + // Generate implied bound cuts + if (settings.implied_bound_cuts != 0) { + f_t cut_start_time = tic(); + generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); + f_t cut_generation_time = toc(cut_start_time); + if (cut_generation_time > 1.0) { + settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); + } + } + + // Build the fractional conflict-graph subgraph once (resolving the async + // clique-table future on the way) so both clique-cut and zero-half cut + // separators consume the same vertex/weight/adjacency tables instead of + // each recomputing them. Done here, after the cut routines that don't + // need the clique table, to give the background clique-table thread as + // much time as possible to finish before we join it. + prepare_fractional_sub_conflict_graph(settings, xstar, start_time); + // Generate Clique cuts (last to give background clique table generation maximum time) if (settings.clique_cuts != 0) { f_t cut_start_time = tic(); @@ -2647,14 +3419,24 @@ bool cut_generation_t::generate_cuts(const lp_problem_t& lp, } } - // Generate implied bound cuts - if (settings.implied_bound_cuts != 0) { + // Generate Zero-half (odd-cycle / odd-wheel) cuts; reuses the clique table built above + if (settings.zero_half_cuts != 0) { + ZERO_HALF_DEBUG("generate_cuts: about to call generate_zero_half_cuts"); f_t cut_start_time = tic(); - generate_implied_bound_cuts(lp, settings, var_types, xstar, start_time); + bool feasible = generate_zero_half_cuts(lp, settings, var_types, xstar, zstar, start_time); + ZERO_HALF_DEBUG("generate_cuts: returned from generate_zero_half_cuts feasible=%d", + static_cast(feasible)); + if (!feasible) { + settings.log.printf("Zero-half cuts proved infeasible\n"); + return false; + } f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { - settings.log.debug("Implied bounds cut generation time %.2f seconds\n", cut_generation_time); + settings.log.debug("Zero-half cut generation time %.2f seconds\n", cut_generation_time); } + } else { + ZERO_HALF_DEBUG("generate_cuts: zero_half_cuts disabled (setting=%d)", + static_cast(settings.zero_half_cuts)); } return true; } @@ -2719,32 +3501,22 @@ bool cut_generation_t::generate_clique_cuts( static_cast(settings.time_limit), static_cast(toc(start_time))); - if (clique_table_ == nullptr) { - CLIQUE_CUTS_DEBUG("generate_clique_cuts signaling background thread and waiting"); - if (signal_extend_) { signal_extend_->store(true, std::memory_order_release); } -#pragma omp taskwait depend(in : *signal_extend_) - if (clique_table_) { - CLIQUE_CUTS_DEBUG("generate_clique_cuts received clique table first=%lld addtl=%lld", - static_cast(clique_table_->first.size()), - static_cast(clique_table_->addtl_cliques.size())); - } - } - - if (clique_table_ == nullptr) { - CLIQUE_CUTS_DEBUG("generate_clique_cuts no clique table available, skipping"); + // The fractional conflict-graph subgraph is built once per cut pass in + // prepare_fractional_sub_conflict_graph() (called from generate_cuts) and shared with + // the zero-half cut separator. Skip if the build was unable to produce a + // useable sub-CG (clique table missing/empty, work/time budget hit, etc.). + if (!sub_cg_.ready) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts sub_cg_ not ready, skipping"); return true; } - CLIQUE_CUTS_DEBUG("generate_clique_cuts using clique table first=%lld addtl=%lld", - static_cast(clique_table_->first.size()), - static_cast(clique_table_->addtl_cliques.size())); - - if (clique_table_->empty()) { - CLIQUE_CUTS_DEBUG("generate_clique_cuts empty clique table, nothing to separate"); + if (sub_cg_.empty_subgraph()) { + CLIQUE_CUTS_DEBUG("generate_clique_cuts no fractional binary vertices"); return true; } - - cuopt_assert(clique_table_->n_variables == num_vars, "Clique table variable count mismatch"); + cuopt_assert(sub_cg_.num_vars == num_vars, "generate_clique_cuts sub_cg_ num_vars mismatch"); cuopt_assert(static_cast(num_vars) <= xstar.size(), "Clique cut xstar size mismatch"); + cuopt_assert(user_problem_.var_types.size() == static_cast(num_vars), + "User problem var_types size mismatch"); const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); const f_t bound_tol = settings.primal_tol; @@ -2754,98 +3526,20 @@ bool cut_generation_t::generate_clique_cuts( f_t work_estimate = 0.0; const f_t max_work_estimate = 1e8; - cuopt_assert(user_problem_.var_types.size() == static_cast(num_vars), - "User problem var_types size mismatch"); - - std::vector vertices; - std::vector weights; - vertices.reserve(num_vars * 2); - weights.reserve(num_vars * 2); + const std::vector& vertices = sub_cg_.vertices; + const std::vector& weights = sub_cg_.weights; + const std::vector>& adj_local = sub_cg_.adj_local; - // create the sub graph induced by fractional binary variables - for (i_t j = 0; j < num_vars; ++j) { - if (user_problem_.var_types[j] == variable_type_t::CONTINUOUS) { continue; } - const f_t lower_bound = user_problem_.lower[j]; - const f_t upper_bound = user_problem_.upper[j]; - if (lower_bound < -bound_tol || upper_bound > 1 + bound_tol) { continue; } - const f_t xj = xstar[j]; - if (std::abs(xj - std::round(xj)) <= settings.integer_tol) { continue; } - vertices.push_back(j); - weights.push_back(xj); - vertices.push_back(j + num_vars); - weights.push_back(1.0 - xj); - } - // Coarse loop estimate: variable scans + selected vertex/weight writes - work_estimate += 4.0 * static_cast(num_vars) + 2.0 * static_cast(vertices.size()); - if (work_estimate > max_work_estimate) { return true; } - - if (vertices.empty()) { - CLIQUE_CUTS_DEBUG("generate_clique_cuts no fractional binary vertices"); - return true; - } CLIQUE_CUTS_DEBUG("generate_clique_cuts fractional subgraph vertices=%lld (literals=%lld)", static_cast(vertices.size() / 2), static_cast(vertices.size())); - std::vector vertex_to_local(2 * num_vars, -1); - std::vector in_subgraph(2 * num_vars, 0); - for (size_t idx = 0; idx < vertices.size(); ++idx) { - if (toc(start_time) >= settings.time_limit) { return true; } - const i_t vertex_idx = vertices[idx]; - vertex_to_local[vertex_idx] = static_cast(idx); - in_subgraph[vertex_idx] = 1; - } - work_estimate += 3.0 * static_cast(vertices.size()); - if (work_estimate > max_work_estimate) { return true; } - - std::vector> adj_local(vertices.size()); - size_t total_adj_entries = 0; - size_t kept_adj_entries = 0; - for (size_t idx = 0; idx < vertices.size(); ++idx) { - if (toc(start_time) >= settings.time_limit) { return true; } - i_t vertex_idx = vertices[idx]; - // returns the complement as well - auto adj_set = clique_table_->get_adj_set_of_var(vertex_idx); - total_adj_entries += adj_set.size(); - auto& adj = adj_local[idx]; - adj.reserve(adj_set.size()); - for (const auto neighbor : adj_set) { - if (toc(start_time) >= settings.time_limit) { return true; } - cuopt_assert(neighbor >= 0 && neighbor < 2 * num_vars, "Neighbor out of range"); - if (!in_subgraph[neighbor]) { continue; } - i_t local_neighbor = vertex_to_local[neighbor]; - cuopt_assert(local_neighbor >= 0, "Local neighbor out of range"); - adj.push_back(local_neighbor); - } - kept_adj_entries += adj.size(); -#ifdef ASSERT_MODE - { - // {k, ~k} as neighbors is legal (vertex_idx is then implicitly fixed to - // 0 by the conflict structure); build_clique_cut handles the resulting - // cliques as fixing cuts or infeasibility signals, so only duplicates - // are a real invariant here. - std::unordered_set adj_global; - adj_global.reserve(adj.size()); - for (const auto neighbor : adj) { - i_t v = vertices[neighbor]; - cuopt_assert(adj_global.insert(v).second, "Duplicate neighbor in adjacency list"); - } - } -#endif - } - work_estimate += static_cast(vertices.size()) + static_cast(total_adj_entries) + - 2.0 * static_cast(kept_adj_entries); - if (work_estimate > max_work_estimate) { return true; } - CLIQUE_CUTS_DEBUG("generate_clique_cuts adjacency raw_entries=%lld kept_entries=%lld", - static_cast(total_adj_entries), - static_cast(kept_adj_entries)); - const size_t words = bitset_words(vertices.size()); std::vector> adj_bitset(vertices.size(), std::vector(words, 0)); size_t local_adj_entries = 0; for (size_t v = 0; v < adj_local.size(); ++v) { local_adj_entries += adj_local[v].size(); - for (const auto neighbor : adj_local[v]) { + for (const i_t neighbor : adj_local[v]) { bitset_set(adj_bitset[v], static_cast(neighbor)); } } @@ -2894,14 +3588,14 @@ bool cut_generation_t::generate_clique_cuts( size_t rejected_cliques = 0; size_t extension_gain = 0; #endif - for (auto& clique_local : ctx.cliques) { + for (std::vector& clique_local : ctx.cliques) { if (toc(start_time) >= settings.time_limit) { return true; } #if DEBUG_CLIQUE_CUTS candidate_cliques++; #endif std::vector clique_vertices; clique_vertices.reserve(clique_local.size()); - for (auto local_idx : clique_local) { + for (const i_t local_idx : clique_local) { clique_vertices.push_back(vertices[local_idx]); } work_estimate += 3.0 * static_cast(clique_local.size()); @@ -2924,18 +3618,18 @@ bool cut_generation_t::generate_clique_cuts( #endif if (work_estimate > max_work_estimate) { return true; } if (toc(start_time) >= settings.time_limit) { return true; } - const auto build_status = build_clique_cut(clique_vertices, - num_vars, - var_types, - user_problem_.lower, - user_problem_.upper, - xstar, - bound_tol, - min_violation, - cut, - cut_rhs, - &work_estimate, - max_work_estimate); + const clique_cut_build_status_t build_status = build_clique_cut(clique_vertices, + num_vars, + var_types, + user_problem_.lower, + user_problem_.upper, + xstar, + bound_tol, + min_violation, + cut, + cut_rhs, + &work_estimate, + max_work_estimate); if (work_estimate > max_work_estimate) { return true; } if (build_status == clique_cut_build_status_t::INFEASIBLE) { settings.log.debug("Detected contradictory variable/complement clique\n"); @@ -2976,6 +3670,194 @@ bool cut_generation_t::generate_clique_cuts( return true; } +template +bool cut_generation_t::generate_zero_half_cuts( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + const std::vector& reduced_costs, + f_t start_time) +{ + if (settings.zero_half_cuts == 0) { return true; } + if (toc(start_time) >= settings.time_limit) { return true; } + + const i_t num_vars = user_problem_.num_cols; + ZERO_HALF_DEBUG( + "generate_zero_half_cuts ENTER num_vars=%lld elapsed=%g time_limit=%g xstar.size=%zu " + "reduced_costs.size=%zu var_types.size=%zu user_problem_.lower.size=%zu " + "user_problem_.upper.size=%zu user_problem_.var_types.size=%zu lp.num_cols=%lld " + "sub_cg_.ready=%d sub_cg_.vertices=%zu", + static_cast(num_vars), + static_cast(toc(start_time)), + static_cast(settings.time_limit), + xstar.size(), + reduced_costs.size(), + var_types.size(), + user_problem_.lower.size(), + user_problem_.upper.size(), + user_problem_.var_types.size(), + static_cast(lp.num_cols), + static_cast(sub_cg_.ready), + sub_cg_.vertices.size()); + + // The fractional conflict-graph subgraph is built once per cut pass in + // prepare_fractional_sub_conflict_graph() (called from generate_cuts) and shared with + // the clique-cut separator. Skip if the build was unable to produce a + // useable sub-CG (clique table missing/empty, work/time budget hit, etc.). + if (!sub_cg_.ready) { + ZERO_HALF_DEBUG("sub_cg_ not ready, skipping"); + return true; + } + if (sub_cg_.empty_subgraph()) { + ZERO_HALF_DEBUG("no fractional binary vertices"); + return true; + } + if (clique_table_ == nullptr) { + ZERO_HALF_DEBUG("no clique table available, skipping"); + return true; + } + cuopt_assert(sub_cg_.num_vars == num_vars, "generate_zero_half_cuts sub_cg_ num_vars mismatch"); + cuopt_assert(clique_table_->n_variables == num_vars, + "Zero-half clique table variable count mismatch"); + cuopt_assert(static_cast(num_vars) <= xstar.size(), "Zero-half xstar size mismatch"); + cuopt_assert(user_problem_.var_types.size() == static_cast(num_vars), + "Zero-half user problem var_types size mismatch"); + + const f_t min_violation = std::max(settings.primal_tol, static_cast(1e-6)); + const f_t bound_tol = settings.primal_tol; + // shortest path of length >= 0.5 - min_violation cannot yield a violated cut + const f_t cutoff = static_cast(0.5) - min_violation; + f_t work_estimate = 0.0; + const f_t max_work_estimate = 1e8; + + const std::vector& vertices = sub_cg_.vertices; + const std::vector& weights = sub_cg_.weights; + const std::vector>& adj_local = sub_cg_.adj_local; + const std::vector& vertex_to_local = sub_cg_.vertex_to_local; + const i_t num_local = sub_cg_.num_local(); + ZERO_HALF_DEBUG("starting separation loop num_local=%lld", static_cast(num_local)); + + sparse_vector_t cut(lp.num_cols, 0); + f_t cut_rhs = 0.0; + std::vector bipartite_path; + std::vector cycle_vertices; + std::vector wheel_centers; + + i_t cycles_found = 0; + i_t cuts_added = 0; + i_t added_per_var = 0; + std::vector already_used(num_local, 0); + dijkstra_scratch_t dijkstra_scratch; + + for (i_t s = 0; s < num_local; ++s) { + if (toc(start_time) >= settings.time_limit) { break; } + if (work_estimate > max_work_estimate) { break; } + if (already_used[s]) { continue; } + ZERO_HALF_DEBUG("separation loop s=%lld / %lld", + static_cast(s), + static_cast(num_local)); + + f_t total_weight = 0; + if (!dijkstra_odd_cycle(s, + adj_local, + weights, + cutoff, + bipartite_path, + total_weight, + &work_estimate, + max_work_estimate, + dijkstra_scratch)) { + continue; + } + if (!path_to_odd_cycle(bipartite_path, + vertices, + num_local, + num_vars, + cycle_vertices, + &work_estimate, + max_work_estimate)) { + continue; + } + cycles_found++; + cuopt_assert(cycle_vertices.size() >= 5 && (cycle_vertices.size() % 2) == 1, + "Zero-half separated cycle must be odd with length >= 5"); + // dijkstra_odd_cycle only returns true when the path stays below the + // half-integer cutoff, the precondition for the cycle to yield a violation. + cuopt_assert(cutoff <= static_cast(0) || total_weight < cutoff, + "Zero-half cycle weight must be below cutoff"); + ZERO_HALF_DEBUG("cycle found s=%lld cycle_vertices.size=%zu", + static_cast(s), + cycle_vertices.size()); + + extend_to_odd_wheel(cycle_vertices, + wheel_centers, + *clique_table_, + reduced_costs, + num_vars, + start_time, + settings.time_limit, + &work_estimate, + max_work_estimate); + + ZERO_HALF_DEBUG("calling build_zero_half_cut cycle=%zu wheel=%zu", + cycle_vertices.size(), + wheel_centers.size()); + const clique_cut_build_status_t build_status = + build_zero_half_cut(cycle_vertices, + wheel_centers, + num_vars, + var_types, + user_problem_.lower, + user_problem_.upper, + xstar, + bound_tol, + min_violation, + cut, + cut_rhs, + &work_estimate, + max_work_estimate); + ZERO_HALF_DEBUG("build_zero_half_cut returned status=%d", static_cast(build_status)); + if (work_estimate > max_work_estimate) { break; } + if (build_status == clique_cut_build_status_t::INFEASIBLE) { + ZERO_HALF_DEBUG("infeasible cycle detected, returning false"); + return false; + } + if (build_status == clique_cut_build_status_t::CUT_ADDED) { + // Only violated cuts are worth pooling; build_zero_half_cut promised a + // violation > min_violation, so re-check it before we commit. + cuopt_assert(cut_rhs - cut.dot(xstar) > min_violation - bound_tol, + "Zero-half cut added to pool must be violated by xstar"); + inequality_t cut_inequality; + cut_inequality.vector = cut; + cut_inequality.rhs = cut_rhs; + ZERO_HALF_DEBUG( + "adding cut to pool nz=%zu rhs=%g", cut.i.size(), static_cast(cut_rhs)); + cut_pool_.add_cut(cut_type_t::ZERO_HALF, cut_inequality); + ZERO_HALF_DEBUG("cut added to pool"); + cuts_added++; + added_per_var++; + // mark all CG vertices that participated so we do not re-derive the same + // cycle from a different source vertex + for (const i_t v : cycle_vertices) { + if (v < 0 || v >= 2 * num_vars) { + ZERO_HALF_DEBUG("mark already_used: cycle v OUT_OF_RANGE v=%lld", + static_cast(v)); + continue; + } + const i_t lv = vertex_to_local[v]; + if (lv >= 0 && lv < num_local) { already_used[lv] = 1; } + } + } + } + + ZERO_HALF_DEBUG("generate_zero_half_cuts EXIT cycles=%lld cuts=%lld work=%g", + static_cast(cycles_found), + static_cast(cuts_added), + static_cast(work_estimate)); + return true; +} + template void cut_generation_t::generate_mir_cuts( const lp_problem_t& lp, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index e8dfba4ccd..ab9d3acca1 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -41,8 +41,9 @@ enum cut_type_t : int8_t { CHVATAL_GOMORY = 3, CLIQUE = 4, IMPLIED_BOUND = 5, - FLOW_COVER = 6, - MAX_CUT_TYPE = 7 + ZERO_HALF = 6, + FLOW_COVER = 7, + MAX_CUT_TYPE = 8 }; template @@ -183,6 +184,7 @@ struct cut_info_t { "Strong CG ", "Clique ", "Implied Bounds", + "Zero-Half ", "Flow Cover "}; std::array num_cuts = {0}; }; @@ -275,6 +277,16 @@ std::vector> find_maximal_cliques_for_test( int max_calls, double time_limit); +// Test-only helper to run the production odd-cycle separator used by zero-half cuts. +// adjacency_list must contain local vertex indices in [0, n_vertices). x_values gives +// the LP value for each vertex. Returns simple odd cycles whose induced edge weight +// sum is < 0.5 - min_violation. +std::vector> find_violated_odd_cycles_for_test( + const std::vector>& adjacency_list, + const std::vector& x_values, + double min_violation, + double time_limit); + template class cut_pool_t { public: @@ -583,6 +595,34 @@ class knapsack_generation_t { template class mixed_integer_rounding_cut_t; +template +class variable_bounds_t; + +template +struct fractional_conflict_subgraph_t { + i_t num_vars{0}; + std::vector vertices; + std::vector weights; + std::vector vertex_to_local; + std::vector in_subgraph; + std::vector> adj_local; + bool ready{false}; + + i_t num_local() const { return static_cast(vertices.size()); } + bool empty_subgraph() const { return vertices.empty(); } + + void clear() + { + num_vars = 0; + vertices.clear(); + weights.clear(); + vertex_to_local.clear(); + in_subgraph.clear(); + adj_local.clear(); + ready = false; + } +}; + template class cut_generation_t { public: @@ -668,6 +708,14 @@ class cut_generation_t { const std::vector& reduced_costs, f_t start_time); + // Generate zero-half (odd-cycle / odd-wheel) cuts from the conflict graph + bool generate_zero_half_cuts(const simplex::lp_problem_t& lp, + const simplex::simplex_solver_settings_t& settings, + const std::vector& var_types, + const std::vector& xstar, + const std::vector& reduced_costs, + f_t start_time); + // Generate implied bounds cuts from probing implications void generate_implied_bound_cuts(const simplex::lp_problem_t& lp, const simplex::simplex_solver_settings_t& settings, @@ -675,6 +723,11 @@ class cut_generation_t { const std::vector& xstar, f_t start_time); + void prepare_fractional_sub_conflict_graph( + const simplex::simplex_solver_settings_t& settings, + const std::vector& xstar, + f_t start_time); + cut_pool_t& cut_pool_; knapsack_generation_t knapsack_generation_; flow_cover_generation_t flow_cover_generation_; @@ -682,6 +735,7 @@ class cut_generation_t { const probing_implied_bound_t& probing_implied_bound_; std::shared_ptr> clique_table_; omp_atomic_t* signal_extend_{nullptr}; + fractional_conflict_subgraph_t sub_cg_; }; template diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index 45168ced8a..5ef4515196 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -88,6 +88,7 @@ struct simplex_solver_settings_t { flow_cover_cuts(-1), implied_bound_cuts(-1), clique_cuts(-1), + zero_half_cuts(-1), strong_chvatal_gomory_cuts(-1), symmetry(-1), reduced_cost_strengthening(-1), @@ -179,6 +180,7 @@ struct simplex_solver_settings_t { i_t flow_cover_cuts; // -1 automatic, 0 to disable, >0 to enable flow cover cuts i_t implied_bound_cuts; // -1 automatic, 0 to disable, >0 to enable implied bound cuts i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique cuts + i_t zero_half_cuts; // -1 automatic, 0 to disable, >0 to enable zero-half cuts i_t strong_chvatal_gomory_cuts; // -1 automatic, 0 to disable, >0 to enable strong Chvatal Gomory // cuts i_t symmetry; // -1 automatic, 0 to disable, >0 to enable different symmetry methods diff --git a/cpp/src/grpc/codegen/field_registry.yaml b/cpp/src/grpc/codegen/field_registry.yaml index e2c26a797a..1981af5091 100644 --- a/cpp/src/grpc/codegen/field_registry.yaml +++ b/cpp/src/grpc/codegen/field_registry.yaml @@ -682,6 +682,10 @@ mip_settings: field_num: 22 type: int32 optional: true + - zero_half_cuts: + field_num: 52 + type: int32 + optional: true - implied_bound_cuts: field_num: 31 type: int32 diff --git a/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto b/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto index 645fee856b..6c3931664e 100644 --- a/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto +++ b/cpp/src/grpc/codegen/generated/cuopt_remote_data.proto @@ -248,6 +248,7 @@ message MIPSolverSettings { optional int32 cycle_detection_length = 49; optional double relaxed_lp_time_limit = 50; optional double related_vars_time_limit = 51; + optional int32 zero_half_cuts = 52; } message PDLPWarmStartData { diff --git a/cpp/src/grpc/codegen/generated/generated_mip_settings_to_proto.inc b/cpp/src/grpc/codegen/generated/generated_mip_settings_to_proto.inc index aeabff36fc..12d95aa6df 100644 --- a/cpp/src/grpc/codegen/generated/generated_mip_settings_to_proto.inc +++ b/cpp/src/grpc/codegen/generated/generated_mip_settings_to_proto.inc @@ -31,6 +31,7 @@ pb_settings->set_mixed_integer_gomory_cuts(settings.mixed_integer_gomory_cuts); pb_settings->set_knapsack_cuts(settings.knapsack_cuts); pb_settings->set_clique_cuts(settings.clique_cuts); + pb_settings->set_zero_half_cuts(settings.zero_half_cuts); pb_settings->set_implied_bound_cuts(settings.implied_bound_cuts); pb_settings->set_strong_chvatal_gomory_cuts(settings.strong_chvatal_gomory_cuts); pb_settings->set_reduced_cost_strengthening(settings.reduced_cost_strengthening); diff --git a/cpp/src/grpc/codegen/generated/generated_proto_to_mip_settings.inc b/cpp/src/grpc/codegen/generated/generated_proto_to_mip_settings.inc index f1d4f40cae..e14ffa1ae0 100644 --- a/cpp/src/grpc/codegen/generated/generated_proto_to_mip_settings.inc +++ b/cpp/src/grpc/codegen/generated/generated_proto_to_mip_settings.inc @@ -73,6 +73,9 @@ if (pb_settings.has_clique_cuts()) { settings.clique_cuts = pb_settings.clique_cuts(); } + if (pb_settings.has_zero_half_cuts()) { + settings.zero_half_cuts = pb_settings.zero_half_cuts(); + } if (pb_settings.has_implied_bound_cuts()) { settings.implied_bound_cuts = pb_settings.implied_bound_cuts(); } diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 611881ebd0..8bfd7cd8ba 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -138,6 +138,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_KNAPSACK_CUTS, &mip_settings.knapsack_cuts, -1, 1, -1}, {CUOPT_MIP_FLOW_COVER_CUTS, &mip_settings.flow_cover_cuts, -1, 1, -1}, {CUOPT_MIP_CLIQUE_CUTS, &mip_settings.clique_cuts, -1, 1, -1}, + {CUOPT_MIP_ZERO_HALF_CUTS, &mip_settings.zero_half_cuts, -1, 1, -1}, {CUOPT_MIP_IMPLIED_BOUND_CUTS, &mip_settings.implied_bound_cuts, -1, 1, -1}, {CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS, &mip_settings.strong_chvatal_gomory_cuts, -1, 1, -1}, {CUOPT_MIP_REDUCED_COST_STRENGTHENING, &mip_settings.reduced_cost_strengthening, -1, std::numeric_limits::max(), -1}, diff --git a/cpp/src/mip_heuristics/diversity/lns/rins.cu b/cpp/src/mip_heuristics/diversity/lns/rins.cu index a670322cb0..55bf3d01aa 100644 --- a/cpp/src/mip_heuristics/diversity/lns/rins.cu +++ b/cpp/src/mip_heuristics/diversity/lns/rins.cu @@ -254,6 +254,7 @@ void rins_t::run_rins() branch_and_bound_settings.reliability_branching = 0; branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.clique_cuts = 0; + branch_and_bound_settings.zero_half_cuts = 0; branch_and_bound_settings.sub_mip = 1; branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; branch_and_bound_settings.log.log = false; diff --git a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh index 0f3625d9aa..004fe041da 100644 --- a/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh +++ b/cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh @@ -111,6 +111,7 @@ class sub_mip_recombiner_t : public recombiner_t { branch_and_bound_settings.reliability_branching = 0; branch_and_bound_settings.max_cut_passes = 0; branch_and_bound_settings.clique_cuts = 0; + branch_and_bound_settings.zero_half_cuts = 0; branch_and_bound_settings.sub_mip = 1; branch_and_bound_settings.strong_branching_simplex_iteration_limit = 200; branch_and_bound_settings.solution_callback = [this](std::vector& solution, diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu index bcf2d1a9dd..923cd5edb8 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cu @@ -667,7 +667,7 @@ void print_clique_table(const clique_table_t& clique_table) template void find_initial_cliques(user_problem_t& problem, typename mip_solver_settings_t::tolerances_t tolerances, - std::shared_ptr>* clique_table_out, + std::shared_ptr>& clique_table_out, cuopt::timer_t& timer, omp_atomic_t* signal_extend) { @@ -698,53 +698,48 @@ void find_initial_cliques(user_problem_t& problem, t_sort = stage_timer.elapsed_time(); #endif clique_config_t clique_config; - std::shared_ptr> clique_table_shared; - clique_table_t clique_table_local(2 * problem.num_cols, - clique_config.min_clique_size, - clique_config.max_clique_size_for_extension); - clique_table_t* clique_table_ptr = &clique_table_local; - if (clique_table_out != nullptr) { - clique_table_shared = - std::make_shared>(2 * problem.num_cols, - clique_config.min_clique_size, - clique_config.max_clique_size_for_extension); - clique_table_ptr = clique_table_shared.get(); - } - clique_table_ptr->tolerances = tolerances; + auto clique_table = + std::make_shared>(2 * problem.num_cols, + clique_config.min_clique_size, + clique_config.max_clique_size_for_extension); + clique_table->tolerances = tolerances; double time_limit_for_additional_cliques = timer.remaining_time() / 2; cuopt::timer_t additional_cliques_timer(time_limit_for_additional_cliques); double find_work_estimate = 0.0; // Always build base cliques in full; signal_extend only gates the extension phase. for (const auto& knapsack_constraint : knapsack_constraints) { if (timer.check_time_limit()) { break; } - find_cliques_from_constraint(knapsack_constraint, *clique_table_ptr, additional_cliques_timer); + find_cliques_from_constraint(knapsack_constraint, *clique_table, additional_cliques_timer); find_work_estimate += knapsack_constraint.entries.size(); } #ifdef DEBUG_CLIQUE_TABLE t_find = stage_timer.elapsed_time(); #endif CUOPT_LOG_DEBUG("Number of cliques: %d, additional cliques: %d, find_work=%.0f", - clique_table_ptr->first.size(), - clique_table_ptr->addtl_cliques.size(), + clique_table->first.size(), + clique_table->addtl_cliques.size(), find_work_estimate); - remove_small_cliques(*clique_table_ptr, timer); + remove_small_cliques(*clique_table, timer); #ifdef DEBUG_CLIQUE_TABLE t_small = stage_timer.elapsed_time(); #endif - fill_var_clique_maps(*clique_table_ptr); + fill_var_clique_maps(*clique_table); #ifdef DEBUG_CLIQUE_TABLE t_maps = stage_timer.elapsed_time(); #endif - if (clique_table_out != nullptr) { *clique_table_out = std::move(clique_table_shared); } + // Publish the base table so cut generation can start using it; the extension + // phase below keeps mutating *clique_table, so the consumer must signal this + // task to stop and join it (taskwait) before reading the table. + clique_table_out = clique_table; double extend_work = 0.0; i_t n_extended_cliques = extend_cliques(knapsack_constraints, - *clique_table_ptr, + *clique_table, timer, &extend_work, clique_config.min_extend_work, clique_config.max_extend_work, signal_extend); - if (n_extended_cliques > 0) { fill_var_clique_maps(*clique_table_ptr); } + if (n_extended_cliques > 0) { fill_var_clique_maps(*clique_table); } #ifdef DEBUG_CLIQUE_TABLE t_extend = stage_timer.elapsed_time(); CUOPT_LOG_DEBUG( @@ -767,7 +762,7 @@ void find_initial_cliques(user_problem_t& problem, template void find_initial_cliques( \ user_problem_t & problem, \ typename mip_solver_settings_t::tolerances_t tolerances, \ - std::shared_ptr> * clique_table_out, \ + std::shared_ptr> & clique_table_out, \ cuopt::timer_t & timer, \ omp_atomic_t * signal_extend); \ template void build_clique_table( \ diff --git a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh index 4abb7bfb52..d9dd0a8920 100644 --- a/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh +++ b/cpp/src/mip_heuristics/presolve/conflict_graph/clique_table.cuh @@ -204,10 +204,16 @@ struct clique_table_t { typename mip_solver_settings_t::tolerances_t tolerances; }; +// Builds the conflict-graph clique table for `problem`. The base cliques are +// published to `clique_table_out` before the (optional, signal-gated) extension +// phase begins, so cut generation can pick up the table while extension keeps +// running concurrently. Consumers MUST set `*signal_extend` and join the +// producing task before reading the table (see prepare_fractional_sub_conflict_graph), since +// the extension phase keeps mutating the same object after it is published. template void find_initial_cliques(simplex::user_problem_t& problem, typename mip_solver_settings_t::tolerances_t tolerances, - std::shared_ptr>* clique_table_out, + std::shared_ptr>& clique_table_out, cuopt::timer_t& timer, omp_atomic_t* signal_extend = nullptr); diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index bac04cb053..d932b1f845 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -24,6 +24,7 @@ #define DETECT_SYMMETRY_AFTER_PRESOLVE #include +#include #include #include @@ -374,6 +375,7 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.flow_cover_cuts = context.settings.flow_cover_cuts; branch_and_bound_settings.implied_bound_cuts = context.settings.implied_bound_cuts; branch_and_bound_settings.clique_cuts = context.settings.clique_cuts; + branch_and_bound_settings.zero_half_cuts = context.settings.zero_half_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; branch_and_bound_settings.cut_change_threshold = context.settings.cut_change_threshold; diff --git a/cpp/src/utilities/omp_helpers.hpp b/cpp/src/utilities/omp_helpers.hpp index cb371fbde4..8080e53fc5 100644 --- a/cpp/src/utilities/omp_helpers.hpp +++ b/cpp/src/utilities/omp_helpers.hpp @@ -39,7 +39,6 @@ class omp_mutex_t { omp_mutex_t(omp_mutex_t&& other) { *this = std::move(other); } - omp_mutex_t(const omp_mutex_t&) = delete; omp_mutex_t& operator=(const omp_mutex_t&) = delete; omp_mutex_t& operator=(omp_mutex_t&& other) diff --git a/cpp/tests/mip/cuts_test.cu b/cpp/tests/mip/cuts_test.cu index 642f69a04e..b4887b6a8f 100644 --- a/cpp/tests/mip/cuts_test.cu +++ b/cpp/tests/mip/cuts_test.cu @@ -63,6 +63,27 @@ End )LP"); } +io::mps_data_model_t create_pairwise_pentagon_set_packing_problem() +{ + return cuopt::test::parse_inline_lp(R"LP( +Minimize + obj: -x0 - x1 - x2 - x3 - x4 +Subject To + c1: x0 + x1 <= 1 + c2: x1 + x2 <= 1 + c3: x2 + x3 <= 1 + c4: x3 + x4 <= 1 + c5: x4 + x0 <= 1 +Binaries + x0 + x1 + x2 + x3 + x4 +End +)LP"); +} + // Same triangle conflicts plus an isolated binary x3 with no conflict rows. io::mps_data_model_t create_pairwise_triangle_with_isolated_variable_problem() { @@ -337,6 +358,18 @@ std::string format_phase2_panic_dump(const io::mps_data_model_t& pr void disable_non_clique_cuts(mip_solver_settings_t& settings) { settings.clique_cuts = 1; + settings.zero_half_cuts = 0; + settings.max_cut_passes = 10; + settings.mixed_integer_gomory_cuts = 0; + settings.knapsack_cuts = 0; + settings.mir_cuts = 0; + settings.strong_chvatal_gomory_cuts = 0; +} + +void disable_non_zero_half_cuts(mip_solver_settings_t& settings) +{ + settings.clique_cuts = 1; + settings.zero_half_cuts = 1; settings.max_cut_passes = 10; settings.mixed_integer_gomory_cuts = 0; settings.knapsack_cuts = 0; @@ -348,6 +381,7 @@ void disable_all_cuts(mip_solver_settings_t& settings) { settings.max_cut_passes = 0; settings.clique_cuts = 0; + settings.zero_half_cuts = 0; settings.mixed_integer_gomory_cuts = 0; settings.knapsack_cuts = 0; settings.mir_cuts = 0; @@ -1311,6 +1345,218 @@ TEST(cuts, clique_neos8_phase4_lp_infeasibility_binary_search) EXPECT_EQ(first_infeasible.value(), injected_index); } +// ---- Zero-half cut tests -------------------------------------------------- + +namespace { + +std::vector> canonicalize_cycles(std::vector> cycles) +{ + for (auto& cycle : cycles) { + if (cycle.empty()) { continue; } + auto min_it = std::min_element(cycle.begin(), cycle.end()); + std::rotate(cycle.begin(), min_it, cycle.end()); + if (cycle.size() >= 3 && cycle[1] > cycle.back()) { + std::reverse(cycle.begin() + 1, cycle.end()); + } + } + std::sort(cycles.begin(), cycles.end()); + cycles.erase(std::unique(cycles.begin(), cycles.end()), cycles.end()); + return cycles; +} + +} // namespace + +TEST(cuts, zero_half_unit_separator_simple_pentagon) +{ + // 5-cycle: 0-1-2-3-4-0. All vertices fractional at 0.5. + std::vector> adj = { + {1, 4}, + {0, 2}, + {1, 3}, + {2, 4}, + {3, 0}, + }; + std::vector x_values(5, 0.5); + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + ASSERT_FALSE(cycles.empty()); + cycles = canonicalize_cycles(std::move(cycles)); + std::vector expected{0, 1, 2, 3, 4}; + bool found = false; + for (const auto& cycle : cycles) { + if (cycle.size() == 5) { + auto sorted = cycle; + std::sort(sorted.begin(), sorted.end()); + if (sorted == expected) { + found = true; + break; + } + } + } + EXPECT_TRUE(found); +} + +TEST(cuts, zero_half_unit_separator_no_cycle_for_4_cycle) +{ + // Even cycle: 0-1-2-3-0 + std::vector> adj = { + {1, 3}, + {0, 2}, + {1, 3}, + {2, 0}, + }; + std::vector x_values(4, 0.5); + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + EXPECT_TRUE(cycles.empty()); +} + +TEST(cuts, zero_half_unit_separator_skips_triangle) +{ + // Triangle 0-1-2-0 ; size-3 cycles must be left to the clique separator. + std::vector> adj = { + {1, 2}, + {0, 2}, + {0, 1}, + }; + std::vector x_values(3, 0.5); + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + for (const auto& cycle : cycles) { + EXPECT_GE(cycle.size(), 5u); + } +} + +TEST(cuts, zero_half_unit_separator_no_cycle_when_integer_solution) +{ + // 5-cycle but x_values are integer feasible: (1, 0, 1, 0, 0) -- no violation. + std::vector> adj = { + {1, 4}, + {0, 2}, + {1, 3}, + {2, 4}, + {3, 0}, + }; + std::vector x_values = {1.0, 0.0, 1.0, 0.0, 0.0}; + // x_v interpreted as conflict-graph vertex weight (here just x_j directly). + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + EXPECT_TRUE(cycles.empty()); +} + +TEST(cuts, zero_half_unit_separator_disjoint_pentagons) +{ + // Two disjoint 5-cycles share no vertices: {0..4} and {5..9}. + std::vector> adj = { + {1, 4}, + {0, 2}, + {1, 3}, + {2, 4}, + {3, 0}, + {6, 9}, + {5, 7}, + {6, 8}, + {7, 9}, + {8, 5}, + }; + std::vector x_values(10, 0.5); + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + ASSERT_GE(cycles.size(), 2u); + cycles = canonicalize_cycles(std::move(cycles)); + bool found_left = false; + bool found_right = false; + for (const auto& cycle : cycles) { + if (cycle.size() != 5) { continue; } + auto sorted = cycle; + std::sort(sorted.begin(), sorted.end()); + if (sorted == std::vector{0, 1, 2, 3, 4}) { found_left = true; } + if (sorted == std::vector{5, 6, 7, 8, 9}) { found_right = true; } + } + EXPECT_TRUE(found_left); + EXPECT_TRUE(found_right); +} + +TEST(cuts, zero_half_unit_separator_overlapping_pentagons) +{ + std::vector> adj = { + {1, 4, 5, 8}, + {0, 2}, + {1, 3}, + {2, 4}, + {3, 0}, + {0, 6}, + {5, 7}, + {6, 8}, + {7, 0}, + }; + std::vector x_values(9, 0.5); + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + cycles = canonicalize_cycles(std::move(cycles)); + + EXPECT_NE(std::find(cycles.begin(), cycles.end(), std::vector{0, 1, 2, 3, 4}), cycles.end()); + EXPECT_NE(std::find(cycles.begin(), cycles.end(), std::vector{0, 5, 6, 7, 8}), cycles.end()); +} + +TEST(cuts, zero_half_end_to_end_pentagon_tightens_lp_relaxation) +{ + const raft::handle_t handle{}; + auto mip_problem = create_pairwise_pentagon_set_packing_problem(); + + // First solve the LP relaxation (no cuts) to confirm the baseline value 2.5. + auto lp_relaxation = mip_problem; + std::vector all_continuous(lp_relaxation.get_n_variables(), 'C'); + lp_relaxation.set_variable_types(all_continuous); + + pdlp_solver_settings_t lp_settings{}; + lp_settings.time_limit = 10.0; + lp_settings.presolver = presolver_t::None; + lp_settings.set_optimality_tolerance(1e-8); + auto lp_solution = solve_lp(&handle, lp_relaxation, lp_settings); + ASSERT_EQ(lp_solution.get_termination_status(), pdlp_termination_status_t::Optimal); + const double lp_obj_no_cuts = lp_solution.get_objective_value(); + EXPECT_NEAR(lp_obj_no_cuts, -2.5, kCliqueTestTol); + + // Optimal IP value is 2 (independent set of size 2), so the LP gap is 0.5. + mip_solver_settings_t settings; + settings.time_limit = 10.0; + settings.presolver = presolver_t::None; + disable_non_zero_half_cuts(settings); + + auto mip_solution = solve_mip(&handle, mip_problem, settings); + ASSERT_EQ(mip_solution.get_termination_status(), mip_termination_status_t::Optimal); + EXPECT_NEAR(mip_solution.get_objective_value(), -2.0, kCliqueTestTol); +} + +TEST(cuts, zero_half_unit_separator_seven_cycle_violated_below_half) +{ + // 7-cycle: 0-1-2-3-4-5-6-0, all weights 0.4. Each edge weight = (1-0.4-0.4)/2 = 0.1 + // total path weight from j1 to j2 of length 7 = 0.7 — not below 0.5, so no cut. + // Make weights slightly higher: 0.45 → edge weight = 0.05, total = 7*0.05 = 0.35 < 0.5. + std::vector> adj = { + {1, 6}, + {0, 2}, + {1, 3}, + {2, 4}, + {3, 5}, + {4, 6}, + {5, 0}, + }; + std::vector x_values(7, 0.45); + auto cycles = mip::find_violated_odd_cycles_for_test( + adj, x_values, 1e-6, std::numeric_limits::infinity()); + ASSERT_FALSE(cycles.empty()); + bool found_seven = false; + for (const auto& cycle : cycles) { + if (cycle.size() == 7) { + found_seven = true; + break; + } + } + EXPECT_TRUE(found_seven); +} + // Minimal 0-1 single-node-flow relaxation for the flow-cover separator. // // y0 + y1 - y2 <= 4 diff --git a/docs/cuopt/source/cuopt-c/mip/mip-c-api.rst b/docs/cuopt/source/cuopt-c/mip/mip-c-api.rst index 9212ce1b3e..85592bb766 100644 --- a/docs/cuopt/source/cuopt-c/mip/mip-c-api.rst +++ b/docs/cuopt/source/cuopt-c/mip/mip-c-api.rst @@ -53,6 +53,7 @@ These constants configure MIP-specific solver behavior. Use them with :c:func:`c .. doxygendefine:: CUOPT_MIP_FLOW_COVER_CUTS .. doxygendefine:: CUOPT_MIP_IMPLIED_BOUND_CUTS .. doxygendefine:: CUOPT_MIP_CLIQUE_CUTS +.. doxygendefine:: CUOPT_MIP_ZERO_HALF_CUTS .. doxygendefine:: CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS .. doxygendefine:: CUOPT_MIP_REDUCED_COST_STRENGTHENING .. doxygendefine:: CUOPT_MIP_OBJECTIVE_STEP diff --git a/docs/cuopt/source/mip-settings.rst b/docs/cuopt/source/mip-settings.rst index f48d632097..025370ea3c 100644 --- a/docs/cuopt/source/mip-settings.rst +++ b/docs/cuopt/source/mip-settings.rst @@ -332,6 +332,16 @@ Set this value to ``1`` to enable clique cuts. .. note:: The default value is ``-1`` (automatic). +Zero-Half Cuts +^^^^^^^^^^^^^^ + +``CUOPT_MIP_ZERO_HALF_CUTS`` controls whether to use zero-half cuts. +Set this value to ``-1`` (automatic) to let the solver decide. +Set this value to ``0`` to disable zero-half cuts. +Set this value to ``1`` to enable zero-half cuts. + +.. note:: The default value is ``-1`` (automatic). + Objective Step ^^^^^^^^^^^^^^^ diff --git a/skills/cuopt-developer/references/contributing.md b/skills/cuopt-developer/references/contributing.md index 34fb75aab1..39f7e9794d 100644 --- a/skills/cuopt-developer/references/contributing.md +++ b/skills/cuopt-developer/references/contributing.md @@ -74,16 +74,36 @@ A few non-YAGNI points worth keeping in mind: When in doubt, mirror how the surrounding cuOpt code handles the same concern. +## Resolving Merge Conflicts + +Don't resolve a conflict by mechanically picking the side that looks like a superset. A small, local conflict (a few changed lines in one function) often sits on top of a larger architectural divergence — one branch refactored a mechanism the other left alone — and the conflict markers only show the tip of it. Picking "the bigger hunk" then strands the rest of that mechanism. + +Before choosing a side, reconstruct what each branch actually did: + +- Diff the conflicting symbols across **both branches and the merge base**, not just the two conflict hunks: `git show :` and `git merge-base A B`. Watch for changes to a member's *type*, an ownership/lifetime model, or a synchronization/threading model (e.g. `std::future` → OpenMP task, `std::atomic` → `omp_atomic_t`). Those changes ripple beyond the conflict region. +- Check how the **already-merged, non-conflicted files** use the symbol. If a caller (constructor call, factory, task spawn) was auto-merged to one branch's signature, the conflicted file must conform to that branch — keeping the other branch's member or wait logic leaves it dead. +- When one branch *removed* a mechanism and the other *built on top of it*, the correct resolution is usually to adopt the removal (the newer baseline) and re-port the feature onto the new mechanism — not to keep both, which yields a member that is never set and a guard that never fires. + +A wrong merge resolution frequently **compiles cleanly and fails silently**: a dead pointer stays `nullptr`, the guard that depended on it never triggers, and a whole feature quietly disables itself with no error. Compilation is not evidence of a correct merge — trace the runtime wiring (who sets this field? who waits on it? is that path still reachable?) before declaring the conflict resolved. + ## Common Tasks ### Adding a Solver Parameter -1. Add to settings struct in `cpp/include/cuopt/` and wire into `set_parameter_from_string()` in `cpp/src/` -2. Expose in Python — if using the string-based interface, the parameter is auto-discovered (no `.pyx` change needed). Add a convenience method in `SolverSettings` if warranted. See [python_bindings.md](python_bindings.md) for the full checklist. -3. Add to server schema (`docs/cuopt/source/cuopt_spec.yaml`) if applicable -4. Add tests at C++ and Python levels -5. Rebuild: `./build.sh libcuopt && ./build.sh cuopt` -6. Update documentation +Internal settings struct fields (e.g. a new MIP cut toggle) are often **already +read by the solver** but not yet exposed to the CLI/string interface, Python, or +the server. To fully expose one, wire it through every layer below — missing any +one silently drops the parameter from that interface. Use a sibling parameter +(e.g. `clique_cuts`) as the template and grep for it across the repo to find +every spot. + +1. **Settings struct** — add the field in `cpp/include/cuopt/.../solver_settings.hpp` (and the internal `simplex_solver_settings.hpp` if the solver reads it there). Default to `-1` (automatic) for cut toggles. +2. **C constant** — add a `#define CUOPT_MIP_ "mip_"` in `cpp/include/cuopt/linear_programming/constants.h`, next to the related parameters. +3. **String-parameter registry** — add a tuple to the `int_parameters` (or `float_parameters`) table in `cpp/src/math_optimization/solver_settings.cu`: `{CUOPT_MIP_, &mip_settings., , , }`. This single registry is what drives CLI parsing, `set_parameter`/`get_parameter`, and Python auto-discovery — no `.pyx` change needed. +4. **gRPC/server path** — add the field under the right message in `cpp/src/grpc/codegen/field_registry.yaml` with the next free `field_num` for that message, then regenerate with `python cpp/src/grpc/codegen/generate_conversions.py` (needs `pyyaml`). The generated `cuopt_remote_data.proto` and `generated_*_to_*.inc` files are **auto-generated — never hand-edit** unless you cannot run the generator, in which case mirror the sibling parameter exactly (proto fields are ordered by `field_num`; the `.inc` files follow registry declaration order). +5. **Docs** — add a `.. doxygendefine:: CUOPT_MIP_` to `docs/cuopt/source/cuopt-c/mip/mip-c-api.rst` and a settings section to `docs/cuopt/source/mip-settings.rst`. Add to the server schema (`docs/cuopt/source/cuopt_spec.yaml`) if applicable. +6. **Tests** — add C++ and Python coverage. +7. **Rebuild**: `./build.sh libcuopt && ./build.sh cuopt` ### Adding a Dependency diff --git a/skills/cuopt-developer/references/conventions.md b/skills/cuopt-developer/references/conventions.md index 3686c900d7..72275f780f 100644 --- a/skills/cuopt-developer/references/conventions.md +++ b/skills/cuopt-developer/references/conventions.md @@ -47,6 +47,29 @@ CUOPT_EXPECTS(condition, "Error message"); CUOPT_FAIL("Unreachable code reached"); ``` +### Assert-only variables + +A variable used only inside `cuopt_assert` (or any assertion that compiles out in +release builds) triggers an unused-variable warning when asserts are disabled. +Mark it `[[maybe_unused]]` at the declaration — do **not** suppress the warning +with `static_cast(var);` (or `(void)var;`) statements after the asserts. + +```cpp +// ❌ WRONG — trailing void-casts to silence the warning +const f_t lower_bound = lower_bounds[var_idx]; +const f_t upper_bound = upper_bounds[var_idx]; +cuopt_assert(lower_bound >= -bound_tol, "..."); +cuopt_assert(upper_bound <= 1 + bound_tol, "..."); +static_cast(lower_bound); +static_cast(upper_bound); + +// ✅ CORRECT — annotate at the declaration +[[maybe_unused]] const f_t lower_bound = lower_bounds[var_idx]; +[[maybe_unused]] const f_t upper_bound = upper_bounds[var_idx]; +cuopt_assert(lower_bound >= -bound_tol, "..."); +cuopt_assert(upper_bound <= 1 + bound_tol, "..."); +``` + ### CUDA Error Checking ```cpp