Implement UnboundedOrInfeasible termination status

cpp/include/cuopt/linear_programming/constants.h

@@ -82,17 +82,18 @@
 #define CUOPT_MODE_DETERMINISTIC 1
 
 /* @brief LP/MIP termination status constants */
-#define CUOPT_TERIMINATION_STATUS_NO_TERMINATION   0
-#define CUOPT_TERIMINATION_STATUS_OPTIMAL          1
-#define CUOPT_TERIMINATION_STATUS_INFEASIBLE       2
-#define CUOPT_TERIMINATION_STATUS_UNBOUNDED        3
-#define CUOPT_TERIMINATION_STATUS_ITERATION_LIMIT  4
-#define CUOPT_TERIMINATION_STATUS_TIME_LIMIT       5
-#define CUOPT_TERIMINATION_STATUS_NUMERICAL_ERROR  6
-#define CUOPT_TERIMINATION_STATUS_PRIMAL_FEASIBLE  7
-#define CUOPT_TERIMINATION_STATUS_FEASIBLE_FOUND   8
-#define CUOPT_TERIMINATION_STATUS_CONCURRENT_LIMIT 9
-#define CUOPT_TERIMINATION_STATUS_WORK_LIMIT       10
+#define CUOPT_TERMINATION_STATUS_NO_TERMINATION          0
+#define CUOPT_TERMINATION_STATUS_OPTIMAL                 1
+#define CUOPT_TERMINATION_STATUS_INFEASIBLE              2
+#define CUOPT_TERMINATION_STATUS_UNBOUNDED               3
+#define CUOPT_TERMINATION_STATUS_ITERATION_LIMIT         4
+#define CUOPT_TERMINATION_STATUS_TIME_LIMIT              5
+#define CUOPT_TERMINATION_STATUS_NUMERICAL_ERROR         6
+#define CUOPT_TERMINATION_STATUS_PRIMAL_FEASIBLE         7
+#define CUOPT_TERMINATION_STATUS_FEASIBLE_FOUND          8
+#define CUOPT_TERMINATION_STATUS_CONCURRENT_LIMIT        9
+#define CUOPT_TERMINATION_STATUS_WORK_LIMIT              10
+#define CUOPT_TERMINATION_STATUS_UNBOUNDED_OR_INFEASIBLE 11
 
 /* @brief The objective sense constants */
 #define CUOPT_MINIMIZE 1

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

