RINS heuristic

cpp/src/dual_simplex/branch_and_bound.cpp

@@ -279,6 +279,7 @@ void branch_and_bound_t<i_t, f_t>::set_new_solution(const std::vector<f_t>& solu
       original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional);
     if (is_feasible) {
       upper_bound_ = obj;
+      incumbent_.set_incumbent_solution(obj, crushed_solution);
     } else {
       attempt_repair         = true;
       constexpr bool verbose = false;
@@ -448,16 +449,25 @@ mip_status_t branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t
     settings_.log.printf("Time limit reached. Stopping the solver...\n");
     mip_status = mip_status_t::TIME_LIMIT;
   }
+  if (status_ == mip_exploration_status_t::NODE_LIMIT) {
+    settings_.log.printf("Node limit reached. Stopping the solver...\n");
+    mip_status = mip_status_t::NODE_LIMIT;
+  }
 
-  f_t upper_bound = get_upper_bound();
-  f_t gap         = upper_bound - lower_bound;
-  f_t obj         = compute_user_objective(original_lp_, upper_bound);
-  f_t user_lower  = compute_user_objective(original_lp_, lower_bound);
-  f_t gap_rel     = user_relative_gap(original_lp_, upper_bound, lower_bound);
+  f_t upper_bound      = get_upper_bound();
+  f_t gap              = upper_bound - lower_bound;
+  f_t obj              = compute_user_objective(original_lp_, upper_bound);
+  f_t user_bound       = compute_user_objective(original_lp_, lower_bound);
+  f_t gap_rel          = user_relative_gap(original_lp_, upper_bound, lower_bound);
+  bool is_maximization = original_lp_.obj_scale < 0.0;
 
   settings_.log.printf(
     "Explored %d nodes in %.2fs.\n", stats_.nodes_explored, toc(stats_.start_time));
-  settings_.log.printf("Absolute Gap %e Objective %.16e Lower Bound %.16e\n", gap, obj, user_lower);
+  settings_.log.printf("Absolute Gap %e Objective %.16e %s Bound %.16e\n",
+                       gap,
+                       obj,
+                       is_maximization ? "Upper" : "Lower",
+                       user_bound);
 
   if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) {
     mip_status = mip_status_t::OPTIMAL;
@@ -483,7 +493,10 @@ mip_status_t branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t
     }
   }
 
-  uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x);
+  if (upper_bound != inf) {
+    assert(incumbent_.has_incumbent);
+    uncrush_primal_solution(original_problem_, original_lp_, incumbent_.x, solution.x);
+  }
   solution.objective          = incumbent_.objective;
   solution.lower_bound        = lower_bound;
   solution.nodes_explored     = stats_.nodes_explored;
