Run FJ heuristics before and during presolve

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)