Merge branch 'main' into reuse-basis-factorization

Author: nguidotti

benchmarks/linear_programming/utils/get_datasets.py

@@ -693,7 +693,7 @@ def extract(file, dir, type):
         raise Exception(f"Unknown file extension found for extraction {file}")
     # download emps and compile
     # Disable emps for now
-    if type == "netlib":
+    if type == "netlib" and False:
         url = MittelmannInstances["emps"]
         file = os.path.join(dir, "emps.c")
         download(url, file)

cpp/src/dual_simplex/barrier.cu

@@ -3575,6 +3575,9 @@ lp_status_t barrier_solver_t<i_t, f_t>::solve(f_t start_time,
   } catch (const raft::cuda_error& e) {
     settings.log.debug("Error in barrier_solver_t: %s\n", e.what());
     return lp_status_t::NUMERICAL_ISSUES;
+  } catch (const rmm::out_of_memory& e) {
+    settings.log.debug("Out of memory in barrier_solver_t: %s\n", e.what());
+    return lp_status_t::NUMERICAL_ISSUES;
   }
 }
 

cpp/src/linear_programming/solve.cu

@@ -28,6 +28,7 @@
 #include <mip/presolve/third_party_presolve.hpp>
 #include <mip/presolve/trivial_presolve.cuh>
 #include <mip/solver.cuh>
+#include <mip/utilities/sort_csr.cuh>
 
 #include <cuopt/linear_programming/pdlp/pdlp_hyper_params.cuh>
 #include <cuopt/linear_programming/pdlp/solver_settings.hpp>
