Merge pull request #534 from aliceb-nv/impl-int-store

aliceb-nv · web-flow · commit 66d652977550 · 2025-10-27T17:30:32.000+01:00
Carry over implied integers from papilo presolve reductions
diff --git a/cpp/src/linear_programming/solve.cu b/cpp/src/linear_programming/solve.cu
@@ -842,19 +842,18 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(optimization_problem_t<i_t, f
       // 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);
     }
diff --git a/cpp/src/mip/presolve/third_party_presolve.cpp b/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>
diff --git a/cpp/src/mip/presolve/third_party_presolve.hpp b/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,
diff --git a/cpp/src/mip/problem/presolve_data.cuh b/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
diff --git a/cpp/src/mip/problem/problem.cu b/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()>>>(
diff --git a/cpp/src/mip/problem/problem.cuh b/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.
diff --git a/cpp/src/mip/solve.cu b/cpp/src/mip/solve.cu
@@ -209,23 +209,27 @@ mip_solution_t<i_t, f_t> solve_mip(optimization_problem_t<i_t, f_t>& op_problem,
       // Note that this is not the presolve time, but the time limit for presolve.
       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,
-                         settings.num_cpu_threads);
-      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::MIP,
+                                     dual_postsolve,
+                                     settings.tolerances.absolute_tolerance,
+                                     settings.tolerances.relative_tolerance,
+                                     presolve_time_limit,
+                                     settings.num_cpu_threads);
+      if (!result.has_value()) {
         return mip_solution_t<i_t, f_t>(mip_termination_status_t::Infeasible,
                                         solver_stats_t<i_t, f_t>{},
                                         op_problem.get_handle_ptr()->get_stream());
       }
 
-      problem       = detail::problem_t<i_t, f_t>(reduced_op_problem);
+      problem = detail::problem_t<i_t, f_t>(result->reduced_problem);
+      problem.set_implied_integers(result->implied_integer_indices);
       presolve_time = timer.elapsed_time();
+      if (result->implied_integer_indices.size() > 0) {
+        CUOPT_LOG_INFO("%d implied integers", result->implied_integer_indices.size());
+      }
+      if (problem.is_objective_integral()) { CUOPT_LOG_INFO("Objective function is integral"); }
       CUOPT_LOG_INFO("Papilo presolve time: %f", presolve_time);
     }
     if (settings.user_problem_file != "") {
diff --git a/cpp/tests/mip/CMakeLists.txt b/cpp/tests/mip/CMakeLists.txt
@@ -43,6 +43,9 @@ ConfigureTest(UNIT_TEST
 ConfigureTest(EMPTY_FIXED_PROBLEMS_TEST
     ${CMAKE_CURRENT_SOURCE_DIR}/empty_fixed_problems_test.cu
 )
+ConfigureTest(PRESOLVE_TEST
+    ${CMAKE_CURRENT_SOURCE_DIR}/presolve_test.cu
+)
 # Disable for now
 # ConfigureTest(FEASIBILITY_JUMP_TEST
 #    ${CMAKE_CURRENT_SOURCE_DIR}/feasibility_jump_tests.cu
diff --git a/cpp/tests/mip/presolve_test.cu b/cpp/tests/mip/presolve_test.cu
@@ -0,0 +1,69 @@
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2024-2025 NVIDIA CORPORATION & AFFILIATES. All rights
+ * reserved. SPDX-License-Identifier: Apache-2.0
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "../linear_programming/utilities/pdlp_test_utilities.cuh"
+
+#include <cuopt/linear_programming/solve.hpp>
+#include <linear_programming/utils.cuh>
+#include <mip/presolve/third_party_presolve.hpp>
+#include <mip/problem/problem.cuh>
+#include <mps_parser/mps_data_model.hpp>
+#include <mps_parser/parser.hpp>
+#include <utilities/common_utils.hpp>
+#include <utilities/copy_helpers.hpp>
+#include <utilities/error.hpp>
+
+#include <raft/core/handle.hpp>
+#include <raft/util/cudart_utils.hpp>
+
+#include <gtest/gtest.h>
+
+#include <cstdint>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace cuopt::linear_programming::test {
+
+TEST(problem, find_implied_integers)
+{
+  const raft::handle_t handle_{};
+
+  auto path           = make_path_absolute("mip/fiball.mps");
+  auto mps_data_model = cuopt::mps_parser::parse_mps<int, double>(path, false);
+  auto op_problem     = mps_data_model_to_optimization_problem(&handle_, mps_data_model);
+  auto presolver      = std::make_unique<detail::third_party_presolve_t<int, double>>();
+  auto result         = presolver->apply(
+    op_problem, cuopt::linear_programming::problem_category_t::MIP, false, 1e-6, 1e-12, 20, 1);
+  ASSERT_TRUE(result.has_value());
+
+  auto problem = detail::problem_t<int, double>(result->reduced_problem);
+  problem.set_implied_integers(result->implied_integer_indices);
+  ASSERT_TRUE(result->implied_integer_indices.size() > 0);
+  auto var_types = host_copy(problem.variable_types);
+  // Find the index of the one continuous variable
+  auto it = std::find_if(var_types.begin(), var_types.end(), [](var_t var_type) {
+    return var_type == var_t::CONTINUOUS;
+  });
+  ASSERT_NE(it, var_types.end());
+  ASSERT_EQ(problem.presolve_data.var_flags.size(), var_types.size());
+  // Ensure it is an implied integer
+  EXPECT_EQ(problem.presolve_data.var_flags.element(it - var_types.begin(), handle_.get_stream()),
+            ((int)detail::problem_t<int, double>::var_flags_t::VAR_IMPLIED_INTEGER));
+}
+
+}  // namespace cuopt::linear_programming::test

Original file line number	Diff line number	Diff line change
`@@ -39,7 +39,8 @@ class presolve_data_t {`
`39`	`39`	`objective_offset(problem.get_objective_offset()),`
`40`	`40`	`objective_scaling_factor(problem.get_objective_scaling_factor()),`
`41`	`41`	`variable_mapping(0, stream),`
`42`		`- fixed_var_assignment(0, stream)`
	`42`	`+ fixed_var_assignment(0, stream),`
	`43`	`+ var_flags(0, stream)`
`43`	`44`	`{`
`44`	`45`	`}`
`45`	`46`
`@@ -50,7 +51,8 @@ class presolve_data_t {`
`50`	`51`	`objective_offset(other.objective_offset),`
`51`	`52`	`objective_scaling_factor(other.objective_scaling_factor),`
`52`	`53`	`variable_mapping(other.variable_mapping, stream),`
`53`		`- fixed_var_assignment(other.fixed_var_assignment, stream)`
	`54`	`+ fixed_var_assignment(other.fixed_var_assignment, stream),`
	`55`	`+ var_flags(other.var_flags, stream)`
`54`	`56`	`{`
`55`	`57`	`}`
`56`	`58`
`@@ -86,6 +88,7 @@ class presolve_data_t {`
`86`	`88`
`87`	`89`	`rmm::device_uvector<i_t> variable_mapping;`
`88`	`90`	`rmm::device_uvector<f_t> fixed_var_assignment;`
	`91`	`+ rmm::device_uvector<i_t> var_flags;`
`89`	`92`	`};`
`90`	`93`
`91`	`94`	`} // namespace linear_programming::detail`