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/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 33a2d983c..dc7d1f4bb 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(); 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 new file mode 100644 index 000000000..44a73b9bd --- /dev/null +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -0,0 +1,441 @@ +/* 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 +#include + +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) + : 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; + 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)); + } + + 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; + 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 be0151665..af595e57e 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 @@ -84,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"); @@ -161,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); @@ -303,6 +309,21 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, callback->template setup(op_problem.get_n_variables()); } + // 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; + 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); + 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); if (settings.mip_scaling != CUOPT_MIP_SCALING_OFF) { detail::mip_scaling_strategy_t scaling(op_problem); @@ -472,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}