Fix inversion of constraint bounds in conditional bounds presolve (#75)

Updating the constraint bounds based on minimal and maximal activity can sometimes cause the update constraint lower bounds to be greater than constraint upper bounds even if the constraint is feasible within tolerance.

This PR fixes the constraint bounds to the same value if this happens.

Authors:
  - Kumar Aatish (https://github.com/kaatish)

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

URL: #75

Author: kaatish

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) -> i_t {
+                     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_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