From 4fe3f943a0ede9a73301ea31e03612862a7196cc Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 14 Apr 2026 18:20:07 -0700 Subject: [PATCH 01/19] Initial symmetry detection using dejavu --- cpp/CMakeLists.txt | 15 ++ cpp/src/branch_and_bound/symmetry.hpp | 317 ++++++++++++++++++++++++++ cpp/src/mip_heuristics/solve.cu | 15 ++ 3 files changed, 347 insertions(+) create mode 100644 cpp/src/branch_and_bound/symmetry.hpp diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index abc6a19ab..5466ee77a 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -279,6 +279,20 @@ FetchContent_MakeAvailable(pslp) set(BUILD_SHARED_LIBS ${BUILD_SHARED_LIBS_SAVED}) +# dejavu - header-only graph automorphism library for MIP symmetry detection +# https://github.com/markusa4/dejavu (header-only, skip its CMakeLists.txt) +FetchContent_Declare( + dejavu + GIT_REPOSITORY "https://github.com/markusa4/dejavu.git" + GIT_TAG "v2.1" + GIT_PROGRESS TRUE + EXCLUDE_FROM_ALL + SYSTEM + SOURCE_SUBDIR _nonexistent +) +FetchContent_MakeAvailable(dejavu) +message(STATUS "dejavu (graph automorphism): ${dejavu_SOURCE_DIR}") + include(${rapids-cmake-dir}/cpm/rapids_logger.cmake) # generate logging macros rapids_cpm_rapids_logger(BUILD_EXPORT_SET cuopt-exports INSTALL_EXPORT_SET cuopt-exports) @@ -469,6 +483,7 @@ target_include_directories(cuopt PRIVATE target_include_directories(cuopt SYSTEM PRIVATE "${pslp_SOURCE_DIR}/include" + "${dejavu_SOURCE_DIR}" ) target_include_directories(cuopt diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp new file mode 100644 index 000000000..5b1befd8f --- /dev/null +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -0,0 +1,317 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#pragma once + +#include +#include +#include +#include +#include + +#include "dejavu.h" + +#include + +namespace cuopt::linear_programming::dual_simplex { + +template +void detect_symmetry(const user_problem_t& user_problem, + const simplex_solver_settings_t& settings) +{ + f_t start_time = tic(); + lp_problem_t problem(user_problem.handle_ptr, 1, 1, 1); + std::vector new_slacks; + dualize_info_t dualize_info; + convert_user_problem(user_problem, settings, problem, new_slacks, dualize_info); + std::vector var_types = user_problem.var_types; + if (problem.num_cols > user_problem.num_cols) { + var_types.resize(problem.num_cols); + for (i_t k = user_problem.num_cols; k < problem.num_cols; k++) { + var_types[k] = variable_type_t::CONTINUOUS; + } + } + + // We now have the problem in the form: + // minimize c^T x + // subject to A * x = b, + // l <= x <= u + // x_j in Z for all j such that var_types[j] == variable_type_t::INTEGER + + // Construct a graph G(V, W, R, E) + // where V is the set of nodes corresponding to the variables + // R is the set of nodes corresponding to the the constraints + // W is the set of nodes corresponding to the nonzero coefficients in the A matrix + // E is the set of edges + // + // Associated with each node in this graph is a color: c_v, c_w, c_r. + + const i_t V_size = problem.num_cols; + const i_t R_size = problem.num_rows; + + const f_t tol = 1e-10; + + // Compute the colors for the variables + std::vector obj_perm(problem.num_cols); + std::iota(obj_perm.begin(), obj_perm.end(), 0); + std::sort(obj_perm.begin(), obj_perm.end(), [&](i_t a, i_t b) { + if (problem.objective[a] != problem.objective[b]) + return problem.objective[a] < problem.objective[b]; + if (problem.lower[a] != problem.lower[b]) return problem.lower[a] < problem.lower[b]; + if (problem.upper[a] != problem.upper[b]) return problem.upper[a] < problem.upper[b]; + return var_types[a] < var_types[b]; + }); + std::vector var_colors(problem.num_cols, -1); + i_t var_color = 0; + f_t last_obj = problem.objective[obj_perm[0]]; + f_t last_lower = problem.lower[obj_perm[0]]; + f_t last_upper = problem.upper[obj_perm[0]]; + variable_type_t last_type = var_types[obj_perm[0]]; + var_colors[obj_perm[0]] = var_color; + for (i_t k = 1; k < problem.num_cols; k++) { + const i_t j = obj_perm[k]; + const f_t obj = problem.objective[j]; + if (obj - last_obj > tol || problem.lower[j] != last_lower || problem.upper[j] != last_upper || + var_types[j] != last_type) { + var_color++; + last_obj = obj; + last_lower = problem.lower[j]; + last_upper = problem.upper[j]; + last_type = var_types[j]; + } + var_colors[j] = var_color; + } + + // Compute the colors for the constraints + std::vector rhs_perm(problem.num_rows); + std::iota(rhs_perm.begin(), rhs_perm.end(), 0); + std::sort(rhs_perm.begin(), rhs_perm.end(), [&](i_t a, i_t b) { + return problem.rhs[a] < problem.rhs[b]; + }); + std::vector rhs_colors(problem.num_rows, -1); + i_t rhs_color = var_color + 1; + f_t last_rhs = problem.rhs[rhs_perm[0]]; + rhs_colors[rhs_perm[0]] = rhs_color; + for (i_t k = 1; k < problem.num_rows; k++) { + const i_t i = rhs_perm[k]; + const f_t rhs = problem.rhs[i]; + if (rhs - last_rhs > tol) { + rhs_color++; + last_rhs = rhs; + } + rhs_colors[i] = rhs_color; + } + + // Calculate the number of colors needed for each nonzero coefficient in the A matrix + const i_t nnz = problem.A.col_start[problem.num_cols]; + // Construct the graph + // We begin by creating the vertex set := V union R union W + + std::vector vertices; + vertices.reserve(V_size + R_size + nnz); + std::vector vertex_colors; + vertex_colors.reserve(V_size + R_size + nnz); + for (i_t j = 0; j < V_size; j++) { + vertices.push_back(j); + vertex_colors.push_back(var_colors[j]); + } + for (i_t i = 0; i < R_size; i++) { + vertices.push_back(V_size + i); + vertex_colors.push_back(rhs_colors[i]); + } + + std::vector edge_in; + std::vector edge_out; + edge_in.reserve(2 * nnz); + edge_out.reserve(2 * nnz); + +#ifdef FULL_GRAPH + // Every nonzero should have an edge between (v, r) with v in V and r in R + // To handle the edge color we create a new node w and an edge (v, w) and (w, r) + // where w is colored according to the edge color. + + std::vector nonzeros = problem.A.x; + std::vector nonzero_perm(nnz); + std::iota(nonzero_perm.begin(), nonzero_perm.end(), 0); + std::sort(nonzero_perm.begin(), nonzero_perm.end(), [&](i_t a, i_t b) { + return nonzeros[a] < nonzeros[b]; + }); + std::vector nonzero_colors(nnz, -1); + i_t edge_color = rhs_color + 1; + f_t last_nz = nonzeros[nonzero_perm[0]]; + nonzero_colors[nonzero_perm[0]] = edge_color; + for (i_t q = 1; q < nnz; q++) { + const i_t p = nonzero_perm[q]; + const f_t val = nonzeros[p]; + if (val - last_nz > tol) { + edge_color++; + last_nz = val; + } + nonzero_colors[p] = edge_color; + } + + for (i_t j = 0; j < problem.num_cols; j++) { + const i_t col_start = problem.A.col_start[j]; + const i_t col_end = problem.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = problem.A.i[p]; + vertices.push_back(V_size + R_size + p); + vertex_colors.push_back(nonzero_colors[p]); + edge_in.push_back(j); + edge_out.push_back(V_size + R_size + p); + edge_in.push_back(V_size + R_size + p); + edge_out.push_back(V_size + i); + } + } +#else + + // Let r_i be the vertex associated with the row i + // Let V_i,c be the set of variables in row i with the same color c + // We create a new vertex w_i,c and edges (v_j, w_i,c) and (w_i,c, r_i) for all v_j in V_i,c + csr_matrix_t A_row(problem.num_rows, problem.num_cols, 0); + problem.A.to_compressed_row(A_row); + + std::vector nonzeros = A_row.x; + std::vector nonzero_perm(nnz); + std::iota(nonzero_perm.begin(), nonzero_perm.end(), 0); + std::sort(nonzero_perm.begin(), nonzero_perm.end(), [&](i_t a, i_t b) { + return nonzeros[a] < nonzeros[b]; + }); + std::vector nonzero_colors(nnz, -1); + i_t edge_color = 0; + f_t last_nz = nonzeros[nonzero_perm[0]]; + nonzero_colors[nonzero_perm[0]] = edge_color; + for (i_t q = 1; q < nnz; q++) { + const i_t p = nonzero_perm[q]; + const f_t val = nonzeros[p]; + if (val - last_nz > tol) { + edge_color++; + last_nz = val; + } + nonzero_colors[p] = edge_color; + } + + i_t num_edge_colors = edge_color + 1; + std::vector edge_color_map(num_edge_colors, 0); + i_t max_nzs_per_row = 0; + for (i_t i = 0; i < problem.num_rows; i++) { + const i_t row_start = A_row.row_start[i]; + const i_t row_end = A_row.row_start[i + 1]; + max_nzs_per_row = std::max(max_nzs_per_row, row_end - row_start); + } + std::vector row_colors; + std::vector sorted_nonzeros_in_row; + row_colors.reserve(max_nzs_per_row); + sorted_nonzeros_in_row.reserve(max_nzs_per_row); + + i_t current_vertex = V_size + R_size; + for (i_t i = 0; i < problem.num_rows; i++) { + const i_t row_start = A_row.row_start[i]; + const i_t row_end = A_row.row_start[i + 1]; + const i_t row_nz = row_end - row_start; + row_colors.clear(); + sorted_nonzeros_in_row.resize(row_nz); + + // Pass 1: Count the number of occurences of each color and the number of unique colors in the + // current row + for (i_t p = row_start; p < row_end; p++) { + const i_t edge_color = nonzero_colors[p]; + edge_color_map[edge_color]++; + if (edge_color_map[edge_color] == 1) { row_colors.push_back(edge_color); } + } + + // Pass 2: Compute the prefix sum directly in edge_color_map + // Note we only touch the colors that are present in the current row + i_t cumulative_sum = 0; + for (i_t k = 0; k < static_cast(row_colors.size()); k++) { + const i_t edge_color = row_colors[k]; + const i_t count = edge_color_map[edge_color]; + edge_color_map[edge_color] = cumulative_sum; + cumulative_sum += count; + } + + // Pass 3: Place the nonzeros in sorted order + for (i_t p = row_start; p < row_end; p++) { + const i_t edge_color = nonzero_colors[p]; + sorted_nonzeros_in_row[edge_color_map[edge_color]++] = p; + } + + // Clear the edge_color_map for the colors we used + for (i_t edge_color : row_colors) { + edge_color_map[edge_color] = 0; + } + + // Pass 4: iterate over the sorted nonzeros, create new vertices and edges + i_t last_color = -1; + for (i_t k = 0; k < row_nz; k++) { + const i_t p = sorted_nonzeros_in_row[k]; + const i_t j = A_row.j[p]; + const i_t edge_color = nonzero_colors[p]; + if (edge_color == last_color) { + // We don't need to create a new vertex + // Add the edge (v_j, w_i_c) + edge_in.push_back(j); + edge_out.push_back(current_vertex - 1); + } else { + last_color = edge_color; + // Create a new vertex w_i_c + vertices.push_back(current_vertex); + vertex_colors.push_back(rhs_color + 1 + edge_color); + // Add the edge (v_j, w_i_c) + edge_in.push_back(j); + edge_out.push_back(current_vertex); + // Add the edge (w_i_c, r_i) + edge_in.push_back(current_vertex); + edge_out.push_back(V_size + i); + current_vertex++; + } + } + } + +#endif + + settings.log.printf("Graph construction time %f\n", toc(start_time)); + f_t dejavu_start_time = tic(); + + // The graph should now be described by: + // vertices, edge_in, edge_out, vertex_colors + + // Dejavu needs the degree of each vertex + std::vector degrees(vertices.size(), 0); + const i_t num_edges = edge_in.size(); + for (i_t i = 0; i < num_edges; i++) { + degrees[edge_in[i]]++; + degrees[edge_out[i]]++; + } + + dejavu::static_graph g; + g.initialize_graph(vertices.size(), edge_in.size()); + + const i_t num_vertices = vertices.size(); + for (i_t i = 0; i < num_vertices; i++) { + g.add_vertex(vertex_colors[i], degrees[i]); + } + for (i_t i = 0; i < num_edges; i++) { + const i_t u = edge_in[i]; + const i_t v = edge_out[i]; + g.add_edge(std::min(u, v), std::max(u, v)); + } + + dejavu::solver d; + i_t num_generators = 0; + dejavu_hook counting_hook = [&](int, const int*, int, const int*) { num_generators++; }; + d.automorphisms(&g, counting_hook); + + std::ostringstream grp_size_str; + grp_size_str << d.get_automorphism_group_size(); + settings.log.printf( + "Automorphism group size %s, %d generators\n", grp_size_str.str().c_str(), num_generators); + settings.log.printf("Dejavu time %f\n", toc(dejavu_start_time)); + settings.log.printf("Total symmetry detection time %f\n", toc(start_time)); +} + +} // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index be0151665..7e1a16784 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -37,6 +37,10 @@ #include #include +#include +#include +#include + #include #include @@ -303,6 +307,17 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, callback->template setup(op_problem.get_n_variables()); } + // Start symmetry detection + { + detail::problem_t problem(op_problem); + dual_simplex::simplex_solver_settings_t simplex_settings; + simplex_settings.set_log(true); + simplex_settings.time_limit = settings.time_limit; + dual_simplex::user_problem_t user_problem = + cuopt_problem_to_simplex_problem(op_problem.get_handle_ptr(), problem); + dual_simplex::detect_symmetry(user_problem, simplex_settings); + } + auto timer = timer_t(time_limit); if (settings.mip_scaling != CUOPT_MIP_SCALING_OFF) { detail::mip_scaling_strategy_t scaling(op_problem); From ada114b3489a045242fc491d0d6dc305a715c42c Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Wed, 15 Apr 2026 19:15:44 -0700 Subject: [PATCH 02/19] First stab at orbital fixing. Using dejavu for graph operations. Fix a bug in branch and bound. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 166 +++++++++++++++++- cpp/src/branch_and_bound/branch_and_bound.hpp | 7 +- cpp/src/branch_and_bound/symmetry.hpp | 134 +++++++++++++- .../diversity/diversity_manager.cu | 2 +- cpp/src/mip_heuristics/solve.cu | 13 +- cpp/src/mip_heuristics/solver.cu | 3 +- cpp/src/mip_heuristics/solver_context.cuh | 6 + cpp/tests/CMakeLists.txt | 4 +- 8 files changed, 320 insertions(+), 15 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 33a2d983c..7db3edcc4 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -248,11 +249,13 @@ branch_and_bound_t::branch_and_bound_t( const simplex_solver_settings_t& solver_settings, f_t start_time, const probing_implied_bound_t& probing_implied_bound, - std::shared_ptr> clique_table) + std::shared_ptr> clique_table, + mip_symmetry_t* symmetry) : original_problem_(user_problem), settings_(solver_settings), probing_implied_bound_(probing_implied_bound), clique_table_(std::move(clique_table)), + symmetry_(symmetry), original_lp_(user_problem.handle_ptr, 1, 1, 1), Arow_(1, 1, 0), incumbent_(1), @@ -1388,6 +1391,163 @@ dual::status_t branch_and_bound_t::solve_node_lp( worker->leaf_edge_norms = edge_norms_; if (feasible) { + + + // Perform orbital fixing + if (symmetry_ != nullptr) { + // First get the set of variables that have been branched down and branched up on + std::vector branched_zero; + std::vector branched_one; + branched_zero.reserve(node_ptr->depth); + branched_one.reserve(node_ptr->depth); + mip_node_t* node = node_ptr; + while (node != nullptr && node->branch_var >= 0) { + if (node->branch_var_upper == 0.0) { + branched_zero.push_back(node->branch_var); + symmetry_->marked_b0[node->branch_var] = 1; + } else if (node->branch_var_lower == 1.0) { + branched_one.push_back(node->branch_var); + symmetry_->marked_b1[node->branch_var] = 1; + } else { + assert(false); // Unexpected non-binary variable. Only binaries supported in symmetry handling. + } + node = node->parent; + } + + { + for (i_t j = 0; j < symmetry_->num_original_vars; j++) { + if (var_types_[j] == variable_type_t::CONTINUOUS) continue; + if (symmetry_->marked_b1[j] == 0 && worker->leaf_problem.lower[j] == 1.0) { + symmetry_->f1.push_back(j); + symmetry_->marked_f1[j] = 1; + } + if (symmetry_->marked_b0[j] == 0 && worker->leaf_problem.upper[j] == 0.0) { + symmetry_->f0.push_back(j); + symmetry_->marked_f0[j] = 1; + } + } + + // Compute Stab(G, B1) and its orbits + std::vector new_base; + new_base.reserve(symmetry_->num_original_vars); + for (i_t j: branched_one) { + new_base.push_back(j); + symmetry_->marked_variables[j] = 1; + } + for (i_t j = 0; j < symmetry_->num_original_vars; j++) { + if (symmetry_->marked_variables[j] == 0) { + new_base.push_back(j); + } + } + for (i_t j: branched_one) { + symmetry_->marked_variables[j] = 0; + } + + symmetry_->schreier->set_base(new_base); + + dejavu::groups::orbit orb; + orb.initialize(symmetry_->domain_size); + symmetry_->schreier->get_stabilizer_orbit(static_cast(branched_one.size()), orb); + + for (i_t v : branched_one) { + symmetry_->orbit_has_b1[orb.find_orbit(v)] = 1; + } + + for (i_t v : branched_zero) { + symmetry_->orbit_has_b0[orb.find_orbit(v)] = 1; + } + + for (i_t v: symmetry_->continuous_variables) { + symmetry_->orbit_has_continuous[orb.find_orbit(v)] = 1; + } + + for (i_t v: symmetry_->f0) { + symmetry_->orbit_has_f0[orb.find_orbit(v)] = 1; + } + + for (i_t v: symmetry_->f1) { + symmetry_->orbit_has_f1[orb.find_orbit(v)] = 1; + } + + std::vector fix_zero; // The set L0 of variables that can be fixed to 0 + std::vector fix_one; // The set L1 of variables that can be fixed to 1 + + for (i_t j = 0; j < symmetry_->num_original_vars; j++) { + i_t o = orb.find_orbit(j); + if (orb.orbit_size(o) < 2) continue; + + if (symmetry_->orbit_has_b1[o] == 1 || symmetry_->orbit_has_continuous[o] == 1) { + // The orbit contains variables in B1 or continuous variables + // So we can't fix any variables in this orbit to 0 + continue; + } + + if (symmetry_->orbit_has_b0[o] == 1 || symmetry_->orbit_has_f0[o] == 1) { + // The orbit of this variable contains variables in B0 or F0 + // So we can fix this variable to zero (provided its not already in B0 or F0) + if (symmetry_->marked_b0[j] == 0 && symmetry_->marked_f0[j] == 0) { + fix_zero.push_back(j); + } + } + + if (symmetry_->orbit_has_f1[o] == 1) { + // The orbit of this variable contains variables in F1 + // So we can fix this variable to one (provided its not already in F1) + if (symmetry_->marked_f1[j] == 0) { + fix_one.push_back(j); + } + } + } + + // Restore the work arrays + for (i_t v: branched_one) { + symmetry_->orbit_has_b1[orb.find_orbit(v)] = 0; + symmetry_->marked_b1[v] = 0; + } + + for (i_t v: branched_zero) { + symmetry_->orbit_has_b0[orb.find_orbit(v)] = 0; + symmetry_->marked_b0[v] = 0; + } + + for (i_t v: symmetry_->continuous_variables) { + symmetry_->orbit_has_continuous[orb.find_orbit(v)] = 0; + } + + for (i_t v: symmetry_->f0) { + symmetry_->orbit_has_f0[orb.find_orbit(v)] = 0; + symmetry_->marked_f0[v] = 0; + } + + for (i_t v: symmetry_->f1) { + symmetry_->orbit_has_f1[orb.find_orbit(v)] = 0; + symmetry_->marked_f1[v] = 0; + } + + symmetry_->f0.clear(); + symmetry_->f1.clear(); + + settings_.log.printf( + "Orbital fixing at node %d: fixing %d variables to 0 and %d variables to 1\n", + node_ptr->node_id, + fix_zero.size(), + fix_one.size()); + // Finally fix the variables in L0 and L1 + for (i_t v: fix_zero) { + settings_.log.printf("Orbital fixing at node %d: fixing variable %d to 0\n", node_ptr->node_id, v); + worker->leaf_problem.lower[v] = 0.0; + worker->leaf_problem.upper[v] = 0.0; + } + for (i_t v: fix_one) { + settings_.log.printf("Orbital fixing at node %d: fixing variable %d to 1\n", node_ptr->node_id, v); + worker->leaf_problem.lower[v] = 1.0; + worker->leaf_problem.upper[v] = 1.0; + } + } + } + + + i_t node_iter = 0; f_t lp_start_time = tic(); @@ -1438,7 +1598,7 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_trecompute_basis = true; worker->recompute_bounds = true; - f_t lower_bound = get_lower_bound(); + f_t lower_bound = std::min(get_lower_bound(), worker->start_node->lower_bound); f_t upper_bound = upper_bound_; f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); @@ -1582,7 +1742,7 @@ void branch_and_bound_t::dive_with(branch_and_bound_worker_t dive_stats.nodes_explored = 0; dive_stats.nodes_unexplored = 1; - f_t lower_bound = get_lower_bound(); + f_t lower_bound = std::min(get_lower_bound(), worker->start_node->lower_bound); f_t upper_bound = upper_bound_; f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index f2917ba93..3954d5296 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -47,6 +47,9 @@ struct clique_table_t; namespace cuopt::linear_programming::dual_simplex { +template +struct mip_symmetry_t; + enum class mip_status_t { OPTIMAL = 0, // The optimal integer solution was found UNBOUNDED = 1, // The problem is unbounded @@ -80,7 +83,8 @@ class branch_and_bound_t { const simplex_solver_settings_t& solver_settings, f_t start_time, const probing_implied_bound_t& probing_implied_bound, - std::shared_ptr> clique_table = nullptr); + std::shared_ptr> clique_table = nullptr, + mip_symmetry_t* symmetry = nullptr); // Set an initial guess based on the user_problem. This should be called before solve. void set_initial_guess(const std::vector& user_guess) { guess_ = user_guess; } @@ -162,6 +166,7 @@ class branch_and_bound_t { const simplex_solver_settings_t settings_; const probing_implied_bound_t& probing_implied_bound_; std::shared_ptr> clique_table_; + mip_symmetry_t* symmetry_; std::future>> clique_table_future_; std::atomic signal_extend_cliques_{false}; diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 5b1befd8f..44a73b9bd 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -15,14 +15,72 @@ #include "dejavu.h" +#include #include namespace cuopt::linear_programming::dual_simplex { template -void detect_symmetry(const user_problem_t& user_problem, - const simplex_solver_settings_t& settings) +struct mip_symmetry_t { + mip_symmetry_t(std::unique_ptr schreier_, + i_t domain_size_, + i_t num_original_vars_, + const std::vector& var_types) + : schreier(std::move(schreier_)), + domain_size(domain_size_), + num_original_vars(num_original_vars_), + marked_variables(num_original_vars_, 0), + orbit_has_b1(domain_size_, 0), + orbit_has_b0(domain_size_, 0), + orbit_has_continuous(domain_size_, 0), + orbit_has_f0(domain_size_, 0), + orbit_has_f1(domain_size_, 0), + marked_b0(num_original_vars_, 0), + marked_b1(num_original_vars_, 0), + marked_f0(num_original_vars_, 0), + marked_f1(num_original_vars_, 0) + { + for (i_t j = 0; j < num_original_vars_; j++) { + if (var_types[j] == variable_type_t::CONTINUOUS) { + continuous_variables.push_back(j); + } + } + f0.reserve(num_original_vars_); + f1.reserve(num_original_vars_); + } + + std::unique_ptr schreier; + i_t domain_size; + i_t num_original_vars; + std::vector marked_variables; + std::vector orbit_has_b1; + std::vector orbit_has_b0; + std::vector orbit_has_continuous; + std::vector orbit_has_f0; + std::vector orbit_has_f1; + std::vector continuous_variables; + std::vector marked_b0; + std::vector marked_b1; + std::vector f0; + std::vector f1; + std::vector marked_f0; + std::vector marked_f1; +}; + +template +std::unique_ptr> detect_symmetry( + const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + bool& has_symmetry) { + has_symmetry = false; + + for (i_t j = 0; j < user_problem.num_cols; j++) { + if (user_problem.var_types[j] != variable_type_t::CONTINUOUS) { + if (user_problem.lower[j] != 0.0 || user_problem.upper[j] != 1.0) { return nullptr; } + } + } + f_t start_time = tic(); lp_problem_t problem(user_problem.handle_ptr, 1, 1, 1); std::vector new_slacks; @@ -301,17 +359,83 @@ void detect_symmetry(const user_problem_t& user_problem, g.add_edge(std::min(u, v), std::max(u, v)); } + auto schreier = std::make_unique(num_vertices); + std::vector base; + base.reserve(user_problem.num_cols); + for (i_t j = 0; j < user_problem.num_cols; j++) { + base.push_back(j); + } + schreier->set_base(base); + + dejavu::hooks::schreier_hook s_hook(*schreier); + + i_t num_generators = 0; + const i_t num_original_vars = user_problem.num_cols; + dejavu_hook counting_hook = [&num_generators, num_original_vars](int, const int*, int nsupp, const int* supp) { + for (int s = 0; s < nsupp; s++) { + if (supp[s] < num_original_vars) { + num_generators++; + return; + } + } + }; + + dejavu::hooks::multi_hook combined; + combined.add_hook(s_hook.get_hook()); + combined.add_hook(&counting_hook); + dejavu::solver d; - i_t num_generators = 0; - dejavu_hook counting_hook = [&](int, const int*, int, const int*) { num_generators++; }; - d.automorphisms(&g, counting_hook); + d.automorphisms(&g, combined); std::ostringstream grp_size_str; grp_size_str << d.get_automorphism_group_size(); settings.log.printf( "Automorphism group size %s, %d generators\n", grp_size_str.str().c_str(), num_generators); settings.log.printf("Dejavu time %f\n", toc(dejavu_start_time)); + + has_symmetry = false; + if (num_generators > 0) { + dejavu::groups::orbit orb; + orb.initialize(num_vertices); + schreier->get_stabilizer_orbit(0, orb); + + i_t num_nontrivial_orbits = 0; + i_t max_orbit_size = 0; + i_t total_vars_in_orbits = 0; + std::vector orbit_histogram(num_original_vars + 1, 0); + for (i_t j = 0; j < num_original_vars; j++) { + if (orb.represents_orbit(j)) { + i_t sz = orb.orbit_size(j); + if (sz >= 2) { + num_nontrivial_orbits++; + max_orbit_size = std::max(max_orbit_size, sz); + total_vars_in_orbits += sz; + orbit_histogram[sz]++; + } + } + } + settings.log.printf("Orbits: %d non-trivial, max size %d, %d/%d (%.1f%%) variables in orbits\n", + num_nontrivial_orbits, max_orbit_size, total_vars_in_orbits, num_original_vars, + 100.0 * total_vars_in_orbits / num_original_vars); + settings.log.printf("Orbit histogram (orbit size: number of orbits):"); + for (i_t sz = 2; sz <= max_orbit_size; sz++) { + if (orbit_histogram[sz] > 0) { + settings.log.printf(" %d:%d", sz, orbit_histogram[sz]); + } + } + settings.log.printf("\n"); + + has_symmetry = (max_orbit_size >= 4) || + (num_nontrivial_orbits >= 3 && max_orbit_size >= 2) || + (total_vars_in_orbits >= 10); + } + settings.log.printf("Total symmetry detection time %f\n", toc(start_time)); + + if (!has_symmetry) { return nullptr; } + + return std::make_unique>( + std::move(schreier), num_vertices, user_problem.num_cols, var_types); } } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/mip_heuristics/diversity/diversity_manager.cu b/cpp/src/mip_heuristics/diversity/diversity_manager.cu index b8dc3d33b..766fd9805 100644 --- a/cpp/src/mip_heuristics/diversity/diversity_manager.cu +++ b/cpp/src/mip_heuristics/diversity/diversity_manager.cu @@ -223,7 +223,7 @@ template bool diversity_manager_t::run_presolve(f_t time_limit, timer_t global_timer) { raft::common::nvtx::range fun_scope("run_presolve"); - CUOPT_LOG_INFO("Running presolve!"); + CUOPT_LOG_INFO("Starting cuOpt presolve"); timer_t presolve_timer(time_limit); auto term_crit = ls.constraint_prop.bounds_update.solve(*problem_ptr); diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 7e1a16784..af595e57e 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -88,7 +88,8 @@ mip_solution_t run_mip(detail::problem_t& problem, mip_solver_settings_t const& settings, timer_t& timer, f_t& initial_upper_bound, - std::vector& initial_incumbent_assignment) + std::vector& initial_incumbent_assignment, + std::unique_ptr> symmetry = nullptr) { try { raft::common::nvtx::range fun_scope("run_mip"); @@ -165,6 +166,7 @@ mip_solution_t run_mip(detail::problem_t& problem, // It will be converted to the target solver-space at each consumption point. solver.context.initial_upper_bound = initial_upper_bound; solver.context.initial_incumbent_assignment = initial_incumbent_assignment; + solver.context.symmetry = std::move(symmetry); if (timer.check_time_limit()) { CUOPT_LOG_INFO("Time limit reached before main solve"); detail::solution_t sol(problem); @@ -308,6 +310,9 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, } // Start symmetry detection + bool has_symmetry = false; + std::unique_ptr> symmetry; + if (1) { detail::problem_t problem(op_problem); dual_simplex::simplex_solver_settings_t simplex_settings; @@ -315,7 +320,8 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, simplex_settings.time_limit = settings.time_limit; dual_simplex::user_problem_t user_problem = cuopt_problem_to_simplex_problem(op_problem.get_handle_ptr(), problem); - dual_simplex::detect_symmetry(user_problem, simplex_settings); + symmetry = dual_simplex::detect_symmetry(user_problem, simplex_settings, has_symmetry); + if (has_symmetry) { settings.presolver = presolver_t::None; } } auto timer = timer_t(time_limit); @@ -487,7 +493,8 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, // early_best_user_obj is in user-space. // run_mip stores it in context.initial_upper_bound and converts to target spaces as needed. - auto sol = run_mip(problem, settings, timer, early_best_user_obj, early_best_user_assignment); + auto sol = run_mip(problem, settings, timer, early_best_user_obj, early_best_user_assignment, + std::move(symmetry)); const f_t cuopt_presolve_time = sol.get_stats().presolve_time; if (run_presolve) { diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index ce6b602fb..fd71e053e 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -392,7 +392,8 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings, timer_.get_tic_start(), probing_implied_bound, - context.problem_ptr->clique_table); + context.problem_ptr->clique_table, + context.symmetry.get()); context.branch_and_bound_ptr = branch_and_bound.get(); // Convert the best external upper bound from user-space to B&B's internal objective space. diff --git a/cpp/src/mip_heuristics/solver_context.cuh b/cpp/src/mip_heuristics/solver_context.cuh index b1bf3fbd7..739f7d130 100644 --- a/cpp/src/mip_heuristics/solver_context.cuh +++ b/cpp/src/mip_heuristics/solver_context.cuh @@ -13,6 +13,7 @@ #include #include +#include #pragma once @@ -22,6 +23,8 @@ template class branch_and_bound_t; } +#include + namespace cuopt::linear_programming::detail { template @@ -71,6 +74,9 @@ struct mip_solver_context_t { // Matching incumbent assignment in original output space from early heuristics. std::vector initial_incumbent_assignment{}; + + // Symmetry information for orbital fixing during B&B. Null if no exploitable symmetry. + std::unique_ptr> symmetry; }; } // namespace cuopt::linear_programming::detail diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index a73a3361c..71fd12be1 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -1,4 +1,4 @@ -# cmake-format: off +# cmake-format: off # SPDX-FileCopyrightText: Copyright (c) 2021-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 # cmake-format: on @@ -17,6 +17,7 @@ if(BUILD_TESTS) PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/../src" "${CMAKE_CURRENT_SOURCE_DIR}" + "${dejavu_SOURCE_DIR}" ) target_link_libraries(cuopttestutils @@ -51,6 +52,7 @@ function(ConfigureTest CMAKE_TEST_NAME) "${papilo_SOURCE_DIR}/src" "${papilo_BINARY_DIR}" "${pslp_SOURCE_DIR}/include" + "${dejavu_SOURCE_DIR}" ) target_link_libraries(${CMAKE_TEST_NAME} From 6425c25de3f1f4a5c6e2c03b09554f710db31c67 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Thu, 16 Apr 2026 16:33:02 +0200 Subject: [PATCH 03/19] fixed incorrect computation of the lower bound when running with a single thread. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 33a2d983c..e69ff7b9a 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -299,10 +299,11 @@ branch_and_bound_t::branch_and_bound_t( template f_t branch_and_bound_t::get_lower_bound() { - f_t lower_bound = lower_bound_ceiling_.load(); - f_t heap_lower_bound = node_queue_.get_lower_bound(); - lower_bound = std::min(heap_lower_bound, lower_bound); - lower_bound = std::min(worker_pool_.get_lower_bound(), lower_bound); + f_t lower_bound = lower_bound_ceiling_.load(); + f_t heap_lower_bound = node_queue_.get_lower_bound(); + f_t worker_lower_bound = worker_pool_.get_lower_bound(); + lower_bound = std::min(heap_lower_bound, lower_bound); + lower_bound = std::min(worker_lower_bound, lower_bound); if (std::isfinite(lower_bound)) { return lower_bound; @@ -1800,7 +1801,8 @@ void branch_and_bound_t::run_scheduler() template void branch_and_bound_t::single_threaded_solve() { - branch_and_bound_worker_t worker(0, original_lp_, Arow_, var_types_, settings_); + worker_pool_.init(1, original_lp_, Arow_, var_types_, settings_); + branch_and_bound_worker_t* worker = worker_pool_.get_idle_worker(); f_t lower_bound = get_lower_bound(); f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound); @@ -1808,7 +1810,6 @@ void branch_and_bound_t::single_threaded_solve() while (solver_status_ == mip_status_t::UNSET && abs_gap > settings_.absolute_mip_gap_tol && rel_gap > settings_.relative_mip_gap_tol && node_queue_.best_first_queue_size() > 0) { - bool launched_any_task = false; repair_heuristic_solutions(); f_t now = toc(exploration_stats_.start_time); @@ -1844,8 +1845,8 @@ void branch_and_bound_t::single_threaded_solve() continue; } - worker.init_best_first(start_node.value(), original_lp_); - plunge_with(&worker); + worker->init_best_first(start_node.value(), original_lp_); + plunge_with(worker); lower_bound = get_lower_bound(); abs_gap = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound); From f447d4c32cbe3f38b73166ad4b41fff8c1c7b53d Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 16 Apr 2026 09:30:07 -0700 Subject: [PATCH 04/19] Revert changes for single-thread MIP hang --- cpp/src/branch_and_bound/branch_and_bound.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 7db3edcc4..dc7d1f4bb 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1598,7 +1598,7 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_trecompute_basis = true; worker->recompute_bounds = true; - f_t lower_bound = std::min(get_lower_bound(), worker->start_node->lower_bound); + f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); @@ -1742,7 +1742,7 @@ void branch_and_bound_t::dive_with(branch_and_bound_worker_t dive_stats.nodes_explored = 0; dive_stats.nodes_unexplored = 1; - f_t lower_bound = std::min(get_lower_bound(), worker->start_node->lower_bound); + f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; f_t rel_gap = user_relative_gap(original_lp_, upper_bound, lower_bound); f_t abs_gap = compute_user_abs_gap(original_lp_, upper_bound, lower_bound); From fcd28d763c6159779bbc7d68f4a6ccd4414a125b Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 16 Apr 2026 16:41:30 -0700 Subject: [PATCH 05/19] Add mip-symmetry setting. Give each worker a unique orbital_fixing_t object. Make orbital fixing work with multiple threads. Allow for general integers (but only fix binary). --- .../cuopt/linear_programming/constants.h | 1 + .../mip/solver_settings.hpp | 1 + cpp/src/branch_and_bound/branch_and_bound.cpp | 162 +--------- .../branch_and_bound_worker.hpp | 8 + cpp/src/branch_and_bound/symmetry.hpp | 298 ++++++++++++++++-- cpp/src/math_optimization/solver_settings.cu | 1 + cpp/src/mip_heuristics/solve.cu | 4 +- 7 files changed, 286 insertions(+), 189 deletions(-) diff --git a/cpp/include/cuopt/linear_programming/constants.h b/cpp/include/cuopt/linear_programming/constants.h index 06eacb340..38ed7705d 100644 --- a/cpp/include/cuopt/linear_programming/constants.h +++ b/cpp/include/cuopt/linear_programming/constants.h @@ -59,6 +59,7 @@ #define CUOPT_MIP_HEURISTICS_ONLY "mip_heuristics_only" #define CUOPT_MIP_SCALING "mip_scaling" #define CUOPT_MIP_PRESOLVE "mip_presolve" +#define CUOPT_MIP_SYMMETRY "mip_symmetry" #define CUOPT_MIP_RELIABILITY_BRANCHING "mip_reliability_branching" #define CUOPT_MIP_CUT_PASSES "mip_cut_passes" #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts" diff --git a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp index 14c4d227b..224fc0a1e 100644 --- a/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp +++ b/cpp/include/cuopt/linear_programming/mip/solver_settings.hpp @@ -90,6 +90,7 @@ class mip_solver_settings_t { bool heuristics_only = false; i_t reliability_branching = -1; i_t num_cpu_threads = -1; // -1 means use default number of threads in branch and bound + i_t symmetry = -1; i_t max_cut_passes = 10; // number of cut passes to make i_t mir_cuts = -1; i_t mixed_integer_gomory_cuts = -1; diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 76a9be9e8..62b816bcc 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1393,162 +1393,12 @@ dual::status_t branch_and_bound_t::solve_node_lp( if (feasible) { - // Perform orbital fixing - if (symmetry_ != nullptr) { - // First get the set of variables that have been branched down and branched up on - std::vector branched_zero; - std::vector branched_one; - branched_zero.reserve(node_ptr->depth); - branched_one.reserve(node_ptr->depth); - mip_node_t* node = node_ptr; - while (node != nullptr && node->branch_var >= 0) { - if (node->branch_var_upper == 0.0) { - branched_zero.push_back(node->branch_var); - symmetry_->marked_b0[node->branch_var] = 1; - } else if (node->branch_var_lower == 1.0) { - branched_one.push_back(node->branch_var); - symmetry_->marked_b1[node->branch_var] = 1; - } else { - assert(false); // Unexpected non-binary variable. Only binaries supported in symmetry handling. - } - node = node->parent; - } - - { - for (i_t j = 0; j < symmetry_->num_original_vars; j++) { - if (var_types_[j] == variable_type_t::CONTINUOUS) continue; - if (symmetry_->marked_b1[j] == 0 && worker->leaf_problem.lower[j] == 1.0) { - symmetry_->f1.push_back(j); - symmetry_->marked_f1[j] = 1; - } - if (symmetry_->marked_b0[j] == 0 && worker->leaf_problem.upper[j] == 0.0) { - symmetry_->f0.push_back(j); - symmetry_->marked_f0[j] = 1; - } - } - - // Compute Stab(G, B1) and its orbits - std::vector new_base; - new_base.reserve(symmetry_->num_original_vars); - for (i_t j: branched_one) { - new_base.push_back(j); - symmetry_->marked_variables[j] = 1; - } - for (i_t j = 0; j < symmetry_->num_original_vars; j++) { - if (symmetry_->marked_variables[j] == 0) { - new_base.push_back(j); - } - } - for (i_t j: branched_one) { - symmetry_->marked_variables[j] = 0; - } - - symmetry_->schreier->set_base(new_base); - - dejavu::groups::orbit orb; - orb.initialize(symmetry_->domain_size); - symmetry_->schreier->get_stabilizer_orbit(static_cast(branched_one.size()), orb); - - for (i_t v : branched_one) { - symmetry_->orbit_has_b1[orb.find_orbit(v)] = 1; - } - - for (i_t v : branched_zero) { - symmetry_->orbit_has_b0[orb.find_orbit(v)] = 1; - } - - for (i_t v: symmetry_->continuous_variables) { - symmetry_->orbit_has_continuous[orb.find_orbit(v)] = 1; - } - - for (i_t v: symmetry_->f0) { - symmetry_->orbit_has_f0[orb.find_orbit(v)] = 1; - } - - for (i_t v: symmetry_->f1) { - symmetry_->orbit_has_f1[orb.find_orbit(v)] = 1; - } - - std::vector fix_zero; // The set L0 of variables that can be fixed to 0 - std::vector fix_one; // The set L1 of variables that can be fixed to 1 - - for (i_t j = 0; j < symmetry_->num_original_vars; j++) { - i_t o = orb.find_orbit(j); - if (orb.orbit_size(o) < 2) continue; - - if (symmetry_->orbit_has_b1[o] == 1 || symmetry_->orbit_has_continuous[o] == 1) { - // The orbit contains variables in B1 or continuous variables - // So we can't fix any variables in this orbit to 0 - continue; - } - - if (symmetry_->orbit_has_b0[o] == 1 || symmetry_->orbit_has_f0[o] == 1) { - // The orbit of this variable contains variables in B0 or F0 - // So we can fix this variable to zero (provided its not already in B0 or F0) - if (symmetry_->marked_b0[j] == 0 && symmetry_->marked_f0[j] == 0) { - fix_zero.push_back(j); - } - } - - if (symmetry_->orbit_has_f1[o] == 1) { - // The orbit of this variable contains variables in F1 - // So we can fix this variable to one (provided its not already in F1) - if (symmetry_->marked_f1[j] == 0) { - fix_one.push_back(j); - } - } - } - - // Restore the work arrays - for (i_t v: branched_one) { - symmetry_->orbit_has_b1[orb.find_orbit(v)] = 0; - symmetry_->marked_b1[v] = 0; - } - - for (i_t v: branched_zero) { - symmetry_->orbit_has_b0[orb.find_orbit(v)] = 0; - symmetry_->marked_b0[v] = 0; - } - - for (i_t v: symmetry_->continuous_variables) { - symmetry_->orbit_has_continuous[orb.find_orbit(v)] = 0; - } - - for (i_t v: symmetry_->f0) { - symmetry_->orbit_has_f0[orb.find_orbit(v)] = 0; - symmetry_->marked_f0[v] = 0; - } - - for (i_t v: symmetry_->f1) { - symmetry_->orbit_has_f1[orb.find_orbit(v)] = 0; - symmetry_->marked_f1[v] = 0; - } - - symmetry_->f0.clear(); - symmetry_->f1.clear(); - - settings_.log.printf( - "Orbital fixing at node %d: fixing %d variables to 0 and %d variables to 1\n", - node_ptr->node_id, - fix_zero.size(), - fix_one.size()); - // Finally fix the variables in L0 and L1 - for (i_t v: fix_zero) { - settings_.log.printf("Orbital fixing at node %d: fixing variable %d to 0\n", node_ptr->node_id, v); - worker->leaf_problem.lower[v] = 0.0; - worker->leaf_problem.upper[v] = 0.0; - } - for (i_t v: fix_one) { - settings_.log.printf("Orbital fixing at node %d: fixing variable %d to 1\n", node_ptr->node_id, v); - worker->leaf_problem.lower[v] = 1.0; - worker->leaf_problem.upper[v] = 1.0; - } - } + auto* of = worker->orbital_fixing.get(); + if (of != nullptr && !of->disabled()) { + of->orbital_fixing(symmetry_, settings_, var_types_, node_ptr, worker->leaf_problem); } - - i_t node_iter = 0; f_t lp_start_time = tic(); @@ -1598,6 +1448,7 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_tstart_node); worker->recompute_basis = true; worker->recompute_bounds = true; + if (worker->orbital_fixing) { worker->orbital_fixing->enable(); } f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; @@ -1723,6 +1574,7 @@ template void branch_and_bound_t::dive_with(branch_and_bound_worker_t* worker) { raft::common::nvtx::range scope("BB::diving_thread"); + if (worker->orbital_fixing) { worker->orbital_fixing->disable(); } logger_t log; log.log = false; @@ -1819,7 +1671,7 @@ void branch_and_bound_t::run_scheduler() std::array max_num_workers_per_type = get_max_workers(num_workers, strategies); - worker_pool_.init(num_workers, original_lp_, Arow_, var_types_, settings_); + worker_pool_.init(num_workers, original_lp_, Arow_, var_types_, symmetry_, settings_); active_workers_per_strategy_.fill(0); #ifdef CUOPT_LOG_DEBUG @@ -1961,7 +1813,7 @@ void branch_and_bound_t::run_scheduler() template void branch_and_bound_t::single_threaded_solve() { - worker_pool_.init(1, original_lp_, Arow_, var_types_, settings_); + worker_pool_.init(1, original_lp_, Arow_, var_types_, symmetry_, settings_); branch_and_bound_worker_t* worker = worker_pool_.get_idle_worker(); f_t lower_bound = get_lower_bound(); diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index 4de2b43ca..6986c2692 100644 --- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp @@ -8,6 +8,7 @@ #pragma once #include +#include #include #include @@ -73,6 +74,8 @@ class branch_and_bound_worker_t { pcgenerator_t rng; + std::unique_ptr> orbital_fixing; + bool recompute_basis = true; bool recompute_bounds = true; @@ -169,6 +172,7 @@ class branch_and_bound_worker_pool_t { const lp_problem_t& original_lp, const csr_matrix_t& Arow, const std::vector& var_type, + mip_symmetry_t* symmetry, const simplex_solver_settings_t& settings) { workers_.resize(num_workers); @@ -176,6 +180,10 @@ class branch_and_bound_worker_pool_t { for (i_t i = 0; i < num_workers; ++i) { workers_[i] = std::make_unique>( i, original_lp, Arow, var_type, settings); + if (symmetry != nullptr) { + workers_[i]->orbital_fixing = + std::make_unique>(*symmetry); + } idle_workers_.push_front(i); } diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 44a73b9bd..da13193a2 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -12,10 +12,13 @@ #include #include #include +#include + #include "dejavu.h" #include +#include #include namespace cuopt::linear_programming::dual_simplex { @@ -25,46 +28,282 @@ struct mip_symmetry_t { mip_symmetry_t(std::unique_ptr schreier_, i_t domain_size_, i_t num_original_vars_, - const std::vector& var_types) + const std::vector& var_types, + const std::vector& lower, + const std::vector& upper) : schreier(std::move(schreier_)), domain_size(domain_size_), num_original_vars(num_original_vars_), - marked_variables(num_original_vars_, 0), - orbit_has_b1(domain_size_, 0), - orbit_has_b0(domain_size_, 0), - orbit_has_continuous(domain_size_, 0), - orbit_has_f0(domain_size_, 0), - orbit_has_f1(domain_size_, 0), - marked_b0(num_original_vars_, 0), - marked_b1(num_original_vars_, 0), - marked_f0(num_original_vars_, 0), - marked_f1(num_original_vars_, 0) + is_binary(num_original_vars_, 0) { for (i_t j = 0; j < num_original_vars_; j++) { if (var_types[j] == variable_type_t::CONTINUOUS) { continuous_variables.push_back(j); + } else if (lower[j] == 0.0 && upper[j] == 1.0) { + is_binary[j] = 1; + } else { + general_integer_variables.push_back(j); } } - f0.reserve(num_original_vars_); - f1.reserve(num_original_vars_); } std::unique_ptr schreier; i_t domain_size; i_t num_original_vars; - std::vector marked_variables; - std::vector orbit_has_b1; - std::vector orbit_has_b0; - std::vector orbit_has_continuous; - std::vector orbit_has_f0; - std::vector orbit_has_f1; std::vector continuous_variables; - std::vector marked_b0; - std::vector marked_b1; - std::vector f0; - std::vector f1; - std::vector marked_f0; - std::vector marked_f1; + std::vector general_integer_variables; + std::vector is_binary; +}; + +template +class orbital_fixing_t { + public: + explicit orbital_fixing_t(mip_symmetry_t& root) + : domain_size_(root.domain_size), + num_original_vars_(root.num_original_vars), + marked_variables_(root.num_original_vars, 0), + orbit_has_b1_(root.domain_size, 0), + orbit_has_b0_(root.domain_size, 0), + orbit_has_continuous_(root.domain_size, 0), + orbit_has_f0_(root.domain_size, 0), + orbit_has_f1_(root.domain_size, 0), + marked_b0_(root.num_original_vars, 0), + marked_b1_(root.num_original_vars, 0), + marked_f0_(root.num_original_vars, 0), + marked_f1_(root.num_original_vars, 0) + { + f0_.reserve(root.num_original_vars); + f1_.reserve(root.num_original_vars); + + // Clone the root Schreier structure into a per-worker copy. The root structure + // is shared and read-only, but set_base/get_stabilizer_orbit mutate internal state, + // so each worker needs its own copy. Since random_schreier is not copyable, we + // reconstruct it by creating a fresh structure with the same base and sifting all + // generators from the root into it. + schreier_ = std::make_unique(root.domain_size); + std::vector base(root.num_original_vars); + std::iota(base.begin(), base.end(), 0); + schreier_->set_base(base); + + int num_gens = root.schreier->get_number_of_generators(); + dejavu::groups::automorphism_workspace ws(root.domain_size); + for (int i = 0; i < num_gens; i++) { + root.schreier->get_generator(i, ws); + schreier_->sift(ws); + } + } + + bool disabled () const { return disabled_; } + void disable() { disabled_ = true; } + void enable() { disabled_ = false; } + + void orbital_fixing(mip_symmetry_t* symmetry, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + mip_node_t* node_ptr, + lp_problem_t& problem) + { + // Collect binary branchings only; skip general integer branchings + std::vector branched_zero; + std::vector branched_one; + branched_zero.reserve(node_ptr->depth); + branched_one.reserve(node_ptr->depth); + mip_node_t* node = node_ptr; + while (node != nullptr && node->branch_var >= 0) { + i_t v = node->branch_var; + bool is_binary = (symmetry->is_binary[v] == 1); + if (is_binary) { + if (node->branch_var_upper == 0.0) { + branched_zero.push_back(v); + marked_b0_[v] = 1; + } else if (node->branch_var_lower == 1.0) { + branched_one.push_back(v); + marked_b1_[v] = 1; + } + } + node = node->parent; + } + + for (i_t j = 0; j < num_original_vars_; j++) { + if (var_types[j] == variable_type_t::CONTINUOUS) continue; + if (marked_b1_[j] == 0 && problem.lower[j] == 1.0) { + f1_.push_back(j); + marked_f1_[j] = 1; + } + if (marked_b0_[j] == 0 && problem.upper[j] == 0.0) { + f0_.push_back(j); + marked_f0_[j] = 1; + } + } + + // Compute Stab(G', B1) where G' = Stab(G, general_integer_variables) + // Base order: [general_int_vars, B1, remaining_vars] + const auto& gen_int = symmetry->general_integer_variables; + std::vector new_base; + new_base.reserve(num_original_vars_); + for (i_t j : gen_int) { + new_base.push_back(j); + marked_variables_[j] = 1; + } + for (i_t j : branched_one) { + new_base.push_back(j); + marked_variables_[j] = 1; + } + for (i_t j = 0; j < num_original_vars_; j++) { + if (marked_variables_[j] == 0) { new_base.push_back(j); } + } + for (i_t j : gen_int) { + marked_variables_[j] = 0; + } + for (i_t j : branched_one) { + marked_variables_[j] = 0; + } + + schreier_->set_base(new_base); + + dejavu::groups::orbit orb; + orb.initialize(domain_size_); + i_t stabilizer_level = static_cast(gen_int.size() + branched_one.size()); + schreier_->get_stabilizer_orbit(static_cast(stabilizer_level), orb); + + // Check if the stabilizer is trivial (all orbits are singletons) + bool trivial_stabilizer = true; + for (i_t j = 0; j < num_original_vars_; j++) { + if (orb.orbit_size(orb.find_orbit(j)) > 1) { + trivial_stabilizer = false; + break; + } + } + + if (trivial_stabilizer) { + disabled_ = true; + for (i_t v : branched_one) { + marked_b1_[v] = 0; + } + for (i_t v : branched_zero) { + marked_b0_[v] = 0; + } + for (i_t v : f0_) { + marked_f0_[v] = 0; + } + for (i_t v : f1_) { + marked_f1_[v] = 0; + } + f0_.clear(); + f1_.clear(); + } else { + for (i_t v : branched_one) { + orbit_has_b1_[orb.find_orbit(v)] = 1; + } + + for (i_t v : branched_zero) { + orbit_has_b0_[orb.find_orbit(v)] = 1; + } + + for (i_t v : symmetry->continuous_variables) { + orbit_has_continuous_[orb.find_orbit(v)] = 1; + } + + for (i_t v : f0_) { + orbit_has_f0_[orb.find_orbit(v)] = 1; + } + + for (i_t v : f1_) { + orbit_has_f1_[orb.find_orbit(v)] = 1; + } + + std::vector fix_zero; // The set L0 of variables that can be fixed to 0 + std::vector fix_one; // The set L1 of variables that can be fixed to 1 + + for (i_t j = 0; j < num_original_vars_; j++) { + i_t o = orb.find_orbit(j); + if (orb.orbit_size(o) < 2) continue; + + if (orbit_has_b1_[o] == 1 || orbit_has_continuous_[o] == 1) { + // The orbit contains variables in B1 or continuous variables + // So we can't fix any variables in this orbit + continue; + } + + if (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1) { + // The orbit of this variable contains variables in B0 or F0 + // So we can fix this variable to zero (provided its not already in B0 or F0) + if (marked_b0_[j] == 0 && marked_f0_[j] == 0) { fix_zero.push_back(j); } + } + + if (orbit_has_f1_[o] == 1) { + // The orbit of this variable contains variables in F1 + // So we can fix this variable to one (provided its not already in F1) + if (marked_f1_[j] == 0) { fix_one.push_back(j); } + } + } + + // Restore the work arrays + for (i_t v : branched_one) { + orbit_has_b1_[orb.find_orbit(v)] = 0; + marked_b1_[v] = 0; + } + + for (i_t v : branched_zero) { + orbit_has_b0_[orb.find_orbit(v)] = 0; + marked_b0_[v] = 0; + } + + for (i_t v : symmetry->continuous_variables) { + orbit_has_continuous_[orb.find_orbit(v)] = 0; + } + + for (i_t v : f0_) { + orbit_has_f0_[orb.find_orbit(v)] = 0; + marked_f0_[v] = 0; + } + + for (i_t v : f1_) { + orbit_has_f1_[orb.find_orbit(v)] = 0; + marked_f1_[v] = 0; + } + + f0_.clear(); + f1_.clear(); + + if (fix_zero.size() > 0 || fix_one.size() > 0) { + settings.log.printf("Orbital fixing at node %d (depth %d): fixing %d to 0 and %d to 1\n", + node_ptr->node_id, + node_ptr->depth, + (int)fix_zero.size(), + (int)fix_one.size()); + } + + // Finally fix the variables in L0 and L1 + for (i_t v : fix_zero) { + problem.lower[v] = 0.0; + problem.upper[v] = 0.0; + } + for (i_t v : fix_one) { + problem.lower[v] = 1.0; + problem.upper[v] = 1.0; + } + } + } + + private: + std::unique_ptr schreier_; + i_t domain_size_; + i_t num_original_vars_; + bool disabled_ = false; + + std::vector marked_variables_; + std::vector orbit_has_b1_; // orbit_has_b1_[o] = 1 if the orbit o contains variables in B1 + std::vector orbit_has_b0_; // orbit_has_b0_[o] = 1 if the orbit o contains variables in B0 + std::vector orbit_has_continuous_; // orbit_has_continuous_[o] = 1 if the orbit o contains continuous variables + std::vector orbit_has_f0_; // orbit_has_f0_[o] = 1 if the orbit o contains variables in F0 + std::vector orbit_has_f1_; // orbit_has_f1_[o] = 1 if the orbit o contains variables in F1 + std::vector marked_b0_; // marked_b0_[v] = 1 if the variable v has been branched down + std::vector marked_b1_; // marked_b1_[v] = 1 if the variable v has been branched up + std::vector f0_; // set of variables that can be fixed to 0 + std::vector f1_; // set of variables that can be fixed to 1 + std::vector marked_f0_; // marked_f0_[v] = 1 if the variable v has been branched down + std::vector marked_f1_; // marked_f1_[v] = 1 if the variable v has been branched up }; template @@ -75,12 +314,6 @@ std::unique_ptr> detect_symmetry( { has_symmetry = false; - for (i_t j = 0; j < user_problem.num_cols; j++) { - if (user_problem.var_types[j] != variable_type_t::CONTINUOUS) { - if (user_problem.lower[j] != 0.0 || user_problem.upper[j] != 1.0) { return nullptr; } - } - } - f_t start_time = tic(); lp_problem_t problem(user_problem.handle_ptr, 1, 1, 1); std::vector new_slacks; @@ -435,7 +668,8 @@ std::unique_ptr> detect_symmetry( if (!has_symmetry) { return nullptr; } return std::make_unique>( - std::move(schreier), num_vertices, user_problem.num_cols, var_types); + std::move(schreier), num_vertices, user_problem.num_cols, var_types, + user_problem.lower, user_problem.upper); } } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index c23b1d27c..78fd1f60e 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -147,6 +147,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -1, std::numeric_limits::max(), -1}, {CUOPT_PDLP_PRECISION, reinterpret_cast(&pdlp_settings.pdlp_precision), CUOPT_PDLP_DEFAULT_PRECISION, CUOPT_PDLP_MIXED_PRECISION, CUOPT_PDLP_DEFAULT_PRECISION}, {CUOPT_MIP_SCALING, &mip_settings.mip_scaling, CUOPT_MIP_SCALING_OFF, CUOPT_MIP_SCALING_NO_OBJECTIVE, CUOPT_MIP_SCALING_ON}, + {CUOPT_MIP_SYMMETRY, &mip_settings.symmetry, -1, 1, -1}, // MIP heuristic hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_HEURISTIC_POPULATION_SIZE, &mip_settings.heuristic_params.population_size, 1, std::numeric_limits::max(), 32, "max solutions in pool"}, {CUOPT_MIP_HYPER_HEURISTIC_NUM_CPUFJ_THREADS, &mip_settings.heuristic_params.num_cpufj_threads, 0, std::numeric_limits::max(), 8, "parallel CPU FJ climbers"}, diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index af595e57e..4db1a5908 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -309,10 +309,11 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, callback->template setup(op_problem.get_n_variables()); } + auto timer = timer_t(time_limit); // Start symmetry detection bool has_symmetry = false; std::unique_ptr> symmetry; - if (1) + if (settings.symmetry != 0) { detail::problem_t problem(op_problem); dual_simplex::simplex_solver_settings_t simplex_settings; @@ -324,7 +325,6 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, if (has_symmetry) { settings.presolver = presolver_t::None; } } - auto timer = timer_t(time_limit); if (settings.mip_scaling != CUOPT_MIP_SCALING_OFF) { detail::mip_scaling_strategy_t scaling(op_problem); scaling.scale_problem(settings.mip_scaling != CUOPT_MIP_SCALING_NO_OBJECTIVE); From 96010da413de38c84343d3848624a15b0654cb57 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 17 Apr 2026 11:25:49 -0700 Subject: [PATCH 06/19] Only do strong branching on representative variables in an orbit. Extract the subgroup on the variables from the symmetry group on the graph --- cpp/src/branch_and_bound/branch_and_bound.cpp | 2 + .../branch_and_bound_worker.hpp | 13 +- cpp/src/branch_and_bound/pseudo_costs.cpp | 122 ++++++++-- cpp/src/branch_and_bound/pseudo_costs.hpp | 4 + cpp/src/branch_and_bound/symmetry.hpp | 221 ++++++++++++------ 5 files changed, 274 insertions(+), 88 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 62b816bcc..f5892ca5f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1448,6 +1448,7 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_tstart_node); worker->recompute_basis = true; worker->recompute_bounds = true; + worker->ensure_orbital_fixing(); if (worker->orbital_fixing) { worker->orbital_fixing->enable(); } f_t lower_bound = get_lower_bound(); @@ -2536,6 +2537,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basic_list, nonbasic_list, basis_update, + symmetry_, pc_); } diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index 6986c2692..3d49fb7b1 100644 --- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp @@ -75,6 +75,14 @@ class branch_and_bound_worker_t { pcgenerator_t rng; std::unique_ptr> orbital_fixing; + mip_symmetry_t* symmetry_ptr = nullptr; + + void ensure_orbital_fixing() + { + if (orbital_fixing == nullptr && symmetry_ptr != nullptr) { + orbital_fixing = std::make_unique>(*symmetry_ptr); + } + } bool recompute_basis = true; bool recompute_bounds = true; @@ -180,10 +188,7 @@ class branch_and_bound_worker_pool_t { for (i_t i = 0; i < num_workers; ++i) { workers_[i] = std::make_unique>( i, original_lp, Arow, var_type, settings); - if (symmetry != nullptr) { - workers_[i]->orbital_fixing = - std::make_unique>(*symmetry); - } + workers_[i]->symmetry_ptr = symmetry; idle_workers_.push_front(i); } diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index c38e98e27..c0d1e7ca9 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -7,6 +7,7 @@ #include #include +#include #include #include @@ -983,21 +984,21 @@ static void batch_pdlp_reliability_branching_task( } template -void strong_branching(const lp_problem_t& original_lp, - const simplex_solver_settings_t& settings, - f_t start_time, - const std::vector& new_slacks, - const std::vector& var_types, - const lp_solution_t& root_solution, - const std::vector& fractional, - f_t root_obj, - f_t upper_bound, - const std::vector& root_vstatus, - const std::vector& edge_norms, - const std::vector& basic_list, - const std::vector& nonbasic_list, - basis_update_mpf_t& basis_factors, - pseudo_costs_t& pc) +void strong_branching_neccesary(const lp_problem_t& original_lp, + const simplex_solver_settings_t& settings, + f_t start_time, + const std::vector& new_slacks, + const std::vector& var_types, + const lp_solution_t& root_solution, + const std::vector& fractional, + f_t root_obj, + f_t upper_bound, + const std::vector& root_vstatus, + const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, + pseudo_costs_t& pc) { constexpr bool verbose = false; @@ -1226,7 +1227,95 @@ void strong_branching(const lp_problem_t& original_lp, } } - pc.update_pseudo_costs_from_strong_branching(fractional, root_solution.x); + +} + +template +void strong_branching(const lp_problem_t& original_lp, + const simplex_solver_settings_t& settings, + f_t start_time, + const std::vector& new_slacks, + const std::vector& var_types, + const lp_solution_t& root_solution, + const std::vector& fractional, + f_t root_obj, + f_t upper_bound, + const std::vector& root_vstatus, + const std::vector& edge_norms, + const std::vector& basic_list, + const std::vector& nonbasic_list, + basis_update_mpf_t& basis_factors, + mip_symmetry_t* symmetry, + pseudo_costs_t& pc) +{ + if (symmetry == nullptr) { + strong_branching_neccesary(original_lp, + settings, + start_time, + new_slacks, + var_types, + root_solution, + fractional, + root_obj, + upper_bound, + root_vstatus, + edge_norms, + basic_list, + nonbasic_list, + basis_factors, + pc); + pc.update_pseudo_costs_from_strong_branching(fractional, root_solution.x); + } else { + // Use precomputed orbit representatives from the full group G. Variables + // in the same orbit have identical branching behavior by symmetry. + // Keep only one fractional variable per orbit. We use the first fractional + // variable encountered as the orbit's delegate (not the orbit representative, + // which might not be fractional itself). + // orbit_delegate[rep] maps each orbit representative to the fractional delegate. + std::vector orbit_delegate(symmetry->num_original_vars, -1); + std::vector necessary_fractional; + for (i_t j : fractional) { + i_t rep = symmetry->orbit_rep[j]; + if (orbit_delegate[rep] == -1) { + orbit_delegate[rep] = j; + necessary_fractional.push_back(j); + } + } + + if (necessary_fractional.size() < fractional.size()) { + settings.log.printf("Strong branching: %d fractional variables reduced to %d by symmetry\n", + (int)fractional.size(), (int)necessary_fractional.size()); + } + + strong_branching_neccesary(original_lp, + settings, + start_time, + new_slacks, + var_types, + root_solution, + necessary_fractional, + root_obj, + upper_bound, + root_vstatus, + edge_norms, + basic_list, + nonbasic_list, + basis_factors, + pc); + + // Copy the pseudocosts from each orbit's delegate to all other + // fractional variables in the same orbit + for (i_t j : fractional) { + i_t delegate = orbit_delegate[symmetry->orbit_rep[j]]; + if (j == delegate) continue; + pc.pseudo_cost_sum_down[j] = pc.pseudo_cost_sum_down[delegate]; + pc.pseudo_cost_sum_up[j] = pc.pseudo_cost_sum_up[delegate]; + pc.pseudo_cost_num_down[j] = pc.pseudo_cost_num_down[delegate]; + pc.pseudo_cost_num_up[j] = pc.pseudo_cost_num_up[delegate]; + } + + pc.update_pseudo_costs_from_strong_branching(fractional, root_solution.x); + } } template @@ -1850,6 +1939,7 @@ template void strong_branching(const lp_problem_t& ori const std::vector& basic_list, const std::vector& nonbasic_list, basis_update_mpf_t& basis_factors, + mip_symmetry_t* symmetry, pseudo_costs_t& pc); #endif diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 009bd8b81..f07739d95 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -27,6 +27,9 @@ namespace cuopt::linear_programming::dual_simplex { +template +struct mip_symmetry_t; + template struct branch_variable_t { i_t variable; @@ -556,6 +559,7 @@ void strong_branching(const lp_problem_t& original_lp, const std::vector& basic_list, const std::vector& nonbasic_list, basis_update_mpf_t& basis_factors, + mip_symmetry_t* symmetry, pseudo_costs_t& pc); } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index da13193a2..b08889c6b 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -18,6 +18,7 @@ #include "dejavu.h" #include +#include #include #include @@ -25,48 +26,32 @@ namespace cuopt::linear_programming::dual_simplex { template struct mip_symmetry_t { - mip_symmetry_t(std::unique_ptr schreier_, - i_t domain_size_, - i_t num_original_vars_, - const std::vector& var_types, - const std::vector& lower, - const std::vector& upper) - : schreier(std::move(schreier_)), - domain_size(domain_size_), - num_original_vars(num_original_vars_), - is_binary(num_original_vars_, 0) - { - for (i_t j = 0; j < num_original_vars_; j++) { - if (var_types[j] == variable_type_t::CONTINUOUS) { - continuous_variables.push_back(j); - } else if (lower[j] == 0.0 && upper[j] == 1.0) { - is_binary[j] = 1; - } else { - general_integer_variables.push_back(j); - } - } - } - + // Schreier projected onto original variables (domain = num_original_vars). + // Generators that permute continuous variables are excluded. std::unique_ptr schreier; - i_t domain_size; i_t num_original_vars; - std::vector continuous_variables; + int num_generators = 0; std::vector general_integer_variables; std::vector is_binary; + + // Precomputed orbit representative for each original variable under the projected group. + // orbit_rep[j] = orbit representative of variable j (for j < num_original_vars). + std::vector orbit_rep; + + // Mutex protecting schreier->get_generator() for concurrent orbital_fixing_t construction. + std::mutex generator_mutex; }; template class orbital_fixing_t { public: explicit orbital_fixing_t(mip_symmetry_t& root) - : domain_size_(root.domain_size), - num_original_vars_(root.num_original_vars), + : num_original_vars_(root.num_original_vars), marked_variables_(root.num_original_vars, 0), - orbit_has_b1_(root.domain_size, 0), - orbit_has_b0_(root.domain_size, 0), - orbit_has_continuous_(root.domain_size, 0), - orbit_has_f0_(root.domain_size, 0), - orbit_has_f1_(root.domain_size, 0), + orbit_has_b1_(root.num_original_vars, 0), + orbit_has_b0_(root.num_original_vars, 0), + orbit_has_f0_(root.num_original_vars, 0), + orbit_has_f1_(root.num_original_vars, 0), marked_b0_(root.num_original_vars, 0), marked_b1_(root.num_original_vars, 0), marked_f0_(root.num_original_vars, 0), @@ -75,20 +60,21 @@ class orbital_fixing_t { f0_.reserve(root.num_original_vars); f1_.reserve(root.num_original_vars); - // Clone the root Schreier structure into a per-worker copy. The root structure - // is shared and read-only, but set_base/get_stabilizer_orbit mutate internal state, - // so each worker needs its own copy. Since random_schreier is not copyable, we - // reconstruct it by creating a fresh structure with the same base and sifting all - // generators from the root into it. - schreier_ = std::make_unique(root.domain_size); - std::vector base(root.num_original_vars); + // Clone the projected Schreier (domain = num_original_vars, only integer-variable + // permutations). The mutex protects get_generator (which may not be thread-safe), + // while the expensive sift calls operate on the per-worker copy and run in parallel. + schreier_ = std::make_unique(num_original_vars_); + std::vector base(num_original_vars_); std::iota(base.begin(), base.end(), 0); schreier_->set_base(base); - int num_gens = root.schreier->get_number_of_generators(); - dejavu::groups::automorphism_workspace ws(root.domain_size); + int num_gens = root.num_generators; + dejavu::groups::automorphism_workspace ws(num_original_vars_); for (int i = 0; i < num_gens; i++) { - root.schreier->get_generator(i, ws); + { + std::lock_guard lock(root.generator_mutex); + root.schreier->get_generator(i, ws); + } schreier_->sift(ws); } } @@ -136,7 +122,7 @@ class orbital_fixing_t { } } - // Compute Stab(G', B1) where G' = Stab(G, general_integer_variables) + // Compute Stab(G', B1) where G' already excludes continuous-variable permutations. // Base order: [general_int_vars, B1, remaining_vars] const auto& gen_int = symmetry->general_integer_variables; std::vector new_base; @@ -162,7 +148,7 @@ class orbital_fixing_t { schreier_->set_base(new_base); dejavu::groups::orbit orb; - orb.initialize(domain_size_); + orb.initialize(num_original_vars_); i_t stabilizer_level = static_cast(gen_int.size() + branched_one.size()); schreier_->get_stabilizer_orbit(static_cast(stabilizer_level), orb); @@ -200,9 +186,8 @@ class orbital_fixing_t { orbit_has_b0_[orb.find_orbit(v)] = 1; } - for (i_t v : symmetry->continuous_variables) { - orbit_has_continuous_[orb.find_orbit(v)] = 1; - } + // No orbit_has_continuous_ needed: generators that permute continuous + // variables were excluded when building the projected schreier. for (i_t v : f0_) { orbit_has_f0_[orb.find_orbit(v)] = 1; @@ -219,8 +204,8 @@ class orbital_fixing_t { i_t o = orb.find_orbit(j); if (orb.orbit_size(o) < 2) continue; - if (orbit_has_b1_[o] == 1 || orbit_has_continuous_[o] == 1) { - // The orbit contains variables in B1 or continuous variables + if (orbit_has_b1_[o] == 1) { + // The orbit contains variables in B1 // So we can't fix any variables in this orbit continue; } @@ -249,10 +234,6 @@ class orbital_fixing_t { marked_b0_[v] = 0; } - for (i_t v : symmetry->continuous_variables) { - orbit_has_continuous_[orb.find_orbit(v)] = 0; - } - for (i_t v : f0_) { orbit_has_f0_[orb.find_orbit(v)] = 0; marked_f0_[v] = 0; @@ -288,22 +269,20 @@ class orbital_fixing_t { private: std::unique_ptr schreier_; - i_t domain_size_; i_t num_original_vars_; bool disabled_ = false; - std::vector marked_variables_; - std::vector orbit_has_b1_; // orbit_has_b1_[o] = 1 if the orbit o contains variables in B1 - std::vector orbit_has_b0_; // orbit_has_b0_[o] = 1 if the orbit o contains variables in B0 - std::vector orbit_has_continuous_; // orbit_has_continuous_[o] = 1 if the orbit o contains continuous variables - std::vector orbit_has_f0_; // orbit_has_f0_[o] = 1 if the orbit o contains variables in F0 - std::vector orbit_has_f1_; // orbit_has_f1_[o] = 1 if the orbit o contains variables in F1 - std::vector marked_b0_; // marked_b0_[v] = 1 if the variable v has been branched down - std::vector marked_b1_; // marked_b1_[v] = 1 if the variable v has been branched up - std::vector f0_; // set of variables that can be fixed to 0 - std::vector f1_; // set of variables that can be fixed to 1 - std::vector marked_f0_; // marked_f0_[v] = 1 if the variable v has been branched down - std::vector marked_f1_; // marked_f1_[v] = 1 if the variable v has been branched up + std::vector marked_variables_; // marked_variables_[j] = 1 if variable j is in the base prefix + std::vector orbit_has_b1_; // orbit_has_b1_[o] = 1 if orbit o contains variables in B1 + std::vector orbit_has_b0_; // orbit_has_b0_[o] = 1 if orbit o contains variables in B0 + std::vector orbit_has_f0_; // orbit_has_f0_[o] = 1 if orbit o contains variables in F0 + std::vector orbit_has_f1_; // orbit_has_f1_[o] = 1 if orbit o contains variables in F1 + std::vector marked_b0_; // marked_b0_[v] = 1 if variable v has been branched down + std::vector marked_b1_; // marked_b1_[v] = 1 if variable v has been branched up + std::vector f0_; // set of variables fixed to 0 by bound propagation + std::vector f1_; // set of variables fixed to 1 by bound propagation + std::vector marked_f0_; // marked_f0_[v] = 1 if variable v is in F0 + std::vector marked_f1_; // marked_f1_[v] = 1 if variable v is in F1 }; template @@ -658,6 +637,41 @@ std::unique_ptr> detect_symmetry( } settings.log.printf("\n"); + // Binary-only orbit statistics: count orbits that consist entirely of binary variables + i_t bin_nontrivial_orbits = 0; + i_t bin_max_orbit_size = 0; + i_t bin_total_vars = 0; + std::vector bin_orbit_histogram(num_original_vars + 1, 0); + i_t num_binary_vars = 0; + for (i_t j = 0; j < num_original_vars; j++) { + bool is_bin = (user_problem.var_types[j] != variable_type_t::CONTINUOUS && + user_problem.lower[j] == 0.0 && user_problem.upper[j] == 1.0); + if (is_bin && orb.represents_orbit(j)) { + i_t sz = orb.orbit_size(j); + if (sz >= 2) { + bin_nontrivial_orbits++; + bin_max_orbit_size = std::max(bin_max_orbit_size, sz); + bin_total_vars += sz; + bin_orbit_histogram[sz]++; + } + } + if (is_bin) { num_binary_vars++; } + } + if (num_binary_vars < num_original_vars) { + settings.log.printf("Binary orbits: %d non-trivial, max size %d, %d/%d (%.1f%%) binary variables in orbits\n", + bin_nontrivial_orbits, bin_max_orbit_size, bin_total_vars, num_binary_vars, + num_binary_vars > 0 ? 100.0 * bin_total_vars / num_binary_vars : 0.0); + if (bin_max_orbit_size >= 2) { + settings.log.printf("Binary orbit histogram (orbit size: number of orbits):"); + for (i_t sz = 2; sz <= bin_max_orbit_size; sz++) { + if (bin_orbit_histogram[sz] > 0) { + settings.log.printf(" %d:%d", sz, bin_orbit_histogram[sz]); + } + } + settings.log.printf("\n"); + } + } + has_symmetry = (max_orbit_size >= 4) || (num_nontrivial_orbits >= 3 && max_orbit_size >= 2) || (total_vars_in_orbits >= 10); @@ -667,9 +681,80 @@ std::unique_ptr> detect_symmetry( if (!has_symmetry) { return nullptr; } - return std::make_unique>( - std::move(schreier), num_vertices, user_problem.num_cols, var_types, - user_problem.lower, user_problem.upper); + auto result = std::make_unique>(); + result->num_original_vars = num_original_vars; + result->is_binary.resize(num_original_vars, 0); + for (i_t j = 0; j < num_original_vars; j++) { + if (var_types[j] != variable_type_t::CONTINUOUS) { + if (user_problem.lower[j] == 0.0 && user_problem.upper[j] == 1.0) { + result->is_binary[j] = 1; + } else { + result->general_integer_variables.push_back(j); + } + } + } + + // Project generators from the full graph (num_vertices) onto original variables + // (num_original_vars). Skip generators that permute any continuous variable. + result->schreier = std::make_unique(num_original_vars); + { + std::vector base(num_original_vars); + std::iota(base.begin(), base.end(), 0); + result->schreier->set_base(base); + } + + int full_num_gens = schreier->get_number_of_generators(); + dejavu::groups::automorphism_workspace full_ws(num_vertices); + dejavu::groups::automorphism_workspace proj_ws(num_original_vars); + int projected_count = 0; + int skipped_continuous = 0; + for (int i = 0; i < full_num_gens; i++) { + schreier->get_generator(i, full_ws); + const int* full_p = full_ws.p(); + + // Project onto original variables and check for continuous variable movement + proj_ws.reset(); + bool moves_continuous = false; + bool is_identity = true; + for (i_t j = 0; j < num_original_vars; j++) { + i_t mapped = full_p[j]; + if (mapped != j) { + if (var_types[j] == variable_type_t::CONTINUOUS) { + moves_continuous = true; + break; + } + proj_ws.write_single_map(j, mapped); + is_identity = false; + } + } + if (moves_continuous) { + skipped_continuous++; + continue; + } + if (!is_identity) { + result->schreier->sift(proj_ws); + projected_count++; + } + } + + result->num_generators = result->schreier->get_number_of_generators(); + settings.log.printf("Projected %d/%d generators onto %d variables (%d skipped continuous), " + "%d stored in projected schreier\n", + projected_count, full_num_gens, (int)num_original_vars, + skipped_continuous, result->num_generators); + + // Precompute orbit representatives from the projected group's orbits. + result->orbit_rep.resize(num_original_vars); + { + dejavu::groups::orbit orb_final; + orb_final.initialize(num_original_vars); + result->schreier->get_stabilizer_orbit(0, orb_final); + for (i_t j = 0; j < num_original_vars; j++) { + result->orbit_rep[j] = orb_final.find_orbit(j); + } + } + + return result; } } // namespace cuopt::linear_programming::dual_simplex From 7ea524fda6be4730bb6d679e8101da39970b54cb Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Fri, 17 Apr 2026 18:18:43 -0700 Subject: [PATCH 07/19] Remove dejavu use for everything but graph automorphism. Dont compute the whole stab(G, B1) but rather a subgroup. Speeds up time for each node. But misses fixings. --- cpp/src/branch_and_bound/branch_and_bound.cpp | 2 +- cpp/src/branch_and_bound/symmetry.hpp | 526 +++++++++++------- 2 files changed, 313 insertions(+), 215 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index f5892ca5f..e7f05d538 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1396,7 +1396,7 @@ dual::status_t branch_and_bound_t::solve_node_lp( // Perform orbital fixing auto* of = worker->orbital_fixing.get(); if (of != nullptr && !of->disabled()) { - of->orbital_fixing(symmetry_, settings_, var_types_, node_ptr, worker->leaf_problem); + of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem); } i_t node_iter = 0; diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index b08889c6b..95b237221 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -18,19 +18,183 @@ #include "dejavu.h" #include -#include #include #include namespace cuopt::linear_programming::dual_simplex { +// permutation_t stores a dense permutation plus its support (non-identity entries). +template +class permutation_t { + public: + // Takes vectors by value so callers can std::move or copy as needed. + permutation_t(std::vector p) : n_(p.size()), p_(std::move(p)) { + for (i_t k = 0; k < n_; k++) { + if (p_[k] != k) { support_.push_back(k); } + } + } + permutation_t(std::vector p, std::vector support) + : n_(p.size()), p_(std::move(p)), support_(std::move(support)) {} + + permutation_t(int n, const int* p, int nsupp, const int* supp) : n_(n) + { + p_.resize(n); + std::iota(p_.begin(), p_.end(), 0); + for (int k = 0; k < nsupp; k++) { + const int i = supp[k]; + p_[i] = p[i]; + support_.push_back(i); + } + } + + i_t size() const { return n_; } + const std::vector& dense_permutation() const { return p_; } + const std::vector& support() const { return support_; } + + private: + i_t n_; + std::vector p_; + std::vector support_; +}; + +// generators_t stores a list of permutations. Can be constructed from a sparse representation. +template +class generators_t { + public: + generators_t() : n_(-1) {} + void add_generator(int n, const int* p, int nsupp, const int* supp) + { + if (n_ == -1) { n_ = n; } + if (n != n_) { + return; + } else { + generators_.emplace_back(n, p, nsupp, supp); + } + } + void add_generator(const permutation_t& p) { + if (n_ == -1) { n_ = p.size(); } + if (p.size() != n_) { return; } + generators_.emplace_back(p); + } + void add_generator(std::vector p, std::vector support) { + const i_t n = static_cast(p.size()); + if (n_ == -1) { n_ = n; } + if (n != n_) { return; } + generators_.emplace_back(std::move(p), std::move(support)); + } + i_t size() const { return n_; } + size_t num_generators() const { return generators_.size(); } + + const permutation_t& get_generator(i_t i) const { return generators_[i]; } + private: + i_t n_; + std::vector> generators_; +}; + +// orbits_t computes the orbits of a set of generators using the union-find algorithm. +template +class orbits_t { + public: + orbits_t(i_t n) : n_(n), parent_(n), size_(n, 1), dirty_(n, 0) + { + std::iota(parent_.begin(), parent_.end(), 0); + dirty_list_.reserve(n); + } + + void compute_orbits(const std::vector& indices, const generators_t& generators) { + for (i_t i : indices) { + const permutation_t& perm = generators.get_generator(i); + const std::vector& p = perm.dense_permutation(); + for (i_t k : perm.support()) { + union_sets(k, p[k]); + } + } + } + + // Incrementally update orbits with a single mapping u -> v. + void add_mapping(i_t u, i_t v) { union_sets(u, v); } + + void compute_orbits(const generators_t& generators) { + std::vector indices(generators.num_generators()); + std::iota(indices.begin(), indices.end(), 0); + compute_orbits(indices, generators); + } + + i_t find_orbit(i_t v) { + return find(v); + } + + bool represents_orbit(i_t v) { + return find(v) == v; + } + + i_t orbit_size(i_t v) { + return size_[find(v)]; + } + + // Reset only the given indices back to identity (parent[j] = j, size[j] = 1). + void reset() { + for (i_t j : dirty_list_) { + parent_[j] = j; + size_[j] = 1; + dirty_[j] = 0; + } + dirty_list_.clear(); + } + +private: + void union_sets(i_t u, i_t v) { + i_t root_u = find(u); + i_t root_v = find(v); + if (root_u == root_v) return; // Already in the same set + mark_dirty(root_u); + mark_dirty(root_v); + if (size_[root_u] < size_[root_v]) { + parent_[root_u] = root_v; + size_[root_v] += size_[root_u]; + } else { + parent_[root_v] = root_u; + size_[root_u] += size_[root_v]; + } + } + + i_t find(i_t v) + { + i_t root = v; + while (parent_[root] != root) { + root = parent_[root]; + } + // Path compression + while (parent_[v] != root) { + i_t next = parent_[v]; + parent_[v] = root; + mark_dirty(v); + v = next; + } + return root; + } + + void mark_dirty(i_t v) { + if (dirty_[v] == 0) { + dirty_list_.push_back(v); + dirty_[v] = 1; + } + } + + + i_t n_; + std::vector parent_; + std::vector size_; + std::vector dirty_; + std::vector dirty_list_; +}; + template struct mip_symmetry_t { - // Schreier projected onto original variables (domain = num_original_vars). - // Generators that permute continuous variables are excluded. - std::unique_ptr schreier; + generators_t generators; i_t num_original_vars; int num_generators = 0; + std::vector binary_variables; std::vector general_integer_variables; std::vector is_binary; @@ -38,8 +202,6 @@ struct mip_symmetry_t { // orbit_rep[j] = orbit representative of variable j (for j < num_original_vars). std::vector orbit_rep; - // Mutex protecting schreier->get_generator() for concurrent orbital_fixing_t construction. - std::mutex generator_mutex; }; template @@ -47,7 +209,7 @@ class orbital_fixing_t { public: explicit orbital_fixing_t(mip_symmetry_t& root) : num_original_vars_(root.num_original_vars), - marked_variables_(root.num_original_vars, 0), + orb_(root.num_original_vars), orbit_has_b1_(root.num_original_vars, 0), orbit_has_b0_(root.num_original_vars, 0), orbit_has_f0_(root.num_original_vars, 0), @@ -59,24 +221,6 @@ class orbital_fixing_t { { f0_.reserve(root.num_original_vars); f1_.reserve(root.num_original_vars); - - // Clone the projected Schreier (domain = num_original_vars, only integer-variable - // permutations). The mutex protects get_generator (which may not be thread-safe), - // while the expensive sift calls operate on the per-worker copy and run in parallel. - schreier_ = std::make_unique(num_original_vars_); - std::vector base(num_original_vars_); - std::iota(base.begin(), base.end(), 0); - schreier_->set_base(base); - - int num_gens = root.num_generators; - dejavu::groups::automorphism_workspace ws(num_original_vars_); - for (int i = 0; i < num_gens; i++) { - { - std::lock_guard lock(root.generator_mutex); - root.schreier->get_generator(i, ws); - } - schreier_->sift(ws); - } } bool disabled () const { return disabled_; } @@ -85,7 +229,6 @@ class orbital_fixing_t { void orbital_fixing(mip_symmetry_t* symmetry, const simplex_solver_settings_t& settings, - const std::vector& var_types, mip_node_t* node_ptr, lp_problem_t& problem) { @@ -110,8 +253,7 @@ class orbital_fixing_t { node = node->parent; } - for (i_t j = 0; j < num_original_vars_; j++) { - if (var_types[j] == variable_type_t::CONTINUOUS) continue; + for (i_t j : symmetry->binary_variables) { if (marked_b1_[j] == 0 && problem.lower[j] == 1.0) { f1_.push_back(j); marked_f1_[j] = 1; @@ -122,40 +264,34 @@ class orbital_fixing_t { } } - // Compute Stab(G', B1) where G' already excludes continuous-variable permutations. - // Base order: [general_int_vars, B1, remaining_vars] - const auto& gen_int = symmetry->general_integer_variables; - std::vector new_base; - new_base.reserve(num_original_vars_); - for (i_t j : gen_int) { - new_base.push_back(j); - marked_variables_[j] = 1; - } - for (i_t j : branched_one) { - new_base.push_back(j); - marked_variables_[j] = 1; - } - for (i_t j = 0; j < num_original_vars_; j++) { - if (marked_variables_[j] == 0) { new_base.push_back(j); } - } - for (i_t j : gen_int) { - marked_variables_[j] = 0; - } - for (i_t j : branched_one) { - marked_variables_[j] = 0; - } + // In true orbital fixing we would compute the group G' = stabilizer(G, B1) + // Instead we compute a subgroup H of G' that is just + // H = { g in generators(G) | g(B1) = B1 } + // This means that we will miss some fixings. But every fixing we perform will be valid. + std::vector surviving_generators; + surviving_generators.reserve(symmetry->generators.num_generators()); + for (size_t k = 0; k < symmetry->generators.num_generators(); k++) { + const permutation_t& perm = symmetry->generators.get_generator(k); + const std::vector& p = perm.dense_permutation(); + bool stabilizes_b1 = true; + + for (i_t j : branched_one) { + const i_t mapped = p[j]; + if (marked_b1_[mapped] == 0) { + stabilizes_b1 = false; + break; + } + } - schreier_->set_base(new_base); + if (stabilizes_b1) { surviving_generators.push_back(static_cast(k)); } + } - dejavu::groups::orbit orb; - orb.initialize(num_original_vars_); - i_t stabilizer_level = static_cast(gen_int.size() + branched_one.size()); - schreier_->get_stabilizer_orbit(static_cast(stabilizer_level), orb); + orb_.compute_orbits(surviving_generators, symmetry->generators); - // Check if the stabilizer is trivial (all orbits are singletons) + // Check if the stabilizer is trivial (all binary orbits are singletons) bool trivial_stabilizer = true; - for (i_t j = 0; j < num_original_vars_; j++) { - if (orb.orbit_size(orb.find_orbit(j)) > 1) { + for (i_t j : symmetry->binary_variables) { + if (orb_.orbit_size(orb_.find_orbit(j)) > 1) { trivial_stabilizer = false; break; } @@ -179,30 +315,27 @@ class orbital_fixing_t { f1_.clear(); } else { for (i_t v : branched_one) { - orbit_has_b1_[orb.find_orbit(v)] = 1; + orbit_has_b1_[orb_.find_orbit(v)] = 1; } for (i_t v : branched_zero) { - orbit_has_b0_[orb.find_orbit(v)] = 1; + orbit_has_b0_[orb_.find_orbit(v)] = 1; } - // No orbit_has_continuous_ needed: generators that permute continuous - // variables were excluded when building the projected schreier. - for (i_t v : f0_) { - orbit_has_f0_[orb.find_orbit(v)] = 1; + orbit_has_f0_[orb_.find_orbit(v)] = 1; } for (i_t v : f1_) { - orbit_has_f1_[orb.find_orbit(v)] = 1; + orbit_has_f1_[orb_.find_orbit(v)] = 1; } std::vector fix_zero; // The set L0 of variables that can be fixed to 0 std::vector fix_one; // The set L1 of variables that can be fixed to 1 - for (i_t j = 0; j < num_original_vars_; j++) { - i_t o = orb.find_orbit(j); - if (orb.orbit_size(o) < 2) continue; + for (i_t j : symmetry->binary_variables) { + i_t o = orb_.find_orbit(j); + if (orb_.orbit_size(o) < 2) continue; if (orbit_has_b1_[o] == 1) { // The orbit contains variables in B1 @@ -225,22 +358,22 @@ class orbital_fixing_t { // Restore the work arrays for (i_t v : branched_one) { - orbit_has_b1_[orb.find_orbit(v)] = 0; + orbit_has_b1_[orb_.find_orbit(v)] = 0; marked_b1_[v] = 0; } for (i_t v : branched_zero) { - orbit_has_b0_[orb.find_orbit(v)] = 0; + orbit_has_b0_[orb_.find_orbit(v)] = 0; marked_b0_[v] = 0; } for (i_t v : f0_) { - orbit_has_f0_[orb.find_orbit(v)] = 0; + orbit_has_f0_[orb_.find_orbit(v)] = 0; marked_f0_[v] = 0; } for (i_t v : f1_) { - orbit_has_f1_[orb.find_orbit(v)] = 0; + orbit_has_f1_[orb_.find_orbit(v)] = 0; marked_f1_[v] = 0; } @@ -265,14 +398,16 @@ class orbital_fixing_t { problem.upper[v] = 1.0; } } + + // Reset orbits for reuse + orb_.reset(); } private: - std::unique_ptr schreier_; i_t num_original_vars_; bool disabled_ = false; + orbits_t orb_; - std::vector marked_variables_; // marked_variables_[j] = 1 if variable j is in the base prefix std::vector orbit_has_b1_; // orbit_has_b1_[o] = 1 if orbit o contains variables in B1 std::vector orbit_has_b0_; // orbit_has_b0_[o] = 1 if orbit o contains variables in B0 std::vector orbit_has_f0_; // orbit_has_f0_[o] = 1 if orbit o contains variables in F0 @@ -571,30 +706,85 @@ std::unique_ptr> detect_symmetry( g.add_edge(std::min(u, v), std::max(u, v)); } - auto schreier = std::make_unique(num_vertices); - std::vector base; - base.reserve(user_problem.num_cols); - for (i_t j = 0; j < user_problem.num_cols; j++) { - base.push_back(j); - } - schreier->set_base(base); + const i_t num_original_vars = user_problem.num_cols; - dejavu::hooks::schreier_hook s_hook(*schreier); + auto result = std::make_unique>(); + result->num_original_vars = num_original_vars; + result->is_binary.resize(num_original_vars, 0); + for (i_t j = 0; j < num_original_vars; j++) { + if (var_types[j] != variable_type_t::CONTINUOUS) { + if (user_problem.lower[j] == 0.0 && user_problem.upper[j] == 1.0) { + result->is_binary[j] = 1; + result->binary_variables.push_back(j); + } else { + result->general_integer_variables.push_back(j); + } + } + } - i_t num_generators = 0; - const i_t num_original_vars = user_problem.num_cols; - dejavu_hook counting_hook = [&num_generators, num_original_vars](int, const int*, int nsupp, const int* supp) { + // Project generators incrementally inside the dejavu callback. + // This avoids storing full-graph dense vectors (size num_vertices per generator). + const size_t max_generators = std::max(size_t{1}, static_cast(64000000 / num_original_vars)); + orbits_t orb(num_original_vars); + int num_dejavu_generators = 0; + int projected_count = 0; + int skipped_non_binary = 0; + std::vector projected_p(num_original_vars); + std::iota(projected_p.begin(), projected_p.end(), 0); + std::vector var_support; + + dejavu_hook generator_hook = [&num_original_vars, &num_dejavu_generators, + &projected_p, &var_support, &orb, &result, &projected_count, + &skipped_non_binary, &max_generators](int n, const int* p, + int nsupp, const int* supp) { + // Check if any support element is an original variable + bool moves_variable = false; for (int s = 0; s < nsupp; s++) { if (supp[s] < num_original_vars) { - num_generators++; - return; + moves_variable = true; + break; } } + if (!moves_variable) { return; } + num_dejavu_generators++; + + // Project onto binary variables only. + // Skip generators that move any non-binary variable (continuous or general integer). + // Only touch support entries; reset them back to identity afterward. + // dejavu guarantees supp is exactly the non-identity entries of p. + bool moves_non_binary = false; + var_support.clear(); + for (int s = 0; s < nsupp; s++) { + const i_t j = supp[s]; + if (j >= num_original_vars) { continue; } + if (result->is_binary[j] == 0) { + moves_non_binary = true; + break; + } + projected_p[j] = p[j]; + var_support.push_back(j); + } + + if (!moves_non_binary && !var_support.empty()) { + for (i_t j : var_support) { + orb.add_mapping(j, projected_p[j]); + } + if (result->generators.num_generators() < max_generators) { + result->generators.add_generator(projected_p, var_support); + } + projected_count++; + } else if (moves_non_binary) { + skipped_non_binary++; + } + + // Reset modified entries back to identity + for (i_t j : var_support) { + projected_p[j] = j; + } }; dejavu::hooks::multi_hook combined; - combined.add_hook(s_hook.get_hook()); - combined.add_hook(&counting_hook); + combined.add_hook(&generator_hook); dejavu::solver d; d.automorphisms(&g, combined); @@ -602,34 +792,44 @@ std::unique_ptr> detect_symmetry( std::ostringstream grp_size_str; grp_size_str << d.get_automorphism_group_size(); settings.log.printf( - "Automorphism group size %s, %d generators\n", grp_size_str.str().c_str(), num_generators); + "Automorphism group size %s, %d dejavu generators (%d move variables)\n", + grp_size_str.str().c_str(), num_dejavu_generators, projected_count); settings.log.printf("Dejavu time %f\n", toc(dejavu_start_time)); + result->num_generators = result->generators.num_generators(); + if (projected_count > static_cast(result->num_generators)) { + settings.log.printf("Generator limit: kept %d/%d projected generators (limit %d)\n", + result->num_generators, projected_count, (int)max_generators); + } + settings.log.printf("Projected %d generators onto %d binary variables (%d skipped non-binary), " + "%d stored\n", + projected_count, (int)num_original_vars, + skipped_non_binary, result->num_generators); + + // Compute orbit statistics from the incrementally built orbits. + // All non-trivial orbits contain only binary variables (non-binary generators were excluded). has_symmetry = false; - if (num_generators > 0) { - dejavu::groups::orbit orb; - orb.initialize(num_vertices); - schreier->get_stabilizer_orbit(0, orb); - - i_t num_nontrivial_orbits = 0; - i_t max_orbit_size = 0; - i_t total_vars_in_orbits = 0; - std::vector orbit_histogram(num_original_vars + 1, 0); - for (i_t j = 0; j < num_original_vars; j++) { - if (orb.represents_orbit(j)) { - i_t sz = orb.orbit_size(j); - if (sz >= 2) { - num_nontrivial_orbits++; - max_orbit_size = std::max(max_orbit_size, sz); - total_vars_in_orbits += sz; - orbit_histogram[sz]++; - } + i_t num_nontrivial_orbits = 0; + i_t max_orbit_size = 0; + i_t total_vars_in_orbits = 0; + std::vector orbit_histogram(num_original_vars + 1, 0); + for (i_t j : result->binary_variables) { + if (orb.represents_orbit(j)) { + i_t sz = orb.orbit_size(j); + if (sz >= 2) { + num_nontrivial_orbits++; + max_orbit_size = std::max(max_orbit_size, sz); + total_vars_in_orbits += sz; + orbit_histogram[sz]++; } } - settings.log.printf("Orbits: %d non-trivial, max size %d, %d/%d (%.1f%%) variables in orbits\n", + } + + if (projected_count > 0) { + settings.log.printf("Binary orbits: %d non-trivial, max size %d, %d/%d (%.1f%%) binary variables in orbits\n", num_nontrivial_orbits, max_orbit_size, total_vars_in_orbits, num_original_vars, 100.0 * total_vars_in_orbits / num_original_vars); - settings.log.printf("Orbit histogram (orbit size: number of orbits):"); + settings.log.printf("Orbit histogram (size: count):"); for (i_t sz = 2; sz <= max_orbit_size; sz++) { if (orbit_histogram[sz] > 0) { settings.log.printf(" %d:%d", sz, orbit_histogram[sz]); @@ -637,41 +837,6 @@ std::unique_ptr> detect_symmetry( } settings.log.printf("\n"); - // Binary-only orbit statistics: count orbits that consist entirely of binary variables - i_t bin_nontrivial_orbits = 0; - i_t bin_max_orbit_size = 0; - i_t bin_total_vars = 0; - std::vector bin_orbit_histogram(num_original_vars + 1, 0); - i_t num_binary_vars = 0; - for (i_t j = 0; j < num_original_vars; j++) { - bool is_bin = (user_problem.var_types[j] != variable_type_t::CONTINUOUS && - user_problem.lower[j] == 0.0 && user_problem.upper[j] == 1.0); - if (is_bin && orb.represents_orbit(j)) { - i_t sz = orb.orbit_size(j); - if (sz >= 2) { - bin_nontrivial_orbits++; - bin_max_orbit_size = std::max(bin_max_orbit_size, sz); - bin_total_vars += sz; - bin_orbit_histogram[sz]++; - } - } - if (is_bin) { num_binary_vars++; } - } - if (num_binary_vars < num_original_vars) { - settings.log.printf("Binary orbits: %d non-trivial, max size %d, %d/%d (%.1f%%) binary variables in orbits\n", - bin_nontrivial_orbits, bin_max_orbit_size, bin_total_vars, num_binary_vars, - num_binary_vars > 0 ? 100.0 * bin_total_vars / num_binary_vars : 0.0); - if (bin_max_orbit_size >= 2) { - settings.log.printf("Binary orbit histogram (orbit size: number of orbits):"); - for (i_t sz = 2; sz <= bin_max_orbit_size; sz++) { - if (bin_orbit_histogram[sz] > 0) { - settings.log.printf(" %d:%d", sz, bin_orbit_histogram[sz]); - } - } - settings.log.printf("\n"); - } - } - has_symmetry = (max_orbit_size >= 4) || (num_nontrivial_orbits >= 3 && max_orbit_size >= 2) || (total_vars_in_orbits >= 10); @@ -681,77 +846,10 @@ std::unique_ptr> detect_symmetry( if (!has_symmetry) { return nullptr; } - auto result = std::make_unique>(); - result->num_original_vars = num_original_vars; - result->is_binary.resize(num_original_vars, 0); - for (i_t j = 0; j < num_original_vars; j++) { - if (var_types[j] != variable_type_t::CONTINUOUS) { - if (user_problem.lower[j] == 0.0 && user_problem.upper[j] == 1.0) { - result->is_binary[j] = 1; - } else { - result->general_integer_variables.push_back(j); - } - } - } - - // Project generators from the full graph (num_vertices) onto original variables - // (num_original_vars). Skip generators that permute any continuous variable. - result->schreier = std::make_unique(num_original_vars); - { - std::vector base(num_original_vars); - std::iota(base.begin(), base.end(), 0); - result->schreier->set_base(base); - } - - int full_num_gens = schreier->get_number_of_generators(); - dejavu::groups::automorphism_workspace full_ws(num_vertices); - dejavu::groups::automorphism_workspace proj_ws(num_original_vars); - int projected_count = 0; - int skipped_continuous = 0; - for (int i = 0; i < full_num_gens; i++) { - schreier->get_generator(i, full_ws); - const int* full_p = full_ws.p(); - - // Project onto original variables and check for continuous variable movement - proj_ws.reset(); - bool moves_continuous = false; - bool is_identity = true; - for (i_t j = 0; j < num_original_vars; j++) { - i_t mapped = full_p[j]; - if (mapped != j) { - if (var_types[j] == variable_type_t::CONTINUOUS) { - moves_continuous = true; - break; - } - proj_ws.write_single_map(j, mapped); - is_identity = false; - } - } - if (moves_continuous) { - skipped_continuous++; - continue; - } - if (!is_identity) { - result->schreier->sift(proj_ws); - projected_count++; - } - } - - result->num_generators = result->schreier->get_number_of_generators(); - settings.log.printf("Projected %d/%d generators onto %d variables (%d skipped continuous), " - "%d stored in projected schreier\n", - projected_count, full_num_gens, (int)num_original_vars, - skipped_continuous, result->num_generators); - // Precompute orbit representatives from the projected group's orbits. result->orbit_rep.resize(num_original_vars); - { - dejavu::groups::orbit orb_final; - orb_final.initialize(num_original_vars); - result->schreier->get_stabilizer_orbit(0, orb_final); - for (i_t j = 0; j < num_original_vars; j++) { - result->orbit_rep[j] = orb_final.find_orbit(j); - } + for (i_t j = 0; j < num_original_vars; j++) { + result->orbit_rep[j] = orb.find_orbit(j); } return result; From 426582a0947df4ad8dc0e8a8095f23a0bc7babf0 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 20 Apr 2026 13:50:23 -0700 Subject: [PATCH 08/19] Use the fact that we are performing a plunge to avoid recomputing the full stabilizer(G, B1) --- cpp/src/branch_and_bound/branch_and_bound.cpp | 4 +- cpp/src/branch_and_bound/symmetry.hpp | 301 ++++++++++-------- 2 files changed, 176 insertions(+), 129 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index e7f05d538..5d964277f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1390,6 +1390,9 @@ dual::status_t branch_and_bound_t::solve_node_lp( bool feasible = worker->set_lp_variable_bounds(node_ptr, settings_); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; worker->leaf_edge_norms = edge_norms_; + if (worker->recompute_bounds && worker->orbital_fixing) { + worker->orbital_fixing->reset(symmetry_, node_ptr); + } if (feasible) { @@ -1449,7 +1452,6 @@ void branch_and_bound_t::plunge_with(branch_and_bound_worker_trecompute_basis = true; worker->recompute_bounds = true; worker->ensure_orbital_fixing(); - if (worker->orbital_fixing) { worker->orbital_fixing->enable(); } f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 95b237221..c8f80fa67 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -209,6 +209,8 @@ class orbital_fixing_t { public: explicit orbital_fixing_t(mip_symmetry_t& root) : num_original_vars_(root.num_original_vars), + max_generators_(root.num_generators), + start_plunge_(true), orb_(root.num_original_vars), orbit_has_b1_(root.num_original_vars, 0), orbit_has_b0_(root.num_original_vars, 0), @@ -219,13 +221,54 @@ class orbital_fixing_t { marked_f0_(root.num_original_vars, 0), marked_f1_(root.num_original_vars, 0) { + branched_zero_.reserve(root.num_original_vars); + branched_one_.reserve(root.num_original_vars); f0_.reserve(root.num_original_vars); f1_.reserve(root.num_original_vars); + fix_zero_.reserve(root.num_original_vars); + fix_one_.reserve(root.num_original_vars); + + surviving_generators_.resize(max_generators_); + std::iota(surviving_generators_.begin(), surviving_generators_.end(), 0); } - bool disabled () const { return disabled_; } - void disable() { disabled_ = true; } - void enable() { disabled_ = false; } + bool disabled () const { return surviving_generators_.empty(); } + void disable() { surviving_generators_.clear(); } + + void reset(mip_symmetry_t* symmetry, mip_node_t* node_ptr) { + for (i_t j : branched_zero_) { + marked_b0_[j] = 0; + orbit_has_b0_[orb_.find_orbit(j)] = 0; + } + for (i_t j : branched_one_) { + marked_b1_[j] = 0; + orbit_has_b1_[orb_.find_orbit(j)] = 0; + } + branched_zero_.clear(); + branched_one_.clear(); + orb_.reset(); + + mip_node_t* node = node_ptr; + while (node != nullptr && node->branch_var >= 0) { + i_t v = node->branch_var; + bool is_binary = (symmetry->is_binary[v] == 1); + if (is_binary) { + if (node->branch_var_upper == 0.0) { + branched_zero_.push_back(v); + marked_b0_[v] = 1; + } else if (node->branch_var_lower == 1.0) { + branched_one_.push_back(v); + marked_b1_[v] = 1; + } + } + node = node->parent; + } + + surviving_generators_.resize(max_generators_); + std::iota(surviving_generators_.begin(), surviving_generators_.end(), 0); + + start_plunge_ = true; + } void orbital_fixing(mip_symmetry_t* symmetry, const simplex_solver_settings_t& settings, @@ -233,24 +276,18 @@ class orbital_fixing_t { lp_problem_t& problem) { // Collect binary branchings only; skip general integer branchings - std::vector branched_zero; - std::vector branched_one; - branched_zero.reserve(node_ptr->depth); - branched_one.reserve(node_ptr->depth); - mip_node_t* node = node_ptr; - while (node != nullptr && node->branch_var >= 0) { - i_t v = node->branch_var; - bool is_binary = (symmetry->is_binary[v] == 1); - if (is_binary) { - if (node->branch_var_upper == 0.0) { - branched_zero.push_back(v); - marked_b0_[v] = 1; - } else if (node->branch_var_lower == 1.0) { - branched_one.push_back(v); - marked_b1_[v] = 1; - } + i_t v = node_ptr->branch_var; + bool is_binary = (symmetry->is_binary[v] == 1); + bool should_recompute = start_plunge_; + if (is_binary) { + if (node_ptr->branch_var_upper == 0.0) { + branched_zero_.push_back(v); + marked_b0_[v] = 1; + } else if (node_ptr->branch_var_lower == 1.0) { + branched_one_.push_back(v); + marked_b1_[v] = 1; + should_recompute = true; } - node = node->parent; } for (i_t j : symmetry->binary_variables) { @@ -268,146 +305,152 @@ class orbital_fixing_t { // Instead we compute a subgroup H of G' that is just // H = { g in generators(G) | g(B1) = B1 } // This means that we will miss some fixings. But every fixing we perform will be valid. - std::vector surviving_generators; - surviving_generators.reserve(symmetry->generators.num_generators()); - for (size_t k = 0; k < symmetry->generators.num_generators(); k++) { - const permutation_t& perm = symmetry->generators.get_generator(k); - const std::vector& p = perm.dense_permutation(); - bool stabilizes_b1 = true; - - for (i_t j : branched_one) { - const i_t mapped = p[j]; - if (marked_b1_[mapped] == 0) { - stabilizes_b1 = false; - break; + // + // We are called within the context of a plunge within the branch and bound tree + // Suppose we computed H for the parent node. + // If we branched on a binary variable x_j then we are in one of two cases: + // Case 1: x_j = 0, then B1 remains unchanged and we do not need to recompute H or its orbits + // Case 2: x_j = 1, then B1 becomes B1 union {x_j} and we only need to check the generators + // used to construct H for the parent node. We call these the "surviving generators". + // + // Note that when we restart the plunge we need to recompute everything. + if (should_recompute) { + for (size_t k = 0; k < surviving_generators_.size(); k++) { + const i_t h = surviving_generators_[k]; + const permutation_t& perm = symmetry->generators.get_generator(h); + const std::vector& p = perm.dense_permutation(); + bool stabilizes_b1 = true; + + if (start_plunge_) { + // At the start of a plunge, surviving_generators_ was reset to all generators. + // We need to check each generator against all of B1. + for (i_t j : branched_one_) { + if (marked_b1_[p[j]] == 0) { + stabilizes_b1 = false; + break; + } + } + } else { + // B1 = v union B1_from_parent. So the surviving generators already stabilize B1_from_parent. + // We only need to check if the generators stabilize v. + if (marked_b1_[p[v]] == 0) { + stabilizes_b1 = false; + } } - } - - if (stabilizes_b1) { surviving_generators.push_back(static_cast(k)); } - } - - orb_.compute_orbits(surviving_generators, symmetry->generators); - // Check if the stabilizer is trivial (all binary orbits are singletons) - bool trivial_stabilizer = true; - for (i_t j : symmetry->binary_variables) { - if (orb_.orbit_size(orb_.find_orbit(j)) > 1) { - trivial_stabilizer = false; - break; + if (!stabilizes_b1) { + std::swap(surviving_generators_[k], surviving_generators_.back()); + surviving_generators_.pop_back(); + k--; + } } - } - if (trivial_stabilizer) { - disabled_ = true; - for (i_t v : branched_one) { - marked_b1_[v] = 0; - } - for (i_t v : branched_zero) { - marked_b0_[v] = 0; - } - for (i_t v : f0_) { - marked_f0_[v] = 0; - } - for (i_t v : f1_) { - marked_f1_[v] = 0; + // Clear old orbit_has values before orbits change + for (i_t v : branched_one_) { + orbit_has_b1_[orb_.find_orbit(v)] = 0; } - f0_.clear(); - f1_.clear(); - } else { - for (i_t v : branched_one) { - orbit_has_b1_[orb_.find_orbit(v)] = 1; + for (i_t v : branched_zero_) { + orbit_has_b0_[orb_.find_orbit(v)] = 0; } - for (i_t v : branched_zero) { - orbit_has_b0_[orb_.find_orbit(v)] = 1; - } + orb_.reset(); + orb_.compute_orbits(surviving_generators_, symmetry->generators); - for (i_t v : f0_) { - orbit_has_f0_[orb_.find_orbit(v)] = 1; + for (i_t v : branched_one_) { + orbit_has_b1_[orb_.find_orbit(v)] = 1; } - - for (i_t v : f1_) { - orbit_has_f1_[orb_.find_orbit(v)] = 1; + for (i_t v : branched_zero_) { + orbit_has_b0_[orb_.find_orbit(v)] = 1; } + } else if (is_binary && node_ptr->branch_var_upper == 0.0) { + orbit_has_b0_[orb_.find_orbit(v)] = 1; + } + start_plunge_ = false; - std::vector fix_zero; // The set L0 of variables that can be fixed to 0 - std::vector fix_one; // The set L1 of variables that can be fixed to 1 - - for (i_t j : symmetry->binary_variables) { - i_t o = orb_.find_orbit(j); - if (orb_.orbit_size(o) < 2) continue; + for (i_t v : f0_) { + orbit_has_f0_[orb_.find_orbit(v)] = 1; + } - if (orbit_has_b1_[o] == 1) { - // The orbit contains variables in B1 - // So we can't fix any variables in this orbit - continue; - } + for (i_t v : f1_) { + orbit_has_f1_[orb_.find_orbit(v)] = 1; + } - if (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1) { - // The orbit of this variable contains variables in B0 or F0 - // So we can fix this variable to zero (provided its not already in B0 or F0) - if (marked_b0_[j] == 0 && marked_f0_[j] == 0) { fix_zero.push_back(j); } - } + fix_zero_.clear(); + fix_one_.clear(); + bool stabilizer_nontrivial = false; + for (i_t j : symmetry->binary_variables) { + i_t o = orb_.find_orbit(j); + if (orb_.orbit_size(o) < 2) continue; - if (orbit_has_f1_[o] == 1) { - // The orbit of this variable contains variables in F1 - // So we can fix this variable to one (provided its not already in F1) - if (marked_f1_[j] == 0) { fix_one.push_back(j); } - } + if (orbit_has_b1_[o] == 1) { + // The orbit contains variables in B1 + // So we can't fix any variables in this orbit + continue; } - // Restore the work arrays - for (i_t v : branched_one) { - orbit_has_b1_[orb_.find_orbit(v)] = 0; - marked_b1_[v] = 0; - } + // A non-B1 orbit of size >= 2 exists, so the stabilizer may produce fixings + // in this node or a descendant. + stabilizer_nontrivial = true; - for (i_t v : branched_zero) { - orbit_has_b0_[orb_.find_orbit(v)] = 0; - marked_b0_[v] = 0; + if (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1) { + // The orbit of this variable contains variables in B0 or F0 + // So we can fix this variable to zero (provided its not already in B0 or F0) + if (marked_b0_[j] == 0 && marked_f0_[j] == 0) { fix_zero_.push_back(j); } } - for (i_t v : f0_) { - orbit_has_f0_[orb_.find_orbit(v)] = 0; - marked_f0_[v] = 0; + if (orbit_has_f1_[o] == 1) { + // The orbit of this variable contains variables in F1 + // So we can fix this variable to one (provided its not already in F1) + if (marked_f1_[j] == 0) { fix_one_.push_back(j); } } + } - for (i_t v : f1_) { - orbit_has_f1_[orb_.find_orbit(v)] = 0; - marked_f1_[v] = 0; - } + if (!stabilizer_nontrivial) { surviving_generators_.clear(); } - f0_.clear(); - f1_.clear(); + // Restore the work arrays + for (i_t v : f0_) { + orbit_has_f0_[orb_.find_orbit(v)] = 0; + marked_f0_[v] = 0; + } - if (fix_zero.size() > 0 || fix_one.size() > 0) { - settings.log.printf("Orbital fixing at node %d (depth %d): fixing %d to 0 and %d to 1\n", - node_ptr->node_id, - node_ptr->depth, - (int)fix_zero.size(), - (int)fix_one.size()); - } + for (i_t v : f1_) { + orbit_has_f1_[orb_.find_orbit(v)] = 0; + marked_f1_[v] = 0; + } - // Finally fix the variables in L0 and L1 - for (i_t v : fix_zero) { - problem.lower[v] = 0.0; - problem.upper[v] = 0.0; - } - for (i_t v : fix_one) { - problem.lower[v] = 1.0; - problem.upper[v] = 1.0; - } + f0_.clear(); + f1_.clear(); + + if (fix_zero_.size() > 0 || fix_one_.size() > 0) { + settings.log.printf("Orbital fixing at node %d (depth %d): fixing %d to 0 and %d to 1\n", + node_ptr->node_id, + node_ptr->depth, + (int)fix_zero_.size(), + (int)fix_one_.size()); } - // Reset orbits for reuse - orb_.reset(); + // Finally fix the variables in L0 and L1 + for (i_t v : fix_zero_) { + problem.lower[v] = 0.0; + problem.upper[v] = 0.0; + } + for (i_t v : fix_one_) { + problem.lower[v] = 1.0; + problem.upper[v] = 1.0; + } } private: i_t num_original_vars_; - bool disabled_ = false; + i_t max_generators_; + bool start_plunge_; orbits_t orb_; + std::vector branched_zero_; + std::vector branched_one_; + + std::vector surviving_generators_; + std::vector orbit_has_b1_; // orbit_has_b1_[o] = 1 if orbit o contains variables in B1 std::vector orbit_has_b0_; // orbit_has_b0_[o] = 1 if orbit o contains variables in B0 std::vector orbit_has_f0_; // orbit_has_f0_[o] = 1 if orbit o contains variables in F0 @@ -418,6 +461,8 @@ class orbital_fixing_t { std::vector f1_; // set of variables fixed to 1 by bound propagation std::vector marked_f0_; // marked_f0_[v] = 1 if variable v is in F0 std::vector marked_f1_; // marked_f1_[v] = 1 if variable v is in F1 + std::vector fix_zero_; // set of variables fixed to 0 by orbital fixing + std::vector fix_one_; // set of variables fixed to 1 by orbital fixing }; template From 07d421be9ea66666e35191fd956c4622387888e8 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 20 Apr 2026 17:40:16 -0700 Subject: [PATCH 09/19] Remove generators that don't perserve bounds after cuts and root bound strengthening. Fixes an issue on lectsched-5-obj --- cpp/src/branch_and_bound/branch_and_bound.cpp | 11 +++ cpp/src/branch_and_bound/symmetry.hpp | 80 ++++++++++++++++--- 2 files changed, 78 insertions(+), 13 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 5d964277f..bf93ddabf 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2609,6 +2609,17 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut node_queue_.push(search_tree_.root.get_down_child()); node_queue_.push(search_tree_.root.get_up_child()); + if (symmetry_ != nullptr) { + i_t removed = symmetry_->generators.template prune_by_bounds( + original_lp_.lower, original_lp_.upper); + if (removed > 0) { + symmetry_->num_generators = static_cast(symmetry_->generators.num_generators()); + settings_.log.printf( + "Pruned %d generators invalidated by root-level bound tightening, %d remain\n", + removed, symmetry_->num_generators); + } + } + settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads); node_concurrent_halt_ = 0; diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index c8f80fa67..1490b2807 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -86,6 +86,37 @@ class generators_t { size_t num_generators() const { return generators_.size(); } const permutation_t& get_generator(i_t i) const { return generators_[i]; } + + // Remove generators that don't preserve the given bounds. + // A generator is invalid if it maps variable j to p[j] where + // lower[j] != lower[p[j]] or upper[j] != upper[p[j]]. + template + i_t prune_by_bounds(const std::vector& lower, const std::vector& upper) + { + i_t num_removed = 0; + size_t k = 0; + while (k < generators_.size()) { + const auto& perm = generators_[k]; + const auto& p = perm.dense_permutation(); + bool valid = true; + for (i_t j : perm.support()) { + i_t pj = p[j]; + if (lower[j] != lower[pj] || upper[j] != upper[pj]) { + valid = false; + break; + } + } + if (!valid) { + std::swap(generators_[k], generators_.back()); + generators_.pop_back(); + num_removed++; + } else { + k++; + } + } + return num_removed; + } + private: i_t n_; std::vector> generators_; @@ -248,7 +279,7 @@ class orbital_fixing_t { branched_one_.clear(); orb_.reset(); - mip_node_t* node = node_ptr; + mip_node_t* node = node_ptr->parent; while (node != nullptr && node->branch_var >= 0) { i_t v = node->branch_var; bool is_binary = (symmetry->is_binary[v] == 1); @@ -392,21 +423,52 @@ class orbital_fixing_t { // in this node or a descendant. stabilizer_nontrivial = true; + // Only fix free variables (not already fixed by branching or propagation) + bool is_free = (marked_b0_[j] == 0 && marked_b1_[j] == 0 && + marked_f0_[j] == 0 && marked_f1_[j] == 0); + if (!is_free) continue; + if (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1) { // The orbit of this variable contains variables in B0 or F0 - // So we can fix this variable to zero (provided its not already in B0 or F0) - if (marked_b0_[j] == 0 && marked_f0_[j] == 0) { fix_zero_.push_back(j); } + // So we can fix this variable to zero + fix_zero_.push_back(j); } if (orbit_has_f1_[o] == 1) { // The orbit of this variable contains variables in F1 - // So we can fix this variable to one (provided its not already in F1) - if (marked_f1_[j] == 0) { fix_one_.push_back(j); } + // So we can fix this variable to one + fix_one_.push_back(j); } } if (!stabilizer_nontrivial) { surviving_generators_.clear(); } + // Detect conflicts: variables in both fix_zero_ and fix_one_ + i_t num_conflicts = 0; + { + std::vector in_fix_zero(num_original_vars_, 0); + for (i_t v : fix_zero_) { in_fix_zero[v] = 1; } + for (i_t v : fix_one_) { + if (in_fix_zero[v]) { num_conflicts++; } + } + } + + if (fix_zero_.size() > 0 || fix_one_.size() > 0 || num_conflicts > 0) { + settings.log.printf("Orbital fixing at node %d (depth %d): " + "B0=%d B1=%d F0=%d F1=%d survivors=%d " + "fixing %d to 0 and %d to 1 (conflicts=%d)\n", + node_ptr->node_id, + node_ptr->depth, + (int)branched_zero_.size(), + (int)branched_one_.size(), + (int)f0_.size(), + (int)f1_.size(), + (int)surviving_generators_.size(), + (int)fix_zero_.size(), + (int)fix_one_.size(), + (int)num_conflicts); + } + // Restore the work arrays for (i_t v : f0_) { orbit_has_f0_[orb_.find_orbit(v)] = 0; @@ -421,14 +483,6 @@ class orbital_fixing_t { f0_.clear(); f1_.clear(); - if (fix_zero_.size() > 0 || fix_one_.size() > 0) { - settings.log.printf("Orbital fixing at node %d (depth %d): fixing %d to 0 and %d to 1\n", - node_ptr->node_id, - node_ptr->depth, - (int)fix_zero_.size(), - (int)fix_one_.size()); - } - // Finally fix the variables in L0 and L1 for (i_t v : fix_zero_) { problem.lower[v] = 0.0; From 24acd6fa7be5455704d62f40559a3ea100a22e5d Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 21 Apr 2026 13:00:42 -0700 Subject: [PATCH 10/19] Fix conflicts coming from root fixings. Fix non-monotonic bound by storing orbital fixing bound changes in mip_node_t --- cpp/src/branch_and_bound/branch_and_bound.cpp | 6 +- cpp/src/branch_and_bound/mip_node.hpp | 8 ++ cpp/src/branch_and_bound/symmetry.hpp | 87 +++++++++++++++++-- 3 files changed, 92 insertions(+), 9 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index bf93ddabf..2d718adf4 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1397,9 +1397,13 @@ dual::status_t branch_and_bound_t::solve_node_lp( if (feasible) { // Perform orbital fixing + i_t orbital_conflicts = 0; auto* of = worker->orbital_fixing.get(); if (of != nullptr && !of->disabled()) { - of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem); + orbital_conflicts = of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem, + worker->start_lower, worker->start_upper); + } else if (of != nullptr) { + of->propagate_cumulative_fixings(node_ptr); } i_t node_iter = 0; diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp index a24f67c3b..260b917fb 100644 --- a/cpp/src/branch_and_bound/mip_node.hpp +++ b/cpp/src/branch_and_bound/mip_node.hpp @@ -271,6 +271,8 @@ class mip_node_t { copy.children[1] = nullptr; copy.status = node_status_t::PENDING; + copy.orbital_fix_zero = orbital_fix_zero; + copy.orbital_fix_one = orbital_fix_one; copy.origin_worker_id = origin_worker_id; copy.creation_seq = creation_seq; return copy; @@ -293,6 +295,12 @@ class mip_node_t { std::vector vstatus; + // Cumulative orbital fixing bound changes from root to this node. + // Stored so that when a child starts a new plunge, the parent's + // orbital fixings can be restored without re-derivation. + std::vector orbital_fix_zero; + std::vector orbital_fix_one; + // Worker-local identification for deterministic ordering: // - origin_worker_id: which worker created this node // - creation_seq: sequence number within that worker (cumulative across horizons, serial) diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 1490b2807..355ade363 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -266,6 +266,15 @@ class orbital_fixing_t { bool disabled () const { return surviving_generators_.empty(); } void disable() { surviving_generators_.clear(); } + // Store the current cumulative orbital fixings into the node without + // running the full orbital fixing computation. Called when orbital + // fixing is disabled (trivial stabilizer) so that children starting + // a new plunge still inherit ancestor fixings. + void propagate_cumulative_fixings(mip_node_t* node_ptr) { + node_ptr->orbital_fix_zero = cumulative_fix_zero_; + node_ptr->orbital_fix_one = cumulative_fix_one_; + } + void reset(mip_symmetry_t* symmetry, mip_node_t* node_ptr) { for (i_t j : branched_zero_) { marked_b0_[j] = 0; @@ -294,18 +303,54 @@ class orbital_fixing_t { } node = node->parent; } + if (node != nullptr && node->parent != nullptr && node->parent->parent != nullptr) { + printf("Orbital fixing at node %d (depth %d): parent is not the root\n", + node_ptr->node_id, + node_ptr->depth); + } surviving_generators_.resize(max_generators_); std::iota(surviving_generators_.begin(), surviving_generators_.end(), 0); + // Seed cumulative fixings from the parent's stored orbital fixings. + // These were accumulated during the parent's plunge and represent all + // orbital fixing bound changes from root to parent. + cumulative_fix_zero_.clear(); + cumulative_fix_one_.clear(); + if (node_ptr->parent != nullptr) { + cumulative_fix_zero_ = node_ptr->parent->orbital_fix_zero; + cumulative_fix_one_ = node_ptr->parent->orbital_fix_one; + } + start_plunge_ = true; } - void orbital_fixing(mip_symmetry_t* symmetry, - const simplex_solver_settings_t& settings, - mip_node_t* node_ptr, - lp_problem_t& problem) + // Returns the number of conflicts (variables in both fix_zero and fix_one). + // When conflicts > 0, no fixings are applied. + i_t orbital_fixing(mip_symmetry_t* symmetry, + const simplex_solver_settings_t& settings, + mip_node_t* node_ptr, + lp_problem_t& problem, + const std::vector& start_lower, + const std::vector& start_upper) { + // At the start of a new plunge, restore the parent's cumulative orbital + // fixings into the problem. These bound changes were derived during the + // parent's plunge and are valid for the entire subtree. Without this, + // orbital fixings from ancestors with larger stabilizers would be lost, + // causing the LP bound to regress and the global lower bound to become + // non-monotonic. + if (start_plunge_) { + for (i_t v : cumulative_fix_zero_) { + problem.lower[v] = 0.0; + problem.upper[v] = 0.0; + } + for (i_t v : cumulative_fix_one_) { + problem.lower[v] = 1.0; + problem.upper[v] = 1.0; + } + } + // Collect binary branchings only; skip general integer branchings i_t v = node_ptr->branch_var; bool is_binary = (symmetry->is_binary[v] == 1); @@ -321,12 +366,15 @@ class orbital_fixing_t { } } + // Collect F0/F1: variables fixed by node-level propagation only. + // Exclude root-level fixings (from reduced cost fixing / bound strengthening + // after cuts) because these can be asymmetric w.r.t. the symmetry group. for (i_t j : symmetry->binary_variables) { - if (marked_b1_[j] == 0 && problem.lower[j] == 1.0) { + if (marked_b1_[j] == 0 && problem.lower[j] == 1.0 && start_lower[j] < 1.0) { f1_.push_back(j); marked_f1_[j] = 1; } - if (marked_b0_[j] == 0 && problem.upper[j] == 0.0) { + if (marked_b0_[j] == 0 && problem.upper[j] == 0.0 && start_upper[j] > 0.0) { f0_.push_back(j); marked_f0_[j] = 1; } @@ -453,7 +501,7 @@ class orbital_fixing_t { } } - if (fix_zero_.size() > 0 || fix_one_.size() > 0 || num_conflicts > 0) { + if (num_conflicts > 0) { settings.log.printf("Orbital fixing at node %d (depth %d): " "B0=%d B1=%d F0=%d F1=%d survivors=%d " "fixing %d to 0 and %d to 1 (conflicts=%d)\n", @@ -483,7 +531,15 @@ class orbital_fixing_t { f0_.clear(); f1_.clear(); - // Finally fix the variables in L0 and L1 + if (num_conflicts > 0) { + // Don't apply any fixings — we want to verify the LP is independently infeasible. + // Still store the cumulative fixings from ancestors (without this node's conflicting fixings). + node_ptr->orbital_fix_zero = cumulative_fix_zero_; + node_ptr->orbital_fix_one = cumulative_fix_one_; + return num_conflicts; + } + + // Apply the fixings for (i_t v : fix_zero_) { problem.lower[v] = 0.0; problem.upper[v] = 0.0; @@ -492,6 +548,15 @@ class orbital_fixing_t { problem.lower[v] = 1.0; problem.upper[v] = 1.0; } + + // Accumulate this node's fixings and store in the node so that + // children starting a new plunge can restore them. + cumulative_fix_zero_.insert(cumulative_fix_zero_.end(), fix_zero_.begin(), fix_zero_.end()); + cumulative_fix_one_.insert(cumulative_fix_one_.end(), fix_one_.begin(), fix_one_.end()); + node_ptr->orbital_fix_zero = cumulative_fix_zero_; + node_ptr->orbital_fix_one = cumulative_fix_one_; + + return 0; } private: @@ -517,6 +582,12 @@ class orbital_fixing_t { std::vector marked_f1_; // marked_f1_[v] = 1 if variable v is in F1 std::vector fix_zero_; // set of variables fixed to 0 by orbital fixing std::vector fix_one_; // set of variables fixed to 1 by orbital fixing + + // Running accumulation of orbital fixings across the current plunge. + // Stored into node_ptr after each successful orbital_fixing() call, + // so children starting a new plunge can restore the parent's fixings. + std::vector cumulative_fix_zero_; + std::vector cumulative_fix_one_; }; template From 659a7c4a1dcebf6604f0f39aa54879e76b111be7 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Tue, 21 Apr 2026 13:10:48 -0700 Subject: [PATCH 11/19] Log orbital fixings --- cpp/src/branch_and_bound/branch_and_bound.cpp | 11 +++++++++++ cpp/src/branch_and_bound/branch_and_bound_worker.hpp | 2 ++ 2 files changed, 13 insertions(+) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 2d718adf4..d24205587 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -733,6 +733,11 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& settings_.log.printf("Explored %d nodes in %.2fs.\n", exploration_stats_.nodes_explored, toc(exploration_stats_.start_time)); + if (exploration_stats_.orbital_fixing_nodes.load() > 0) { + settings_.log.printf("Orbital fixing applied at %lld nodes, %lld total variable fixings\n", + (long long)exploration_stats_.orbital_fixing_nodes.load(), + (long long)exploration_stats_.orbital_fixings_applied.load()); + } settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n", gap, obj, @@ -1400,8 +1405,14 @@ dual::status_t branch_and_bound_t::solve_node_lp( i_t orbital_conflicts = 0; auto* of = worker->orbital_fixing.get(); if (of != nullptr && !of->disabled()) { + i_t prev_fix = node_ptr->orbital_fix_zero.size() + node_ptr->orbital_fix_one.size(); orbital_conflicts = of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem, worker->start_lower, worker->start_upper); + i_t new_fix = node_ptr->orbital_fix_zero.size() + node_ptr->orbital_fix_one.size(); + if (new_fix > prev_fix) { + ++stats.orbital_fixing_nodes; + stats.orbital_fixings_applied += (new_fix - prev_fix); + } } else if (of != nullptr) { of->propagate_cumulative_fixings(node_ptr); } diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index 3d49fb7b1..7dee668bb 100644 --- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp @@ -47,6 +47,8 @@ struct branch_and_bound_stats_t { omp_atomic_t total_lp_iters = 0; omp_atomic_t nodes_since_last_log = 0; omp_atomic_t last_log = 0.0; + omp_atomic_t orbital_fixing_nodes = 0; + omp_atomic_t orbital_fixings_applied = 0; }; template From 7bd9569634c51089720430fc230ea3df0951e623 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Wed, 22 Apr 2026 11:30:20 -0700 Subject: [PATCH 12/19] Skip orbits with conflicting fixings. Fix bug where orbital fixing was accidently being enabled in diving --- cpp/src/branch_and_bound/branch_and_bound.cpp | 3 +- cpp/src/branch_and_bound/symmetry.hpp | 44 +++++++------------ 2 files changed, 19 insertions(+), 28 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index d24205587..8c2897519 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1395,7 +1395,8 @@ dual::status_t branch_and_bound_t::solve_node_lp( bool feasible = worker->set_lp_variable_bounds(node_ptr, settings_); dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED; worker->leaf_edge_norms = edge_norms_; - if (worker->recompute_bounds && worker->orbital_fixing) { + if (worker->recompute_bounds && worker->orbital_fixing && + worker->search_strategy == BEST_FIRST) { worker->orbital_fixing->reset(symmetry_, node_ptr); } diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 355ade363..7a7c018a7 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -303,11 +303,6 @@ class orbital_fixing_t { } node = node->parent; } - if (node != nullptr && node->parent != nullptr && node->parent->parent != nullptr) { - printf("Orbital fixing at node %d (depth %d): parent is not the root\n", - node_ptr->node_id, - node_ptr->depth); - } surviving_generators_.resize(max_generators_); std::iota(surviving_generators_.begin(), surviving_generators_.end(), 0); @@ -456,6 +451,7 @@ class orbital_fixing_t { fix_zero_.clear(); fix_one_.clear(); + i_t num_conflicts = 0; bool stabilizer_nontrivial = false; for (i_t j : symmetry->binary_variables) { i_t o = orb_.find_orbit(j); @@ -476,13 +472,24 @@ class orbital_fixing_t { marked_f0_[j] == 0 && marked_f1_[j] == 0); if (!is_free) continue; - if (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1) { + bool has_zero_source = (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1); + bool has_one_source = (orbit_has_f1_[o] == 1); + + // An orbit with both fix-to-0 and fix-to-1 sources is conflicting: + // branching (B0) broke the symmetry within this orbit, so propagation + // produced asymmetric fixings. Skip it — no valid fixing can be derived. + if (has_zero_source && has_one_source) { + num_conflicts++; + continue; + } + + if (has_zero_source) { // The orbit of this variable contains variables in B0 or F0 // So we can fix this variable to zero fix_zero_.push_back(j); } - if (orbit_has_f1_[o] == 1) { + if (has_one_source) { // The orbit of this variable contains variables in F1 // So we can fix this variable to one fix_one_.push_back(j); @@ -491,20 +498,10 @@ class orbital_fixing_t { if (!stabilizer_nontrivial) { surviving_generators_.clear(); } - // Detect conflicts: variables in both fix_zero_ and fix_one_ - i_t num_conflicts = 0; - { - std::vector in_fix_zero(num_original_vars_, 0); - for (i_t v : fix_zero_) { in_fix_zero[v] = 1; } - for (i_t v : fix_one_) { - if (in_fix_zero[v]) { num_conflicts++; } - } - } - if (num_conflicts > 0) { settings.log.printf("Orbital fixing at node %d (depth %d): " "B0=%d B1=%d F0=%d F1=%d survivors=%d " - "fixing %d to 0 and %d to 1 (conflicts=%d)\n", + "fixing %d to 0 and %d to 1 (skipped %d vars in conflicting orbits)\n", node_ptr->node_id, node_ptr->depth, (int)branched_zero_.size(), @@ -531,15 +528,7 @@ class orbital_fixing_t { f0_.clear(); f1_.clear(); - if (num_conflicts > 0) { - // Don't apply any fixings — we want to verify the LP is independently infeasible. - // Still store the cumulative fixings from ancestors (without this node's conflicting fixings). - node_ptr->orbital_fix_zero = cumulative_fix_zero_; - node_ptr->orbital_fix_one = cumulative_fix_one_; - return num_conflicts; - } - - // Apply the fixings + // Apply the fixings from non-conflicting orbits for (i_t v : fix_zero_) { problem.lower[v] = 0.0; problem.upper[v] = 0.0; @@ -957,6 +946,7 @@ std::unique_ptr> detect_symmetry( combined.add_hook(&generator_hook); dejavu::solver d; + d.set_print(false); d.automorphisms(&g, combined); std::ostringstream grp_size_str; From 2aba39ffe1ae0fb32adaec611427816ad81a6aa4 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Wed, 22 Apr 2026 16:28:06 -0700 Subject: [PATCH 13/19] More fixes --- cpp/src/branch_and_bound/branch_and_bound.cpp | 15 +-- .../branch_and_bound_worker.hpp | 1 + cpp/src/branch_and_bound/symmetry.hpp | 93 ++++++++++--------- cpp/src/mip_heuristics/solve.cu | 11 +++ 4 files changed, 69 insertions(+), 51 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 8c2897519..5d00e79bb 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -733,10 +733,13 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& settings_.log.printf("Explored %d nodes in %.2fs.\n", exploration_stats_.nodes_explored, toc(exploration_stats_.start_time)); - if (exploration_stats_.orbital_fixing_nodes.load() > 0) { - settings_.log.printf("Orbital fixing applied at %lld nodes, %lld total variable fixings\n", + if (exploration_stats_.orbital_fixing_nodes.load() > 0 || + exploration_stats_.orbital_conflict_nodes.load() > 0) { + settings_.log.printf("Orbital fixing applied at %lld nodes, %lld total variable fixings, " + "%lld nodes with conflicting orbits\n", (long long)exploration_stats_.orbital_fixing_nodes.load(), - (long long)exploration_stats_.orbital_fixings_applied.load()); + (long long)exploration_stats_.orbital_fixings_applied.load(), + (long long)exploration_stats_.orbital_conflict_nodes.load()); } settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n", gap, @@ -1403,17 +1406,17 @@ dual::status_t branch_and_bound_t::solve_node_lp( if (feasible) { // Perform orbital fixing - i_t orbital_conflicts = 0; auto* of = worker->orbital_fixing.get(); if (of != nullptr && !of->disabled()) { i_t prev_fix = node_ptr->orbital_fix_zero.size() + node_ptr->orbital_fix_one.size(); - orbital_conflicts = of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem, - worker->start_lower, worker->start_upper); + i_t conflicts = of->orbital_fixing(symmetry_, settings_, node_ptr, worker->leaf_problem, + worker->start_lower, worker->start_upper); i_t new_fix = node_ptr->orbital_fix_zero.size() + node_ptr->orbital_fix_one.size(); if (new_fix > prev_fix) { ++stats.orbital_fixing_nodes; stats.orbital_fixings_applied += (new_fix - prev_fix); } + if (conflicts > 0) { ++stats.orbital_conflict_nodes; } } else if (of != nullptr) { of->propagate_cumulative_fixings(node_ptr); } diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index 7dee668bb..6cced4d47 100644 --- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp @@ -49,6 +49,7 @@ struct branch_and_bound_stats_t { omp_atomic_t last_log = 0.0; omp_atomic_t orbital_fixing_nodes = 0; omp_atomic_t orbital_fixings_applied = 0; + omp_atomic_t orbital_conflict_nodes = 0; }; template diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 7a7c018a7..1eaea3025 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -320,8 +320,8 @@ class orbital_fixing_t { start_plunge_ = true; } - // Returns the number of conflicts (variables in both fix_zero and fix_one). - // When conflicts > 0, no fixings are applied. + // Returns the number of free variables in conflicting orbits (orbits with + // both zero and one sources). i_t orbital_fixing(mip_symmetry_t* symmetry, const simplex_solver_settings_t& settings, mip_node_t* node_ptr, @@ -354,6 +354,9 @@ class orbital_fixing_t { if (node_ptr->branch_var_upper == 0.0) { branched_zero_.push_back(v); marked_b0_[v] = 1; + // Need to recompute: propagation at this node may have fixed new variables + // to 1, invalidating generators that map B0 variables to those newly-fixed vars. + should_recompute = true; } else if (node_ptr->branch_var_lower == 1.0) { branched_one_.push_back(v); marked_b1_[v] = 1; @@ -376,43 +379,61 @@ class orbital_fixing_t { } // In true orbital fixing we would compute the group G' = stabilizer(G, B1) - // Instead we compute a subgroup H of G' that is just - // H = { g in generators(G) | g(B1) = B1 } - // This means that we will miss some fixings. But every fixing we perform will be valid. - // - // We are called within the context of a plunge within the branch and bound tree - // Suppose we computed H for the parent node. - // If we branched on a binary variable x_j then we are in one of two cases: - // Case 1: x_j = 0, then B1 remains unchanged and we do not need to recompute H or its orbits - // Case 2: x_j = 1, then B1 becomes B1 union {x_j} and we only need to check the generators - // used to construct H for the parent node. We call these the "surviving generators". + // Instead we compute a subgroup H of G' using generator filtering, additionally + // checking that for each variable in B0, its image is not fixed to 1. + // This is necessary because propagation using asymmetric cuts can fix variables + // to 1 that the symmetry group maps to 0-branched variables, creating false conflicts. + // See van Doornmalen & Hojny (2022), "A Unified Framework for Symmetry Handling". // - // Note that when we restart the plunge we need to recompute everything. + // We always recompute because propagation at each node can invalidate generators. if (should_recompute) { for (size_t k = 0; k < surviving_generators_.size(); k++) { const i_t h = surviving_generators_[k]; const permutation_t& perm = symmetry->generators.get_generator(h); const std::vector& p = perm.dense_permutation(); - bool stabilizes_b1 = true; + bool stabilizes = true; if (start_plunge_) { // At the start of a plunge, surviving_generators_ was reset to all generators. - // We need to check each generator against all of B1. + // We need to check each generator against all of B1 and B0. for (i_t j : branched_one_) { if (marked_b1_[p[j]] == 0) { - stabilizes_b1 = false; + stabilizes = false; break; } } + // Check B0: for each variable branched to 0, its image must not be fixed to 1. + if (stabilizes) { + for (i_t j : branched_zero_) { + if (problem.lower[p[j]] >= 1.0) { + stabilizes = false; + break; + } + } + } } else { - // B1 = v union B1_from_parent. So the surviving generators already stabilize B1_from_parent. - // We only need to check if the generators stabilize v. - if (marked_b1_[p[v]] == 0) { - stabilizes_b1 = false; + // Incremental check. The surviving generators already stabilize the parent's B1. + // Check the new branching variable. + if (is_binary && node_ptr->branch_var_lower == 1.0) { + // Branched to 1: check g(v) is also in B1 + if (marked_b1_[p[v]] == 0) { + stabilizes = false; + } + } + // Check B0: propagation at this node may have newly fixed variables to 1. + // A generator that was valid at the parent may now map a B0 variable to + // something that propagation fixed to 1. + if (stabilizes) { + for (i_t j : branched_zero_) { + if (problem.lower[p[j]] >= 1.0) { + stabilizes = false; + break; + } + } } } - if (!stabilizes_b1) { + if (!stabilizes) { std::swap(surviving_generators_[k], surviving_generators_.back()); surviving_generators_.pop_back(); k--; @@ -436,8 +457,9 @@ class orbital_fixing_t { for (i_t v : branched_zero_) { orbit_has_b0_[orb_.find_orbit(v)] = 1; } - } else if (is_binary && node_ptr->branch_var_upper == 0.0) { - orbit_has_b0_[orb_.find_orbit(v)] = 1; + } else if (!is_binary) { + // Non-binary variable branching: no symmetry recomputation needed. + // (Non-binary variables are excluded from orbital fixing.) } start_plunge_ = false; @@ -475,9 +497,8 @@ class orbital_fixing_t { bool has_zero_source = (orbit_has_b0_[o] == 1 || orbit_has_f0_[o] == 1); bool has_one_source = (orbit_has_f1_[o] == 1); - // An orbit with both fix-to-0 and fix-to-1 sources is conflicting: - // branching (B0) broke the symmetry within this orbit, so propagation - // produced asymmetric fixings. Skip it — no valid fixing can be derived. + // An orbit with both zero and one sources is conflicting. + // Skip the orbit — no valid additional fixing can be derived from it. if (has_zero_source && has_one_source) { num_conflicts++; continue; @@ -498,33 +519,15 @@ class orbital_fixing_t { if (!stabilizer_nontrivial) { surviving_generators_.clear(); } - if (num_conflicts > 0) { - settings.log.printf("Orbital fixing at node %d (depth %d): " - "B0=%d B1=%d F0=%d F1=%d survivors=%d " - "fixing %d to 0 and %d to 1 (skipped %d vars in conflicting orbits)\n", - node_ptr->node_id, - node_ptr->depth, - (int)branched_zero_.size(), - (int)branched_one_.size(), - (int)f0_.size(), - (int)f1_.size(), - (int)surviving_generators_.size(), - (int)fix_zero_.size(), - (int)fix_one_.size(), - (int)num_conflicts); - } - // Restore the work arrays for (i_t v : f0_) { orbit_has_f0_[orb_.find_orbit(v)] = 0; marked_f0_[v] = 0; } - for (i_t v : f1_) { orbit_has_f1_[orb_.find_orbit(v)] = 0; marked_f1_[v] = 0; } - f0_.clear(); f1_.clear(); @@ -545,7 +548,7 @@ class orbital_fixing_t { node_ptr->orbital_fix_zero = cumulative_fix_zero_; node_ptr->orbital_fix_one = cumulative_fix_one_; - return 0; + return num_conflicts; } private: diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 4db1a5908..4784b7718 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -159,7 +159,18 @@ mip_solution_t run_mip(detail::problem_t& problem, // only call preprocess on scaled problem, so we can compute feasibility on the original problem scaled_problem.preprocess_problem(); scaled_problem.related_vars_time_limit = settings.heuristic_params.related_vars_time_limit; + const i_t n_vars_before = scaled_problem.n_variables; detail::trivial_presolve(scaled_problem); + // Trivial presolve may remove unused variables and renumber the remaining ones. + // When that happens the symmetry generators and binary_variables reference the + // original (pre-trivial-presolve) column indices which are now invalid. + // Discard symmetry to avoid out-of-bounds accesses in orbital fixing. + if (symmetry != nullptr && scaled_problem.n_variables != n_vars_before) { + CUOPT_LOG_INFO("Trivial presolve changed variable count (%d -> %d); " + "disabling orbital fixing to avoid index mismatch", + n_vars_before, scaled_problem.n_variables); + symmetry.reset(); + } detail::mip_solver_t solver(scaled_problem, settings, timer); // initial_upper_bound is in user-space (representation-invariant). From d4658ed1df05e6ac52235d9eea4243df0977a6b4 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 23 Apr 2026 18:48:40 -0700 Subject: [PATCH 14/19] Add lexical reduction --- cpp/src/branch_and_bound/branch_and_bound.cpp | 80 +++++++----- .../branch_and_bound_worker.hpp | 6 + cpp/src/branch_and_bound/symmetry.hpp | 116 ++++++++++++++++++ .../dual_simplex/simplex_solver_settings.hpp | 2 + cpp/src/math_optimization/solver_settings.cu | 2 +- cpp/src/mip_heuristics/solver.cu | 1 + 6 files changed, 175 insertions(+), 32 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 5d00e79bb..d0f389f02 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -741,6 +741,11 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& (long long)exploration_stats_.orbital_fixings_applied.load(), (long long)exploration_stats_.orbital_conflict_nodes.load()); } + if (exploration_stats_.lexical_reduction_nodes.load() > 0) { + settings_.log.printf("Lexical reduction applied at %lld nodes, %lld total variable fixings\n", + (long long)exploration_stats_.lexical_reduction_nodes.load(), + (long long)exploration_stats_.lexical_reduction_fixings_applied.load()); + } settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n", gap, obj, @@ -1421,40 +1426,53 @@ dual::status_t branch_and_bound_t::solve_node_lp( of->propagate_cumulative_fixings(node_ptr); } - i_t node_iter = 0; - f_t lp_start_time = tic(); - - lp_status = dual_phase2_with_advanced_basis(2, - 0, - worker->recompute_basis, - lp_start_time, - worker->leaf_problem, - lp_settings, - leaf_vstatus, - worker->basis_factors, - worker->basic_list, - worker->nonbasic_list, - worker->leaf_solution, - node_iter, - worker->leaf_edge_norms); + if (settings_.symmetry == 2 && worker->lexical_reduction != nullptr) { + i_t lexical_reductions_info = + worker->lexical_reduction->lexical_reduce(symmetry_, node_ptr, worker->leaf_problem); + if (lexical_reductions_info > 0) { + stats.lexical_reduction_nodes++; + stats.lexical_reduction_fixings_applied += lexical_reductions_info; + } + if (lexical_reductions_info == -1) { feasible = false; } + } - if (lp_status == dual::status_t::NUMERICAL) { - log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); - lp_status_t second_status = solve_linear_program_with_advanced_basis(worker->leaf_problem, - lp_start_time, - lp_settings, - worker->leaf_solution, - worker->basis_factors, - worker->basic_list, - worker->nonbasic_list, - leaf_vstatus, - worker->leaf_edge_norms); + if (feasible) { + i_t node_iter = 0; + f_t lp_start_time = tic(); + + lp_status = dual_phase2_with_advanced_basis(2, + 0, + worker->recompute_basis, + lp_start_time, + worker->leaf_problem, + lp_settings, + leaf_vstatus, + worker->basis_factors, + worker->basic_list, + worker->nonbasic_list, + worker->leaf_solution, + node_iter, + worker->leaf_edge_norms); + + if (lp_status == dual::status_t::NUMERICAL) { + log.debug("Numerical issue node %d. Resolving from scratch.\n", node_ptr->node_id); + lp_status_t second_status = + solve_linear_program_with_advanced_basis(worker->leaf_problem, + lp_start_time, + lp_settings, + worker->leaf_solution, + worker->basis_factors, + worker->basic_list, + worker->nonbasic_list, + leaf_vstatus, + worker->leaf_edge_norms); + + lp_status = convert_lp_status_to_dual_status(second_status); + } - lp_status = convert_lp_status_to_dual_status(second_status); + stats.total_lp_solve_time += toc(lp_start_time); + stats.total_lp_iters += node_iter; } - - stats.total_lp_solve_time += toc(lp_start_time); - stats.total_lp_iters += node_iter; } #ifdef LOG_NODE_SIMPLEX diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index 6cced4d47..fdce29014 100644 --- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp @@ -50,6 +50,8 @@ struct branch_and_bound_stats_t { omp_atomic_t orbital_fixing_nodes = 0; omp_atomic_t orbital_fixings_applied = 0; omp_atomic_t orbital_conflict_nodes = 0; + omp_atomic_t lexical_reduction_nodes = 0; + omp_atomic_t lexical_reduction_fixings_applied = 0; }; template @@ -78,6 +80,7 @@ class branch_and_bound_worker_t { pcgenerator_t rng; std::unique_ptr> orbital_fixing; + std::unique_ptr> lexical_reduction; mip_symmetry_t* symmetry_ptr = nullptr; void ensure_orbital_fixing() @@ -85,6 +88,9 @@ class branch_and_bound_worker_t { if (orbital_fixing == nullptr && symmetry_ptr != nullptr) { orbital_fixing = std::make_unique>(*symmetry_ptr); } + if (lexical_reduction == nullptr && symmetry_ptr != nullptr) { + lexical_reduction = std::make_unique>(symmetry_ptr->num_original_vars); + } } bool recompute_basis = true; diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 1eaea3025..7b327e557 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -51,6 +51,15 @@ class permutation_t { const std::vector& dense_permutation() const { return p_; } const std::vector& support() const { return support_; } + permutation_t inverse() const { + std::vector pinv(n_); + for (i_t k = 0; k < n_; k++) { + pinv[p_[k]] = k; + } + return permutation_t(std::move(pinv), support_); + } + + private: i_t n_; std::vector p_; @@ -233,6 +242,8 @@ struct mip_symmetry_t { // orbit_rep[j] = orbit representative of variable j (for j < num_original_vars). std::vector orbit_rep; + generators_t inverse_generators; + }; template @@ -582,6 +593,107 @@ class orbital_fixing_t { std::vector cumulative_fix_one_; }; +template +class lexical_reduction_t { + public: + lexical_reduction_t(i_t num_original_vars) + : branched_lower_(num_original_vars, 0.0), branched_upper_(num_original_vars, 1.0) + { + fixings_.reserve(num_original_vars); + reverse_branched_variables_.reserve(num_original_vars); + } + // Return -1 to prune the node, otherwise return the number of fixings applied. + i_t lexical_reduce(mip_symmetry_t* symmetry, + mip_node_t* node_ptr, + lp_problem_t& problem) + { + reverse_branched_variables_.clear(); + mip_node_t* node = node_ptr; + while (node != nullptr && node->branch_var >= 0) { + i_t v = node->branch_var; + if (symmetry->is_binary[v] == 1) { + reverse_branched_variables_.push_back(v); + branched_lower_[v] = node->branch_var_lower; + branched_upper_[v] = node->branch_var_upper; + } + node = node->parent; + } + + fixings_.clear(); + + bool prune = false; + i_t num_fixings = 0; + for (size_t k = 0; k < symmetry->inverse_generators.num_generators(); k++) { + const permutation_t& perm = symmetry->inverse_generators.get_generator(k); + const std::vector& p = perm.dense_permutation(); + const size_t reverse_size = reverse_branched_variables_.size(); + for (size_t h = reverse_size; h > 0; --h) { + i_t j = reverse_branched_variables_[h-1]; // This orders the variables from the root down to the current node + const i_t p_j = p[j]; + if (p_j == j) continue; + // Compare x[j] with x[p[j] + // x[j] = 1, x[p[j]] = 1, continue to next variable + // x[j] = 1, x[p[j]] = 0, strict greater. stop (continue to next generator), constraint is + // satisfied x[j] = 1, x[p[j]] = free, stop (continue to next generator) x[j] = 0, x[p[j]] + // = 0, continue to next variable x[j] = 0, x[p[j]] = 1, violated. Prune the node x[j] = + // 0, x[p[j]] = free, fix x[p[j]] to 0, continue to the next pair x[j] = free, x[p[j]] = + // any, stop (continue to next generator) + i_t val_j = -1; + if (branched_lower_[j] == branched_upper_[j]) { val_j = branched_lower_[j]; } + i_t val_p_j = -1; + if (branched_lower_[p_j] == branched_upper_[p_j]) { val_p_j = branched_lower_[p_j]; } + if (val_j == -1) { // free. continue to next generator + break; + } + if (val_j == 1 && val_p_j == 1) { + continue; // continue to next variable + } + if (val_j == 1 && val_p_j == 0) { + break; // stop. continue to next generator. Lex constraint is satisfied + } + if (val_j == 1 && val_p_j == -1) { + break; // stop. continue to the next generator. + } + if (val_j == 0 && val_p_j == 0) { + continue; // continue to next variable + } + if (val_j == 0 && val_p_j == 1) { + prune = true; // violated. Prune the node + break; + } + if (val_j == 0 && val_p_j == -1) { + problem.lower[p_j] = 0.0; + problem.upper[p_j] = 0.0; + branched_lower_[p_j] = 0.0; + branched_upper_[p_j] = 0.0; + fixings_.push_back(p_j); + num_fixings++; + continue; // continue to the next pair + } + } + if (prune) break; + } + + for (i_t v : reverse_branched_variables_) { + branched_lower_[v] = 0.0; + branched_upper_[v] = 1.0; + } + + for (i_t v : fixings_) { + branched_lower_[v] = 0.0; + branched_upper_[v] = 1.0; + } + + return prune ? -1 : num_fixings; + } + + private: + std::vector branched_lower_; + std::vector branched_upper_; + std::vector fixings_; + std::vector reverse_branched_variables_; +}; + template std::unique_ptr> detect_symmetry( const user_problem_t& user_problem, @@ -1015,6 +1127,10 @@ std::unique_ptr> detect_symmetry( result->orbit_rep[j] = orb.find_orbit(j); } + for (size_t i = 0; i < result->generators.num_generators(); i++) { + result->inverse_generators.add_generator(result->generators.get_generator(i).inverse()); + } + return result; } diff --git a/cpp/src/dual_simplex/simplex_solver_settings.hpp b/cpp/src/dual_simplex/simplex_solver_settings.hpp index cfc120e47..83f8f7637 100644 --- a/cpp/src/dual_simplex/simplex_solver_settings.hpp +++ b/cpp/src/dual_simplex/simplex_solver_settings.hpp @@ -103,6 +103,7 @@ struct simplex_solver_settings_t { implied_bound_cuts(-1), clique_cuts(-1), strong_chvatal_gomory_cuts(-1), + symmetry(-1), reduced_cost_strengthening(-1), cut_change_threshold(1e-3), cut_min_orthogonality(0.5), @@ -188,6 +189,7 @@ struct simplex_solver_settings_t { i_t clique_cuts; // -1 automatic, 0 to disable, >0 to enable clique 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 i_t reduced_cost_strengthening; // -1 automatic, 0 to disable, >0 to enable reduced cost // strengthening f_t cut_change_threshold; // threshold for cut change diff --git a/cpp/src/math_optimization/solver_settings.cu b/cpp/src/math_optimization/solver_settings.cu index 78fd1f60e..e25e9fee0 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -147,7 +147,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_MIP_RELIABILITY_BRANCHING, &mip_settings.reliability_branching, -1, std::numeric_limits::max(), -1}, {CUOPT_PDLP_PRECISION, reinterpret_cast(&pdlp_settings.pdlp_precision), CUOPT_PDLP_DEFAULT_PRECISION, CUOPT_PDLP_MIXED_PRECISION, CUOPT_PDLP_DEFAULT_PRECISION}, {CUOPT_MIP_SCALING, &mip_settings.mip_scaling, CUOPT_MIP_SCALING_OFF, CUOPT_MIP_SCALING_NO_OBJECTIVE, CUOPT_MIP_SCALING_ON}, - {CUOPT_MIP_SYMMETRY, &mip_settings.symmetry, -1, 1, -1}, + {CUOPT_MIP_SYMMETRY, &mip_settings.symmetry, -1, 2, -1}, // MIP heuristic hyper-parameters (hidden from default --help: name contains "hyper_") {CUOPT_MIP_HYPER_HEURISTIC_POPULATION_SIZE, &mip_settings.heuristic_params.population_size, 1, std::numeric_limits::max(), 32, "max solutions in pool"}, {CUOPT_MIP_HYPER_HEURISTIC_NUM_CPUFJ_THREADS, &mip_settings.heuristic_params.num_cpufj_threads, 0, std::numeric_limits::max(), 8, "parallel CPU FJ climbers"}, diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index fd71e053e..091435ffb 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -355,6 +355,7 @@ solution_t mip_solver_t::run_solver() context.settings.reduced_cost_strengthening == -1 ? 2 : context.settings.reduced_cost_strengthening; + branch_and_bound_settings.symmetry = context.settings.symmetry; if (context.settings.num_cpu_threads < 0) { branch_and_bound_settings.num_threads = std::max(1, omp_get_max_threads() - 1); From 6cdc089a8384d5d2816b498515599f8ef1d73d89 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 23 Apr 2026 19:09:20 -0700 Subject: [PATCH 15/19] Fix malformed comment --- cpp/src/branch_and_bound/symmetry.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 7b327e557..00fcd08a5 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -632,12 +632,13 @@ class lexical_reduction_t { const i_t p_j = p[j]; if (p_j == j) continue; // Compare x[j] with x[p[j] - // x[j] = 1, x[p[j]] = 1, continue to next variable - // x[j] = 1, x[p[j]] = 0, strict greater. stop (continue to next generator), constraint is - // satisfied x[j] = 1, x[p[j]] = free, stop (continue to next generator) x[j] = 0, x[p[j]] - // = 0, continue to next variable x[j] = 0, x[p[j]] = 1, violated. Prune the node x[j] = - // 0, x[p[j]] = free, fix x[p[j]] to 0, continue to the next pair x[j] = free, x[p[j]] = - // any, stop (continue to next generator) + // x[j] = 1, x[p[j]] = 1, continue to next variable + // x[j] = 1, x[p[j]] = 0, strict greater. stop (continue to next generator), constraint is satisfied + // x[j] = 1, x[p[j]] = free, stop (continue to next generator) + // x[j] = 0, x[p[j]] = 0, continue to next variable + // x[j] = 0, x[p[j]] = 1, violated. Prune the node + // x[j] = 0, x[p[j]] = free, fix x[p[j]] to 0, continue to the next variable + // x[j] = free, x[p[j]] = any, stop (continue to next generator) i_t val_j = -1; if (branched_lower_[j] == branched_upper_[j]) { val_j = branched_lower_[j]; } i_t val_p_j = -1; From 1dfa7843b7ead209a8f7640a107bd09dd8e61336 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 23 Apr 2026 19:10:50 -0700 Subject: [PATCH 16/19] Add clang-format off/on to prevent comment from being mangled again --- cpp/src/branch_and_bound/symmetry.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 00fcd08a5..69e3e409f 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -631,6 +631,7 @@ class lexical_reduction_t { i_t j = reverse_branched_variables_[h-1]; // This orders the variables from the root down to the current node const i_t p_j = p[j]; if (p_j == j) continue; + // clang-format off // Compare x[j] with x[p[j] // x[j] = 1, x[p[j]] = 1, continue to next variable // x[j] = 1, x[p[j]] = 0, strict greater. stop (continue to next generator), constraint is satisfied @@ -639,6 +640,7 @@ class lexical_reduction_t { // x[j] = 0, x[p[j]] = 1, violated. Prune the node // x[j] = 0, x[p[j]] = free, fix x[p[j]] to 0, continue to the next variable // x[j] = free, x[p[j]] = any, stop (continue to next generator) + // clang-format on i_t val_j = -1; if (branched_lower_[j] == branched_upper_[j]) { val_j = branched_lower_[j]; } i_t val_p_j = -1; From 9c41f40e70de79cdec78a597de82d133a981541f Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Thu, 23 Apr 2026 21:44:15 -0700 Subject: [PATCH 17/19] Repeat symmetry detection if trivial presolve changes the problems. Use fixings in lexical reduction. Track pruned nodes in lexical reduction --- cpp/src/branch_and_bound/branch_and_bound.cpp | 10 ++++--- .../branch_and_bound_worker.hpp | 1 + cpp/src/branch_and_bound/symmetry.hpp | 26 ++----------------- cpp/src/mip_heuristics/solve.cu | 11 +++++++- 4 files changed, 20 insertions(+), 28 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index d0f389f02..3761831e3 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -742,9 +742,10 @@ void branch_and_bound_t::set_final_solution(mip_solution_t& (long long)exploration_stats_.orbital_conflict_nodes.load()); } if (exploration_stats_.lexical_reduction_nodes.load() > 0) { - settings_.log.printf("Lexical reduction applied at %lld nodes, %lld total variable fixings\n", + settings_.log.printf("Lexical reduction applied at %lld nodes, %lld total variable fixings, %lld nodes pruned\n", (long long)exploration_stats_.lexical_reduction_nodes.load(), - (long long)exploration_stats_.lexical_reduction_fixings_applied.load()); + (long long)exploration_stats_.lexical_reduction_fixings_applied.load(), + (long long)exploration_stats_.lexical_reduction_pruned_nodes.load()); } settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n", gap, @@ -1433,7 +1434,10 @@ dual::status_t branch_and_bound_t::solve_node_lp( stats.lexical_reduction_nodes++; stats.lexical_reduction_fixings_applied += lexical_reductions_info; } - if (lexical_reductions_info == -1) { feasible = false; } + if (lexical_reductions_info == -1) { + feasible = false; + stats.lexical_reduction_pruned_nodes++; + } } if (feasible) { diff --git a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index fdce29014..39d98d0cb 100644 --- a/cpp/src/branch_and_bound/branch_and_bound_worker.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp @@ -52,6 +52,7 @@ struct branch_and_bound_stats_t { omp_atomic_t orbital_conflict_nodes = 0; omp_atomic_t lexical_reduction_nodes = 0; omp_atomic_t lexical_reduction_fixings_applied = 0; + omp_atomic_t lexical_reduction_pruned_nodes = 0; }; template diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 69e3e409f..0c080ab84 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -597,9 +597,7 @@ template class lexical_reduction_t { public: lexical_reduction_t(i_t num_original_vars) - : branched_lower_(num_original_vars, 0.0), branched_upper_(num_original_vars, 1.0) { - fixings_.reserve(num_original_vars); reverse_branched_variables_.reserve(num_original_vars); } // Return -1 to prune the node, otherwise return the number of fixings applied. @@ -613,14 +611,10 @@ class lexical_reduction_t { i_t v = node->branch_var; if (symmetry->is_binary[v] == 1) { reverse_branched_variables_.push_back(v); - branched_lower_[v] = node->branch_var_lower; - branched_upper_[v] = node->branch_var_upper; } node = node->parent; } - fixings_.clear(); - bool prune = false; i_t num_fixings = 0; for (size_t k = 0; k < symmetry->inverse_generators.num_generators(); k++) { @@ -642,9 +636,9 @@ class lexical_reduction_t { // x[j] = free, x[p[j]] = any, stop (continue to next generator) // clang-format on i_t val_j = -1; - if (branched_lower_[j] == branched_upper_[j]) { val_j = branched_lower_[j]; } + if (problem.lower[j] == problem.upper[j]) { val_j = static_cast(problem.lower[j]); } i_t val_p_j = -1; - if (branched_lower_[p_j] == branched_upper_[p_j]) { val_p_j = branched_lower_[p_j]; } + if (problem.lower[p_j] == problem.upper[p_j]) { val_p_j = static_cast(problem.lower[p_j]); } if (val_j == -1) { // free. continue to next generator break; } @@ -667,9 +661,6 @@ class lexical_reduction_t { if (val_j == 0 && val_p_j == -1) { problem.lower[p_j] = 0.0; problem.upper[p_j] = 0.0; - branched_lower_[p_j] = 0.0; - branched_upper_[p_j] = 0.0; - fixings_.push_back(p_j); num_fixings++; continue; // continue to the next pair } @@ -677,23 +668,10 @@ class lexical_reduction_t { if (prune) break; } - for (i_t v : reverse_branched_variables_) { - branched_lower_[v] = 0.0; - branched_upper_[v] = 1.0; - } - - for (i_t v : fixings_) { - branched_lower_[v] = 0.0; - branched_upper_[v] = 1.0; - } - return prune ? -1 : num_fixings; } private: - std::vector branched_lower_; - std::vector branched_upper_; - std::vector fixings_; std::vector reverse_branched_variables_; }; diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 4784b7718..703260e48 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -167,9 +167,18 @@ mip_solution_t run_mip(detail::problem_t& problem, // Discard symmetry to avoid out-of-bounds accesses in orbital fixing. if (symmetry != nullptr && scaled_problem.n_variables != n_vars_before) { CUOPT_LOG_INFO("Trivial presolve changed variable count (%d -> %d); " - "disabling orbital fixing to avoid index mismatch", + "re-detecting symmetry on reduced problem", n_vars_before, scaled_problem.n_variables); symmetry.reset(); + if (settings.symmetry != 0) { + dual_simplex::simplex_solver_settings_t simplex_settings; + simplex_settings.set_log(true); + simplex_settings.time_limit = settings.time_limit; + dual_simplex::user_problem_t reduced_user_problem = + cuopt_problem_to_simplex_problem(scaled_problem.original_problem_ptr->get_handle_ptr(), scaled_problem); + bool has_symmetry_reduced = false; + symmetry = dual_simplex::detect_symmetry(reduced_user_problem, simplex_settings, has_symmetry_reduced); + } } detail::mip_solver_t solver(scaled_problem, settings, timer); From c607017bb579741dd55ef4b7993df01b44bb0efd Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 27 Apr 2026 11:01:17 -0700 Subject: [PATCH 18/19] Detect symmetry after presolve --- cpp/src/branch_and_bound/symmetry.hpp | 6 +++++- cpp/src/mip_heuristics/solve.cu | 29 ++++++++++++++++++++++++++- 2 files changed, 33 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/symmetry.hpp b/cpp/src/branch_and_bound/symmetry.hpp index 0c080ab84..45366744d 100644 --- a/cpp/src/branch_and_bound/symmetry.hpp +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -1100,7 +1100,11 @@ std::unique_ptr> detect_symmetry( settings.log.printf("Total symmetry detection time %f\n", toc(start_time)); - if (!has_symmetry) { return nullptr; } + if (!has_symmetry) { + settings.log.printf("No exploitable symmetry found (%d generators, %d non-trivial orbits, max orbit size %d)\n", + projected_count, num_nontrivial_orbits, max_orbit_size); + return nullptr; + } // Precompute orbit representatives from the projected group's orbits. result->orbit_rep.resize(num_original_vars); diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index f23cb19f9..f73dc3149 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -41,6 +41,13 @@ #include #include +// Choose when to detect symmetry: +// DETECT_SYMMETRY_BEFORE_PRESOLVE: detect on original problem, disable presolve if symmetry found. +// Finds maximum symmetry but loses presolve benefits. +// DETECT_SYMMETRY_AFTER_PRESOLVE: detect after PaPILO + trivial presolve on the reduced problem. +// Presolve runs at full power; symmetry detection on whatever structure remains. +#define DETECT_SYMMETRY_AFTER_PRESOLVE + #include #include @@ -161,10 +168,12 @@ mip_solution_t run_mip(detail::problem_t& problem, scaled_problem.related_vars_time_limit = settings.heuristic_params.related_vars_time_limit; const i_t n_vars_before = scaled_problem.n_variables; detail::trivial_presolve(scaled_problem); + +#ifdef DETECT_SYMMETRY_BEFORE_PRESOLVE // Trivial presolve may remove unused variables and renumber the remaining ones. // When that happens the symmetry generators and binary_variables reference the // original (pre-trivial-presolve) column indices which are now invalid. - // Discard symmetry to avoid out-of-bounds accesses in orbital fixing. + // Re-detect symmetry on the reduced problem. if (symmetry != nullptr && scaled_problem.n_variables != n_vars_before) { CUOPT_LOG_INFO("Trivial presolve changed variable count (%d -> %d); " "re-detecting symmetry on reduced problem", @@ -180,6 +189,21 @@ mip_solution_t run_mip(detail::problem_t& problem, symmetry = dual_simplex::detect_symmetry(reduced_user_problem, simplex_settings, has_symmetry_reduced); } } +#endif + +#ifdef DETECT_SYMMETRY_AFTER_PRESOLVE + // Detect symmetry on the final problem after both PaPILO and trivial presolve. + // Generators will have correct indices for the reduced problem. + if (settings.symmetry != 0) { + dual_simplex::simplex_solver_settings_t simplex_settings; + simplex_settings.set_log(true); + simplex_settings.time_limit = settings.time_limit; + dual_simplex::user_problem_t post_presolve_problem = + cuopt_problem_to_simplex_problem(scaled_problem.original_problem_ptr->get_handle_ptr(), scaled_problem); + bool has_symmetry_post = false; + symmetry = dual_simplex::detect_symmetry(post_presolve_problem, simplex_settings, has_symmetry_post); + } +#endif detail::mip_solver_t solver(scaled_problem, settings, timer); // initial_upper_bound is in user-space (representation-invariant). @@ -333,6 +357,8 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, // Start symmetry detection bool has_symmetry = false; std::unique_ptr> symmetry; + +#ifdef DETECT_SYMMETRY_BEFORE_PRESOLVE if (settings.symmetry != 0) { detail::problem_t problem(op_problem); @@ -344,6 +370,7 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, symmetry = dual_simplex::detect_symmetry(user_problem, simplex_settings, has_symmetry); if (has_symmetry) { settings.presolver = presolver_t::None; } } +#endif if (settings.mip_scaling != CUOPT_MIP_SCALING_OFF) { detail::mip_scaling_strategy_t scaling(op_problem); From 98cef5af3bc0d0ebab5cbe2257e24be2cbe2a141 Mon Sep 17 00:00:00 2001 From: Christopher Maes Date: Mon, 27 Apr 2026 11:42:11 -0700 Subject: [PATCH 19/19] Move detect symmetry back until all of cuOpt presolve is complete --- cpp/src/mip_heuristics/solve.cu | 15 ++------------- cpp/src/mip_heuristics/solver.cu | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+), 13 deletions(-) diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index f73dc3149..a03ba980b 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -191,19 +191,8 @@ mip_solution_t run_mip(detail::problem_t& problem, } #endif -#ifdef DETECT_SYMMETRY_AFTER_PRESOLVE - // Detect symmetry on the final problem after both PaPILO and trivial presolve. - // Generators will have correct indices for the reduced problem. - if (settings.symmetry != 0) { - dual_simplex::simplex_solver_settings_t simplex_settings; - simplex_settings.set_log(true); - simplex_settings.time_limit = settings.time_limit; - dual_simplex::user_problem_t post_presolve_problem = - cuopt_problem_to_simplex_problem(scaled_problem.original_problem_ptr->get_handle_ptr(), scaled_problem); - bool has_symmetry_post = false; - symmetry = dual_simplex::detect_symmetry(post_presolve_problem, simplex_settings, has_symmetry_post); - } -#endif + // Note: DETECT_SYMMETRY_AFTER_PRESOLVE detection is done in solver.cu::run_solver() + // after cuOpt's presolve (probing cache, bounds propagation, trivial presolve) completes. detail::mip_solver_t solver(scaled_problem, settings, timer); // initial_upper_bound is in user-space (representation-invariant). diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 091435ffb..8567c2568 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -15,8 +15,14 @@ #include #include +#include #include #include +#include + +// Must match the setting in solve.cu +#define DETECT_SYMMETRY_AFTER_PRESOLVE + #include #include @@ -286,6 +292,22 @@ solution_t mip_solver_t::run_solver() } context.work_unit_scheduler_.register_context(context.gpu_heur_loop); +#ifdef DETECT_SYMMETRY_AFTER_PRESOLVE + // Detect symmetry after all presolve steps (PaPILO, cuOpt probing, bounds, trivial presolve). + // context.problem_ptr is the final reduced problem with correct variable indices. + if (context.settings.symmetry != 0 && !context.problem_ptr->empty) { + cuopt::linear_programming::dual_simplex::simplex_solver_settings_t simplex_settings; + simplex_settings.set_log(true); + simplex_settings.time_limit = context.settings.time_limit; + cuopt::linear_programming::dual_simplex::user_problem_t post_presolve_problem = + cuopt_problem_to_simplex_problem( + context.problem_ptr->handle_ptr, *context.problem_ptr); + bool has_symmetry_post = false; + context.symmetry = cuopt::linear_programming::dual_simplex::detect_symmetry( + post_presolve_problem, simplex_settings, has_symmetry_post); + } +#endif + namespace dual_simplex = cuopt::linear_programming::dual_simplex; std::future branch_and_bound_status_future; dual_simplex::user_problem_t branch_and_bound_problem(context.problem_ptr->handle_ptr);