Heuristic for arc-flow formulations

cpp/src/mip_heuristics/CMakeLists.txt

@@ -48,7 +48,9 @@ set(MIP_NON_LP_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump_kernels.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/fj_cpu.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_cpufj.cu
-  ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_gpufj.cu)
+  ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_gpufj.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/structural/early_structural.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/structural/arc_flow.cu)
 
 # Choose which files to include based on build mode
 if(BUILD_LP_ONLY)

cpp/src/mip_heuristics/early_heuristic.cuh

@@ -74,7 +74,9 @@ class early_heuristic_t {
 
   // NOT thread-safe. solver_obj is in solver-space (always minimization).
   // Uses a private CUDA stream to avoid racing with the FJ solver's stream.
-  void try_update_best(f_t solver_obj, const std::vector<f_t>& assignment)
+  void try_update_best(f_t solver_obj,
+                       const std::vector<f_t>& assignment,
+                       const char* heuristic_name = Derived::name())
   {
     if (solver_obj >= best_objective_) { return; }
     best_objective_ = solver_obj;
@@ -92,7 +94,7 @@ class early_heuristic_t {
     // Log and callback are deferred to the shared incumbent_callback_ which enforces
     // global monotonicity across all early heuristic instances.
     if (incumbent_callback_) {
-      incumbent_callback_(solver_obj, user_obj, user_assignment, Derived::name());
+      incumbent_callback_(solver_obj, user_obj, user_assignment, heuristic_name);
     }
   }
 

cpp/src/mip_heuristics/mip_constants.hpp

@@ -15,11 +15,13 @@
 #define PDLP_INSTANTIATE_FLOAT 1
 
 /* @brief Minimimum number of threads to enable each part of the MIP Solver */
-#define CUOPT_MIP_FJ_REQUIRED_THREAD_COUNT          8
-#define CUOPT_MIP_EARLY_GPUFJ_REQUIRED_THREAD_COUNT 3
-#define CUOPT_MIP_EARLY_CPUFJ_REQUIRED_THREAD_COUNT 2
-#define CUOPT_MIP_BATCH_PDLP_REQUIRED_THREAD_COUNT  3
-#define CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT 3
+#define CUOPT_MIP_FJ_REQUIRED_THREAD_COUNT               8
+#define CUOPT_MIP_EARLY_GPUFJ_REQUIRED_THREAD_COUNT      3
+#define CUOPT_MIP_EARLY_CPUFJ_REQUIRED_THREAD_COUNT      2
+#define CUOPT_MIP_EARLY_STRUCTURAL_REQUIRED_THREAD_COUNT 2
+#define CUOPT_MIP_ROOT_STRUCTURAL_REQUIRED_THREAD_COUNT  3
+#define CUOPT_MIP_BATCH_PDLP_REQUIRED_THREAD_COUNT       3
+#define CUOPT_MIP_CLIQUE_CUTS_REQUIRED_THREAD_COUNT      3
 
 // MIP-only gate: skip the concurrent barrier when fewer threads are available than this
 // (1 PDLP + 1 dual simplex + 1 barrier). Stand-alone LP always runs all three.

cpp/src/mip_heuristics/solve.cu

@@ -19,6 +19,7 @@
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
 #include <mip_heuristics/presolve/trivial_presolve.cuh>
 #include <mip_heuristics/solver.cuh>
+#include <mip_heuristics/structural/early_structural.cuh>
 #include <mip_heuristics/utils.cuh>
 
 #include <pdlp/pdlp.cuh>
@@ -65,6 +66,7 @@
 #include <omp.h>
 
 #include <cmath>
+#include <mutex>
 #include <sstream>
 
 namespace cuopt::mathematical_optimization {
@@ -249,7 +251,13 @@ mip_solution_t<i_t, f_t> run_mip_solver(
     // optimization_problem_t). Its solver-space differs from both the first-pass FJ (original
     // problem) and B&B (post-trivial- presolve), so initial_upper_bound (user-space) is converted
     // via problem.get_solver_obj_from_user_obj.
+
+    // Must outlive early_cpufj/early_structural below, whose destructors join the tasks.
+    std::mutex papilo_callback_mutex;
+    f_t papilo_best_solver_obj = std::numeric_limits<f_t>::infinity();
+
     std::unique_ptr<mip::early_cpufj_t<i_t, f_t>> early_cpufj;
+    std::unique_ptr<mip::early_structural_t<i_t, f_t>> early_structural;
     bool run_early_cpufj = problem.has_papilo_presolve_data() &&
                            settings.determinism_mode != CUOPT_MODE_DETERMINISTIC &&
                            problem.original_problem_ptr->get_n_integers() > 0;
@@ -268,13 +276,21 @@ mip_solution_t<i_t, f_t> run_mip_solver(
          semi_continuous_original_num_variables =
            mip_solver_settings_accessor<i_t, f_t>::get_semi_continuous_original_num_variables(
              settings),
-         ctx_ptr = &solver.context,
+         ctx_ptr                  = &solver.context,
+         papilo_num_original_vars = problem.get_papilo_original_num_variables(),
+         &papilo_callback_mutex,
+         &papilo_best_solver_obj,
          early_fj_start](f_t solver_obj,
                          f_t user_obj,
                          const std::vector<f_t>& assignment,
                          const char* heuristic_name) {
+          std::lock_guard<std::mutex> lock(papilo_callback_mutex);
+          if (solver_obj >= papilo_best_solver_obj) { return; }
+          papilo_best_solver_obj = solver_obj;
+
           std::vector<f_t> user_assignment;
           presolver_ptr->uncrush_primal_solution(assignment, user_assignment);
+          cuopt_assert(user_assignment.size() == (size_t)papilo_num_original_vars, "Size mismatch");
           ctx_ptr->initial_incumbent_assignment = user_assignment;
           ctx_ptr->initial_upper_bound          = user_obj;
           double elapsed =
@@ -302,6 +318,17 @@ mip_solution_t<i_t, f_t> run_mip_solver(
       early_cpufj->start();
       solver.context.early_cpufj_ptr = early_cpufj.get();
       CUOPT_LOG_DEBUG("Started early CPUFJ on papilo-presolved problem during cuOpt presolve");
+
+      early_structural = mip::early_structural_t<i_t, f_t>::create(
+        *problem.original_problem_ptr, settings.get_tolerances(), incumbent_callback);
+      if (early_structural) {
+        if (std::isfinite(initial_upper_bound)) {
+          early_structural->set_best_objective(
+            problem.get_solver_obj_from_user_obj(initial_upper_bound));
+        }
+        early_structural->start();
+        solver.context.early_structural_ptr = early_structural.get();
+      }
     }
 
     auto presolved_sol            = solver.run_solver();
@@ -498,6 +525,7 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
 
     std::unique_ptr<mip::early_cpufj_t<i_t, f_t>> early_cpufj;
     std::unique_ptr<mip::early_gpufj_t<i_t, f_t>> early_gpufj;
+    std::unique_ptr<mip::early_structural_t<i_t, f_t>> early_structural;
 
     bool run_early_fj = run_presolve && settings.determinism_mode != CUOPT_MODE_DETERMINISTIC &&
                         op_problem.get_n_integers() > 0 && op_problem.get_n_constraints() > 0;
@@ -556,6 +584,9 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
         std::make_unique<mip::early_gpufj_t<i_t, f_t>>(op_problem, settings, early_fj_callback);
       early_gpufj->start();
       CUOPT_LOG_DEBUG("Started early GPUFJ during presolve");
+      early_structural = mip::early_structural_t<i_t, f_t>::create(
+        op_problem, settings.get_tolerances(), early_fj_callback);
+      if (early_structural) { early_structural->start(); }
     }
 
     auto constexpr const dual_postsolve = false;
@@ -650,6 +681,17 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
       early_cpufj.reset();
     }
 
+    if (early_structural) {
+      early_structural->stop();
+      if (early_structural->solution_found()) {
+        CUOPT_LOG_DEBUG(
+          "Early structural heuristic (original) found incumbent with objective %.6e "
+          "during presolve",
+          early_structural->get_best_objective());
+      }
+      early_structural.reset();
+    }
+
     // Add early-heuristic incumbents (original-space) to initial_solutions.
     // PaPILO crushing + validation happens downstream in add_user_given_solutions().
     if (!early_incumbent_pool.empty()) {

cpp/src/mip_heuristics/solver.cu

@@ -25,6 +25,7 @@
 
 #include <mip_heuristics/feasibility_jump/early_cpufj.cuh>
 #include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
+#include <mip_heuristics/structural/early_structural.cuh>
 
 #include <raft/sparse/detail/cusparse_wrappers.h>
 #include <raft/core/cusparse_macros.hpp>
@@ -227,6 +228,16 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
     }
   }
 
+  if (context.early_structural_ptr) {
+    context.early_structural_ptr->stop();
+    if (context.early_structural_ptr->solution_found()) {
+      CUOPT_LOG_DEBUG(
+        "Early structural heuristic found incumbent with user-space objective %g "
+        "during presolve",
+        context.early_structural_ptr->get_best_user_objective());
+    }
+  }
+
   if (!presolve_success) {
     CUOPT_LOG_INFO("Problem proven infeasible in presolve");
     sol.set_problem_fully_reduced();
@@ -480,6 +491,20 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
     }
   }
 
+  std::unique_ptr<mip::root_structural_t<i_t, f_t>> root_structural;
+  if (num_threads >= CUOPT_MIP_ROOT_STRUCTURAL_REQUIRED_THREAD_COUNT &&
+      context.settings.determinism_mode != CUOPT_MODE_DETERMINISTIC &&
+      !context.settings.heuristics_only) {
+    root_structural = std::make_unique<mip::root_structural_t<i_t, f_t>>(
+      *context.problem_ptr,
+      context.settings.get_tolerances(),
+      context.preempt_heuristic_solver_,
+      [&dm](const std::vector<f_t>& assignment, f_t objective) {
+        dm.population.add_external_solution(assignment, objective, solution_origin_t::EXTERNAL);
+      });
+    if (!root_structural->recognized()) { root_structural.reset(); }
+  }
+
 #pragma omp taskgroup
   {
     if (!context.settings.heuristics_only) {
@@ -489,11 +514,25 @@ solution_t<i_t, f_t> mip_solver_t<i_t, f_t>::run_solver()
       }
     }
 
+    if (root_structural) {
+#pragma omp task default(shared) priority(CUOPT_DEFAULT_TASK_PRIORITY)
+      {
+        root_structural->run();
+      }
+    }
+
     // Start the primal heuristics
     context.diversity_manager_ptr = &dm;
     sol                           = dm.run_solver();
   }  // implicit barrier for all tasks created in B&B and heuristics
 
