Add Implied bound cuts. Find more knapsack cuts. Lift knapsack cuts

cpp/include/cuopt/linear_programming/constants.h

@@ -64,6 +64,7 @@
 #define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts"
 #define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS   "mip_mixed_integer_gomory_cuts"
 #define CUOPT_MIP_KNAPSACK_CUTS               "mip_knapsack_cuts"
+#define CUOPT_MIP_IMPLIED_BOUND_CUTS          "mip_implied_bound_cuts"
 #define CUOPT_MIP_CLIQUE_CUTS                 "mip_clique_cuts"
 #define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS  "mip_strong_chvatal_gomory_cuts"
 #define CUOPT_MIP_REDUCED_COST_STRENGTHENING  "mip_reduced_cost_strengthening"

cpp/include/cuopt/linear_programming/mip/solver_settings.hpp

@@ -94,6 +94,7 @@ class mip_solver_settings_t {
   i_t mixed_integer_gomory_cuts = -1;
   i_t knapsack_cuts             = -1;
   i_t clique_cuts               = -1;
+  i_t implied_bound_cuts        = -1;
   i_t strong_chvatal_gomory_cuts      = -1;
   i_t reduced_cost_strengthening      = -1;
   f_t cut_change_threshold            = -1.0;

cpp/src/branch_and_bound/branch_and_bound.cpp

@@ -205,7 +205,11 @@ std::string user_mip_gap(f_t obj_value, f_t lower_bound)
   } else {
     constexpr int BUFFER_LEN = 32;
     char buffer[BUFFER_LEN];
-    snprintf(buffer, BUFFER_LEN - 1, "%5.1f%%", user_mip_gap * 100);
+    if (user_mip_gap > 1e-3) {
+      snprintf(buffer, BUFFER_LEN - 1, "%5.1f%%", user_mip_gap * 100);
+    } else {
+      snprintf(buffer, BUFFER_LEN - 1, "%5.2f%%", user_mip_gap * 100);
+    }
     return std::string(buffer);
   }
 }
@@ -243,9 +247,11 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
   const user_problem_t<i_t, f_t>& user_problem,
   const simplex_solver_settings_t<i_t, f_t>& solver_settings,
   f_t start_time,
+  const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
   std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table)
   : original_problem_(user_problem),
     settings_(solver_settings),
+    probing_implied_bound_(probing_implied_bound),
     clique_table_(std::move(clique_table)),
     original_lp_(user_problem.handle_ptr, 1, 1, 1),
     Arow_(1, 1, 0),
