Integrate CliqueMerging presolver and tune presolver settings

ci/test_wheel_cuopt.sh

@@ -73,10 +73,8 @@ timeout 10m bash ./python/libcuopt/libcuopt/tests/test_cli.sh
 # Run Python tests
 RAPIDS_DATASET_ROOT_DIR=./datasets timeout 30m python -m pytest --verbose --capture=no ./python/cuopt/cuopt/tests/
 
-# run cvxpy integration tests
-./ci/thirdparty-testing/run_cvxpy_tests.sh
-
-# run jump tests for only nightly builds
+# run jump tests and cvxpy integration tests for only nightly builds
 if [[ "${RAPIDS_BUILD_TYPE}" == "nightly" ]]; then
     ./ci/thirdparty-testing/run_jump_tests.sh
+    ./ci/thirdparty-testing/run_cvxpy_tests.sh
 fi

ci/utils/install_boost_tbb.sh

@@ -24,7 +24,8 @@ if [ -f /etc/os-release ]; then
         echo "Detected Rocky Linux. Installing Boost and TBB via dnf..."
         dnf clean all
         dnf -y update
-        dnf install -y boost-devel tbb-devel
+        dnf install -y epel-release
+        dnf install -y boost1.78-devel tbb-devel
         if [[ "$(uname -m)" == "x86_64" ]]; then
             dnf install -y gcc-toolset-14-libquadmath-devel
         fi

cpp/CMakeLists.txt

@@ -183,7 +183,14 @@ endif()
 FetchContent_Declare(
   papilo
   GIT_REPOSITORY "https://github.com/scipopt/papilo.git"
-  GIT_TAG "v2.4.3"
+  # We would want to get the main branch. However, the main branch
+  # does not have some of the presolvers and settings that we need
+  # Mainly, probing and clique merging.
+  # This is the reason we are using the development branch
+  # commit from Sep 26, 2025. Once these changes are merged into the main branch,
+  #we can switch to the main branch.
+  GIT_TAG "34a40781fa14f8870cb6368cffb6c0eda2f47511"
+  GIT_PROGRESS TRUE
   SYSTEM
 )
 
@@ -422,7 +429,6 @@ if(NOT BUILD_LP_ONLY)
     cuopt
     OpenMP::OpenMP_CXX
     PRIVATE
-    papilo-core
   )
   set_property(TARGET cuopt_cli PROPERTY INSTALL_RPATH "$ORIGIN/../${lib_dir}")
 
@@ -446,7 +452,6 @@ if(BUILD_MIP_BENCHMARKS AND NOT BUILD_LP_ONLY)
     cuopt
     OpenMP::OpenMP_CXX
     PRIVATE
-    papilo-core
   )
 endif()
 
@@ -462,7 +467,6 @@ if(BUILD_LP_BENCHMARKS)
     cuopt
     OpenMP::OpenMP_CXX
     PRIVATE
-    papilo-core
   )
 endif()
 

cpp/include/cuopt/linear_programming/constants.h

@@ -51,6 +51,7 @@
 #define CUOPT_LOG_TO_CONSOLE              "log_to_console"
 #define CUOPT_CROSSOVER                   "crossover"
 #define CUOPT_PRESOLVE                    "presolve"
+#define CUOPT_DUAL_POSTSOLVE              "dual_postsolve"
 #define CUOPT_MIP_ABSOLUTE_TOLERANCE      "mip_absolute_tolerance"
 #define CUOPT_MIP_RELATIVE_TOLERANCE      "mip_relative_tolerance"
 #define CUOPT_MIP_INTEGRALITY_TOLERANCE   "mip_integrality_tolerance"

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

@@ -211,6 +211,7 @@ class pdlp_solver_settings_t {
   bool save_best_primal_so_far{false};
   bool first_primal_feasible{false};
   bool presolve{false};
+  bool dual_postsolve{true};
   method_t method{method_t::Concurrent};
   // For concurrent termination
   std::atomic<i_t>* concurrent_halt;

cpp/src/dual_simplex/branch_and_bound.cpp

@@ -567,8 +567,6 @@ dual::status_t branch_and_bound_t<i_t, f_t>::node_dual_simplex(
         leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms);
       lp_status = convert_lp_status_to_dual_status(second_status);
     }
-  } else {
-    log.printf("Infeasible after bounds strengthening. Fathoming node %d.\n", leaf_id);
   }
 
   mutex_stats_.lock();