+  dm.population.add_external_solutions_to_population();
+  if (dm.population.is_feasible() &&
+      (!sol.get_feasible() ||
+       dm.population.best_feasible().get_objective() < sol.get_objective())) {
+    sol = solution_t<i_t, f_t>(dm.population.best_feasible());
+  }
+
   if (!context.settings.heuristics_only && branch_and_bound->has_solver_space_incumbent()) {
     solution_t<i_t, f_t> branch_and_bound_sol(*context.problem_ptr);
     branch_and_bound_sol.copy_new_assignment(branch_and_bound_solution.x);

cpp/src/mip_heuristics/solver_context.cuh

@@ -33,6 +33,9 @@ class diversity_manager_t;
 template <typename i_t, typename f_t>
 class early_cpufj_t;
 
+template <typename i_t, typename f_t>
+class early_structural_t;
+
 // Aggregate structure containing the global context of the solving process for convenience:
 // The current problem, user settings, raft handle and statistics objects
 template <typename i_t, typename f_t>
@@ -66,6 +69,7 @@ struct mip_solver_context_t {
   work_unit_scheduler_t work_unit_scheduler_{5.0};
 
   early_cpufj_t<i_t, f_t>* early_cpufj_ptr{nullptr};
+  early_structural_t<i_t, f_t>* early_structural_ptr{nullptr};
   // Best upper bound from early heuristics, in user-space.
   // Must be converted to the target solver-space before use:
   //   - B&B: problem_ptr->get_solver_obj_from_user_obj(initial_upper_bound)