Run FJ heuristics before and during presolve (NVIDIA#899)

This presolve lets FJ (both CPU and GPU) run before presolving is completed, in order to provide incumbent callbacks as early as possible.
While we cannot yet crush these incumbents into the fully-presolved problem representation, their primal values can be used to prune the B&B and provide stronger reduced cost fixings.

With this PR, we now find 191 feasible in the first 15s (145 before), and the primal integral in 10min is now 18.7% (was 19.8% before).


## Summary by CodeRabbit

* **New Features**
  * Early feasibility-jump heuristics (CPU & GPU) during presolve that report incumbents via callbacks and seed the main solve.
  * New early-incumbent heuristic framework with lifecycle, timing, and callback support.

* **Bug Fixes**
  * Null-pointer guards for scaling and safer Python callback invocation with GIL handling.
  * Improvement callbacks now fire only for strictly better solutions.

* **Chores**
  * Centralized cutoff handling and propagation of presolve incumbents; support for providing an initial cutoff.

Authors:
  - Alice Boucher (https://github.com/aliceb-nv)

Approvers:
  - Nicolas L. Guidotti (https://github.com/nguidotti)
  - Rajesh Gandham (https://github.com/rg20)
  - Ramakrishnap (https://github.com/rgsl888prabhu)

URL: NVIDIA#899

Author: aliceb-nv

ci/validate_wheel.sh

@@ -26,7 +26,7 @@ if [[ "${package_dir}" == "python/libcuopt" ]]; then
         )
     else
         PYDISTCHECK_ARGS+=(
-            --max-allowed-size-compressed '485Mi'
+            --max-allowed-size-compressed '490Mi'
         )
     fi
 elif [[ "${package_dir}" != "python/cuopt" ]] && \

cpp/CMakeLists.txt

@@ -78,6 +78,28 @@ if(CMAKE_COMPILER_IS_GNUCXX)
   list(APPEND CUOPT_CXX_FLAGS -Werror -Wno-error=deprecated-declarations)
 endif(CMAKE_COMPILER_IS_GNUCXX)
 
+# Papilo pulls in Boost.Multiprecision float128 support, which expects quadmath.h from the GCC
+# toolchain internals. Conda clang ships libquadmath, but does not surface the matching GCC
+# internal include directory by default. Add it late in the search order so clang still prefers its
+# own builtin intrinsic headers.
+if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
+  execute_process(
+    COMMAND ${CMAKE_CXX_COMPILER} --print-file-name=libquadmath.a
+    OUTPUT_VARIABLE CUOPT_QUADMATH_LIB
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+  )
+
+  if(IS_ABSOLUTE "${CUOPT_QUADMATH_LIB}")
+    get_filename_component(CUOPT_QUADMATH_LIBDIR "${CUOPT_QUADMATH_LIB}" DIRECTORY)
+    set(CUOPT_QUADMATH_INCLUDEDIR "${CUOPT_QUADMATH_LIBDIR}/include")
+
+    if(EXISTS "${CUOPT_QUADMATH_INCLUDEDIR}/quadmath.h")
+      message(STATUS "Adding clang fallback include for quadmath: ${CUOPT_QUADMATH_INCLUDEDIR}")
+      add_compile_options("$<$<COMPILE_LANGUAGE:CXX>:-idirafter${CUOPT_QUADMATH_INCLUDEDIR}>")
+    endif()
+  endif()
+endif()
+
 # To use sanitizer with cuda runtime, one must follow a few steps:
 # 1. Run the binary with env var set: LD_PRELOAD="$(gcc -print-file-name=libasan.so)" ASAN_OPTIONS='protect_shadow_gap=0:replace_intrin=0'
 # 2. (Optional) To run with a debugger (gdb or cuda-gdb) use the additional ASAN option alloc_dealloc_mismatch=0
@@ -204,6 +226,12 @@ if (FETCH_RAPIDS)
   include(cmake/thirdparty/get_cccl.cmake)
   include(cmake/thirdparty/get_rmm.cmake)
   include(cmake/thirdparty/get_raft.cmake)
+  # Source-built RMM can hide out-of-line utility symbols such as
+  # rmm::align_up / rmm::get_current_cuda_device when built with hidden visibility on clang.
+  # Force default visibility on the fetched rmm target until this is fixed upstream/figured out.
+  if(TARGET rmm)
+    set_target_properties(rmm PROPERTIES CXX_VISIBILITY_PRESET default)
+  endif()
 else()
   find_package(CCCL REQUIRED)
   find_package(RMM REQUIRED)
@@ -415,9 +443,13 @@ add_library(cuopt::cuopt ALIAS cuopt)
 # - include paths ---------------------------------------------------------------------------------
 message(STATUS "target include directories CUDSS_INCLUDES = ${CUDSS_INCLUDE}")
 
-target_include_directories(cuopt SYSTEM PRIVATE
+# Adding Papilo as a system include messes up clang's include resolution if papilo is already installed as a conda package
+target_include_directories(cuopt PRIVATE
   "${papilo_SOURCE_DIR}/src"
   "${papilo_BINARY_DIR}"
+)
+
+target_include_directories(cuopt SYSTEM PRIVATE
   "${pslp_SOURCE_DIR}/include"
 )
 

cpp/include/cuopt/linear_programming/cpu_optimization_problem_solution.hpp

@@ -132,9 +132,9 @@ class cpu_lp_solution_t : public lp_solution_interface_t<i_t, f_t> {
 
   i_t get_reduced_cost_size() const override { return reduced_cost_.size(); }
 
-  f_t get_objective_value(i_t = 0) const override { return primal_objective_; }
+  f_t get_objective_value(i_t) const override { return primal_objective_; }
 
-  f_t get_dual_objective_value(i_t = 0) const override { return dual_objective_; }
+  f_t get_dual_objective_value(i_t) const override { return dual_objective_; }
 
   pdlp_termination_status_t get_termination_status(i_t = 0) const override
   {

cpp/include/cuopt/linear_programming/optimization_problem_solution.hpp

@@ -90,9 +90,9 @@ class gpu_lp_solution_t : public lp_solution_interface_t<i_t, f_t> {
     return result;
   }
 
-  f_t get_objective_value(i_t id = 0) const override { return solution_.get_objective_value(id); }
+  f_t get_objective_value(i_t id) const override { return solution_.get_objective_value(id); }
 
-  f_t get_dual_objective_value(i_t id = 0) const override
+  f_t get_dual_objective_value(i_t id) const override
   {
     return solution_.get_dual_objective_value(id);
   }

cpp/include/cuopt/linear_programming/optimization_problem_solution_interface.hpp

@@ -237,7 +237,7 @@ class lp_solution_interface_t : public optimization_problem_solution_interface_t
    * @brief Get solve time
    * @return Total solve time in seconds
    */
-  virtual f_t get_solve_time() const = 0;
+  virtual f_t get_solve_time() const override = 0;
 
   /**
    * @brief Get primal objective value
@@ -391,31 +391,31 @@ class mip_solution_interface_t : public optimization_problem_solution_interface_
    * @brief Get solution as host vector
    * @return Host vector of solution
    */
-  virtual std::vector<f_t> get_solution_host() const = 0;
+  virtual std::vector<f_t> get_solution_host() const override = 0;
 
   /**
    * @brief Get objective value
    * @return Objective value
    */
-  virtual f_t get_objective_value() const = 0;
+  virtual f_t get_objective_value() const override = 0;
 
   /**
    * @brief Get solve time
    * @return Total solve time in seconds
    */
-  virtual f_t get_solve_time() const = 0;
+  virtual f_t get_solve_time() const override = 0;
 
   /**
    * @brief Get MIP gap
    * @return MIP gap
    */
-  virtual f_t get_mip_gap() const = 0;
+  virtual f_t get_mip_gap() const override = 0;
 
   /**
    * @brief Get solution bound
    * @return Solution bound
    */
-  virtual f_t get_solution_bound() const = 0;
+  virtual f_t get_solution_bound() const override = 0;
 
   /**
    * @brief Get termination status

cpp/include/cuopt/linear_programming/utilities/callbacks_implems.hpp

@@ -32,11 +32,12 @@ class default_get_solution_callback_t : public get_solution_callback_t {
                     void* solution_bound,
                     void* user_data) override
   {
-    PyObject* numpy_matrix = get_numpy_array(data, n_variables);
-    PyObject* numpy_array  = get_numpy_array(objective_value, 1);
-    PyObject* numpy_bound  = get_numpy_array(solution_bound, 1);
-    PyObject* py_user_data = user_data == nullptr ? Py_None : static_cast<PyObject*>(user_data);
-    PyObject* res          = PyObject_CallMethod(this->pyCallbackClass,
+    PyGILState_STATE gstate = PyGILState_Ensure();
+    PyObject* numpy_matrix  = get_numpy_array(data, n_variables);
+    PyObject* numpy_array   = get_numpy_array(objective_value, 1);
+    PyObject* numpy_bound   = get_numpy_array(solution_bound, 1);
+    PyObject* py_user_data  = user_data == nullptr ? Py_None : static_cast<PyObject*>(user_data);
+    PyObject* res           = PyObject_CallMethod(this->pyCallbackClass,
                                         "get_solution",
                                         "(OOOO)",
                                         numpy_matrix,
@@ -47,6 +48,7 @@ class default_get_solution_callback_t : public get_solution_callback_t {
     Py_DECREF(numpy_array);
     Py_DECREF(numpy_bound);
     if (res != nullptr) { Py_DECREF(res); }
+    PyGILState_Release(gstate);
   }
 
   PyObject* pyCallbackClass;
@@ -69,11 +71,12 @@ class default_set_solution_callback_t : public set_solution_callback_t {
                     void* solution_bound,
                     void* user_data) override
   {
-    PyObject* numpy_matrix = get_numpy_array(data, n_variables);
-    PyObject* numpy_array  = get_numpy_array(objective_value, 1);
-    PyObject* numpy_bound  = get_numpy_array(solution_bound, 1);
-    PyObject* py_user_data = user_data == nullptr ? Py_None : static_cast<PyObject*>(user_data);
-    PyObject* res          = PyObject_CallMethod(this->pyCallbackClass,
+    PyGILState_STATE gstate = PyGILState_Ensure();
+    PyObject* numpy_matrix  = get_numpy_array(data, n_variables);
+    PyObject* numpy_array   = get_numpy_array(objective_value, 1);
+    PyObject* numpy_bound   = get_numpy_array(solution_bound, 1);
+    PyObject* py_user_data  = user_data == nullptr ? Py_None : static_cast<PyObject*>(user_data);
+    PyObject* res           = PyObject_CallMethod(this->pyCallbackClass,
                                         "set_solution",
                                         "(OOOO)",
                                         numpy_matrix,
@@ -84,6 +87,7 @@ class default_set_solution_callback_t : public set_solution_callback_t {
     Py_DECREF(numpy_array);
     Py_DECREF(numpy_bound);
     if (res != nullptr) { Py_DECREF(res); }
+    PyGILState_Release(gstate);
   }
 
   PyObject* pyCallbackClass;

cpp/src/branch_and_bound/branch_and_bound.cpp

@@ -919,7 +919,7 @@ struct nondeterministic_policy_t : tree_update_policy_t<i_t, f_t> {
   {
   }
 
-  f_t upper_bound() const override { return bnb.upper_bound_.load(); }
+  f_t upper_bound() const override { return bnb.get_cutoff(); }
 
   void update_pseudo_costs(mip_node_t<i_t, f_t>* node, f_t leaf_obj) override
   {
@@ -1337,10 +1337,11 @@ dual::status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
 
   simplex_solver_settings_t lp_settings = settings_;
   lp_settings.set_log(false);
+  f_t cutoff = get_cutoff();
   if (original_lp_.objective_is_integral) {
-    lp_settings.cut_off = std::ceil(upper_bound_ - settings_.integer_tol) + settings_.dual_tol;
+    lp_settings.cut_off = std::ceil(cutoff - settings_.integer_tol) + settings_.dual_tol;
   } else {
-    lp_settings.cut_off = upper_bound_ + settings_.dual_tol;
+    lp_settings.cut_off = cutoff + settings_.dual_tol;
   }
   lp_settings.inside_mip    = 2;
   lp_settings.time_limit    = settings_.time_limit - toc(exploration_stats_.start_time);
@@ -1447,7 +1448,7 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
     // - The lower bound of the parent is lower or equal to its children
     worker->lower_bound = lower_bound;
 
-    if (lower_bound > upper_bound) {
+    if (lower_bound > get_cutoff()) {
       search_tree_.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound);
       search_tree_.update(node_ptr, node_status_t::FATHOMED);
       worker->recompute_basis  = true;
@@ -1557,7 +1558,7 @@ void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>
     f_t rel_gap         = user_relative_gap(original_lp_, upper_bound, lower_bound);
     worker->lower_bound = lower_bound;
 
-    if (node_ptr->lower_bound > upper_bound) {
+    if (node_ptr->lower_bound > get_cutoff()) {
       worker->recompute_basis  = true;
       worker->recompute_bounds = true;
       continue;
@@ -1696,7 +1697,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
         std::optional<mip_node_t<i_t, f_t>*> start_node = node_queue_.pop_best_first();
 
         if (!start_node.has_value()) { continue; }
-        if (upper_bound_ < start_node.value()->lower_bound) {
+        if (get_cutoff() < start_node.value()->lower_bound) {
           // This node was put on the heap earlier but its lower bound is now greater than the
           // current upper bound
           search_tree_.graphviz_node(
@@ -1720,7 +1721,7 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
         std::optional<mip_node_t<i_t, f_t>*> start_node = node_queue_.pop_diving();
 
         if (!start_node.has_value()) { continue; }
-        if (upper_bound_ < start_node.value()->lower_bound ||
+        if (get_cutoff() < start_node.value()->lower_bound ||
             start_node.value()->depth < diving_settings.min_node_depth) {
           continue;
         }
@@ -1788,7 +1789,7 @@ void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
     std::optional<mip_node_t<i_t, f_t>*> start_node = node_queue_.pop_best_first();
 
     if (!start_node.has_value()) { continue; }
-    if (upper_bound_ < start_node.value()->lower_bound) {
+    if (get_cutoff() < start_node.value()->lower_bound) {
       // This node was put on the heap earlier but its lower bound is now greater than the
       // current upper bound
       search_tree_.graphviz_node(
@@ -2275,12 +2276,12 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
         return mip_status_t::NUMERICAL;
       }
 
-      if (settings_.reduced_cost_strengthening >= 1 && upper_bound_.load() < last_upper_bound) {
+      if (settings_.reduced_cost_strengthening >= 1 && get_cutoff() < last_upper_bound) {
         mutex_upper_.lock();
-        last_upper_bound = upper_bound_.load();
+        last_upper_bound = get_cutoff();
         std::vector<f_t> lower_bounds;
         std::vector<f_t> upper_bounds;
-        find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds);
+        find_reduced_cost_fixings(get_cutoff(), lower_bounds, upper_bounds);
         mutex_upper_.unlock();
         mutex_original_lp_.lock();
         original_lp_.lower = lower_bounds;
@@ -2467,10 +2468,10 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     return solver_status_;
   }
 
-  if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) {
+  if (settings_.reduced_cost_strengthening >= 2 && get_cutoff() < last_upper_bound) {
     std::vector<f_t> lower_bounds;
     std::vector<f_t> upper_bounds;
-    i_t num_fixed = find_reduced_cost_fixings(upper_bound_.load(), lower_bounds, upper_bounds);
+    i_t num_fixed = find_reduced_cost_fixings(get_cutoff(), lower_bounds, upper_bounds);
     if (num_fixed > 0) {
       std::vector<bool> bounds_changed(original_lp_.num_cols, true);
       std::vector<char> row_sense;
@@ -2574,7 +2575,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
         std::optional<mip_node_t<i_t, f_t>*> start_node = node_queue_.pop_best_first();
 
         if (!start_node.has_value()) { continue; }
-        if (upper_bound_ < start_node.value()->lower_bound) {
+        if (get_cutoff() < start_node.value()->lower_bound) {
           // This node was put on the heap earlier but its lower bound is now greater than the
           // current upper bound
           search_tree_.graphviz_node(
@@ -3416,7 +3417,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_sort_replay_events(
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::deterministic_prune_worker_nodes_vs_incumbent()
 {
-  f_t upper_bound = upper_bound_.load();
+  f_t upper_bound = get_cutoff();
 
   for (auto& worker : *deterministic_workers_) {
     // Check nodes in plunge stack - filter in place
@@ -3552,14 +3553,14 @@ void branch_and_bound_t<i_t, f_t>::deterministic_populate_diving_heap()
   const int num_diving                  = deterministic_diving_workers_->size();
   constexpr int target_nodes_per_worker = 10;
   const int target_total                = num_diving * target_nodes_per_worker;
-  f_t upper_bound                       = upper_bound_.load();
+  f_t cutoff                            = get_cutoff();
 
   // Collect candidate nodes from BFS worker backlog heaps
   std::vector<std::pair<mip_node_t<i_t, f_t>*, f_t>> candidates;
 
   for (auto& worker : *deterministic_workers_) {
     for (auto* node : worker.backlog.data()) {
-      if (node->lower_bound < upper_bound) {
+      if (node->lower_bound < cutoff) {
         f_t score = node->objective_estimate;
         if (score >= inf) { score = node->lower_bound; }
         candidates.push_back({node, score});

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -115,6 +115,20 @@ class branch_and_bound_t {
 
   void set_concurrent_lp_root_solve(bool enable) { enable_concurrent_lp_root_solve_ = enable; }
 
+  // Set a cutoff bound from an external source (e.g., early FJ during presolve).
+  // Used for node pruning and reduced cost strengthening but NOT for gap computation.
+  // Unlike upper_bound_, this does not imply a verified incumbent solution exists.
+  //
+  // IMPORTANT: `bound` must be in B&B's internal objective space, i.e. the space of
+  // original_lp_ where:  user_obj = obj_scale * (internal_obj + obj_constant).
+  // The caller (solver.cu) converts from user-space via
+  //   problem_ptr->get_solver_obj_from_user_obj(user_cutoff)
+  // which accounts for both the presolve objective offset and maximization.
+  void set_initial_cutoff(f_t bound) { initial_cutoff_ = bound; }
+
+  // Effective cutoff for node pruning: min of verified incumbent and external cutoff.
+  f_t get_cutoff() const { return std::min(upper_bound_.load(), initial_cutoff_); }
+
   // Repair a low-quality solution from the heuristics.
   bool repair_solution(const std::vector<f_t>& leaf_edge_norms,
                        const std::vector<f_t>& potential_solution,
@@ -179,9 +193,13 @@ class branch_and_bound_t {
   // Mutex for upper bound
   omp_mutex_t mutex_upper_;
 
-  // Global variable for upper bound
+  // Verified incumbent bound (only set when B&B has an actual integer-feasible solution).
   omp_atomic_t<f_t> upper_bound_;
 
+  // External cutoff from early heuristics (for pruning only, no verified solution).
+  // Must be in B&B internal objective space (see set_initial_cutoff).
+  f_t initial_cutoff_{std::numeric_limits<f_t>::infinity()};
+
   // Global variable for incumbent. The incumbent should be updated with the upper bound
   mip_solution_t<i_t, f_t> incumbent_;
 

cpp/src/branch_and_bound/mip_node.hpp

@@ -36,6 +36,27 @@ bool inactive_status(node_status_t status);
 template <typename i_t, typename f_t>
 class mip_node_t {
  public:
+  ~mip_node_t()
+  {
+    // Iterative teardown to avoid stack overflow on deep trees.
+    // Detach all descendants breadth-first, then destroy them as leaves.
+    std::vector<std::unique_ptr<mip_node_t>> nodes;
+    for (auto& c : children) {
+      if (c) { nodes.push_back(std::move(c)); }
+    }
+    // nodes.size() grows so that this loop only terminates when only leaves remain
+    for (size_t i = 0; i < nodes.size(); ++i) {
+      for (auto& c : nodes[i]->children) {
+        if (c) { nodes.push_back(std::move(c)); }
+      }
+    }
+
+    // scope-exit ensure destruction of all detached leaves
+  }
+
+  mip_node_t(mip_node_t&&)            = default;
+  mip_node_t& operator=(mip_node_t&&) = default;
+
   mip_node_t()
     : status(node_status_t::PENDING),
       lower_bound(-std::numeric_limits<f_t>::infinity()),

cpp/src/mip_heuristics/CMakeLists.txt

@@ -41,7 +41,9 @@ set(MIP_NON_LP_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/presolve/conflict_graph/clique_table.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/feasibility_jump.cu
   ${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/fj_cpu.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_cpufj.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump/early_gpufj.cu)
 
 # Choose which files to include based on build mode
 if(BUILD_LP_ONLY)