cpp/src/linear_programming/solve.cu

@@ -700,11 +700,12 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(optimization_problem_t<i_t, f
     if (run_presolve) {
       // allocate no more than 10% of the time limit to presolve.
       // Note that this is not the presolve time, but the time limit for presolve.
-      const double presolve_time_limit = 0.1 * lp_timer.remaining_time();
+      const double presolve_time_limit = std::min(0.1 * lp_timer.remaining_time(), 60.0);
       presolver = std::make_unique<detail::third_party_presolve_t<i_t, f_t>>();
       auto [reduced_problem, feasible] =
         presolver->apply(op_problem,
                          cuopt::linear_programming::problem_category_t::LP,
+                         settings.dual_postsolve,
                          settings.tolerances.absolute_primal_tolerance,
                          settings.tolerances.relative_primal_tolerance,
                          presolve_time_limit);
@@ -714,7 +715,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(optimization_problem_t<i_t, f
       }
       problem       = detail::problem_t<i_t, f_t>(reduced_problem);
       presolve_time = lp_timer.elapsed_time();
-      CUOPT_LOG_INFO("Third party presolve time: %f", presolve_time);
+      CUOPT_LOG_INFO("Papilo presolve time: %f", presolve_time);
     }
 
     CUOPT_LOG_INFO(

cpp/src/math_optimization/solver_settings.cu

@@ -107,7 +107,8 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_LOG_TO_CONSOLE, &mip_settings.log_to_console, true},
     {CUOPT_CROSSOVER, &pdlp_settings.crossover, false},
     {CUOPT_PRESOLVE, &pdlp_settings.presolve, false},
-    {CUOPT_PRESOLVE, &mip_settings.presolve, true}
+    {CUOPT_PRESOLVE, &mip_settings.presolve, true},
+    {CUOPT_DUAL_POSTSOLVE, &pdlp_settings.dual_postsolve, true}
   };
   // String parameters
   string_parameters = {

cpp/src/mip/presolve/third_party_presolve.cpp

@@ -34,7 +34,8 @@ static papilo::PostsolveStorage<double> post_solve_storage_;
 static bool maximize_ = false;
 
 template <typename i_t, typename f_t>
-papilo::Problem<f_t> build_papilo_problem(const optimization_problem_t<i_t, f_t>& op_problem)
+papilo::Problem<f_t> build_papilo_problem(const optimization_problem_t<i_t, f_t>& op_problem,
+                                          problem_category_t category)
 {
   // Build papilo problem from optimization problem
   papilo::ProblemBuilder<f_t> builder;
@@ -167,7 +168,11 @@ papilo::Problem<f_t> build_papilo_problem(const optimization_problem_t<i_t, f_t>
 
   if (h_entries.size()) {
     auto constexpr const sorted_entries = true;
-    auto csr_storage = papilo::SparseStorage<f_t>(h_entries, num_rows, num_cols, sorted_entries);
+    // MIP reductions like clique merging and substituition require more fillin
+    const double spare_ratio      = category == problem_category_t::MIP ? 4.0 : 2.0;
+    const int min_inter_row_space = category == problem_category_t::MIP ? 30 : 4;
+    auto csr_storage              = papilo::SparseStorage<f_t>(
+      h_entries, num_rows, num_cols, sorted_entries, spare_ratio, min_inter_row_space);
     problem.setConstraintMatrix(csr_storage, h_constr_lb, h_constr_ub, h_row_flags);
 
     papilo::ConstraintMatrix<f_t>& matrix = problem.getConstraintMatrix();
@@ -304,14 +309,16 @@ void check_postsolve_status(const papilo::PostsolveStatus& status)
 }
 
 template <typename f_t>
-void set_presolve_methods(papilo::Presolve<f_t>& presolver, problem_category_t category)
+void set_presolve_methods(papilo::Presolve<f_t>& presolver,
+                          problem_category_t category,
+                          bool dual_postsolve)
 {
   using uptr = std::unique_ptr<papilo::PresolveMethod<f_t>>;
 
-  // cuopt custom presolvers
-  if (category == problem_category_t::MIP)
+  if (category == problem_category_t::MIP) {
+    // cuOpt custom GF2 presolver
     presolver.addPresolveMethod(uptr(new cuopt::linear_programming::detail::GF2Presolve<f_t>()));
-
+  }
   // fast presolvers
   presolver.addPresolveMethod(uptr(new papilo::SingletonCols<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::CoefficientStrengthening<f_t>()));
@@ -326,16 +333,21 @@ void set_presolve_methods(papilo::Presolve<f_t>& presolver, problem_category_t c
   presolver.addPresolveMethod(uptr(new papilo::SingletonStuffing<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::DualFix<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::SimplifyInequalities<f_t>()));
+  presolver.addPresolveMethod(uptr(new papilo::CliqueMerging<f_t>()));
 
   // exhaustive presolvers
   presolver.addPresolveMethod(uptr(new papilo::ImplIntDetection<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::DominatedCols<f_t>()));
   presolver.addPresolveMethod(uptr(new papilo::Probing<f_t>()));
 
-  presolver.addPresolveMethod(uptr(new papilo::DualInfer<f_t>));
-  presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution<f_t>()));
-  presolver.addPresolveMethod(uptr(new papilo::Sparsify<f_t>()));
-  presolver.addPresolveMethod(uptr(new papilo::Substitution<f_t>()));
+  if (!dual_postsolve) {
+    presolver.addPresolveMethod(uptr(new papilo::DualInfer<f_t>()));
+    presolver.addPresolveMethod(uptr(new papilo::SimpleSubstitution<f_t>()));
+    presolver.addPresolveMethod(uptr(new papilo::Sparsify<f_t>()));
+    presolver.addPresolveMethod(uptr(new papilo::Substitution<f_t>()));
+  } else {
+    CUOPT_LOG_INFO("Disabling the presolver methods that do not support dual postsolve");
+  }
 }
 
 template <typename i_t, typename f_t>
@@ -351,26 +363,51 @@ void set_presolve_options(papilo::Presolve<f_t>& presolver,
   presolver.getPresolveOptions().feastol = 1e-5;
 }
 
+template <typename f_t>
+void set_presolve_parameters(papilo::Presolve<f_t>& presolver,
+                             problem_category_t category,
+                             int nrows,
+                             int ncols)
+{
+  // It looks like a copy. But this copy has the pointers to relevant variables in papilo
+  auto params = presolver.getParameters();
+  if (category == problem_category_t::MIP) {
+    // Papilo has work unit measurements for probing. Because of this when the first batch fails to
+    // produce any reductions, the algorithm stops. To avoid stopping the algorithm, we set a
+    // minimum badge size to a huge value. The time limit makes sure that we exit if it takes too
+    // long
+    int min_badgesize = std::max(ncols / 2, 32);
+    params.setParameter("probing.minbadgesize", min_badgesize);
+    params.setParameter("cliquemerging.enabled", true);
+    params.setParameter("cliquemerging.maxcalls", 50);
+  }
+}
+
 template <typename i_t, typename f_t>
 std::pair<optimization_problem_t<i_t, f_t>, bool> third_party_presolve_t<i_t, f_t>::apply(
   optimization_problem_t<i_t, f_t> const& op_problem,
   problem_category_t category,
+  bool dual_postsolve,
   f_t absolute_tolerance,
   f_t relative_tolerance,
   double time_limit,
   i_t num_cpu_threads)
 {
-  papilo::Problem<f_t> papilo_problem = build_papilo_problem(op_problem);
+  papilo::Problem<f_t> papilo_problem = build_papilo_problem(op_problem, category);
 
   CUOPT_LOG_INFO("Unpresolved problem:: %d constraints, %d variables, %d nonzeros",
                  papilo_problem.getNRows(),
                  papilo_problem.getNCols(),
                  papilo_problem.getConstraintMatrix().getNnz());
 
+  CUOPT_LOG_INFO("Calling Papilo presolver");
+  if (category == problem_category_t::MIP) { dual_postsolve = false; }
   papilo::Presolve<f_t> presolver;
-  set_presolve_methods<f_t>(presolver, category);
+  set_presolve_methods<f_t>(presolver, category, dual_postsolve);
   set_presolve_options<i_t, f_t>(
     presolver, category, absolute_tolerance, relative_tolerance, time_limit, num_cpu_threads);
+  set_presolve_parameters<f_t>(
+    presolver, category, op_problem.get_n_constraints(), op_problem.get_n_variables());
 
   // Disable papilo logs
   presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet);
@@ -423,9 +460,12 @@ void third_party_presolve_t<i_t, f_t>::undo(rmm::device_uvector<f_t>& primal_sol
   check_postsolve_status(status);
 
   primal_solution.resize(full_sol.primal.size(), stream_view);
-  dual_solution.resize(full_sol.primal.size(), stream_view);
-  reduced_costs.resize(full_sol.primal.size(), stream_view);
+  dual_solution.resize(full_sol.dual.size(), stream_view);
+  reduced_costs.resize(full_sol.reducedCosts.size(), stream_view);
   raft::copy(primal_solution.data(), full_sol.primal.data(), full_sol.primal.size(), stream_view);
+  raft::copy(dual_solution.data(), full_sol.dual.data(), full_sol.dual.size(), stream_view);
+  raft::copy(
+    reduced_costs.data(), full_sol.reducedCosts.data(), full_sol.reducedCosts.size(), stream_view);
 }
 
 #if MIP_INSTANTIATE_FLOAT

cpp/src/mip/presolve/third_party_presolve.hpp

@@ -29,6 +29,7 @@ class third_party_presolve_t {
   std::pair<optimization_problem_t<i_t, f_t>, bool> apply(
     optimization_problem_t<i_t, f_t> const& op_problem,
     problem_category_t category,
+    bool dual_postsolve,
     f_t absolute_tolerance,
     f_t relative_tolerance,
     double time_limit,

cpp/src/mip/solve.cu

@@ -202,11 +202,13 @@ mip_solution_t<i_t, f_t> solve_mip(optimization_problem_t<i_t, f_t>& op_problem,
     if (run_presolve) {
       // allocate not more than 10% of the time limit to presolve.
       // Note that this is not the presolve time, but the time limit for presolve.
-      const double presolve_time_limit = 0.1 * time_limit;
+      const double presolve_time_limit = std::min(0.1 * time_limit, 60.0);
+      const bool dual_postsolve        = false;
       presolver = std::make_unique<detail::third_party_presolve_t<i_t, f_t>>();
       auto [reduced_op_problem, feasible] =
         presolver->apply(op_problem,
                          cuopt::linear_programming::problem_category_t::MIP,
+                         dual_postsolve,
                          settings.tolerances.absolute_tolerance,
                          settings.tolerances.relative_tolerance,
                          presolve_time_limit,
@@ -219,7 +221,7 @@ mip_solution_t<i_t, f_t> solve_mip(optimization_problem_t<i_t, f_t>& op_problem,
 
       problem       = detail::problem_t<i_t, f_t>(reduced_op_problem);
       presolve_time = timer.elapsed_time();
-      CUOPT_LOG_INFO("Third party presolve time: %f", presolve_time);
+      CUOPT_LOG_INFO("Papilo presolve time: %f", presolve_time);
     }
     if (settings.user_problem_file != "") {
       CUOPT_LOG_INFO("Writing user problem to file: %s", settings.user_problem_file.c_str());

cpp/tests/CMakeLists.txt

@@ -34,7 +34,6 @@ if(BUILD_TESTS)
       cuopt
       GTest::gmock
       GTest::gtest
-      papilo-core
   )
 endif()
 

cpp/tests/examples/routing/CMakeLists.txt

@@ -48,6 +48,5 @@ foreach(target
       cuopt
       cuopttestutils
       OpenMP::OpenMP_CXX
-      papilo-core
   )
 endforeach()

docs/cuopt/source/lp-features.rst

@@ -78,6 +78,13 @@ Crossover
 Crossover allows you to obtain a high-quality basic solution from the results of a PDLP solve. More details can be found :ref:`here <crossover>`.
 
 
+Presolve
+--------
+
+Presolve procedure is applied to the problem before the solver is called. It can be used to reduce the problem size and improve solve time. It is enabled by default for MIP problems, and disabled by default for LP problems.
+Furthermore, for LP problems, when the dual solution is not needed, additional presolve procedures can be applied to further improve solve times. This is achived by turned off dual postsolve.
+
+
 Logging
 -------
 

docs/cuopt/source/lp-milp-settings.rst

@@ -65,6 +65,11 @@ Presolve
 ^^^^^^^^
 ``CUOPT_PRESOLVE`` controls whether presolve is enabled. Presolve can reduce problem size and improve solve time. Enabled by default for MIP, disabled by default for LP.
 
+Dual Postsolve
+^^^^^^^^^^^^^^
+``CUOPT_DUAL_POSTSOLVE`` controls whether dual postsolve is enabled. Disabling dual postsolve can improve solve time at the expense of not having
+access to the dual solution. Enabled by default for LP when presolve is enabled. This is not relevant for MIP problems.
+
 Linear Programming
 ------------------
 

python/cuopt/cuopt/linear_programming/solver/solver_parameters.pyx

@@ -60,6 +60,7 @@ cdef extern from "cuopt/linear_programming/constants.h": # noqa
     cdef const char* c_CUOPT_LOG_TO_CONSOLE "CUOPT_LOG_TO_CONSOLE" # noqa
     cdef const char* c_CUOPT_CROSSOVER "CUOPT_CROSSOVER" # noqa
     cdef const char* c_CUOPT_PRESOLVE "CUOPT_PRESOLVE" # noqa
+    cdef const char* c_CUOPT_DUAL_POSTSOLVE "CUOPT_DUAL_POSTSOLVE" # noqa
     cdef const char* c_CUOPT_MIP_ABSOLUTE_TOLERANCE "CUOPT_MIP_ABSOLUTE_TOLERANCE" # noqa
     cdef const char* c_CUOPT_MIP_RELATIVE_TOLERANCE "CUOPT_MIP_RELATIVE_TOLERANCE" # noqa
     cdef const char* c_CUOPT_MIP_INTEGRALITY_TOLERANCE "CUOPT_MIP_INTEGRALITY_TOLERANCE" # noqa
@@ -94,6 +95,7 @@ CUOPT_LOG_FILE = c_CUOPT_LOG_FILE.decode('utf-8') # noqa
 CUOPT_LOG_TO_CONSOLE = c_CUOPT_LOG_TO_CONSOLE.decode('utf-8') # noqa
 CUOPT_CROSSOVER = c_CUOPT_CROSSOVER.decode('utf-8') # noqa
 CUOPT_PRESOLVE = c_CUOPT_PRESOLVE.decode('utf-8') # noqa
+CUOPT_DUAL_POSTSOLVE = c_CUOPT_DUAL_POSTSOLVE.decode('utf-8') # noqa
 CUOPT_MIP_ABSOLUTE_TOLERANCE = c_CUOPT_MIP_ABSOLUTE_TOLERANCE.decode('utf-8') # noqa
 CUOPT_MIP_RELATIVE_TOLERANCE = c_CUOPT_MIP_RELATIVE_TOLERANCE.decode('utf-8') # noqa
 CUOPT_MIP_INTEGRALITY_TOLERANCE = c_CUOPT_MIP_INTEGRALITY_TOLERANCE.decode('utf-8') # noqa

python/cuopt/cuopt/linear_programming/solver_settings/solver_settings.py

@@ -21,6 +21,7 @@
     CUOPT_ABSOLUTE_PRIMAL_TOLERANCE,
     CUOPT_CROSSOVER,
     CUOPT_DUAL_INFEASIBLE_TOLERANCE,
+    CUOPT_DUAL_POSTSOLVE,
     CUOPT_FIRST_PRIMAL_FEASIBLE,
     CUOPT_INFEASIBILITY_DETECTION,
     CUOPT_ITERATION_LIMIT,
@@ -37,6 +38,7 @@
     CUOPT_NUM_CPU_THREADS,
     CUOPT_PDLP_SOLVER_MODE,
     CUOPT_PER_CONSTRAINT_RESIDUAL,
+    CUOPT_PRESOLVE,
     CUOPT_PRIMAL_INFEASIBLE_TOLERANCE,
     CUOPT_RELATIVE_DUAL_TOLERANCE,
     CUOPT_RELATIVE_GAP_TOLERANCE,
@@ -370,6 +372,8 @@ def toDict(self):
             "iteration_limit": self.get_parameter(CUOPT_ITERATION_LIMIT),
             "pdlp_solver_mode": self.get_parameter(CUOPT_PDLP_SOLVER_MODE),
             "method": self.get_parameter(CUOPT_METHOD),
+            "presolve": self.get_parameter(CUOPT_PRESOLVE),
+            "dual_postsolve": self.get_parameter(CUOPT_DUAL_POSTSOLVE),
             "mip_scaling": self.get_parameter(CUOPT_MIP_SCALING),
             "mip_heuristics_only": self.get_parameter(
                 CUOPT_MIP_HEURISTICS_ONLY