Root Heuristics (RINS, CPU FJ)

cpp/include/cuopt/mathematical_optimization/mip/submip_hyper_params.hpp

@@ -41,6 +41,11 @@ struct mip_submip_hyper_params_t {
   // number of simplex iteration from the parent B&B.
   f_t iteration_limit_ratio = 0.8;
 
+  // If there is not enough variables fixed or we already found an improving solution,
+  // perform a short DFS to quickly find a feasible solution. This setting controls
+  // the maximum number of nodes allow for backtracking.
+  i_t dfs_max_backtrack = 5;
+
   // Run CPU FJ over the sub-MIP
   bool enable_cpufj = true;
 };

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -34,12 +34,14 @@
 #include <cuopt/mathematical_optimization/pdlp/solver_settings.hpp>
 
 #include <mip_heuristics/presolve/third_party_presolve.hpp>
+#include <mip_heuristics/root_heuristics.hpp>
 
 #include <omp.h>
 
 #include <atomic>
 #include <functional>
 #include <future>
+#include <list>
 #include <memory>
 #include <vector>
 
@@ -370,22 +372,26 @@ class branch_and_bound_t {
   bool launch_rins_worker(const std::vector<f_t>& sol);
   void set_solution_from_submip(const std::vector<f_t>& solution,
                                 const third_party_presolve_t<i_t, f_t>& presolver,
-                                f_t fixrate,
-                                f_t obj);
+                                f_t fixrate);
 
   // Solve the RINS sub-MIP
   void solve_submip(diving_worker_t<i_t, f_t>* worker,
                     const std::vector<f_t>& current_incumbent,
                     i_t num_var_fixed,
                     i_t num_integers,
                     i_t submip_level,
-                    std::string_view log_prefix);
+                    std::string_view log_prefix,
+                    bool is_root_heuristic);
 
   // Creates and solves the RINS sub-MIP
-  void rins(diving_worker_t<i_t, f_t>* rins_worker, const std::vector<f_t>& node_solution);
-
-  // Get the simplex settings for solving the LP of a single node
-  simplex::simplex_solver_settings_t<i_t, f_t> get_node_lp_settings();
+  void rins(diving_worker_t<i_t, f_t>* worker,
+            const std::vector<f_t>& current_incumbent,
+            bool is_root_heuristic);
+
+  void launch_root_heuristics(const simplex::lp_problem_t<i_t, f_t>& lp,
+                              const std::vector<f_t>& sol,
+                              std::list<root_heuristics_t<i_t, f_t>>& heuristics,
+                              omp_atomic_t<i_t>* worker_count);
 
   // Solve the LP relaxation of a leaf node
   simplex::dual_status_t solve_node_lp(mip_node_t<i_t, f_t>* node_ptr,

cpp/src/branch_and_bound/worker.hpp

@@ -243,6 +243,8 @@ class diving_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
   // The best-first worker that is associated with this diving worker. Used for controlling the
   // number of active diving workers.
   bfs_worker_t<i_t, f_t>* bfs_worker{nullptr};
+
+  std::atomic<int> halt = false;
 };
 
 struct submip_stats_t {

cpp/src/dual_simplex/solve.cpp

@@ -104,6 +104,12 @@ f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj)
   return user_obj;
 }
 
+template <typename i_t, typename f_t>
+f_t compute_solver_objective(const lp_problem_t<i_t, f_t>& lp, f_t user_obj)
+{
+  return user_obj / lp.obj_scale - lp.obj_constant;
+}
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original_lp,
                                           const f_t start_time,
