Fix inversion of constraint bounds in conditional bounds presolve

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

@@ -70,11 +70,12 @@ class mip_solver_settings_t {
   bool has_initial_solution() const;
 
   struct tolerances_t {
-    f_t absolute_tolerance    = 1.0e-4;
-    f_t relative_tolerance    = 1.0e-6;
-    f_t integrality_tolerance = 1.0e-5;
-    f_t absolute_mip_gap      = 1.0e-10;
-    f_t relative_mip_gap      = 1.0e-4;
+    f_t presolve_absolute_tolerance = 1.0e-6;
+    f_t absolute_tolerance          = 1.0e-4;
+    f_t relative_tolerance          = 1.0e-6;
+    f_t integrality_tolerance       = 1.0e-5;
+    f_t absolute_mip_gap            = 1.0e-10;
+    f_t relative_mip_gap            = 1.0e-4;
   };
 
   /**

cpp/src/mip/diversity/diversity_manager.cu

@@ -246,15 +246,15 @@ bool diversity_manager_t<i_t, f_t>::run_presolve(f_t time_limit)
   if (termination_criterion_t::NO_UPDATE != term_crit) {
     ls.constraint_prop.bounds_update.set_updated_bounds(*problem_ptr);
     trivial_presolve(*problem_ptr);
-    if (!problem_ptr->empty) { check_bounds_sanity(*problem_ptr); }
+    if (!problem_ptr->empty && !check_bounds_sanity(*problem_ptr)) { return false; }
   }
   if (!problem_ptr->empty) {
     // do the resizing no-matter what, bounds presolve might not change the bounds but initial
     // trivial presolve might have
     ls.constraint_prop.bounds_update.resize(*problem_ptr);
     ls.constraint_prop.conditional_bounds_update.update_constraint_bounds(
       *problem_ptr, ls.constraint_prop.bounds_update);
-    check_bounds_sanity(*problem_ptr);
+    if (!check_bounds_sanity(*problem_ptr)) { return false; }
   }
   stats.presolve_time = presolve_timer.elapsed_time();
   return true;

cpp/src/mip/presolve/bounds_presolve.cu

@@ -20,7 +20,9 @@
 
 #include <thrust/count.h>
 #include <thrust/extrema.h>
+#include <thrust/functional.h>
 #include <thrust/iterator/zip_iterator.h>
+#include <thrust/transform_reduce.h>
 #include <thrust/tuple.h>
 #include <utilities/copy_helpers.hpp>
 #include <utilities/device_utils.cuh>
@@ -350,18 +352,29 @@ void bound_presolve_t<i_t, f_t>::calc_and_set_updated_constraint_bounds(problem_
 {
   calculate_activity_on_problem_bounds(pb);
 
-  thrust::transform(pb.handle_ptr->get_thrust_policy(),
-                    upd.max_activity.begin(),
-                    upd.max_activity.end(),
-                    pb.constraint_upper_bounds.begin(),
-                    pb.constraint_upper_bounds.begin(),
-                    thrust::minimum<f_t>());
-  thrust::transform(pb.handle_ptr->get_thrust_policy(),
-                    upd.min_activity.begin(),
-                    upd.min_activity.end(),
-                    pb.constraint_lower_bounds.begin(),
-                    pb.constraint_lower_bounds.begin(),
-                    thrust::maximum<f_t>());
+  thrust::for_each(pb.handle_ptr->get_thrust_policy(),
+                   thrust::make_counting_iterator(0),
+                   thrust::make_counting_iterator(pb.n_constraints),
+                   [pb      = pb.view(),
+                    min_act = make_span(upd.min_activity),
+                    max_act = make_span(upd.max_activity),
+                    cnst_lb = make_span(pb.constraint_lower_bounds),
+                    cnst_ub = make_span(pb.constraint_upper_bounds)] __device__(i_t idx) {
+                     auto min_a    = min_act[idx];
+                     auto max_a    = max_act[idx];
+                     auto c_lb     = cnst_lb[idx];
+                     auto c_ub     = cnst_ub[idx];
+                     auto new_c_lb = max(c_lb, min_a);
+                     auto new_c_ub = min(c_ub, max_a);
+                     i_t infeas    = check_infeasibility<i_t, f_t>(
+                       min_a, max_a, new_c_lb, new_c_ub, pb.tolerances.presolve_absolute_tolerance);
+                     if (!infeas && (new_c_lb > new_c_ub)) {
+                       new_c_lb = (new_c_lb + new_c_ub) / 2;
+                       new_c_ub = new_c_lb;
+                     }
+                     cnst_lb[idx] = new_c_lb;
+                     cnst_ub[idx] = new_c_ub;
+                   });
 }
 
 #if MIP_INSTANTIATE_FLOAT

cpp/src/mip/presolve/bounds_update_helpers.cuh

@@ -142,6 +142,12 @@ __global__ void calc_activity_kernel(typename problem_t<i_t, f_t>::view_t pb,
 
 // Update bounds
 
+template <typename i_t, typename f_t>
+inline __device__ bool check_infeasibility(f_t min_a, f_t max_a, f_t cnst_lb, f_t cnst_ub, f_t eps)
+{
+  return (min_a > cnst_ub + eps) || (max_a < cnst_lb - eps);
+}
+
 template <typename i_t, typename f_t>
 inline __device__ bool check_infeasibility(
   f_t min_a, f_t max_a, f_t cnst_lb, f_t cnst_ub, f_t abs_tol, f_t rel_tol)

cpp/src/mip/problem/problem.cu

@@ -388,23 +388,25 @@ void problem_t<i_t, f_t>::check_problem_representation(bool check_transposed,
                "Sizes for vectors related to the constraints are not the same.");
 
   // Check the validity of bounds
-  cuopt_assert(
+  cuopt_expects(
     thrust::all_of(handle_ptr->get_thrust_policy(),
                    thrust::make_counting_iterator<i_t>(0),
                    thrust::make_counting_iterator<i_t>(n_variables),
                    [variable_lower_bounds = variable_lower_bounds.data(),
                     variable_upper_bounds = variable_upper_bounds.data()] __device__(i_t idx) {
                      return variable_lower_bounds[idx] <= variable_upper_bounds[idx];
                    }),
+    error_type_t::ValidationError,
     "Variable bounds are invalid");
-  cuopt_assert(
+  cuopt_expects(
     thrust::all_of(handle_ptr->get_thrust_policy(),
                    thrust::make_counting_iterator<i_t>(0),
                    thrust::make_counting_iterator<i_t>(n_constraints),
                    [constraint_lower_bounds = constraint_lower_bounds.data(),
                     constraint_upper_bounds = constraint_upper_bounds.data()] __device__(i_t idx) {
                      return constraint_lower_bounds[idx] <= constraint_upper_bounds[idx];
                    }),
+    error_type_t::ValidationError,
     "Constraints bounds are invalid");
 
   if (check_mip_related_data) {

cpp/src/mip/problem/problem_helpers.cuh

@@ -248,42 +248,40 @@ static void check_csr_representation([[maybe_unused]] const rmm::device_uvector<
 }
 
 template <typename i_t, typename f_t>
-static void check_var_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
+static bool check_var_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
 {
   bool crossing_bounds_detected =
     thrust::any_of(problem.handle_ptr->get_thrust_policy(),
                    thrust::counting_iterator(0),
                    thrust::counting_iterator((i_t)problem.variable_lower_bounds.size()),
-                   [lb = make_span(problem.variable_lower_bounds),
-                    ub = make_span(problem.variable_upper_bounds)] __device__(i_t index) {
-                     return lb[index] > ub[index];
+                   [tolerance = problem.tolerances.presolve_absolute_tolerance,
+                    lb        = make_span(problem.variable_lower_bounds),
+                    ub        = make_span(problem.variable_upper_bounds)] __device__(i_t index) {
+                     return (lb[index] > ub[index] + tolerance);
                    });
-  cuopt_expects(!crossing_bounds_detected,
-                error_type_t::ValidationError,
-                "There shouldn't be any crossing bounds in variable bounds.");
+  return !crossing_bounds_detected;
 }
 
 template <typename i_t, typename f_t>
-static void check_constraint_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
+static bool check_constraint_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
 {
   bool crossing_bounds_detected =
     thrust::any_of(problem.handle_ptr->get_thrust_policy(),
                    thrust::counting_iterator(0),
                    thrust::counting_iterator((i_t)problem.constraint_lower_bounds.size()),
-                   [lb = make_span(problem.constraint_lower_bounds),
-                    ub = make_span(problem.constraint_upper_bounds)] __device__(i_t index) {
-                     return lb[index] > ub[index];
+                   [tolerance = problem.tolerances.presolve_absolute_tolerance,
+                    lb        = make_span(problem.constraint_lower_bounds),
+                    ub        = make_span(problem.constraint_upper_bounds)] __device__(i_t index) {
+                     return (lb[index] > ub[index] + tolerance);
                    });
-  cuopt_expects(!crossing_bounds_detected,
-                error_type_t::ValidationError,
-                "There shouldn't be any crossing bounds in constraints bounds.");
+  return !crossing_bounds_detected;
 }
 
 template <typename i_t, typename f_t>
-static void check_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
+static bool check_bounds_sanity(const detail::problem_t<i_t, f_t>& problem)
 {
-  check_var_bounds_sanity<i_t, f_t>(problem);
-  check_constraint_bounds_sanity<i_t, f_t>(problem);
+  return check_var_bounds_sanity<i_t, f_t>(problem) &&
+         check_constraint_bounds_sanity<i_t, f_t>(problem);
 }
 
 }  // namespace cuopt::linear_programming::detail