@@ -444,7 +445,7 @@ optimization_problem_solution_t<i_t, f_t> run_barrier(
 {
   // Convert data structures to dual simplex format and back
   dual_simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
-    cuopt_problem_to_simplex_problem<i_t, f_t>(problem);
+    cuopt_problem_to_simplex_problem<i_t, f_t>(problem.handle_ptr, problem);
   auto sol_dual_simplex = run_barrier(dual_simplex_problem, settings, timer);
   return convert_dual_simplex_sol(problem,
                                   std::get<0>(sol_dual_simplex),
@@ -515,7 +516,7 @@ optimization_problem_solution_t<i_t, f_t> run_dual_simplex(
 {
   // Convert data structures to dual simplex format and back
   dual_simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
-    cuopt_problem_to_simplex_problem<i_t, f_t>(problem);
+    cuopt_problem_to_simplex_problem<i_t, f_t>(problem.handle_ptr, problem);
   auto sol_dual_simplex = run_dual_simplex(dual_simplex_problem, settings, timer);
   return convert_dual_simplex_sol(problem,
                                   std::get<0>(sol_dual_simplex),
@@ -671,16 +672,14 @@ optimization_problem_solution_t<i_t, f_t> run_concurrent(
   // Initialize the dual simplex structures before we run PDLP.
   // Otherwise, CUDA API calls to the problem stream may occur in both threads and throw graph
   // capture off
-  auto barrier_handle = raft::handle_t(*op_problem.get_handle_ptr());
-  detail::problem_t<i_t, f_t> d_barrier_problem(problem);
+  auto barrier_handle                  = raft::handle_t(*op_problem.get_handle_ptr());
   rmm::cuda_stream_view barrier_stream = rmm::cuda_stream_per_thread;
-  d_barrier_problem.handle_ptr         = &barrier_handle;
   raft::resource::set_cuda_stream(barrier_handle, barrier_stream);
   // Make sure allocations are done on the original stream
   problem.handle_ptr->sync_stream();
 
   dual_simplex::user_problem_t<i_t, f_t> dual_simplex_problem =
-    cuopt_problem_to_simplex_problem<i_t, f_t>(d_barrier_problem);
+    cuopt_problem_to_simplex_problem<i_t, f_t>(&barrier_handle, problem);
   // Create a thread for dual simplex
   std::unique_ptr<
     std::tuple<dual_simplex::lp_solution_t<i_t, f_t>, dual_simplex::lp_status_t, f_t, f_t, f_t>>
@@ -837,24 +836,24 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(optimization_problem_t<i_t, f
     if (!run_presolve) { CUOPT_LOG_INFO("Third-party presolve is disabled, skipping"); }
 
     if (run_presolve) {
+      detail::sort_csr(op_problem);
       // 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.
       // But no less than 1 second, to avoid early timeout triggering known crashes
       const double presolve_time_limit =
         std::max(1.0, 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);
-      if (!feasible) {
+      presolver   = std::make_unique<detail::third_party_presolve_t<i_t, f_t>>();
+      auto result = 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);
+      if (!result.has_value()) {
         return optimization_problem_solution_t<i_t, f_t>(
           pdlp_termination_status_t::PrimalInfeasible, op_problem.get_handle_ptr()->get_stream());
       }
-      problem       = detail::problem_t<i_t, f_t>(reduced_problem);
+      problem       = detail::problem_t<i_t, f_t>(result->reduced_problem);
       presolve_time = lp_timer.elapsed_time();
       CUOPT_LOG_INFO("Papilo presolve time: %f", presolve_time);
     }

cpp/src/linear_programming/translate.hpp

@@ -28,21 +28,21 @@ namespace cuopt::linear_programming {
 
 template <typename i_t, typename f_t>
 static dual_simplex::user_problem_t<i_t, f_t> cuopt_problem_to_simplex_problem(
-  detail::problem_t<i_t, f_t>& model)
+  raft::handle_t const* handle_ptr, detail::problem_t<i_t, f_t>& model)
 {
-  dual_simplex::user_problem_t<i_t, f_t> user_problem(model.handle_ptr);
+  dual_simplex::user_problem_t<i_t, f_t> user_problem(handle_ptr);
 
   int m                  = model.n_constraints;
   int n                  = model.n_variables;
   int nz                 = model.nnz;
   user_problem.num_rows  = m;
   user_problem.num_cols  = n;
-  user_problem.objective = cuopt::host_copy(model.objective_coefficients);
+  user_problem.objective = cuopt::host_copy(model.objective_coefficients, handle_ptr->get_stream());
 
   dual_simplex::csr_matrix_t<i_t, f_t> csr_A(m, n, nz);
-  csr_A.x         = cuopt::host_copy(model.coefficients);
-  csr_A.j         = cuopt::host_copy(model.variables);
-  csr_A.row_start = cuopt::host_copy(model.offsets);
+  csr_A.x         = cuopt::host_copy(model.coefficients, handle_ptr->get_stream());
+  csr_A.j         = cuopt::host_copy(model.variables, handle_ptr->get_stream());
+  csr_A.row_start = cuopt::host_copy(model.offsets, handle_ptr->get_stream());
 
   csr_A.to_compressed_col(user_problem.A);
 
@@ -51,8 +51,10 @@ static dual_simplex::user_problem_t<i_t, f_t> cuopt_problem_to_simplex_problem(
   user_problem.range_rows.clear();
   user_problem.range_value.clear();
 
-  auto model_constraint_lower_bounds = cuopt::host_copy(model.constraint_lower_bounds);
-  auto model_constraint_upper_bounds = cuopt::host_copy(model.constraint_upper_bounds);
+  auto model_constraint_lower_bounds =
+    cuopt::host_copy(model.constraint_lower_bounds, handle_ptr->get_stream());
+  auto model_constraint_upper_bounds =
+    cuopt::host_copy(model.constraint_upper_bounds, handle_ptr->get_stream());
 
   // All constraints have lower and upper bounds
   // lr <= a_i^T x <= ur
@@ -79,7 +81,7 @@ static dual_simplex::user_problem_t<i_t, f_t> cuopt_problem_to_simplex_problem(
   }
   user_problem.num_range_rows = user_problem.range_rows.size();
   std::tie(user_problem.lower, user_problem.upper) =
-    extract_host_bounds<f_t>(model.variable_bounds, model.handle_ptr);
+    extract_host_bounds<f_t>(model.variable_bounds, handle_ptr);
   user_problem.problem_name = model.original_problem_ptr->get_problem_name();
   if (model.row_names.size() > 0) {
     user_problem.row_names.resize(m);
@@ -97,7 +99,7 @@ static dual_simplex::user_problem_t<i_t, f_t> cuopt_problem_to_simplex_problem(
   user_problem.obj_scale    = model.presolve_data.objective_scaling_factor;
   user_problem.var_types.resize(n);
 
-  auto model_variable_types = cuopt::host_copy(model.variable_types);
+  auto model_variable_types = cuopt::host_copy(model.variable_types, handle_ptr->get_stream());
   for (int j = 0; j < n; ++j) {
     user_problem.var_types[j] =
       model_variable_types[j] == var_t::CONTINUOUS

cpp/src/mip/presolve/third_party_presolve.cpp

@@ -391,7 +391,7 @@ void set_presolve_parameters(papilo::Presolve<f_t>& presolver,
 }
 
 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(
+std::optional<third_party_presolve_result_t<i_t, f_t>> 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,
@@ -423,7 +423,7 @@ std::pair<optimization_problem_t<i_t, f_t>, bool> third_party_presolve_t<i_t, f_
   check_presolve_status(result.status);
   if (result.status == papilo::PresolveStatus::kInfeasible ||
       result.status == papilo::PresolveStatus::kUnbndOrInfeas) {
-    return std::make_pair(optimization_problem_t<i_t, f_t>(op_problem.get_handle_ptr()), false);
+    return std::nullopt;
   }
   post_solve_storage_ = result.postsolve;
   CUOPT_LOG_INFO("Presolve removed: %d constraints, %d variables, %d nonzeros",
@@ -435,8 +435,16 @@ std::pair<optimization_problem_t<i_t, f_t>, bool> third_party_presolve_t<i_t, f_
                  papilo_problem.getNCols(),
                  papilo_problem.getConstraintMatrix().getNnz());
 
-  return std::make_pair(
-    build_optimization_problem<i_t, f_t>(papilo_problem, op_problem.get_handle_ptr()), true);
+  auto opt_problem =
+    build_optimization_problem<i_t, f_t>(papilo_problem, op_problem.get_handle_ptr());
+  auto col_flags = papilo_problem.getColFlags();
+  std::vector<i_t> implied_integer_indices;
+  for (size_t i = 0; i < col_flags.size(); i++) {
+    if (col_flags[i].test(papilo::ColFlag::kImplInt)) implied_integer_indices.push_back(i);
+  }
+
+  return std::make_optional(
+    third_party_presolve_result_t<i_t, f_t>{opt_problem, implied_integer_indices});
 }
 
 template <typename i_t, typename f_t>

cpp/src/mip/presolve/third_party_presolve.hpp

@@ -17,16 +17,25 @@
 
 #pragma once
 
+#include <optional>
+
 #include <cuopt/linear_programming/optimization_problem.hpp>
 
 namespace cuopt::linear_programming::detail {
 
+template <typename i_t, typename f_t>
+struct third_party_presolve_result_t {
+  optimization_problem_t<i_t, f_t> reduced_problem;
+  std::vector<i_t> implied_integer_indices;
+  // clique info, etc...
+};
+
 template <typename i_t, typename f_t>
 class third_party_presolve_t {
  public:
   third_party_presolve_t() = default;
 
-  std::pair<optimization_problem_t<i_t, f_t>, bool> apply(
+  std::optional<third_party_presolve_result_t<i_t, f_t>> apply(
     optimization_problem_t<i_t, f_t> const& op_problem,
     problem_category_t category,
     bool dual_postsolve,

cpp/src/mip/problem/presolve_data.cuh

@@ -39,7 +39,8 @@ class presolve_data_t {
       objective_offset(problem.get_objective_offset()),
       objective_scaling_factor(problem.get_objective_scaling_factor()),
       variable_mapping(0, stream),
-      fixed_var_assignment(0, stream)
+      fixed_var_assignment(0, stream),
+      var_flags(0, stream)
   {
   }
 
@@ -50,7 +51,8 @@ class presolve_data_t {
       objective_offset(other.objective_offset),
       objective_scaling_factor(other.objective_scaling_factor),
       variable_mapping(other.variable_mapping, stream),
-      fixed_var_assignment(other.fixed_var_assignment, stream)
+      fixed_var_assignment(other.fixed_var_assignment, stream),
+      var_flags(other.var_flags, stream)
   {
   }
 
@@ -86,6 +88,7 @@ class presolve_data_t {
 
   rmm::device_uvector<i_t> variable_mapping;
   rmm::device_uvector<f_t> fixed_var_assignment;
+  rmm::device_uvector<i_t> var_flags;
 };
 
 }  // namespace linear_programming::detail

cpp/src/mip/problem/problem.cu

@@ -89,6 +89,11 @@ void problem_t<i_t, f_t>::op_problem_cstr_body(const optimization_problem_t<i_t,
   // If maximization problem, convert the problem
   if (maximize) convert_to_maximization_problem(*this);
   if (is_mip) {
+    presolve_data.var_flags.resize(n_variables, handle_ptr->get_stream());
+    thrust::fill(handle_ptr->get_thrust_policy(),
+                 presolve_data.var_flags.begin(),
+                 presolve_data.var_flags.end(),
+                 0);
     integer_indices.resize(n_variables, handle_ptr->get_stream());
     is_binary_variable.resize(n_variables, handle_ptr->get_stream());
     compute_n_integer_vars();
@@ -190,6 +195,7 @@ problem_t<i_t, f_t>::problem_t(const problem_t<i_t, f_t>& problem_)
     objective_name(problem_.objective_name),
     is_scaled_(problem_.is_scaled_),
     preprocess_called(problem_.preprocess_called),
+    objective_is_integral(problem_.objective_is_integral),
     lp_state(problem_.lp_state),
     fixing_helpers(problem_.fixing_helpers, handle_ptr),
     vars_with_objective_coeffs(problem_.vars_with_objective_coeffs),
@@ -284,6 +290,7 @@ problem_t<i_t, f_t>::problem_t(const problem_t<i_t, f_t>& problem_, bool no_deep
     objective_name(problem_.objective_name),
     is_scaled_(problem_.is_scaled_),
     preprocess_called(problem_.preprocess_called),
+    objective_is_integral(problem_.objective_is_integral),
     lp_state(problem_.lp_state),
     fixing_helpers(problem_.fixing_helpers, handle_ptr),
     vars_with_objective_coeffs(problem_.vars_with_objective_coeffs),
@@ -933,6 +940,8 @@ typename problem_t<i_t, f_t>::view_t problem_t<i_t, f_t>::view()
   v.variable_types = raft::device_span<var_t>{variable_types.data(), variable_types.size()};
   v.is_binary_variable =
     raft::device_span<i_t>{is_binary_variable.data(), is_binary_variable.size()};
+  v.var_flags =
+    raft::device_span<i_t>{presolve_data.var_flags.data(), presolve_data.var_flags.size()};
   v.related_variables = raft::device_span<i_t>{related_variables.data(), related_variables.size()};
   v.related_variables_offsets =
     raft::device_span<i_t>{related_variables_offsets.data(), related_variables_offsets.size()};
@@ -953,6 +962,7 @@ void problem_t<i_t, f_t>::resize_variables(size_t size)
   variable_types.resize(size, handle_ptr->get_stream());
   objective_coefficients.resize(size, handle_ptr->get_stream());
   is_binary_variable.resize(size, handle_ptr->get_stream());
+  presolve_data.var_flags.resize(size, handle_ptr->get_stream());  // 0 is default - no flag
   related_variables_offsets.resize(size, handle_ptr->get_stream());
 }
 
@@ -1060,6 +1070,32 @@ void problem_t<i_t, f_t>::insert_constraints(constraints_delta_t<i_t, f_t>& h_co
   combine_constraint_bounds<i_t, f_t>(*this, combined_bounds);
 }
 
+template <typename i_t, typename f_t>
+void problem_t<i_t, f_t>::set_implied_integers(const std::vector<i_t>& implied_integer_indices)
+{
+  raft::common::nvtx::range fun_scope("set_implied_integers");
+  auto d_indices = cuopt::device_copy(implied_integer_indices, handle_ptr->get_stream());
+  print("implied integer indices", d_indices);
+  thrust::for_each(handle_ptr->get_thrust_policy(),
+                   d_indices.begin(),
+                   d_indices.end(),
+                   [var_flags = make_span(presolve_data.var_flags),
+                    var_types = make_span(variable_types)] __device__(i_t idx) {
+                     cuopt_assert(idx < var_flags.size(), "Index out of bounds");
+                     cuopt_assert(var_types[idx] == var_t::CONTINUOUS, "Variable is integer");
+                     var_flags[idx] |= (i_t)VAR_IMPLIED_INTEGER;
+                   });
+  objective_is_integral = thrust::all_of(handle_ptr->get_thrust_policy(),
+                                         thrust::make_counting_iterator(0),
+                                         thrust::make_counting_iterator(n_variables),
+                                         [v = view()] __device__(i_t var_idx) {
+                                           if (v.objective_coefficients[var_idx] == 0) return true;
+                                           return v.is_integer(v.objective_coefficients[var_idx]) &&
+                                                  (v.variable_types[var_idx] == var_t::INTEGER ||
+                                                   (v.var_flags[var_idx] & VAR_IMPLIED_INTEGER));
+                                         });
+}
+
 template <typename i_t, typename f_t>
 void problem_t<i_t, f_t>::fix_given_variables(problem_t<i_t, f_t>& original_problem,
                                               rmm::device_uvector<f_t>& assignment,
@@ -1249,6 +1285,13 @@ void problem_t<i_t, f_t>::remove_given_variables(problem_t<i_t, f_t>& original_p
                  original_problem.variable_types.begin(),
                  variable_types.begin());
   variable_types.resize(variable_map.size(), handle_ptr->get_stream());
+  // keep implied-integer and other flags consistent with new variable set
+  thrust::gather(handle_ptr->get_thrust_policy(),
+                 variable_map.begin(),
+                 variable_map.end(),
+                 original_problem.presolve_data.var_flags.begin(),
+                 presolve_data.var_flags.begin());
+  presolve_data.var_flags.resize(variable_map.size(), handle_ptr->get_stream());
   const i_t TPB = 64;
   // compute new offsets
   compute_new_offsets<i_t, f_t><<<variable_map.size(), TPB, 0, handle_ptr->get_stream()>>>(

cpp/src/mip/problem/problem.cuh

@@ -91,6 +91,7 @@ class problem_t {
 
   void insert_variables(variables_delta_t<i_t, f_t>& h_vars);
   void insert_constraints(constraints_delta_t<i_t, f_t>& h_constraints);
+  void set_implied_integers(const std::vector<i_t>& implied_integer_indices);
   void resize_variables(size_t size);
   void resize_constraints(size_t matrix_size, size_t constraint_size, size_t var_size);
   void preprocess_problem();
@@ -100,6 +101,7 @@ class problem_t {
   void post_process_solution(solution_t<i_t, f_t>& solution);
   void compute_transpose_of_problem();
   f_t get_user_obj_from_solver_obj(f_t solver_obj) const;
+  bool is_objective_integral() const { return objective_is_integral; }
   void compute_integer_fixed_problem();
   void fill_integer_fixed_problem(rmm::device_uvector<f_t>& assignment,
                                   const raft::handle_t* handle_ptr);
@@ -116,6 +118,10 @@ class problem_t {
   void compute_vars_with_objective_coeffs();
   void test_problem_fixing_time();
 
+  enum var_flags_t : i_t {
+    VAR_IMPLIED_INTEGER = 1 << 0,
+  };
+
   struct view_t {
     HDI std::pair<i_t, i_t> reverse_range_for_var(i_t v) const
     {
@@ -193,6 +199,7 @@ class problem_t {
     raft::device_span<f_t> constraint_upper_bounds;
     raft::device_span<var_t> variable_types;
     raft::device_span<i_t> is_binary_variable;
+    raft::device_span<i_t> var_flags;
     raft::device_span<i_t> integer_indices;
     raft::device_span<i_t> binary_indices;
     raft::device_span<i_t> nonbinary_indices;
@@ -269,6 +276,7 @@ class problem_t {
   f_t objective_offset;
   bool is_scaled_{false};
   bool preprocess_called{false};
+  bool objective_is_integral{false};
   // this LP state keeps the warm start data of some solution of
   // 1. Original problem: it is unchanged and part of it is used
   // to warm start slightly modified problems.