@@ -813,6 +819,8 @@ template double compute_user_objective<int, double>(const lp_problem_t<int, doub
 
 template double compute_user_objective(const lp_problem_t<int, double>& lp, double obj);
 
+template double compute_solver_objective(const lp_problem_t<int, double>& lp, double user_obj);
+
 template lp_status_t solve_linear_program_advanced(
   const lp_problem_t<int, double>& original_lp,
   const double start_time,

cpp/src/dual_simplex/solve.hpp

@@ -63,6 +63,9 @@ f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, const std::vector<f
 template <typename i_t, typename f_t>
 f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj);
 
+template <typename i_t, typename f_t>
+f_t compute_solver_objective(const lp_problem_t<i_t, f_t>& lp, f_t user_obj);
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_advanced(const lp_problem_t<i_t, f_t>& original_lp,
                                           const f_t start_time,

cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cu

@@ -1036,9 +1036,20 @@ void fj_t<i_t, f_t>::resize_vectors(const raft::handle_t* handle_ptr)
   climbers[0]->grid_delta_buf.resize(update_weights_launch_dims.first.x, handle_ptr->get_stream());
 
   // FJ related vars
-  cstr_weights.resize(pb_ptr->n_constraints, handle_ptr->get_stream());
-  cstr_right_weights.resize(pb_ptr->n_constraints, handle_ptr->get_stream());
-  cstr_left_weights.resize(pb_ptr->n_constraints, handle_ptr->get_stream());
+  // the problem can gain constraints between two runs (e.g. the objective cutting plane added by
+  // the feasibility pump), and resize leaves the new elements uninitialized: give them the default
+  // weight
+  auto resize_weights = [&](rmm::device_uvector<f_t>& weights) {
+    const auto old_size = weights.size();
+    weights.resize(pb_ptr->n_constraints, handle_ptr->get_stream());
+    if (old_size < weights.size()) {
+      thrust::uninitialized_fill(
+        handle_ptr->get_thrust_policy(), weights.begin() + old_size, weights.end(), 1.);
+    }
+  };
+  resize_weights(cstr_weights);
+  resize_weights(cstr_right_weights);
+  resize_weights(cstr_left_weights);
   constraint_lower_bounds_csr.resize(pb_ptr->coefficients.size(), handle_ptr->get_stream());
   constraint_upper_bounds_csr.resize(pb_ptr->coefficients.size(), handle_ptr->get_stream());
   cstr_coeff_reciprocal.resize(pb_ptr->coefficients.size(), handle_ptr->get_stream());

cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu

@@ -1814,13 +1814,19 @@ void fj_cpu_worker_t<i_t, f_t>::create_worker(
 }
 
 template <typename i_t, typename f_t>
-void fj_cpu_worker_t<i_t, f_t>::run_async(f_t time_limit, double work_unit_limit)
+void fj_cpu_worker_t<i_t, f_t>::run_async(f_t time_limit,
+                                          double work_unit_limit,
+                                          omp_atomic_t<i_t>* worker_count)
 {
   if (!fj_cpu) return;
 
-#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit) \
+  if (worker_count) ++(*worker_count);
+#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit, worker_count) \
   priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : *fj_cpu)
-  cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit);
+  {
+    cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit);
+    if (worker_count) --(*worker_count);
+  }
 }
 
 template <typename i_t, typename f_t>

cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh

