Root heuristics

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>
 
@@ -368,24 +370,31 @@ class branch_and_bound_t {
 
   // Launch a new RINS worker
   bool launch_rins_worker(const std::vector<f_t>& sol);
-  void set_solution_from_submip(const std::vector<f_t>& solution,
+  void set_solution_from_submip(const simplex::lp_problem_t<i_t, f_t>& lp,
+                                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,
+                    const std::vector<simplex::variable_type_t>& var_types,
                     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 = false);
 
   // 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,
+            const std::vector<simplex::variable_type_t>& var_types,
+            bool is_root_heuristic = false);
+
+  void launch_root_heuristics(const simplex::lp_problem_t<i_t, f_t>& lp,
+                              const std::vector<f_t>& sol,
+                              i_t cut_pass,
+                              root_heuristics_t<i_t, f_t>& root_heuristics);
 
   // 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/constants.hpp

@@ -27,15 +27,31 @@ enum class heuristics_origin_t {
 // doi: 10.1007/s10107-004-0518-7.
 enum class search_strategy_t : int {
   BEST_FIRST           = 0,  // Best-First + Plunging.
-  PSEUDOCOST_DIVING    = 1,  // Pseudocost diving (9.2.5)
-  LINE_SEARCH_DIVING   = 2,  // Line search diving (9.2.4)
-  GUIDED_DIVING        = 3,  // Guided diving (9.2.3).
-  COEFFICIENT_DIVING   = 4,  // Coefficient diving (9.2.1)
+  PSEUDOCOST_DIVING    = 1,  // Pseudocost diving [1, Section 9.2.5]
+  LINE_SEARCH_DIVING   = 2,  // Line search diving [1, Section 9.2.4]
+  GUIDED_DIVING        = 3,  // Guided diving. [1, Section 9.2.3]
+  COEFFICIENT_DIVING   = 4,  // Coefficient diving [1, Section 9.2.1]
   FARKAS_DIVING        = 5,  // Farkas Diving (see [2])
-  VECTOR_LENGTH_DIVING = 6,  // Vector Length Diving (9.2.6)
-  SUBMIP               = 7   // RINS (see [3])
+  VECTOR_LENGTH_DIVING = 6,  // Vector Length Diving [1, Section 9.2.6]
+  RINS                 = 7,  // RINS (see [3])
 };
 
 enum class branch_direction_t { NONE = -1, DOWN = 0, UP = 1 };
 
+inline const char* search_strategy_to_string(search_strategy_t search_strategy)
+{
+  switch (search_strategy) {
+    case search_strategy_t::BEST_FIRST: return "BEST_FIRST";
+    case search_strategy_t::PSEUDOCOST_DIVING: return "PSEUDOCOST_DIVING";
+    case search_strategy_t::LINE_SEARCH_DIVING: return "LINE_SEARCH_DIVING";
+    case search_strategy_t::GUIDED_DIVING: return "GUIDED_DIVING";
+    case search_strategy_t::COEFFICIENT_DIVING: return "COEFFICIENT_DIVING";
+    case search_strategy_t::FARKAS_DIVING: return "FARKAS_DIVING";
+    case search_strategy_t::VECTOR_LENGTH_DIVING: return "VECTOR_LENGTH_DIVING";
+    case search_strategy_t::RINS: return "RINS";
+  }
+
+  return "UNKNOWN";
+}
+
 }  // namespace cuopt::mathematical_optimization::mip

cpp/src/branch_and_bound/deterministic_workers.hpp

@@ -89,8 +89,10 @@ class deterministic_worker_base_t : public branch_and_bound_worker_t<i_t, f_t> {
                               const csr_matrix_t<i_t, f_t>& Arow,
                               const std::vector<simplex::variable_type_t>& var_types,
                               const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+                              const std::vector<f_t>& root_solution,
+                              const std::vector<f_t>& root_edge_norm,
                               const std::string& context_name)
-    : base_t(id, original_lp, Arow, var_types, settings),
+    : base_t(id, original_lp, Arow, var_types, settings, root_solution, root_edge_norm),
       work_context(context_name),
       pc_snapshot(1, settings)
   {
@@ -140,8 +142,17 @@ class deterministic_bfs_worker_t
                                       const simplex::lp_problem_t<i_t, f_t>& original_lp,
                                       const csr_matrix_t<i_t, f_t>& Arow,
                                       const std::vector<simplex::variable_type_t>& var_types,
-                                      const simplex::simplex_solver_settings_t<i_t, f_t>& settings)
-    : base_t(id, original_lp, Arow, var_types, settings, "BB_Worker_" + std::to_string(id))
+                                      const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+                                      const std::vector<f_t>& root_solution,
+                                      const std::vector<f_t>& root_edge_norm)
+    : base_t(id,
+             original_lp,
+             Arow,
+             var_types,
+             settings,
+             root_solution,
+             root_edge_norm,
+             "BB_Worker_" + std::to_string(id))
   {
   }
 
