Diving on cut passes

cpp/include/cuopt/mathematical_optimization/constants.h

@@ -131,7 +131,6 @@
 #define CUOPT_MIP_HYPER_DIVING_FARKAS        "mip_hyper_diving_farkas"
 #define CUOPT_MIP_HYPER_DIVING_VECTOR_LENGTH "mip_hyper_diving_vector_length"
 /* @brief Diving heuristic limits */
-#define CUOPT_MIP_HYPER_DIVING_MIN_NODE_DEPTH         "mip_hyper_diving_min_node_depth"
 #define CUOPT_MIP_HYPER_DIVING_NODE_LIMIT             "mip_hyper_diving_node_limit"
 #define CUOPT_MIP_HYPER_DIVING_ITERATION_LIMIT_FACTOR "mip_hyper_diving_iteration_limit_factor"
 #define CUOPT_MIP_HYPER_DIVING_BACKTRACK_LIMIT        "mip_hyper_diving_backtrack_limit"

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

@@ -6,6 +6,7 @@
 /* clang-format on */
 
 #pragma once
+#include <cstdint>
 
 namespace cuopt::mathematical_optimization {
 
@@ -22,16 +23,13 @@ struct mip_diving_hyper_params_t {
   i_t farkas_diving        = -1;
   i_t vector_length_diving = -1;
 
-  // The minimum depth to start diving from.
-  i_t min_node_depth = 10;
-
   // The maximum number of nodes when performing a dive.
   i_t node_limit = 500;
 
   // The maximum number of dual simplex iteration allowed
-  // in a single dive. This set in terms of the total number of
-  // iterations in the best-first threads.
-  f_t iteration_limit_factor = 0.05;
+  // in a single dive.
+  f_t iteration_limit_factor     = 0.05;
+  int64_t iteration_limit_offset = 10000;
 
   // The maximum backtracking allowed.
   i_t backtrack_limit = 5;

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -208,12 +208,6 @@ class branch_and_bound_t {
   std::vector<i_t> new_slacks_;
   std::vector<simplex::variable_type_t> var_types_;
 
-  // Variable locks (see definition 3.3 from T. Achterberg, “Constraint Integer Programming,”
-  // PhD, Technischen Universität Berlin, Berlin, 2007. doi: 10.14279/depositonce-1634).
-  // Here we assume that the constraints are in the form `Ax = b, l <= x <= u`.
-  std::vector<i_t> var_up_locks_;
-  std::vector<i_t> var_down_locks_;
-
   // Mutex for the original LP
   // The heuristics threads look at the original LP. But the main thread modifies the
   // size of the original LP by adding slacks for cuts. Heuristic threads should lock
@@ -296,7 +290,8 @@ class branch_and_bound_t {
 
   void print_table_header();
   void report_heuristic(f_t obj, heuristics_origin_t origin);
-  void report(char symbol,
+  void report(const simplex::lp_problem_t<i_t, f_t>& lp,
+              char symbol,
               f_t obj,
               f_t lower_bound,
               i_t node_depth,
@@ -333,14 +328,15 @@ class branch_and_bound_t {
   // Set the solution when found at the root node
   void set_solution_at_root(simplex::mip_solution_t<i_t, f_t>& solution,
                             const cut_info_t<i_t, f_t>& cut_info);
-  void update_user_bound(f_t lower_bound);
+  void update_user_bound(const simplex::lp_problem_t<i_t, f_t>& lp, f_t lower_bound);
 
   // Set the final solution.
   void set_final_solution(simplex::mip_solution_t<i_t, f_t>& solution, f_t lower_bound);
 
   // Update the incumbent solution with the new feasible solution
   // found during branch and bound.
-  void add_feasible_solution(f_t leaf_objective,
+  void add_feasible_solution(const simplex::lp_problem_t<i_t, f_t>& lp,
+                             f_t leaf_objective,
                              const std::vector<f_t>& leaf_solution,
                              i_t leaf_depth,
                              search_strategy_t thread_type);
@@ -351,7 +347,9 @@ class branch_and_bound_t {
   // Launch a new diving worker from a given best-first worker.
   bool launch_diving_worker(bfs_worker_t<i_t, f_t>* bfs_worker);
 
-  void snap_to_lattice(mip_node_t<i_t, f_t>* node_ptr, f_t leaf_obj);
+  void snap_to_lattice(const simplex::lp_problem_t<i_t, f_t>& lp,
+                       mip_node_t<i_t, f_t>* node_ptr,
+                       f_t leaf_obj);
 
   // Launch a new best-first worker from a given bfs worker.
   void launch_bfs_worker(bfs_worker_t<i_t, f_t>* worker);
@@ -369,7 +367,8 @@ class branch_and_bound_t {
 
   // Perform a deep dive in the subtree determined by the `start_node` in order
   // to find integer feasible solutions.
-  void dive_with(diving_worker_t<i_t, f_t>* worker, i_t backtrack_limit);
+  void dive_with(diving_worker_t<i_t, f_t>* worker,
+                 const simplex::simplex_solver_settings_t<i_t, f_t>& settings);
 
   // Launch a new RINS worker
   bool launch_submip_worker(const std::vector<f_t>& sol);
@@ -382,21 +381,21 @@ class branch_and_bound_t {
 
   // 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,
                     submip_stats_t& submip_stats,
                     f_t fixrate,
                     i_t simplex_iter_used,
                     simplex::simplex_solver_settings_t<i_t, f_t> submip_settings);
 
   // Creates and solves the RINS/RENS sub-MIP.
   void recursive_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,
                         simplex::simplex_solver_settings_t<i_t, f_t> submip_settings);
 
   void launch_root_heuristics(const simplex::lp_problem_t<i_t, f_t>& lp,
-                              const std::vector<f_t>& sol,
+                              const simplex::lp_solution_t<i_t, f_t>& lp_solution,
+                              const std::vector<i_t>& fractional,
+                              const std::vector<i_t>& basic_list,
+                              const std::vector<i_t>& nonbasic_list,
+                              simplex::basis_update_mpf_t<i_t, f_t>& basis_factor,
                               i_t cut_pass,
                               root_heuristics_t<i_t, f_t>& root_heuristics);
 
@@ -405,7 +404,7 @@ class branch_and_bound_t {
                                        branch_and_bound_worker_t<i_t, f_t>* worker,
                                        branch_and_bound_stats_t<i_t, f_t>& stats,
                                        simplex::logger_t& log,
-                                       i_t iter_limit = std::numeric_limits<i_t>::max());
+                                       int64_t iter_limit = std::numeric_limits<i_t>::max());
 
   // Apply symmetry-based bound reductions (orbital fixing and, when
   // settings_.symmetry == 2, lexical reduction) to the current node.

cpp/src/branch_and_bound/deterministic_workers.hpp

@@ -89,10 +89,11 @@ 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,
+                              pseudo_costs_t<i_t, f_t>& pc,
                               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, root_solution, root_edge_norm),
+    : base_t(id, original_lp, Arow, var_types, settings, pc, root_solution, root_edge_norm),
       work_context(context_name),
       pc_snapshot(1, settings)
   {
@@ -143,13 +144,15 @@ class deterministic_bfs_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,
+                                      pseudo_costs_t<i_t, f_t>& pc,
                                       const std::vector<f_t>& root_solution,
                                       const std::vector<f_t>& root_edge_norm)
     : base_t(id,
              original_lp,
              Arow,
              var_types,
              settings,
+             pc,
              root_solution,
              root_edge_norm,
              "BB_Worker_" + std::to_string(id))
@@ -308,13 +311,15 @@ 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,
+    pseudo_costs_t<i_t, f_t>& pc,
     const std::vector<f_t>& root_solution,
     const std::vector<f_t>& root_edge_norm)
     : base_t(id,
              original_lp,
              Arow,
              var_types,
              settings,
+             pc,
              root_solution,
              root_edge_norm,
              "Diving_Worker_" + std::to_string(id)),
@@ -423,13 +428,14 @@ class deterministic_bfs_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,
+                                  pseudo_costs_t<i_t, f_t>& pc,
                                   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, root_solution, root_edge_norm);
+        i, original_lp, Arow, var_types, settings, pc, root_solution, root_edge_norm);
     }
   }
 
@@ -461,14 +467,15 @@ 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,
+                                     pseudo_costs_t<i_t, f_t>& pc,
                                      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, root_edge_norm);
+        i, type, original_lp, Arow, var_types, settings, pc, root_solution, root_edge_norm);
     }
   }
 