@@ -47,8 +47,9 @@ struct fj_cpu_worker_t {
 
   // Run the worker asynchronously (i.e., launch an openmp task and then continue the
   // execution). Call `stop()` for stopping the worker
-  void run_async(f_t time_limit         = std::numeric_limits<f_t>::infinity(),
-                 double work_unit_limit = std::numeric_limits<double>::infinity());
+  void run_async(f_t time_limit                  = std::numeric_limits<f_t>::infinity(),
+                 double work_unit_limit          = std::numeric_limits<double>::infinity(),
+                 omp_atomic_t<i_t>* worker_count = nullptr);
 
   // Run the CPU FJ synchronously (i.e., wait for it to finish before proceeding)
   void run_sync(f_t time_limit         = std::numeric_limits<f_t>::infinity(),

cpp/src/mip_heuristics/presolve/third_party_presolve.cpp

@@ -1223,9 +1223,9 @@ void third_party_presolve_t<i_t, f_t>::undo(std::vector<f_t>& primal_solution,
 
 template <typename i_t, typename f_t>
 void third_party_presolve_t<i_t, f_t>::uncrush_primal_solution(
-  const std::vector<f_t>& reduced_primal, std::vector<f_t>& full_primal) const
+  const std::vector<f_t>& reduced_primal, std::vector<f_t>& full_primal, bool check_postsolve) const
 {
-  if (presolver_ == cuopt::mathematical_optimization::presolver_t::PSLP) {
+  if (presolver_ == PSLP) {
     cuopt_expects(false,
                   error_type_t::RuntimeError,
                   "This code path should be never called, as this is meant for callbacks and they "
@@ -1241,7 +1241,7 @@ void third_party_presolve_t<i_t, f_t>::uncrush_primal_solution(
 
   bool is_optimal = false;
   auto status = post_solver.undo(reduced_sol, full_sol, *papilo_post_solve_storage_, is_optimal);
-  check_postsolve_status(status);
+  if (check_postsolve) check_postsolve_status(status);
   full_primal = std::move(full_sol.primal);
 }
 

cpp/src/mip_heuristics/presolve/third_party_presolve.hpp

@@ -123,7 +123,8 @@ class third_party_presolve_t {
             bool dual_postsolve);
 
   void uncrush_primal_solution(const std::vector<f_t>& reduced_primal,
-                               std::vector<f_t>& full_primal) const;
+                               std::vector<f_t>& full_primal,
+                               bool check_postsolve = true) const;
 
   void crush_primal_solution(const optimization_problem_t<i_t, f_t>& reduced_problem,
                              const std::vector<f_t>& original_primal,

cpp/src/mip_heuristics/root_heuristics.hpp

@@ -0,0 +1,57 @@
+/* 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 <branch_and_bound/worker.hpp>
+#include <dual_simplex/user_problem.hpp>
+#include "feasibility_jump/fj_cpu_worker.cuh"
+
+namespace cuopt::mathematical_optimization::mip {
+
+template <typename i_t, typename f_t>
+struct root_heuristics_t {
+  std::unique_ptr<diving_worker_t<i_t, f_t>> submip_worker_;
+  fj_cpu_worker_t<i_t, f_t> fj_cpu_worker_;
+  std::vector<simplex::variable_type_t> var_types_;
+  csr_matrix_t<i_t, f_t> Arow_;
+
+  root_heuristics_t(const csr_matrix_t<i_t, f_t>& Arow,
+                    const std::vector<simplex::variable_type_t>& var_types)
+    : submip_worker_(nullptr), var_types_(var_types), Arow_(Arow) {};
+
+  ~root_heuristics_t() { stop(); }
+
+  void stop()
+  {
+    fj_cpu_worker_.stop();
+    if (submip_worker_) {
+      submip_worker_->halt              = true;
+      diving_worker_t<i_t, f_t>* worker = submip_worker_.get();
+#pragma omp taskwait depend(in : *worker)
+    }
+  }
+
+  void create_submip_worker(i_t id,
+                            const simplex::lp_problem_t<i_t, f_t>& lp,
+                            const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+                            f_t root_obj,
+                            const std::vector<simplex::variable_status_t>& root_vstatus,
+                            const std::vector<f_t>& sol)
+  {
+    submip_worker_ =
+      std::make_unique<diving_worker_t<i_t, f_t>>(id, lp, Arow_, var_types_, settings);
+    submip_worker_->start_node       = mip_node_t<i_t, f_t>(root_obj, root_vstatus);
+    submip_worker_->leaf_vstatus     = root_vstatus;
+    submip_worker_->leaf_solution.x  = sol;
+    submip_worker_->recompute_bounds = false;
+    submip_worker_->recompute_basis  = true;
+    submip_worker_->search_strategy  = search_strategy_t::SUBMIP;
+    submip_worker_->set_active();
+  }
+};
+}  // namespace cuopt::mathematical_optimization::mip