@@ -24,13 +24,14 @@
 namespace cuopt::linear_programming {
 
 enum class mip_termination_status_t : int8_t {
-  NoTermination = CUOPT_TERIMINATION_STATUS_NO_TERMINATION,
-  Optimal       = CUOPT_TERIMINATION_STATUS_OPTIMAL,
-  FeasibleFound = CUOPT_TERIMINATION_STATUS_FEASIBLE_FOUND,
-  Infeasible    = CUOPT_TERIMINATION_STATUS_INFEASIBLE,
-  Unbounded     = CUOPT_TERIMINATION_STATUS_UNBOUNDED,
-  TimeLimit     = CUOPT_TERIMINATION_STATUS_TIME_LIMIT,
-  WorkLimit     = CUOPT_TERIMINATION_STATUS_WORK_LIMIT,
+  NoTermination         = CUOPT_TERMINATION_STATUS_NO_TERMINATION,
+  Optimal               = CUOPT_TERMINATION_STATUS_OPTIMAL,
+  FeasibleFound         = CUOPT_TERMINATION_STATUS_FEASIBLE_FOUND,
+  Infeasible            = CUOPT_TERMINATION_STATUS_INFEASIBLE,
+  Unbounded             = CUOPT_TERMINATION_STATUS_UNBOUNDED,
+  TimeLimit             = CUOPT_TERMINATION_STATUS_TIME_LIMIT,
+  WorkLimit             = CUOPT_TERMINATION_STATUS_WORK_LIMIT,
+  UnboundedOrInfeasible = CUOPT_TERMINATION_STATUS_UNBOUNDED_OR_INFEASIBLE,
 };
 
 template <typename i_t, typename f_t>

cpp/include/cuopt/linear_programming/pdlp/solver_solution.hpp

@@ -25,15 +25,16 @@ namespace cuopt::linear_programming {
 
 // Possible reasons for terminating
 enum class pdlp_termination_status_t : int8_t {
-  NoTermination    = CUOPT_TERIMINATION_STATUS_NO_TERMINATION,
-  NumericalError   = CUOPT_TERIMINATION_STATUS_NUMERICAL_ERROR,
-  Optimal          = CUOPT_TERIMINATION_STATUS_OPTIMAL,
-  PrimalInfeasible = CUOPT_TERIMINATION_STATUS_INFEASIBLE,
-  DualInfeasible   = CUOPT_TERIMINATION_STATUS_UNBOUNDED,
-  IterationLimit   = CUOPT_TERIMINATION_STATUS_ITERATION_LIMIT,
-  TimeLimit        = CUOPT_TERIMINATION_STATUS_TIME_LIMIT,
-  PrimalFeasible   = CUOPT_TERIMINATION_STATUS_PRIMAL_FEASIBLE,
-  ConcurrentLimit  = CUOPT_TERIMINATION_STATUS_CONCURRENT_LIMIT
+  NoTermination         = CUOPT_TERMINATION_STATUS_NO_TERMINATION,
+  NumericalError        = CUOPT_TERMINATION_STATUS_NUMERICAL_ERROR,
+  Optimal               = CUOPT_TERMINATION_STATUS_OPTIMAL,
+  PrimalInfeasible      = CUOPT_TERMINATION_STATUS_INFEASIBLE,
+  DualInfeasible        = CUOPT_TERMINATION_STATUS_UNBOUNDED,
+  IterationLimit        = CUOPT_TERMINATION_STATUS_ITERATION_LIMIT,
+  TimeLimit             = CUOPT_TERMINATION_STATUS_TIME_LIMIT,
+  PrimalFeasible        = CUOPT_TERMINATION_STATUS_PRIMAL_FEASIBLE,
+  ConcurrentLimit       = CUOPT_TERMINATION_STATUS_CONCURRENT_LIMIT,
+  UnboundedOrInfeasible = CUOPT_TERMINATION_STATUS_UNBOUNDED_OR_INFEASIBLE
 };
 
 /**

cpp/src/dual_simplex/solve.cpp

@@ -212,11 +212,13 @@ lp_status_t solve_linear_program_with_advanced_basis(
                                 edge_norms,
                                 work_unit_context);
   }
-  if (phase1_status == dual::status_t::NUMERICAL ||
-      phase1_status == dual::status_t::DUAL_UNBOUNDED) {
+  if (phase1_status == dual::status_t::NUMERICAL) {
     settings.log.printf("Failed in Phase 1\n");
     return lp_status_t::NUMERICAL_ISSUES;
   }
+  if (phase1_status == dual::status_t::DUAL_UNBOUNDED) {
+    return lp_status_t::UNBOUNDED_OR_INFEASIBLE;
+  }
   if (phase1_status == dual::status_t::TIME_LIMIT) { return lp_status_t::TIME_LIMIT; }
   if (phase1_status == dual::status_t::WORK_LIMIT) { return lp_status_t::WORK_LIMIT; }
   if (phase1_status == dual::status_t::ITERATION_LIMIT) { return lp_status_t::ITERATION_LIMIT; }
@@ -312,7 +314,7 @@ lp_status_t solve_linear_program_with_advanced_basis(
     if (status == dual::status_t::CUTOFF) { lp_status = lp_status_t::CUTOFF; }
     original_solution.iterations = iter;
   } else {
-    // Dual infeasible -> Primal unbounded
+    // Dual infeasible -> Primal unbounded or infeasible
     settings.log.printf("Dual infeasible\n");
     original_solution.objective = -inf;
     if (lp.obj_scale == 1.0) {
@@ -323,7 +325,7 @@ lp_status_t solve_linear_program_with_advanced_basis(
       original_solution.user_objective = inf;
     }
     original_solution.iterations = iter;
-    return lp_status_t::UNBOUNDED;
+    return lp_status_t::UNBOUNDED_OR_INFEASIBLE;
   }
   return lp_status;
 }

cpp/src/dual_simplex/solve.hpp

@@ -23,16 +23,17 @@ template <typename i_t, typename f_t>
 bool is_mip(const user_problem_t<i_t, f_t>& problem);
 
 enum class lp_status_t {
-  OPTIMAL          = 0,
-  INFEASIBLE       = 1,
-  UNBOUNDED        = 2,
-  ITERATION_LIMIT  = 3,
-  TIME_LIMIT       = 4,
-  NUMERICAL_ISSUES = 5,
-  CUTOFF           = 6,
-  CONCURRENT_LIMIT = 7,
-  WORK_LIMIT       = 8,
-  UNSET            = 9
+  OPTIMAL                 = 0,
+  INFEASIBLE              = 1,
+  UNBOUNDED               = 2,
+  ITERATION_LIMIT         = 3,
+  TIME_LIMIT              = 4,
+  NUMERICAL_ISSUES        = 5,
+  CUTOFF                  = 6,
+  CONCURRENT_LIMIT        = 7,
+  WORK_LIMIT              = 8,
+  UNBOUNDED_OR_INFEASIBLE = 9,
+  UNSET                   = 10
 };
 
 static std::string lp_status_to_string(lp_status_t status)
@@ -47,6 +48,7 @@ static std::string lp_status_to_string(lp_status_t status)
     case lp_status_t::CUTOFF: return "CUTOFF";
     case lp_status_t::CONCURRENT_LIMIT: return "CONCURRENT_LIMIT";
     case lp_status_t::WORK_LIMIT: return "WORK_LIMIT";
+    case lp_status_t::UNBOUNDED_OR_INFEASIBLE: return "UNBOUNDED_OR_INFEASIBLE";
     case lp_status_t::UNSET: return "UNSET";
   }
   return "UNKNOWN";

cpp/src/mip_heuristics/presolve/third_party_presolve.cpp

@@ -477,6 +477,36 @@ void check_presolve_status(const papilo::PresolveStatus& status)
   }
 }
 
+third_party_presolve_status_t convert_papilo_presolve_status_to_third_party_presolve_status(
+  const papilo::PresolveStatus& status)
+{
+  switch (status) {
+    case papilo::PresolveStatus::kUnchanged: return third_party_presolve_status_t::UNCHANGED;
+    case papilo::PresolveStatus::kReduced: return third_party_presolve_status_t::REDUCED;
+    case papilo::PresolveStatus::kUnbndOrInfeas:
+      return third_party_presolve_status_t::UNBNDORINFEAS;
+    case papilo::PresolveStatus::kInfeasible: return third_party_presolve_status_t::INFEASIBLE;
+    case papilo::PresolveStatus::kUnbounded:
+      return third_party_presolve_status_t::UNBOUNDED;
+      // Do not implement default case to trigger compile time error if new enum is added
+  }
+  return third_party_presolve_status_t::UNCHANGED;
+}
+
+third_party_presolve_status_t convert_pslp_presolve_status_to_third_party_presolve_status(
+  const PresolveStatus& status)
+{
+  switch (status) {
+    case PresolveStatus_::UNCHANGED: return third_party_presolve_status_t::UNCHANGED;
+    case PresolveStatus_::REDUCED: return third_party_presolve_status_t::REDUCED;
+    case PresolveStatus_::INFEASIBLE: return third_party_presolve_status_t::INFEASIBLE;
+    case PresolveStatus_::UNBNDORINFEAS:
+      return third_party_presolve_status_t::UNBNDORINFEAS;
+      // Do not implement default case to trigger compile time error if new enum is added
+  }
+  return third_party_presolve_status_t::UNCHANGED;
+}
+
 void check_postsolve_status(const papilo::PostsolveStatus& status)
 {
   switch (status) {
@@ -570,7 +600,7 @@ void set_presolve_parameters(papilo::Presolve<f_t>& presolver,
 }
 
 template <typename i_t, typename f_t>
-std::optional<third_party_presolve_result_t<i_t, f_t>> third_party_presolve_t<i_t, f_t>::apply_pslp(
+third_party_presolve_result_t<i_t, f_t> third_party_presolve_t<i_t, f_t>::apply_pslp(
   optimization_problem_t<i_t, f_t> const& op_problem, const double time_limit)
 {
   if constexpr (std::is_same_v<f_t, double>) {
@@ -584,23 +614,29 @@ std::optional<third_party_presolve_result_t<i_t, f_t>> third_party_presolve_t<i_
     pslp_presolver_ = ctx.presolver;
     pslp_stgs_      = ctx.settings;
 
+    auto status = convert_pslp_presolve_status_to_third_party_presolve_status(ctx.status);
     if (ctx.status == PresolveStatus_::INFEASIBLE || ctx.status == PresolveStatus_::UNBNDORINFEAS) {
-      return std::nullopt;
+      optimization_problem_t<i_t, f_t> empty_problem(op_problem.get_handle_ptr());
+      return third_party_presolve_result_t<i_t, f_t>{status, std::move(empty_problem), {}, {}, {}};
     }
 
     auto opt_problem = build_optimization_problem_from_pslp<i_t, f_t>(
       pslp_presolver_, op_problem.get_handle_ptr(), maximize_, original_obj_offset);
-    opt_problem.set_problem_name(op_problem.get_problem_name());
-    return std::make_optional(third_party_presolve_result_t<i_t, f_t>{opt_problem, {}});
+    return third_party_presolve_result_t<i_t, f_t>{status, std::move(opt_problem), {}, {}, {}};
   } else {
     cuopt_expects(
       false, error_type_t::ValidationError, "PSLP presolver only supports double precision");
-    return std::nullopt;
+    return third_party_presolve_result_t<i_t, f_t>{
+      third_party_presolve_status_t::UNCHANGED,
+      optimization_problem_t<i_t, f_t>(op_problem.get_handle_ptr()),
+      {},
+      {},
+      {}};
   }
 }
 
 template <typename i_t, typename f_t>
-std::optional<third_party_presolve_result_t<i_t, f_t>> third_party_presolve_t<i_t, f_t>::apply(
+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,
   cuopt::linear_programming::presolver_t presolver,
@@ -648,9 +684,12 @@ std::optional<third_party_presolve_result_t<i_t, f_t>> third_party_presolve_t<i_
 
   auto result = papilo_presolver.apply(papilo_problem);
   check_presolve_status(result.status);
+  auto status = convert_papilo_presolve_status_to_third_party_presolve_status(result.status);
   if (result.status == papilo::PresolveStatus::kInfeasible ||
-      result.status == papilo::PresolveStatus::kUnbndOrInfeas) {
-    return std::nullopt;
+      result.status == papilo::PresolveStatus::kUnbndOrInfeas ||
+      result.status == papilo::PresolveStatus::kUnbounded) {
+    optimization_problem_t<i_t, f_t> empty_problem(op_problem.get_handle_ptr());
+    return third_party_presolve_result_t<i_t, f_t>{status, std::move(empty_problem), {}, {}, {}};
   }
   papilo_post_solve_storage_.reset(new papilo::PostsolveStorage<f_t>(result.postsolve));
   CUOPT_LOG_INFO("Presolve removed: %d constraints, %d variables, %d nonzeros",
@@ -686,8 +725,11 @@ std::optional<third_party_presolve_result_t<i_t, f_t>> third_party_presolve_t<i_
     }
   }
 
-  return std::make_optional(third_party_presolve_result_t<i_t, f_t>{
-    opt_problem, implied_integer_indices, reduced_to_original_map_, original_to_reduced_map_});
+  return third_party_presolve_result_t<i_t, f_t>{status,
+                                                 std::move(opt_problem),
+                                                 implied_integer_indices,
+                                                 reduced_to_original_map_,
+                                                 original_to_reduced_map_};
 }
 
 template <typename i_t, typename f_t>

cpp/src/mip_heuristics/presolve/third_party_presolve.hpp

@@ -27,8 +27,18 @@ struct papilo_postsolve_deleter {
   void operator()(papilo::PostsolveStorage<f_t>* ptr) const;
 };
 
+enum class third_party_presolve_status_t {
+  INFEASIBLE,
+  UNBNDORINFEAS,
+  UNBOUNDED,
+  OPTIMAL,
+  REDUCED,
+  UNCHANGED,
+};
+
 template <typename i_t, typename f_t>
 struct third_party_presolve_result_t {
+  third_party_presolve_status_t status;
   optimization_problem_t<i_t, f_t> reduced_problem;
   std::vector<i_t> implied_integer_indices;
   std::vector<i_t> reduced_to_original_map;
@@ -48,15 +58,14 @@ class third_party_presolve_t {
   third_party_presolve_t(third_party_presolve_t&&)                 = delete;
   third_party_presolve_t& operator=(third_party_presolve_t&&)      = delete;
 
-  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,
-    cuopt::linear_programming::presolver_t presolver,
-    bool dual_postsolve,
-    f_t absolute_tolerance,
-    f_t relative_tolerance,
-    double time_limit,
-    i_t num_cpu_threads = 0);
+  third_party_presolve_result_t<i_t, f_t> apply(optimization_problem_t<i_t, f_t> const& op_problem,
+                                                problem_category_t category,
+                                                cuopt::linear_programming::presolver_t presolver,
+                                                bool dual_postsolve,
+                                                f_t absolute_tolerance,
+                                                f_t relative_tolerance,
+                                                double time_limit,
+                                                i_t num_cpu_threads = 0);
 
   void undo(rmm::device_uvector<f_t>& primal_solution,
             rmm::device_uvector<f_t>& dual_solution,
@@ -74,7 +83,7 @@ class third_party_presolve_t {
   ~third_party_presolve_t();
 
  private:
-  std::optional<third_party_presolve_result_t<i_t, f_t>> apply_pslp(
+  third_party_presolve_result_t<i_t, f_t> apply_pslp(
     optimization_problem_t<i_t, f_t> const& op_problem, const double time_limit);
 
   void undo_pslp(rmm::device_uvector<f_t>& primal_solution,

cpp/src/mip_heuristics/solve.cu

@@ -244,7 +244,7 @@ mip_solution_t<i_t, f_t> solve_mip(optimization_problem_t<i_t, f_t>& op_problem,
 
     double presolve_time = 0.0;
     std::unique_ptr<detail::third_party_presolve_t<i_t, f_t>> presolver;
-    std::optional<detail::third_party_presolve_result_t<i_t, f_t>> presolve_result;
+    std::optional<detail::third_party_presolve_result_t<i_t, f_t>> presolve_result_opt;
     detail::problem_t<i_t, f_t> problem(
       op_problem, settings.get_tolerances(), settings.determinism_mode == CUOPT_MODE_DETERMINISTIC);
 
@@ -283,22 +283,32 @@ mip_solution_t<i_t, f_t> solve_mip(optimization_problem_t<i_t, f_t>& op_problem,
                                      settings.tolerances.relative_tolerance,
                                      presolve_time_limit,
                                      settings.num_cpu_threads);
-      if (!result.has_value()) {
+      if (result.status == detail::third_party_presolve_status_t::INFEASIBLE) {
         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());
       }
-      presolve_result.emplace(std::move(*result));
+      if (result.status == detail::third_party_presolve_status_t::UNBNDORINFEAS) {
+        return mip_solution_t<i_t, f_t>(mip_termination_status_t::UnboundedOrInfeasible,
+                                        solver_stats_t<i_t, f_t>{},
+                                        op_problem.get_handle_ptr()->get_stream());
+      }
+      if (result.status == detail::third_party_presolve_status_t::UNBOUNDED) {
+        return mip_solution_t<i_t, f_t>(mip_termination_status_t::Unbounded,
+                                        solver_stats_t<i_t, f_t>{},
+                                        op_problem.get_handle_ptr()->get_stream());
+      }
+      presolve_result_opt.emplace(std::move(result));
 
-      problem = detail::problem_t<i_t, f_t>(presolve_result->reduced_problem);
+      problem = detail::problem_t<i_t, f_t>(presolve_result_opt->reduced_problem);
       problem.set_papilo_presolve_data(presolver.get(),
-                                       presolve_result->reduced_to_original_map,
-                                       presolve_result->original_to_reduced_map,
+                                       presolve_result_opt->reduced_to_original_map,
+                                       presolve_result_opt->original_to_reduced_map,
                                        op_problem.get_n_variables());
-      problem.set_implied_integers(presolve_result->implied_integer_indices);
+      problem.set_implied_integers(presolve_result_opt->implied_integer_indices);
       presolve_time = timer.elapsed_time();
-      if (presolve_result->implied_integer_indices.size() > 0) {
-        CUOPT_LOG_INFO("%d implied integers", presolve_result->implied_integer_indices.size());
+      if (presolve_result_opt->implied_integer_indices.size() > 0) {
+        CUOPT_LOG_INFO("%d implied integers", presolve_result_opt->implied_integer_indices.size());
       }
       CUOPT_LOG_INFO("Papilo presolve time: %.2f", presolve_time);
     }

cpp/src/mip_heuristics/solver_solution.cu

@@ -137,8 +137,9 @@ std::string mip_solution_t<i_t, f_t>::get_termination_status_string(
     case mip_termination_status_t::Infeasible: return "Infeasible";
     case mip_termination_status_t::TimeLimit: return "TimeLimit";
     case mip_termination_status_t::WorkLimit: return "WorkLimit";
-    case mip_termination_status_t::Unbounded:
-      return "Unbounded";
+    case mip_termination_status_t::Unbounded: return "Unbounded";
+    case mip_termination_status_t::UnboundedOrInfeasible:
+      return "UnboundedOrInfeasible";
       // Do not implement default case to trigger compile time error if new enum is added
   }
   return std::string();