diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 9fdcd3ada..8c596dc4f 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -276,6 +276,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) @@ -471,7 +485,8 @@ target_include_directories(cuopt PRIVATE ) target_include_directories(cuopt SYSTEM PRIVATE - "${pslp_SOURCE_DIR}/include" + "${pslp_SOURCE_DIR}/include" + "${dejavu_SOURCE_DIR}" ) target_include_directories(cuopt 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 e69ff7b9a..3761831e3 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), @@ -730,6 +733,20 @@ 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 || + 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_conflict_nodes.load()); + } + if (exploration_stats_.lexical_reduction_nodes.load() > 0) { + 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_pruned_nodes.load()); + } settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n", gap, obj, @@ -1387,42 +1404,79 @@ 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->search_strategy == BEST_FIRST) { + worker->orbital_fixing->reset(symmetry_, node_ptr); + } 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); + // Perform orbital fixing + 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(); + 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); + } - lp_status = convert_lp_status_to_dual_status(second_status); + 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; + stats.lexical_reduction_pruned_nodes++; + } } - stats.total_lp_solve_time += toc(lp_start_time); - stats.total_lp_iters += node_iter; + 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); + } + + stats.total_lp_solve_time += toc(lp_start_time); + stats.total_lp_iters += node_iter; + } } #ifdef LOG_NODE_SIMPLEX @@ -1438,6 +1492,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(); f_t lower_bound = get_lower_bound(); f_t upper_bound = upper_bound_; @@ -1563,6 +1618,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; @@ -1659,7 +1715,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 @@ -1801,7 +1857,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(); @@ -2524,6 +2580,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basic_list, nonbasic_list, basis_update, + symmetry_, pc_); } @@ -2593,6 +2650,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/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/branch_and_bound_worker.hpp b/cpp/src/branch_and_bound/branch_and_bound_worker.hpp index 4de2b43ca..39d98d0cb 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 @@ -46,6 +47,12 @@ 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; + 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 @@ -73,6 +80,20 @@ 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() + { + 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; bool recompute_bounds = true; @@ -169,6 +190,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 +198,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); + workers_[i]->symmetry_ptr = symmetry; idle_workers_.push_front(i); } 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/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 new file mode 100644 index 000000000..45366744d --- /dev/null +++ b/cpp/src/branch_and_bound/symmetry.hpp @@ -0,0 +1,1122 @@ +/* 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 + + +#include "dejavu.h" + +#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_; } + + 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_; + 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]; } + + // 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_; +}; + +// 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 { + 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; + + // 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; + + generators_t inverse_generators; + +}; + +template +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), + 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), + 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 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; + 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->parent; + 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); + + // 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; + } + + // 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, + 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); + 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; + // 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; + should_recompute = true; + } + } + + // 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 && start_lower[j] < 1.0) { + f1_.push_back(j); + marked_f1_[j] = 1; + } + if (marked_b0_[j] == 0 && problem.upper[j] == 0.0 && start_upper[j] > 0.0) { + f0_.push_back(j); + marked_f0_[j] = 1; + } + } + + // In true orbital fixing we would compute the group G' = stabilizer(G, B1) + // 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". + // + // 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 = 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 and B0. + for (i_t j : branched_one_) { + if (marked_b1_[p[j]] == 0) { + 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 { + // 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) { + std::swap(surviving_generators_[k], surviving_generators_.back()); + surviving_generators_.pop_back(); + k--; + } + } + + // Clear old orbit_has values before orbits change + for (i_t v : branched_one_) { + orbit_has_b1_[orb_.find_orbit(v)] = 0; + } + for (i_t v : branched_zero_) { + orbit_has_b0_[orb_.find_orbit(v)] = 0; + } + + orb_.reset(); + orb_.compute_orbits(surviving_generators_, symmetry->generators); + + 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; + } + } else if (!is_binary) { + // Non-binary variable branching: no symmetry recomputation needed. + // (Non-binary variables are excluded from orbital fixing.) + } + start_plunge_ = false; + + 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; + } + + 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); + if (orb_.orbit_size(o) < 2) continue; + + if (orbit_has_b1_[o] == 1) { + // The orbit contains variables in B1 + // So we can't fix any variables in this orbit + continue; + } + + // A non-B1 orbit of size >= 2 exists, so the stabilizer may produce fixings + // 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; + + 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 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; + } + + 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 (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); + } + } + + if (!stabilizer_nontrivial) { surviving_generators_.clear(); } + + // 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(); + + // Apply the fixings from non-conflicting orbits + 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; + } + + // 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 num_conflicts; + } + + private: + i_t num_original_vars_; + 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 + 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 + 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 +class lexical_reduction_t { + public: + lexical_reduction_t(i_t 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); + } + node = node->parent; + } + + 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; + // 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 + // 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) + // clang-format on + i_t val_j = -1; + if (problem.lower[j] == problem.upper[j]) { val_j = static_cast(problem.lower[j]); } + i_t val_p_j = -1; + 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; + } + 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; + num_fixings++; + continue; // continue to the next pair + } + } + if (prune) break; + } + + return prune ? -1 : num_fixings; + } + + private: + std::vector reverse_branched_variables_; +}; + +template +std::unique_ptr> detect_symmetry( + const user_problem_t& user_problem, + const simplex_solver_settings_t& settings, + bool& has_symmetry) +{ + has_symmetry = false; + + 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)); + } + + const i_t num_original_vars = user_problem.num_cols; + + 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); + } + } + } + + // 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) { + 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(&generator_hook); + + dejavu::solver d; + d.set_print(false); + 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 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; + 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]++; + } + } + } + + 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 (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]); + } + } + 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) { + 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); + for (i_t j = 0; j < num_original_vars; j++) { + 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; +} + +} // namespace cuopt::linear_programming::dual_simplex 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 39bcdeb8f..a553891c6 100644 --- a/cpp/src/math_optimization/solver_settings.cu +++ b/cpp/src/math_optimization/solver_settings.cu @@ -146,6 +146,7 @@ solver_settings_t::solver_settings_t() : pdlp_settings(), mip_settings {CUOPT_RANDOM_SEED, &mip_settings.seed, -1, std::numeric_limits::max(), -1}, {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_SYMMETRY, &mip_settings.symmetry, -1, 2, -1}, {CUOPT_MIP_SCALING, &mip_settings.mip_scaling, CUOPT_MIP_SCALING_OFF, CUOPT_MIP_SCALING_NO_OBJECTIVE, CUOPT_MIP_SCALING_NO_OBJECTIVE}, // 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"}, 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 978871a4e..a03ba980b 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -37,6 +37,17 @@ #include #include +#include +#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 @@ -84,7 +95,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"); @@ -154,13 +166,40 @@ 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); +#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. + // 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", + 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); + } + } +#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). // 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); @@ -304,6 +343,24 @@ mip_solution_t solve_mip(optimization_problem_t& op_problem, } auto timer = timer_t(time_limit); + // 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); + 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; } + } +#endif + 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); @@ -472,7 +529,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..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); @@ -355,6 +377,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); @@ -392,7 +415,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 4e225270b..dc6cb3c6b 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -17,6 +17,7 @@ if(BUILD_TESTS) PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/../src" "${CMAKE_CURRENT_SOURCE_DIR}" + "${dejavu_SOURCE_DIR}" ) target_link_libraries(cuopttestutils @@ -55,6 +56,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}