cpp/src/branch_and_bound/pseudo_costs.cpp

@@ -1369,6 +1369,69 @@ void strong_branching(const lp_problem_t<i_t, f_t>& original_lp,
   }
 }
 
+template <typename i_t, typename f_t>
+void pseudo_costs_t<i_t, f_t>::initialize_with_estimate(
+  const lp_problem_t<i_t, f_t>& lp,
+  const std::vector<simplex::variable_status_t>& vstatus,
+  const std::vector<i_t>& fractional,
+  const lp_solution_t<i_t, f_t>& lp_solution,
+  const std::vector<i_t>& basic_list,
+  const std::vector<i_t>& nonbasic_list,
+  basis_update_mpf_t<i_t, f_t>& basis_factors)
+{
+  i_t m = lp.num_rows;
+  i_t n = lp.num_cols;
+
+  std::vector<f_t> delta_z(n, 0);
+  std::vector<i_t> workspace(n, 0);
+
+  f_t work_estimate = 0;
+
+  std::vector<i_t> basic_map(n, -1);
+  for (i_t i = 0; i < m; i++) {
+    basic_map[basic_list[i]] = i;
+  }
+
+  // compute_initial_nonbasic_end permutes columns in place; copy so pc.Arow is unchanged
+  csr_matrix_t<i_t, f_t> local_Arow = Arow;
+
+  std::vector<i_t> nonbasic_end(m);
+  compute_initial_nonbasic_end(basic_map, local_Arow, nonbasic_end);
+
+  for (i_t k = 0; k < fractional.size(); k++) {
+    const i_t j = fractional[k];
+    assert(j >= 0);
+
+    if (pseudo_cost_num_up[j] == 0 || pseudo_cost_num_down[j] == 0) {
+      objective_change_estimate_t<f_t> estimate =
+        single_pivot_objective_change_estimate(lp,
+                                               settings,
+                                               local_Arow,
+                                               vstatus,
+                                               j,
+                                               basic_map[j],
+                                               lp_solution,
+                                               basic_list,
+                                               nonbasic_list,
+                                               nonbasic_end,
+                                               basis_factors,
+                                               workspace,
+                                               delta_z,
+                                               work_estimate);
+
+      if (pseudo_cost_num_down[j] == 0) {
+        pseudo_cost_sum_down[j] += estimate.down_obj_change;
+        ++pseudo_cost_num_down[j];
+      }
+
+      if (pseudo_cost_num_up[j] == 0) {
+        pseudo_cost_sum_up[j] += estimate.up_obj_change;
+        ++pseudo_cost_num_up[j];
+      }
+    }
+  }
+}
+
 template <typename i_t, typename f_t>
 inline f_t pseudo_costs_t<i_t, f_t>::compute_pseudocost_average_down()
 {

cpp/src/branch_and_bound/pseudo_costs.hpp

@@ -216,6 +216,14 @@ class pseudo_costs_t {
   f_t compute_pseudocost_average_down();
   f_t compute_pseudocost_average_up();
 
+  void initialize_with_estimate(const simplex::lp_problem_t<i_t, f_t>& lp,
+                                const std::vector<simplex::variable_status_t>& vstatus,
+                                const std::vector<i_t>& fractional,
+                                const simplex::lp_solution_t<i_t, f_t>& lp_solution,
+                                const std::vector<i_t>& basic_list,
+                                const std::vector<i_t>& nonbasic_list,
+                                simplex::basis_update_mpf_t<i_t, f_t>& basis_factors);
+
   f_t obj_estimate(const std::vector<i_t>& fractional,
                    const std::vector<f_t>& solution,
                    f_t lower_bound);