Detect MIP symmetry using dejavu; Perform orbital fixing and lexical reduction

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

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"

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;

cpp/src/branch_and_bound/branch_and_bound.cpp

@@ -8,6 +8,7 @@
 #include <branch_and_bound/branch_and_bound.hpp>
 #include <branch_and_bound/mip_node.hpp>
 #include <branch_and_bound/pseudo_costs.hpp>
+#include <branch_and_bound/symmetry.hpp>
 
 #include <cuts/cuts.hpp>
 #include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
@@ -248,11 +249,13 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
   const simplex_solver_settings_t<i_t, f_t>& solver_settings,
   f_t start_time,
   const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
-  std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table)
+  std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table,
+  mip_symmetry_t<i_t, f_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),
@@ -299,10 +302,11 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
 template <typename i_t, typename f_t>
 f_t branch_and_bound_t<i_t, f_t>::get_lower_bound()
 {
-  f_t lower_bound      = lower_bound_ceiling_.load();
-  f_t heap_lower_bound = node_queue_.get_lower_bound();
-  lower_bound          = std::min(heap_lower_bound, lower_bound);
-  lower_bound          = std::min(worker_pool_.get_lower_bound(), lower_bound);
+  f_t lower_bound        = lower_bound_ceiling_.load();
+  f_t heap_lower_bound   = node_queue_.get_lower_bound();
+  f_t worker_lower_bound = worker_pool_.get_lower_bound();
+  lower_bound            = std::min(heap_lower_bound, lower_bound);
+  lower_bound            = std::min(worker_lower_bound, lower_bound);
 
   if (std::isfinite(lower_bound)) {
     return lower_bound;
@@ -729,6 +733,14 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_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());
+  }
   settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n",
                        gap,
                        obj,
@@ -1386,8 +1398,29 @@ dual::status_t branch_and_bound_t<i_t, f_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) {
+
+    // 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);
+    }
+
     i_t node_iter     = 0;
     f_t lp_start_time = tic();
 
@@ -1437,6 +1470,7 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
   stack.push_front(worker->start_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_;
@@ -1562,6 +1596,7 @@ template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>* worker)
 {
   raft::common::nvtx::range scope("BB::diving_thread");
+  if (worker->orbital_fixing) { worker->orbital_fixing->disable(); }
   logger_t log;
   log.log = false;
 
@@ -1658,7 +1693,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
   std::array<i_t, num_search_strategies> 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
@@ -1800,15 +1835,15 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
 {
-  branch_and_bound_worker_t<i_t, f_t> worker(0, original_lp_, Arow_, var_types_, settings_);
+  worker_pool_.init(1, original_lp_, Arow_, var_types_, symmetry_, settings_);
+  branch_and_bound_worker_t<i_t, f_t>* worker = worker_pool_.get_idle_worker();
 
   f_t lower_bound = get_lower_bound();
   f_t abs_gap     = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound);
   f_t rel_gap     = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
 
   while (solver_status_ == mip_status_t::UNSET && abs_gap > settings_.absolute_mip_gap_tol &&
          rel_gap > settings_.relative_mip_gap_tol && node_queue_.best_first_queue_size() > 0) {
-    bool launched_any_task = false;
     repair_heuristic_solutions();
 
     f_t now = toc(exploration_stats_.start_time);
@@ -1844,8 +1879,8 @@ void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
       continue;
     }
 
-    worker.init_best_first(start_node.value(), original_lp_);
-    plunge_with(&worker);
+    worker->init_best_first(start_node.value(), original_lp_);
+    plunge_with(worker);
 
     lower_bound = get_lower_bound();
     abs_gap     = compute_user_abs_gap(original_lp_, upper_bound_.load(), lower_bound);
@@ -2523,6 +2558,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
                                basic_list,
                                nonbasic_list,
                                basis_update,
+                               symmetry_,
                                pc_);
   }
 
@@ -2592,6 +2628,17 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_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<f_t>(
+      original_lp_.lower, original_lp_.upper);
+    if (removed > 0) {
+      symmetry_->num_generators = static_cast<int>(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;
 

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -47,6 +47,9 @@ struct clique_table_t;
 
 namespace cuopt::linear_programming::dual_simplex {
 
+template <typename i_t, typename f_t>
+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<i_t, f_t>& solver_settings,
                      f_t start_time,
                      const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
-                     std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table = nullptr);
+                     std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table = nullptr,
+                     mip_symmetry_t<i_t, f_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<f_t>& user_guess) { guess_ = user_guess; }
@@ -162,6 +166,7 @@ class branch_and_bound_t {
   const simplex_solver_settings_t<i_t, f_t> settings_;
   const probing_implied_bound_t<i_t, f_t>& probing_implied_bound_;
   std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table_;
+  mip_symmetry_t<i_t, f_t>* symmetry_;
   std::future<std::shared_ptr<detail::clique_table_t<i_t, f_t>>> clique_table_future_;
   std::atomic<bool> signal_extend_cliques_{false};
 

cpp/src/branch_and_bound/branch_and_bound_worker.hpp

@@ -8,6 +8,7 @@
 #pragma once
 
 #include <branch_and_bound/mip_node.hpp>
+#include <branch_and_bound/symmetry.hpp>
 
 #include <dual_simplex/basis_updates.hpp>
 #include <dual_simplex/bounds_strengthening.hpp>
@@ -46,6 +47,9 @@ struct branch_and_bound_stats_t {
   omp_atomic_t<int64_t> total_lp_iters   = 0;
   omp_atomic_t<i_t> nodes_since_last_log = 0;
   omp_atomic_t<f_t> last_log             = 0.0;
+  omp_atomic_t<int64_t> orbital_fixing_nodes = 0;
+  omp_atomic_t<int64_t> orbital_fixings_applied = 0;
+  omp_atomic_t<int64_t> orbital_conflict_nodes = 0;
 };
 
 template <typename i_t, typename f_t>
@@ -73,6 +77,16 @@ class branch_and_bound_worker_t {
 
   pcgenerator_t rng;
 
+  std::unique_ptr<orbital_fixing_t<i_t, f_t>> orbital_fixing;
+  mip_symmetry_t<i_t, f_t>* symmetry_ptr = nullptr;
+
+  void ensure_orbital_fixing()
+  {
+    if (orbital_fixing == nullptr && symmetry_ptr != nullptr) {
+      orbital_fixing = std::make_unique<orbital_fixing_t<i_t, f_t>>(*symmetry_ptr);
+    }
+  }
+
   bool recompute_basis  = true;
   bool recompute_bounds = true;
 
@@ -169,13 +183,15 @@ class branch_and_bound_worker_pool_t {
             const lp_problem_t<i_t, f_t>& original_lp,
             const csr_matrix_t<i_t, f_t>& Arow,
             const std::vector<variable_type_t>& var_type,
+            mip_symmetry_t<i_t, f_t>* symmetry,
             const simplex_solver_settings_t<i_t, f_t>& settings)
   {
     workers_.resize(num_workers);
     num_idle_workers_ = num_workers;
     for (i_t i = 0; i < num_workers; ++i) {
       workers_[i] = std::make_unique<branch_and_bound_worker_t<i_t, f_t>>(
         i, original_lp, Arow, var_type, settings);
+      workers_[i]->symmetry_ptr = symmetry;
       idle_workers_.push_front(i);
     }
 

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<variable_status_t> 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<i_t> orbital_fix_zero;
+  std::vector<i_t> 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)