@@ -556,7 +562,7 @@ template <typename i_t, typename f_t>
 bool branch_and_bound_t<i_t, f_t>::repair_solution(const std::vector<f_t>& edge_norms,
                                                    const std::vector<f_t>& potential_solution,
                                                    f_t& repaired_obj,
-                                                   std::vector<f_t>& repaired_solution) const
+                                                   std::vector<f_t>& repaired_solution)
 {
   bool feasible = false;
   repaired_obj  = std::numeric_limits<f_t>::quiet_NaN();
@@ -579,9 +585,10 @@ bool branch_and_bound_t<i_t, f_t>::repair_solution(const std::vector<f_t>& edge_
   i_t iter                               = 0;
   f_t lp_start_time                      = tic();
   simplex_solver_settings_t lp_settings  = settings_;
+  lp_settings.concurrent_halt            = &node_concurrent_halt_;
   std::vector<variable_status_t> vstatus = root_vstatus_;
   lp_settings.set_log(false);
-  lp_settings.inside_mip           = true;
+  lp_settings.inside_mip           = 2;
   std::vector<f_t> leaf_edge_norms = edge_norms;
   // should probably set the cut off here lp_settings.cut_off
   dual::status_t lp_status = dual_phase2(
@@ -724,12 +731,6 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
                        obj,
                        is_maximization ? "Upper" : "Lower",
                        user_bound);
-  {
-    const f_t root_lp_obj = root_lp_current_lower_bound_.load();
-    if (std::isfinite(root_lp_obj)) {
-      settings_.log.printf("Root LP dual objective (last): %.16e\n", root_lp_obj);
-    }
-  }
 
   if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) {
     solver_status_ = mip_status_t::OPTIMAL;
@@ -1336,6 +1337,7 @@ dual::status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
   assert(leaf_vstatus.size() == worker->leaf_problem.num_cols);
 
   simplex_solver_settings_t lp_settings = settings_;
+  lp_settings.concurrent_halt           = &node_concurrent_halt_;
   lp_settings.set_log(false);
   f_t cutoff = get_cutoff();
   if (original_lp_.objective_is_integral) {
@@ -1432,23 +1434,25 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
   worker->recompute_basis  = true;
   worker->recompute_bounds = true;
 
-  while (stack.size() > 0 && solver_status_ == mip_status_t::UNSET) {
+  f_t lower_bound = get_lower_bound();
+  f_t upper_bound = upper_bound_;
+  f_t rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
+  f_t abs_gap     = upper_bound - lower_bound;
+
+  while (stack.size() > 0 && (solver_status_ == mip_status_t::UNSET && is_running_) &&
+         rel_gap > settings_.relative_mip_gap_tol && abs_gap > settings_.absolute_mip_gap_tol) {
     mip_node_t<i_t, f_t>* node_ptr = stack.front();
     stack.pop_front();
 
-    f_t lower_bound = node_ptr->lower_bound;
-    f_t upper_bound = upper_bound_;
-    f_t rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
-
     // This is based on three assumptions:
     // - The stack only contains sibling nodes, i.e., the current node and it sibling (if
     // applicable)
     // - The current node and its siblings uses the lower bound of the parent before solving the LP
     // relaxation
     // - The lower bound of the parent is lower or equal to its children
-    worker->lower_bound = lower_bound;
+    worker->lower_bound = node_ptr->lower_bound;
 
-    if (lower_bound > get_cutoff()) {
+    if (node_ptr->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;
@@ -1472,6 +1476,9 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
     if (lp_status == dual::status_t::TIME_LIMIT) {
       solver_status_ = mip_status_t::TIME_LIMIT;
       break;
+    } else if (lp_status == dual::status_t::CONCURRENT_LIMIT) {
+      stack.push_front(node_ptr);
+      break;
     } else if (lp_status == dual::status_t::ITERATION_LIMIT) {
       break;
     }
@@ -1517,6 +1524,28 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(branch_and_bound_worker_t<i_t, f_
         stack.push_front(node_ptr->get_down_child());
       }
     }
+
+    lower_bound = get_lower_bound();
+    upper_bound = upper_bound_;
+    rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
+    abs_gap     = upper_bound - lower_bound;
+  }
+
+  lower_bound = get_lower_bound();
+  upper_bound = upper_bound_;
+  rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
+  abs_gap     = upper_bound - lower_bound;
+
+  if (stack.size() > 0 &&
+      (rel_gap <= settings_.relative_mip_gap_tol || abs_gap <= settings_.absolute_mip_gap_tol)) {
+    // If the solver converged according to the gap rules, but we still have nodes to explore
+    // in the stack, then we should add all the pending nodes back to the heap so the lower
+    // bound of the solver is set to the correct value.
+    while (!stack.empty()) {
+      auto node = stack.front();
+      stack.pop_front();
+      node_queue_.push(node);
+    }
   }
 
   if (settings_.num_threads > 1) {
@@ -1549,14 +1578,17 @@ void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>
   dive_stats.nodes_explored      = 0;
   dive_stats.nodes_unexplored    = 1;
 
-  while (stack.size() > 0 && solver_status_ == mip_status_t::UNSET && is_running_) {
+  f_t lower_bound = get_lower_bound();
+  f_t upper_bound = upper_bound_;
+  f_t rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
+  f_t abs_gap     = upper_bound - lower_bound;
+
+  while (stack.size() > 0 && (solver_status_ == mip_status_t::UNSET && is_running_) &&
+         rel_gap > settings_.relative_mip_gap_tol && abs_gap > settings_.absolute_mip_gap_tol) {
     mip_node_t<i_t, f_t>* node_ptr = stack.front();
     stack.pop_front();
 
-    f_t lower_bound     = node_ptr->lower_bound;
-    f_t upper_bound     = upper_bound_;
-    f_t rel_gap         = user_relative_gap(original_lp_, upper_bound, lower_bound);
-    worker->lower_bound = lower_bound;
+    worker->lower_bound = node_ptr->lower_bound;
 
     if (node_ptr->lower_bound > get_cutoff()) {
       worker->recompute_basis  = true;
@@ -1572,6 +1604,8 @@ void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>
     if (lp_status == dual::status_t::TIME_LIMIT) {
       solver_status_ = mip_status_t::TIME_LIMIT;
       break;
+    } else if (lp_status == dual::status_t::CONCURRENT_LIMIT) {
+      break;
     } else if (lp_status == dual::status_t::ITERATION_LIMIT) {
       break;
     }
@@ -1598,6 +1632,11 @@ void branch_and_bound_t<i_t, f_t>::dive_with(branch_and_bound_worker_t<i_t, f_t>
     if (stack.size() > 1 && stack.front()->depth - stack.back()->depth > diving_backtrack_limit) {
       stack.pop_back();
     }
+
+    lower_bound = get_lower_bound();
+    upper_bound = upper_bound_;
+    rel_gap     = user_relative_gap(original_lp_, upper_bound, lower_bound);
+    abs_gap     = upper_bound - lower_bound;
   }
 
   worker_pool_.return_worker_to_pool(worker);
@@ -1637,9 +1676,6 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
          rel_gap > settings_.relative_mip_gap_tol &&
          (active_workers_per_strategy_[0] > 0 || node_queue_.best_first_queue_size() > 0)) {
     bool launched_any_task = false;
-    lower_bound            = get_lower_bound();
-    abs_gap                = upper_bound_ - lower_bound;
-    rel_gap                = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
 
     repair_heuristic_solutions();
 
@@ -1740,6 +1776,16 @@ void branch_and_bound_t<i_t, f_t>::run_scheduler()
       }
     }
 
+    lower_bound = get_lower_bound();
+    abs_gap     = upper_bound_ - lower_bound;
+    rel_gap     = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
+
+    if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) {
+      node_concurrent_halt_ = 1;
+      solver_status_        = mip_status_t::OPTIMAL;
+      break;
+    }
+
     // If no new task was launched in this iteration, suspend temporarily the
     // execution of the scheduler. As of 8/Jan/2026, GCC does not
     // implement taskyield, but LLVM does.
@@ -1759,10 +1805,6 @@ void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
   while (solver_status_ == mip_status_t::UNSET && abs_gap > settings_.absolute_mip_gap_tol &&
          rel_gap > settings_.relative_mip_gap_tol && node_queue_.best_first_queue_size() > 0) {
     bool launched_any_task = false;
-    lower_bound            = get_lower_bound();
-    abs_gap                = upper_bound_ - lower_bound;
-    rel_gap                = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
-
     repair_heuristic_solutions();
 
     f_t now = toc(exploration_stats_.start_time);
@@ -1800,6 +1842,15 @@ void branch_and_bound_t<i_t, f_t>::single_threaded_solve()
 
     worker.init_best_first(start_node.value(), original_lp_);
     plunge_with(&worker);
+
+    lower_bound = get_lower_bound();
+    abs_gap     = upper_bound_ - lower_bound;
+    rel_gap     = user_relative_gap(original_lp_, upper_bound_.load(), lower_bound);
+
+    if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) {
+      solver_status_ = mip_status_t::OPTIMAL;
+      break;
+    }
   }
 }
 
@@ -2151,6 +2202,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
                                             new_slacks_,
                                             var_types_,
                                             original_problem_,
+                                            probing_implied_bound_,
                                             clique_table_,
                                             &clique_table_future_,
                                             &signal_extend_cliques_);
@@ -2529,6 +2581,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   node_queue_.push(search_tree_.root.get_up_child());
 
   settings_.log.printf("Exploring the B&B tree using %d threads\n\n", settings_.num_threads);
+  node_concurrent_halt_ = 0;
 
   exploration_stats_.nodes_explored       = 0;
   exploration_stats_.nodes_unexplored     = 2;

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -77,6 +77,7 @@ class branch_and_bound_t {
   branch_and_bound_t(const user_problem_t<i_t, f_t>& user_problem,
                      const simplex_solver_settings_t<i_t, f_t>& solver_settings,
                      f_t start_time,
+                     const probing_implied_bound_t<i_t, f_t>& probing_implied_bound,
                      std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table = nullptr);
 
   // Set an initial guess based on the user_problem. This should be called before solve.
@@ -133,7 +134,7 @@ class branch_and_bound_t {
   bool repair_solution(const std::vector<f_t>& leaf_edge_norms,
                        const std::vector<f_t>& potential_solution,
                        f_t& repaired_obj,
-                       std::vector<f_t>& repaired_solution) const;
+                       std::vector<f_t>& repaired_solution);
 
   f_t get_lower_bound();
   bool enable_concurrent_lp_root_solve() const { return enable_concurrent_lp_root_solve_; }
@@ -162,6 +163,7 @@ class branch_and_bound_t {
  private:
   const user_problem_t<i_t, f_t>& original_problem_;
   const simplex_solver_settings_t<i_t, f_t> settings_;
+  const probing_implied_bound_t<i_t, f_t>& probing_implied_bound_;
   std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table_;
   std::future<std::shared_ptr<detail::clique_table_t<i_t, f_t>>> clique_table_future_;
   std::atomic<bool> signal_extend_cliques_{false};
@@ -222,6 +224,7 @@ class branch_and_bound_t {
   omp_atomic_t<bool> solving_root_relaxation_{false};
   bool enable_concurrent_lp_root_solve_{false};
   std::atomic<int> root_concurrent_halt_{0};
+  std::atomic<int> node_concurrent_halt_{0};
   bool is_root_solution_set{false};
 
   // Pseudocosts