@@ -637,6 +650,12 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node(search_tree_t<i_t, f_t>&
     search_tree.graphviz_node(log, node_ptr, "lower bound", leaf_objective);
     pc_.update_pseudo_costs(node_ptr, leaf_objective);
 
+    if (settings_.node_processed_callback != nullptr) {
+      std::vector<f_t> original_x;
+      uncrush_primal_solution(original_problem_, original_lp_, leaf_solution.x, original_x);
+      settings_.node_processed_callback(original_x, leaf_objective);
+    }
+
     if (leaf_num_fractional == 0) {
       // Found a integer feasible solution
       add_feasible_solution(leaf_objective, leaf_solution.x, node_ptr->depth, thread_type);
@@ -845,6 +864,10 @@ void branch_and_bound_t<i_t, f_t>::explore_subtree(i_t task_id,
       status_ = mip_exploration_status_t::TIME_LIMIT;
       return;
     }
+    if (stats_.nodes_explored >= settings_.node_limit) {
+      status_ = mip_exploration_status_t::NODE_LIMIT;
+      return;
+    }
 
     // Set the correct bounds for the leaf problem
     leaf_problem.lower = original_lp_.lower;
@@ -1018,13 +1041,15 @@ void branch_and_bound_t<i_t, f_t>::diving_thread(lp_problem_t<i_t, f_t>& leaf_pr
 }
 
 template <typename i_t, typename f_t>
-mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solution)
+mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solution,
+                                                 std::string log_prefix)
 {
   logger_t log;
-  log.log                 = false;
-  status_                 = mip_exploration_status_t::UNSET;
-  stats_.nodes_unexplored = 0;
-  stats_.nodes_explored   = 0;
+  log.log        = false;
+  log.log_prefix = settings_.log.log_prefix = log_prefix;
+  status_                                   = mip_exploration_status_t::UNSET;
+  stats_.nodes_unexplored                   = 0;
+  stats_.nodes_explored                     = 0;
 
   if (guess_.size() != 0) {
     std::vector<f_t> crushed_guess;

cpp/src/dual_simplex/branch_and_bound.hpp

@@ -142,7 +142,7 @@ class branch_and_bound_t {
   i_t get_heap_size();
 
   // The main entry routine. Returns the solver status and populates solution with the incumbent.
-  mip_status_t solve(mip_solution_t<i_t, f_t>& solution);
+  mip_status_t solve(mip_solution_t<i_t, f_t>& solution, std::string log_prefix = "");
 
  private:
   const user_problem_t<i_t, f_t>& original_problem_;

cpp/src/dual_simplex/logger.hpp

@@ -71,7 +71,7 @@ class logger_t {
 
         size_t len = strlen(buffer);
         if (len > 0 && buffer[len - 1] == '\n') { buffer[len - 1] = '\0'; }
-        CUOPT_LOG_INFO(buffer);
+        CUOPT_LOG_INFO("%s%s", log_prefix.c_str(), buffer);
       }
 #else
       if (log_to_console) {
@@ -105,7 +105,7 @@ class logger_t {
 
         size_t len = strlen(buffer);
         if (len > 0 && buffer[len - 1] == '\n') { buffer[len - 1] = '\0'; }
-        CUOPT_LOG_TRACE(buffer);
+        CUOPT_LOG_TRACE("%s%s", log_prefix.c_str(), buffer);
       }
 #else
       if (log_to_console) {
@@ -128,6 +128,7 @@ class logger_t {
 
   bool log;
   bool log_to_console;
+  std::string log_prefix;
 
  private:
   bool log_to_file;

cpp/src/dual_simplex/presolve.cpp

@@ -1381,7 +1381,10 @@ void uncrush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
 {
   user_solution.resize(user_problem.num_cols);
   assert(problem.num_cols >= user_problem.num_cols);
-  std::copy(solution.begin(), solution.begin() + user_problem.num_cols, user_solution.data());
+  assert(solution.size() >= user_problem.num_cols);
+  std::copy(solution.begin(),
+            solution.begin() + std::min((i_t)solution.size(), user_problem.num_cols),
+            user_solution.data());
 }
 
 template <typename i_t, typename f_t>

cpp/src/dual_simplex/simplex_solver_settings.hpp

@@ -33,6 +33,7 @@ struct simplex_solver_settings_t {
  public:
   simplex_solver_settings_t()
     : iteration_limit(std::numeric_limits<i_t>::max()),
+      node_limit(std::numeric_limits<i_t>::max()),
       time_limit(std::numeric_limits<f_t>::infinity()),
       absolute_mip_gap_tol(0.0),
       relative_mip_gap_tol(1e-3),
@@ -91,6 +92,7 @@ struct simplex_solver_settings_t {
   void set_log_filename(const std::string& log_filename) { log.set_log_file(log_filename); }
   void close_log_file() { log.close_log_file(); }
   i_t iteration_limit;
+  i_t node_limit;
   f_t time_limit;
   f_t absolute_mip_gap_tol;  // Tolerance on mip gap to declare optimal
   f_t relative_mip_gap_tol;  // Tolerance on mip gap to declare optimal
@@ -144,6 +146,7 @@ struct simplex_solver_settings_t {
   i_t num_diving_threads;          // number of threads dedicated to diving
   i_t inside_mip;  // 0 if outside MIP, 1 if inside MIP at root node, 2 if inside MIP at leaf node
   std::function<void(std::vector<f_t>&, f_t)> solution_callback;
+  std::function<void(const std::vector<f_t>&, f_t)> node_processed_callback;
   std::function<void()> heuristic_preemption_callback;
   std::function<void(std::vector<f_t>&, std::vector<f_t>&, f_t)> set_simplex_solution_callback;
   mutable logger_t log;

cpp/src/dual_simplex/solution.hpp

@@ -61,16 +61,21 @@ class lp_solution_t {
 template <typename i_t, typename f_t>
 class mip_solution_t {
  public:
-  mip_solution_t(i_t n) : x(n), objective(std::numeric_limits<f_t>::quiet_NaN()), lower_bound(-inf)
+  mip_solution_t(i_t n)
+    : x(n),
+      objective(std::numeric_limits<f_t>::quiet_NaN()),
+      lower_bound(-inf),
+      has_incumbent(false)
   {
   }
 
   void resize(i_t n) { x.resize(n); }
 
   void set_incumbent_solution(f_t primal_objective, const std::vector<f_t>& primal_solution)
   {
-    x         = primal_solution;
-    objective = primal_objective;
+    x             = primal_solution;
+    objective     = primal_objective;
+    has_incumbent = true;
   }
 
   // Primal solution vector
@@ -79,6 +84,7 @@ class mip_solution_t {
   f_t lower_bound;
   i_t nodes_explored;
   i_t simplex_iterations;
+  bool has_incumbent;
 };
 
 }  // namespace cuopt::linear_programming::dual_simplex

cpp/src/mip/CMakeLists.txt

@@ -33,6 +33,7 @@ set(MIP_NON_LP_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/diversity/diversity_manager.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/diversity/multi_armed_bandit.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/diversity/population.cu
+  ${CMAKE_CURRENT_SOURCE_DIR}/diversity/lns/rins.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/relaxed_lp/relaxed_lp.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/local_search/local_search.cu
   ${CMAKE_CURRENT_SOURCE_DIR}/local_search/rounding/bounds_repair.cu

cpp/src/mip/diversity/diversity_manager.cu

@@ -57,6 +57,7 @@ diversity_manager_t<i_t, f_t>::diversity_manager_t(mip_solver_context_t<i_t, f_t
     lp_dual_optimal_solution(context.problem_ptr->n_constraints,
                              context.problem_ptr->handle_ptr->get_stream()),
     ls(context, lp_optimal_solution),
+    rins(context, *this),
     timer(diversity_config.default_time_limit),
     bound_prop_recombiner(context,
                           context.problem_ptr->n_variables,
@@ -417,6 +418,8 @@ solution_t<i_t, f_t> diversity_manager_t<i_t, f_t>::run_solver()
     run_fj_alone(sol);
     return sol;
   }
+  rins.enable();
+
   generate_solution(timer.remaining_time(), false);
   if (timer.check_time_limit()) {
     population.add_external_solutions_to_population();
@@ -731,7 +734,7 @@ void diversity_manager_t<i_t, f_t>::set_simplex_solution(const std::vector<f_t>&
 {
   CUOPT_LOG_DEBUG("Setting simplex solution with objective %f", objective);
   using sol_t = solution_t<i_t, f_t>;
-  cudaSetDevice(context.handle_ptr->get_device());
+  RAFT_CUDA_TRY(cudaSetDevice(context.handle_ptr->get_device()));
   context.handle_ptr->sync_stream();
   cuopt_func_call(sol_t new_sol(*problem_ptr));
   cuopt_assert(new_sol.assignment.size() == solution.size(), "Assignment size mismatch");

cpp/src/mip/diversity/diversity_manager.cuh

@@ -31,6 +31,7 @@
 #include <cuopt/linear_programming/mip/solver_settings.hpp>
 #include <cuopt/linear_programming/mip/solver_stats.hpp>
 
+#include <mip/diversity/lns/rins.cuh>
 #include <mip/local_search/local_search.cuh>
 #include <mip/solution/solution.cuh>
 #include <mip/solver.cuh>
@@ -105,6 +106,8 @@ class diversity_manager_t {
   // atomic for signalling pdlp to stop
   volatile int global_concurrent_halt{0};
 
+  rins_t<i_t, f_t> rins;
+
   bool run_only_ls_recombiner{false};
   bool run_only_bp_recombiner{false};
   bool run_only_fp_recombiner{false};