@@ -282,9 +293,6 @@ class deterministic_diving_worker_t
   std::vector<f_t> dive_lower;
   std::vector<f_t> dive_upper;
 
-  // Root LP relaxation solution (constant, set once at construction)
-  const std::vector<f_t>* root_solution{nullptr};
-
   // Diving state
   bool recompute_bounds_and_basis{true};
 
@@ -300,10 +308,17 @@ class deterministic_diving_worker_t
     const csr_matrix_t<i_t, f_t>& Arow,
     const std::vector<simplex::variable_type_t>& var_types,
     const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-    const std::vector<f_t>* root_sol)
-    : base_t(id, original_lp, Arow, var_types, settings, "Diving_Worker_" + std::to_string(id)),
-      diving_type(type),
-      root_solution(root_sol)
+    const std::vector<f_t>& root_solution,
+    const std::vector<f_t>& root_edge_norm)
+    : base_t(id,
+             original_lp,
+             Arow,
+             var_types,
+             settings,
+             root_solution,
+             root_edge_norm,
+             "Diving_Worker_" + std::to_string(id)),
+      diving_type(type)
   {
     dive_lower = original_lp.lower;
     dive_upper = original_lp.upper;
@@ -407,11 +422,14 @@ class deterministic_bfs_worker_pool_t
                                   const simplex::lp_problem_t<i_t, f_t>& original_lp,
                                   const csr_matrix_t<i_t, f_t>& Arow,
                                   const std::vector<simplex::variable_type_t>& var_types,
-                                  const simplex::simplex_solver_settings_t<i_t, f_t>& settings)
+                                  const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+                                  const std::vector<f_t>& root_solution,
+                                  const std::vector<f_t>& root_edge_norm)
   {
     this->workers_.reserve(num_workers);
     for (int i = 0; i < num_workers; ++i) {
-      this->workers_.emplace_back(i, original_lp, Arow, var_types, settings);
+      this->workers_.emplace_back(
+        i, original_lp, Arow, var_types, settings, root_solution, root_edge_norm);
     }
   }
 
@@ -443,12 +461,14 @@ class deterministic_diving_worker_pool_t
                                      const csr_matrix_t<i_t, f_t>& Arow,
                                      const std::vector<simplex::variable_type_t>& var_types,
                                      const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
-                                     const std::vector<f_t>* root_solution)
+                                     const std::vector<f_t>& root_solution,
+                                     const std::vector<f_t>& root_edge_norm)
   {
     this->workers_.reserve(num_workers);
     for (int i = 0; i < num_workers; ++i) {
       search_strategy_t type = diving_types[i % diving_types.size()];
-      this->workers_.emplace_back(i, type, original_lp, Arow, var_types, settings, root_solution);
+      this->workers_.emplace_back(
+        i, type, original_lp, Arow, var_types, settings, root_solution, root_edge_norm);
     }
   }
 

cpp/src/branch_and_bound/worker.hpp

