Add clique cuts (#937)

This PR adds clique cuts to the set of cuts we have. It uses the ideas from Preprocessing and Cutting Planes with Conflict Graphs, Brito et.al. 

In a previous PR, I have added the clique table generation as a preprocessing step (currently doesn't change the problem). This PR builds on that by using the extended cliques during cut passes. At each cut pass, a subgraph induced by the fractinoal variables on the current relaxation solution is used to generate clique cuts. We want to generate the set of violated clique inequalities that have at least min violation (1e-6). The maximal cliques in the fractional subgraph is found by Bron Kerbosh dynamic programming algorithm. The found maximal cliques are extended with the suitable integer valued variables on the original conflict graph. The extension is guided by reduced cost of the variables at the current relaxation optimal. We try to extend the cliques with variables that have lower reduced cost. The reason is that it disturbs dual simplex less (causes fewer refactors and iterations for resolve). Also the variables with lower reduced costs are likely to be active in the next resolve. This also results in more dominant cliques preventing multiple smaller sized cliques in the next iterations.

Benchmark results:
MIP GAP wins/losses defined as at least 10% difference:
main wins (6):
  Instance                                   dir1 gap   dir2 gap       diff
  app1-2                                       0.8077     1.1818    +0.3741
  bab2                                         0.4148     0.6148    +0.2001
  glass4                                       0.3333     0.4483    +0.1149
  map16715-04                                  0.9269     1.4085    +0.4816
  neos-4413714-turia                           0.7027     0.8454    +0.1427
  rail01                                       0.5476     0.7934    +0.2458

clique cuts wins (11):
  Instance                                   dir1 gap   dir2 gap       diff
  30n20b8                                      0.7688     0.6667    +0.1021
  bab6                                         0.9255     0.3611    +0.5643
  buildingenergy                               0.1706     0.0036    +0.1670
  co-100                                       0.8410     0.7089    +0.1321
  neos-3754480-nidda                          20.3512    19.6495    +0.7017
  neos-5188808-nattai                          0.2032     0.0000    +0.2032
  netdiversion                                 0.3710     0.1429    +0.2281
  physiciansched3-3                            0.2444     0.0269    +0.2175
  rocII-5-11                                   1.0893     0.7762    +0.3131
  satellites2-60-fs                           13.5000     1.4167   +12.0833
  sorrell3                                     9.5625     9.4375    +0.1250

Optimality wins/loses:
Wins

 graphdraw-domain  (443.68s)
 neos-5188808-nattai  (462.84s)
 supportcase26  (537.66s)

Losses

dano3_3  (121.30s)


Time to optimality geomean (shifted geomean by 1.0):
on 50 common optimal instances
Common optimal instances: 50
Geometric mean of ratios (main/clique): 1.1925

Authors:
  - Akif ÇÖRDÜK (https://github.com/akifcorduk)
  - Alice Boucher (https://github.com/aliceb-nv)

Approvers:
  - Chris Maes (https://github.com/chris-maes)

URL: #937

Author: akifcorduk

benchmarks/linear_programming/cuopt/run_mip.cpp

@@ -183,6 +183,10 @@ int run_single_file(std::string file_path,
     CUOPT_LOG_ERROR("Parsing MPS failed exiting!");
     return -1;
   }
+  // Use the benchmark filename for downstream instance-level reporting.
+  // This keeps per-instance metrics aligned with the run list even if the MPS NAME card differs.
+  mps_data_model.set_problem_name(base_filename);
+
   if (initial_solution_dir.has_value()) {
     auto initial_solutions = read_solution_from_dir(
       initial_solution_dir.value(), base_filename, mps_data_model.get_variable_names());
@@ -209,6 +213,7 @@ int run_single_file(std::string file_path,
   settings.tolerances.absolute_tolerance = 1e-6;
   settings.presolver                     = cuopt::linear_programming::presolver_t::Default;
   settings.reliability_branching         = reliability_branching;
+  settings.clique_cuts                   = -1;
   settings.seed                          = 42;
   cuopt::linear_programming::benchmark_info_t benchmark_info;
   settings.benchmark_info_ptr = &benchmark_info;
@@ -413,7 +418,16 @@ int main(int argc, char* argv[])
   int reliability_branching = program.get<int>("--reliability-branching");
   bool deterministic        = program.get<bool>("--determinism");
 
-  if (num_cpu_threads < 0) { num_cpu_threads = omp_get_max_threads() / n_gpus; }
+  if (num_cpu_threads < 0) {
+    num_cpu_threads = omp_get_max_threads() / n_gpus;
+    // std::ifstream smt_file("/sys/devices/system/cpu/smt/active");
+    // if (smt_file.is_open()) {
+    //   int smt_active = 0;
+    //   smt_file >> smt_active;
+    //   if (smt_active) { num_cpu_threads /= 2; }
+    // }
+    num_cpu_threads = std::max(num_cpu_threads, 1);
+  }
 
   if (program.is_used("--out-dir")) {
     out_dir     = program.get<std::string>("--out-dir");

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_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"
 #define CUOPT_MIP_CUT_CHANGE_THRESHOLD        "mip_cut_change_threshold"

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

@@ -93,6 +93,7 @@ class mip_solver_settings_t {
   i_t mir_cuts                  = -1;
   i_t mixed_integer_gomory_cuts = -1;
   i_t knapsack_cuts             = -1;
+  i_t clique_cuts               = -1;
   i_t strong_chvatal_gomory_cuts      = -1;
   i_t reduced_cost_strengthening      = -1;
   f_t cut_change_threshold            = 1e-3;

cpp/src/branch_and_bound/branch_and_bound.cpp

@@ -10,6 +10,7 @@
 #include <branch_and_bound/pseudo_costs.hpp>
 
 #include <cuts/cuts.hpp>
+#include <mip_heuristics/presolve/conflict_graph/clique_table.cuh>
 
 #include <dual_simplex/basis_solves.hpp>
 #include <dual_simplex/bounds_strengthening.hpp>
@@ -241,9 +242,11 @@ template <typename i_t, typename f_t>
 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)
+  f_t start_time,
+  std::shared_ptr<detail::clique_table_t<i_t, f_t>> clique_table)
   : original_problem_(user_problem),
     settings_(solver_settings),
+    clique_table_(std::move(clique_table)),
     original_lp_(user_problem.handle_ptr, 1, 1, 1),
     Arow_(1, 1, 0),
     incumbent_(1),
@@ -1967,6 +1970,31 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
 
   root_relax_soln_.resize(original_lp_.num_rows, original_lp_.num_cols);
 
+  if (settings_.clique_cuts != 0 && clique_table_ == nullptr) {
+    signal_extend_cliques_.store(false, std::memory_order_release);
+    typename ::cuopt::linear_programming::mip_solver_settings_t<i_t, f_t>::tolerances_t
+      tolerances_for_clique{};
+    tolerances_for_clique.presolve_absolute_tolerance = settings_.primal_tol;
+    tolerances_for_clique.absolute_tolerance          = settings_.primal_tol;
+    tolerances_for_clique.relative_tolerance          = settings_.zero_tol;
+    tolerances_for_clique.integrality_tolerance       = settings_.integer_tol;
+    tolerances_for_clique.absolute_mip_gap            = settings_.absolute_mip_gap_tol;
+    tolerances_for_clique.relative_mip_gap            = settings_.relative_mip_gap_tol;
+    auto* signal_ptr                                  = &signal_extend_cliques_;
+    clique_table_future_ =
+      std::async(std::launch::async,
+                 [this,
+                  tolerances_for_clique,
+                  signal_ptr]() -> std::shared_ptr<detail::clique_table_t<i_t, f_t>> {
+                   user_problem_t<i_t, f_t> problem_copy = original_problem_;
+                   cuopt::timer_t timer(std::numeric_limits<double>::infinity());
+                   std::shared_ptr<detail::clique_table_t<i_t, f_t>> table;
+                   detail::find_initial_cliques(
+                     problem_copy, tolerances_for_clique, &table, timer, false, signal_ptr);
+                   return table;
+                 });
+  }
+
   i_t original_rows                     = original_lp_.num_rows;
   simplex_solver_settings_t lp_settings = settings_;
   lp_settings.inside_mip                = 1;
@@ -2002,38 +2030,44 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   exploration_stats_.total_lp_iters      = root_relax_soln_.iterations;
   exploration_stats_.total_lp_solve_time = toc(exploration_stats_.start_time);
 
+  auto finish_clique_thread = [this]() {
+    if (clique_table_future_.valid()) {
+      signal_extend_cliques_.store(true, std::memory_order_release);
+      clique_table_ = clique_table_future_.get();
+    }
+  };
+
   if (root_status == lp_status_t::INFEASIBLE) {
     settings_.log.printf("MIP Infeasible\n");
-    // FIXME: rarely dual simplex detects infeasible whereas it is feasible.
-    // to add a small safety net, check if there is a primal solution already.
-    // Uncomment this if the issue with cost266-UUE is resolved
-    // if (settings.heuristic_preemption_callback != nullptr) {
-    //   settings.heuristic_preemption_callback();
-    // }
+    finish_clique_thread();
     return mip_status_t::INFEASIBLE;
   }
   if (root_status == lp_status_t::UNBOUNDED) {
     settings_.log.printf("MIP Unbounded\n");
     if (settings_.heuristic_preemption_callback != nullptr) {
       settings_.heuristic_preemption_callback();
     }
+    finish_clique_thread();
     return mip_status_t::UNBOUNDED;
   }
   if (root_status == lp_status_t::TIME_LIMIT) {
     solver_status_ = mip_status_t::TIME_LIMIT;
     set_final_solution(solution, -inf);
+    finish_clique_thread();
     return solver_status_;
   }
 
   if (root_status == lp_status_t::WORK_LIMIT) {
     solver_status_ = mip_status_t::WORK_LIMIT;
     set_final_solution(solution, -inf);
+    finish_clique_thread();
     return solver_status_;
   }
 
   if (root_status == lp_status_t::NUMERICAL_ISSUES) {
     solver_status_ = mip_status_t::NUMERICAL;
     set_final_solution(solution, -inf);
+    finish_clique_thread();
     return solver_status_;
   }
 
@@ -2064,6 +2098,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
 
   if (num_fractional == 0) {
     set_solution_at_root(solution, cut_info);
+    finish_clique_thread();
     return mip_status_t::OPTIMAL;
   }
 
@@ -2078,8 +2113,16 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   cut_pool_t<i_t, f_t> cut_pool(original_lp_.num_cols, settings_);
-  cut_generation_t<i_t, f_t> cut_generation(
-    cut_pool, original_lp_, settings_, Arow_, new_slacks_, var_types_);
+  cut_generation_t<i_t, f_t> cut_generation(cut_pool,
+                                            original_lp_,
+                                            settings_,
+                                            Arow_,
+                                            new_slacks_,
+                                            var_types_,
+                                            original_problem_,
+                                            clique_table_,
+                                            &clique_table_future_,
+                                            &signal_extend_cliques_);
 
   std::vector<f_t> saved_solution;
 #ifdef CHECK_CUTS_AGAINST_SAVED_SOLUTION
@@ -2090,7 +2133,8 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   f_t last_objective       = root_objective_;
   f_t root_relax_objective = root_objective_;
 
-  i_t cut_pool_size = 0;
+  f_t cut_generation_start_time = tic();
+  i_t cut_pool_size             = 0;
   for (i_t cut_pass = 0; cut_pass < settings_.max_cut_passes; cut_pass++) {
     if (num_fractional == 0) {
       set_solution_at_root(solution, cut_info);
@@ -2109,16 +2153,25 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
 #endif
 
       // Generate cuts and add them to the cut pool
-      f_t cut_start_time = tic();
-      cut_generation.generate_cuts(original_lp_,
-                                   settings_,
-                                   Arow_,
-                                   new_slacks_,
-                                   var_types_,
-                                   basis_update,
-                                   root_relax_soln_.x,
-                                   basic_list,
-                                   nonbasic_list);
+      f_t cut_start_time    = tic();
+      bool problem_feasible = cut_generation.generate_cuts(original_lp_,
+                                                           settings_,
+                                                           Arow_,
+                                                           new_slacks_,
+                                                           var_types_,
+                                                           basis_update,
+                                                           root_relax_soln_.x,
+                                                           root_relax_soln_.z,
+                                                           basic_list,
+                                                           nonbasic_list,
+                                                           exploration_stats_.start_time);
+      if (!problem_feasible) {
+        if (settings_.heuristic_preemption_callback != nullptr) {
+          settings_.heuristic_preemption_callback();
+        }
+        finish_clique_thread();
+        return mip_status_t::INFEASIBLE;
+      }
       f_t cut_generation_time = toc(cut_start_time);
       if (cut_generation_time > 1.0) {
         settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time);
@@ -2339,8 +2392,9 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   print_cut_info(settings_, cut_info);
-
+  f_t cut_generation_time = toc(cut_generation_start_time);
   if (cut_info.has_cuts()) {
+    settings_.log.printf("Cut generation time: %.2f seconds\n", cut_generation_time);
     settings_.log.printf("Cut pool size  : %d\n", cut_pool_size);
     settings_.log.printf("Size with cuts : %d constraints, %d variables, %d nonzeros\n",
                          original_lp_.num_rows,

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -32,9 +32,17 @@
 
 #include <omp.h>
 
+#include <atomic>
 #include <functional>
+#include <future>
+#include <memory>
 #include <vector>
 
+namespace cuopt::linear_programming::detail {
+template <typename i_t, typename f_t>
+struct clique_table_t;
+}
+
 namespace cuopt::linear_programming::dual_simplex {
 
 enum class mip_status_t {
@@ -68,7 +76,8 @@ class branch_and_bound_t {
  public:
   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);
+                     f_t start_time,
+                     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.
   void set_initial_guess(const std::vector<f_t>& user_guess) { guess_ = user_guess; }
@@ -106,8 +115,6 @@ class branch_and_bound_t {
 
   void set_concurrent_lp_root_solve(bool enable) { enable_concurrent_lp_root_solve_ = enable; }
 
-  bool stop_for_time_limit(mip_solution_t<i_t, f_t>& solution);
-
   // Repair a low-quality solution from the heuristics.
   bool repair_solution(const std::vector<f_t>& leaf_edge_norms,
                        const std::vector<f_t>& potential_solution,
@@ -141,6 +148,9 @@ 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_;
+  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};
 
   work_limit_context_t work_unit_context_{"B&B"};