Rounding and crash fixes (#670)

This PR fixes errouneous logic in the simple rounding and edge cases causing supposedly integer solutions to contain fractional variables. RINS fixes are also included (wrong problem pointer being used for some checks). 

Also, it seems that the cuSparse CSR-to-CSC routine occasionally produces an incorrect output transpose. As a workaround, a manual transpose implementation is used instead. Further work will be required to identify a more easily reproducible example and identify the underlying issue (and/or file a cuSparse bug)


## Summary by CodeRabbit

* **Bug Fixes**
  * Improved numerical robustness with tighter clamping, added bounds checks, and reduced spurious runtime assertions.
* **Performance Improvements**
  * Faster matrix transpose and improved GPU rounding kernels for better throughput.
* **Stability**
  * Consistent use of copied problem state during heuristics and safer handling of bound/rounding edge cases.
* **Behavioral**
  * Rounding and feasibility flows now enforce post-clamp validation and per-variable locking/decision logic.

<sub>✏️ Tip: You can customize this high-level summary in your review settings.</sub>

Authors:
  - Alice Boucher (https://github.com/aliceb-nv)

Approvers:
  - Rajesh Gandham (https://github.com/rg20)

URL: #670

Author: aliceb-nv

cpp/src/mip/diversity/lns/rins.cu

@@ -99,13 +99,16 @@ void rins_t<i_t, f_t>::run_rins()
 
   if (!dm.population.is_feasible()) return;
 
-  cuopt_assert(lp_optimal_solution.size() == problem_ptr->n_variables, "Assignment size mismatch");
+  cuopt_assert(lp_optimal_solution.size() == problem_copy->n_variables, "Assignment size mismatch");
   cuopt_assert(problem_copy->handle_ptr == &rins_handle, "Handle mismatch");
-  cuopt_assert(problem_copy->n_variables == problem_ptr->n_variables, "Problem size mismatch");
-  cuopt_assert(problem_copy->n_constraints == problem_ptr->n_constraints, "Problem size mismatch");
-  cuopt_assert(problem_copy->n_integer_vars == problem_ptr->n_integer_vars,
-               "Problem size mismatch");
-  cuopt_assert(problem_copy->n_binary_vars == problem_ptr->n_binary_vars, "Problem size mismatch");
+  // Do not make assertions based on problem_ptr. The original problem may have been modified within
+  // the FP loop relaxing integers cuopt_assert(problem_copy->n_variables ==
+  // problem_ptr->n_variables, "Problem size mismatch"); cuopt_assert(problem_copy->n_constraints ==
+  // problem_ptr->n_constraints, "Problem size mismatch"); cuopt_assert(problem_copy->n_integer_vars
+  // == problem_ptr->n_integer_vars,
+  //              "Problem size mismatch");
+  // cuopt_assert(problem_copy->n_binary_vars == problem_ptr->n_binary_vars, "Problem size
+  // mismatch");
   cuopt_assert(dm.population.current_size() > 0, "No solutions in population");
 
   solution_t<i_t, f_t> best_sol(*problem_copy);
@@ -128,21 +131,21 @@ void rins_t<i_t, f_t>::run_rins()
 
   i_t sol_size_before_rins = best_sol.assignment.size();
   auto lp_opt_device = cuopt::device_copy(this->lp_optimal_solution, rins_handle.get_stream());
-  cuopt_assert(lp_opt_device.size() == problem_ptr->n_variables, "Assignment size mismatch");
-  cuopt_assert(best_sol.assignment.size() == problem_ptr->n_variables, "Assignment size mismatch");
+  cuopt_assert(lp_opt_device.size() == problem_copy->n_variables, "Assignment size mismatch");
+  cuopt_assert(best_sol.assignment.size() == problem_copy->n_variables, "Assignment size mismatch");
 
-  rmm::device_uvector<i_t> vars_to_fix(problem_ptr->n_integer_vars, rins_handle.get_stream());
+  rmm::device_uvector<i_t> vars_to_fix(problem_copy->n_integer_vars, rins_handle.get_stream());
   auto end = thrust::copy_if(rins_handle.get_thrust_policy(),
-                             problem_ptr->integer_indices.begin(),
-                             problem_ptr->integer_indices.end(),
+                             problem_copy->integer_indices.begin(),
+                             problem_copy->integer_indices.end(),
                              vars_to_fix.begin(),
                              [lpopt     = lp_opt_device.data(),
-                              pb        = problem_ptr->view(),
+                              pb        = problem_copy->view(),
                               incumbent = best_sol.assignment.data()] __device__(i_t var_idx) {
                                return pb.integer_equal(lpopt[var_idx], incumbent[var_idx]);
                              });
   vars_to_fix.resize(end - vars_to_fix.begin(), rins_handle.get_stream());
-  f_t fractional_ratio = (f_t)(vars_to_fix.size()) / (f_t)problem_ptr->n_integer_vars;
+  f_t fractional_ratio = (f_t)(vars_to_fix.size()) / (f_t)problem_copy->n_integer_vars;
 
   // abort if the fractional ratio is too low
   if (fractional_ratio < settings.min_fractional_ratio) {
@@ -164,7 +167,7 @@ void rins_t<i_t, f_t>::run_rins()
   cuopt_assert(thrust::all_of(rins_handle.get_thrust_policy(),
                               vars_to_fix.begin(),
                               vars_to_fix.end(),
-                              [pb = problem_ptr->view()] __device__(i_t var_idx) {
+                              [pb = problem_copy->view()] __device__(i_t var_idx) {
                                 return pb.is_integer_var(var_idx);
                               }),
                "All variables to fix must be integer variables");
@@ -180,15 +183,15 @@ void rins_t<i_t, f_t>::run_rins()
   CUOPT_LOG_DEBUG("Running RINS on solution with objective %g, fixing %d/%d",
                   best_sol.get_user_objective(),
                   vars_to_fix.size(),
-                  problem_ptr->n_integer_vars);
+                  problem_copy->n_integer_vars);
   CUOPT_LOG_DEBUG("RINS fixrate %g time limit %g", fixrate, time_limit);
   CUOPT_LOG_DEBUG("RINS fractional ratio %g%%", fractional_ratio * 100);
 
   f_t prev_obj = best_sol.get_user_objective();
 
   auto [fixed_problem, fixed_assignment, variable_map] = best_sol.fix_variables(vars_to_fix);
   CUOPT_LOG_DEBUG(
-    "new var count %d var_count %d", fixed_problem.n_variables, problem_ptr->n_integer_vars);
+    "new var count %d var_count %d", fixed_problem.n_variables, problem_copy->n_integer_vars);
 
   // should probably just do an spmv to get the objective instead. ugly mess of copies
   solution_t<i_t, f_t> best_sol_fixed_space(fixed_problem);
@@ -271,7 +274,8 @@ void rins_t<i_t, f_t>::run_rins()
     CUOPT_LOG_DEBUG("RINS submip solution found. Objective %.16e. Status %d",
                     branch_and_bound_solution.objective,
                     int(branch_and_bound_status));
-    cuopt_assert(rins_solution_queue.size() > 0, "RINS solution queue is unexpectedly empty");
+    // RINS submip may have just proved the initial guess is the optimal, therefore the queue might
+    // be empty in that case
   }
   if (branch_and_bound_status == dual_simplex::mip_status_t::OPTIMAL) {
     CUOPT_LOG_DEBUG("RINS submip optimal");
@@ -328,7 +332,7 @@ void rins_t<i_t, f_t>::run_rins()
       cuopt_assert(best_sol.test_number_all_integer(), "All must be integers after RINS");
       if (best_sol.get_user_objective() < prev_obj) { improvement_found = true; }
       cuopt_assert(best_sol.assignment.size() == sol_size_before_rins, "Assignment size mismatch");
-      cuopt_assert(best_sol.assignment.size() == problem_ptr->n_variables,
+      cuopt_assert(best_sol.assignment.size() == problem_copy->n_variables,
                    "Assignment size mismatch");
       dm.population.add_external_solution(
         best_sol.get_host_assignment(), best_sol.get_objective(), solution_origin_t::RINS);

cpp/src/mip/diversity/recombiners/bound_prop_recombiner.cuh

@@ -204,8 +204,10 @@ class bound_prop_recombiner_t : public recombiner_t<i_t, f_t> {
       offspring.handle_ptr->sync_stream();
       offspring.unfix_variables(fixed_assignment, variable_map);
       cuopt_func_call(bool feasible_after_unfix = offspring.get_feasible());
-      cuopt_assert(feasible_after_unfix == feasible_after_bounds_prop,
-                   "Feasible after unfix should be same as feasible after bounds prop!");
+      // May be triggered due to numerical issues
+      // TODO: investigate further
+      // cuopt_assert(feasible_after_unfix == feasible_after_bounds_prop,
+      //              "Feasible after unfix should be same as feasible after bounds prop!");
       a.handle_ptr->sync_stream();
     } else {
       timer_t timer(bp_recombiner_config_t::bounds_prop_time_limit);

cpp/src/mip/diversity/recombiners/sub_mip.cuh

@@ -141,6 +141,8 @@ class sub_mip_recombiner_t : public recombiner_t<i_t, f_t> {
       scaling.unscale_solutions(fixed_assignment, dummy);
       // unfix the assignment on given result no matter if it is feasible
       offspring.unfix_variables(fixed_assignment, variable_map);
+      offspring
+        .clamp_within_bounds();  // Scaling might bring some very slight variable bound violations
     } else {
       offspring.round_nearest();
     }
@@ -171,6 +173,7 @@ class sub_mip_recombiner_t : public recombiner_t<i_t, f_t> {
       rmm::device_uvector<f_t> dummy(0, offspring.handle_ptr->get_stream());
       scaling.unscale_solutions(fixed_assignment, dummy);
       sol.unfix_variables(fixed_assignment, variable_map);
+      sol.clamp_within_bounds();  // Scaling might bring some very slight variable bound violations
       sol.compute_feasibility();
       cuopt_func_call(sol.test_variable_bounds());
       population.add_solution(std::move(sol));

cpp/src/mip/feasibility_jump/fj_cpu.cu

@@ -710,10 +710,12 @@ static void perturb(fj_cpu_climber_t<i_t, f_t>& fj_cpu)
   raft::random::PCGenerator rng(fj_cpu.settings.seed + fj_cpu.iterations, 0, 0);
 
   for (auto var_idx : sampled_vars) {
-    f_t lb  = ceil(std::max(get_lower(fj_cpu.h_var_bounds[var_idx]), -1e7));
-    f_t ub  = floor(std::min(get_upper(fj_cpu.h_var_bounds[var_idx]), 1e7));
+    f_t lb  = std::max(get_lower(fj_cpu.h_var_bounds[var_idx]), -1e7);
+    f_t ub  = std::min(get_upper(fj_cpu.h_var_bounds[var_idx]), 1e7);
     f_t val = lb + (ub - lb) * rng.next_double();
     if (fj_cpu.view.pb.is_integer_var(var_idx)) {
+      lb  = std::ceil(lb);
+      ub  = std::floor(ub);
       val = std::round(val);
       val = std::min(std::max(val, lb), ub);
     }
@@ -860,6 +862,13 @@ static void init_fj_cpu(fj_cpu_climber_t<i_t, f_t>& fj_cpu,
 template <typename i_t, typename f_t>
 static void sanity_checks(fj_cpu_climber_t<i_t, f_t>& fj_cpu)
 {
+  // Check that each variable is within its bounds
+  for (i_t var_idx = 0; var_idx < fj_cpu.view.pb.n_variables; ++var_idx) {
+    f_t val = fj_cpu.h_assignment[var_idx];
+    cuopt_assert(fj_cpu.view.pb.check_variable_within_bounds(var_idx, val),
+                 "Variable is out of bounds");
+  }
+
   // Check that each violated constraint is actually violated and not present in
   // satisfied_constraints
   for (const auto& cstr_idx : fj_cpu.violated_constraints) {

cpp/src/mip/local_search/rounding/constraint_prop.cu

@@ -933,6 +933,7 @@ bool constraint_prop_t<i_t, f_t>::find_integer(
     update_host_assignment(sol);
     if (max_timer.check_time_limit()) {
       CUOPT_LOG_DEBUG("Second time limit is reached returning nearest rounding!");
+      collapse_crossing_bounds(*sol.problem_ptr, *orig_sol.problem_ptr, sol.handle_ptr);
       sol.round_nearest();
       timeout_happened = true;
       break;

cpp/src/mip/local_search/rounding/simple_rounding_kernels.cuh

@@ -19,51 +19,56 @@ template <typename i_t, typename f_t>
 __global__ void simple_rounding_kernel(typename solution_t<i_t, f_t>::view_t solution,
                                        bool* successful)
 {
-  if (TH_ID_X >= solution.problem.n_integer_vars) { return; }
-  i_t var_id   = solution.problem.integer_indices[TH_ID_X];
-  f_t curr_val = solution.assignment[var_id];
-  if (solution.problem.is_integer(curr_val)) { return; }
-
-  i_t up_locks   = 0;
-  i_t down_locks = 0;
-
-  auto [offset_begin, offset_end] = solution.problem.reverse_range_for_var(var_id);
-  for (i_t i = offset_begin; i < offset_end; i += 1) {
-    auto cstr_idx   = solution.problem.reverse_constraints[i];
-    auto cstr_coeff = solution.problem.reverse_coefficients[i];
+  for (i_t idx = TH_ID_X; idx < solution.problem.n_integer_vars; idx += gridDim.x * blockDim.x) {
+    i_t var_id   = solution.problem.integer_indices[idx];
+    f_t curr_val = solution.assignment[var_id];
+    if (solution.problem.is_integer(curr_val)) { continue; }
 
-    // boxed constraint. can't be rounded safely
-    if (std::isfinite(solution.problem.constraint_lower_bounds[cstr_idx]) &&
-        std::isfinite(solution.problem.constraint_upper_bounds[cstr_idx])) {
-      up_locks += 1;
-      down_locks += 1;
-      continue;
+    i_t up_locks   = 0;
+    i_t down_locks = 0;
+
+    auto [offset_begin, offset_end] = solution.problem.reverse_range_for_var(var_id);
+    for (i_t i = offset_begin; i < offset_end; i += 1) {
+      auto cstr_idx   = solution.problem.reverse_constraints[i];
+      auto cstr_coeff = solution.problem.reverse_coefficients[i];
+
+      // boxed constraint. can't be rounded safely
+      if (std::isfinite(solution.problem.constraint_lower_bounds[cstr_idx]) &&
+          std::isfinite(solution.problem.constraint_upper_bounds[cstr_idx])) {
+        up_locks += 1;
+        down_locks += 1;
+        continue;
+      }
+
+      f_t sign = std::isfinite(solution.problem.constraint_upper_bounds[cstr_idx]) ? 1 : -1;
+
+      if (cstr_coeff * sign > 0) {
+        up_locks += 1;
+      } else {
+        down_locks += 1;
+      }
     }
 
-    f_t sign = std::isfinite(solution.problem.constraint_upper_bounds[cstr_idx]) ? 1 : -1;
+    auto var_bnd = solution.problem.variable_bounds[var_id];
+    f_t lb       = get_lower(var_bnd);
+    f_t ub       = get_upper(var_bnd);
 
-    if (cstr_coeff * sign > 0) {
-      up_locks += 1;
-    } else {
-      down_locks += 1;
-    }
-  }
+    bool can_round_up   = up_locks == 0 && ceil(curr_val) <= floor(ub);
+    bool can_round_down = down_locks == 0 && floor(curr_val) >= ceil(lb);
 
-  bool can_round_up   = up_locks == 0;
-  bool can_round_down = down_locks == 0;
-
-  if (can_round_up && can_round_down) {
-    if (solution.problem.objective_coefficients[var_id] > 0) {
+    if (can_round_up && can_round_down) {
+      if (solution.problem.objective_coefficients[var_id] > 0) {
+        solution.assignment[var_id] = floor(curr_val);
+      } else {
+        solution.assignment[var_id] = ceil(curr_val);
+      }
+    } else if (can_round_up) {
+      solution.assignment[var_id] = ceil(curr_val);
+    } else if (can_round_down) {
       solution.assignment[var_id] = floor(curr_val);
     } else {
-      solution.assignment[var_id] = ceil(curr_val);
+      *successful = false;
     }
-  } else if (can_round_up) {
-    solution.assignment[var_id] = ceil(curr_val);
-  } else if (can_round_down) {
-    solution.assignment[var_id] = floor(curr_val);
-  } else {
-    *successful = false;
   }
 }