@@ -78,6 +78,9 @@ class branch_and_bound_worker_t {
   bool recompute_basis  = true;
   bool recompute_bounds = true;
 
+  const std::vector<f_t>& root_solution;
+  const std::vector<f_t>& root_edge_norm;
+
   void ensure_orbital_fixing()
   {
     if (orbital_fixing == nullptr && symmetry_ptr != nullptr) {
@@ -94,6 +97,8 @@ class branch_and_bound_worker_t {
                             const csr_matrix_t<i_t, f_t>& Arow,
                             const std::vector<simplex::variable_type_t>& var_type,
                             const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+                            const std::vector<f_t>& root_solution,
+                            const std::vector<f_t>& root_edge_norm,
                             uint64_t rng_offset = 0)
     : worker_id(worker_id),
       search_strategy(search_strategy_t::BEST_FIRST),
@@ -108,7 +113,9 @@ class branch_and_bound_worker_t {
       node_presolver(leaf_problem, Arow, {}, var_type),
       bounds_changed(original_lp.num_cols, false),
       rng(settings.random_seed + pcgenerator_t::default_seed + rng_offset + worker_id,
-          pcgenerator_t::default_stream ^ (worker_id + rng_offset))
+          pcgenerator_t::default_stream ^ (worker_id + rng_offset)),
+      root_solution(root_solution),
+      root_edge_norm(root_edge_norm)
   {
   }
 
@@ -146,8 +153,11 @@ class bfs_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
                const csr_matrix_t<i_t, f_t>& Arow,
                const std::vector<simplex::variable_type_t>& var_type,
                const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+               const std::vector<f_t>& root_solution,
+               const std::vector<f_t>& root_edge_norm,
                uint64_t rng_offset = 0)
-    : Base(worker_id, original_lp, Arow, var_type, settings, rng_offset)
+    : Base(
+        worker_id, original_lp, Arow, var_type, settings, root_solution, root_edge_norm, rng_offset)
   {
     this->start_lower     = original_lp.lower;
     this->start_upper     = original_lp.upper;
@@ -243,6 +253,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/branch_and_bound/worker_pool.hpp

@@ -24,6 +24,8 @@ class worker_pool_t {
             const std::vector<simplex::variable_type_t>& var_type,
             mip_symmetry_t<i_t, f_t>* symmetry,
             const simplex::simplex_solver_settings_t<i_t, f_t>& settings,
+            const std::vector<f_t>& root_solution,
+            const std::vector<f_t>& root_edge_norm,
             const uint64_t rng_offset = 0)
   {
     assert(!is_initialized_);
@@ -33,8 +35,8 @@ class worker_pool_t {
     num_idle_workers_ = num_workers;
     idle_workers_.clear_resize(num_workers);
     for (i_t i = 0; i < num_workers; ++i) {
-      workers_[i] =
-        std::make_unique<WorkerType>(i, original_lp, Arow, var_type, settings, rng_offset);
+      workers_[i] = std::make_unique<WorkerType>(
+        i, original_lp, Arow, var_type, settings, root_solution, root_edge_norm, rng_offset);
       idle_workers_.push_back(i);
       // Propagate the (possibly null) symmetry pointer; workers lazily build
       // their orbital_fixing/lexical_reduction state via ensure_orbital_fixing().
@@ -61,14 +63,17 @@ class worker_pool_t {
 
   void return_worker_to_pool(WorkerType* worker)
   {
-    std::lock_guard lock(mutex_);
     assert(worker != nullptr);
+    worker->set_inactive();
+    assert(!worker->is_active.load());
+
+    if (!is_initialized_) return;
+
+    std::lock_guard lock(mutex_);
     assert(workers_[worker->worker_id].get() == worker);
     assert(static_cast<size_t>(num_idle_workers_.load()) == idle_workers_.size());
     assert(idle_workers_.size() <= workers_.size());
 
-    worker->set_inactive();
-    assert(!worker->is_active.load());
     idle_workers_.push_back(worker->worker_id);
     num_idle_workers_++;
   }

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_presolved_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_presolved_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_presolved_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/fj_cpu.cu

@@ -2109,37 +2109,52 @@ void fj_cpu_worker_t<i_t, f_t>::create_worker(
   fj_cpu.reset(new_climber.release());
   fj_cpu->log_prefix           = std::move(log_prefix);
   fj_cpu->improvement_callback = improvement_callback;
+  fj_cpu->halted               = false;
+  preemption_flag              = false;
+  is_initialized               = true;
 }
 
 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)
 {
-  if (!fj_cpu) return;
+  if (!is_initialized) return;
 
-#pragma omp task shared(fj_cpu) firstprivate(time_limit, work_unit_limit) \
-  priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : *fj_cpu)
-  cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit);
+  auto& fj_ptr = fj_cpu;
+#pragma omp task shared(fj_cpu, is_initialized, fj_ptr) firstprivate(time_limit, work_unit_limit) \
+  priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) depend(out : fj_ptr)
+  {
+    if (is_initialized) { cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit); }
+  }
 }
 
 template <typename i_t, typename f_t>
 void fj_cpu_worker_t<i_t, f_t>::run_sync(f_t time_limit, double work_unit_limit)
 {
-  if (!fj_cpu) return;
+  if (!is_initialized) return;
   cpufj_solve(fj_cpu.get(), time_limit, work_unit_limit);
+  is_initialized = false;
   fj_cpu.reset();
 }
 
 template <typename i_t, typename f_t>
 void fj_cpu_worker_t<i_t, f_t>::stop()
 {
-  if (!fj_cpu) return;
+  if (!is_initialized) return;
+
+  preemption_flag = true;
 
-  fj_cpu->preemption_flag = true;
-  fj_cpu->halted          = true;
-#pragma omp taskwait depend(in : *fj_cpu)
+  auto& fj_ptr = fj_cpu;
+#pragma omp taskwait depend(in : fj_ptr)
+  is_initialized = false;
   fj_cpu.reset();
 }
 
+template <typename i_t, typename f_t>
+void fj_cpu_worker_t<i_t, f_t>::send_stop_signal()
+{
+  preemption_flag = true;
+}
+
 #if MIP_INSTANTIATE_FLOAT
 template class fj_t<int, float>;
 template struct fj_cpu_worker_t<int, float>;