Fixing FJ bugs, Termination status and MIP GAP parameters

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

@@ -224,8 +224,8 @@ class mip_solver_settings_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      = 0.;
-    f_t relative_mip_gap      = 1.0e-3;
+    f_t absolute_mip_gap      = 1.0e-10;
+    f_t relative_mip_gap      = 1.0e-4;
   };
 
   /**

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

@@ -53,6 +53,8 @@ class mip_solution_t : public base_solution_t {
                  f_t max_constraint_violation,
                  f_t max_int_violation,
                  f_t max_variable_bound_violation,
+                 i_t num_nodes,
+                 i_t num_simplex_iterations,
                  std::vector<rmm::device_uvector<f_t>> solution_pool = {});
 
   mip_solution_t(mip_termination_status_t termination_status, rmm::cuda_stream_view stream_view);
@@ -71,6 +73,8 @@ class mip_solution_t : public base_solution_t {
   f_t get_max_constraint_violation() const;
   f_t get_max_int_violation() const;
   f_t get_max_variable_bound_violation() const;
+  i_t get_num_nodes() const;
+  i_t get_num_simplex_iterations() const;
   const std::vector<std::string>& get_variable_names() const;
   const std::vector<rmm::device_uvector<f_t>>& get_solution_pool() const;
   void write_to_sol_file(std::string_view filename, rmm::cuda_stream_view stream_view) const;
@@ -87,6 +91,8 @@ class mip_solution_t : public base_solution_t {
   f_t max_constraint_violation_;
   f_t max_int_violation_;
   f_t max_variable_bound_violation_;
+  i_t num_nodes_;
+  i_t num_simplex_iterations_;
   std::vector<rmm::device_uvector<f_t>> solution_pool_;
 };
 

cpp/include/cuopt/linear_programming/optimization_problem.hpp

@@ -38,14 +38,12 @@ enum class var_t { CONTINUOUS = 0, INTEGER };
 enum class problem_category_t : int8_t { LP = 0, MIP = 1, IP = 2 };
 
 /**
- * @brief Standard form representation of a Linear Program with a twist from the standard form for
- * Quadratic Programs.
+ * @brief A representation of a linear programming (LP) optimization problem
  *
  * @tparam f_t  Data type of the variables and their weights in the equations
  *
- * Standard form representation follows the description from the wiki article here:
- * https://en.wikipedia.org/wiki/Linear_programming#Standard_form. In other words,
- * this structure stores all information used to represent the following LP equation:
+ * This structure stores all the information necessary to represent the following LP:
+ *
  * <pre>
  * Minimize:
  *   dot(c, x)
@@ -63,7 +61,7 @@ enum class problem_category_t : int8_t { LP = 0, MIP = 1, IP = 2 };
  *
  * Objective value can be scaled and offset accordingly:
  * objective_scaling_factor * (dot(c, x) + objective_offset)
- * please refeto to the `set_objective_scaling_factor()` and `set_objective_offset()` method.
+ * please refer to the `set_objective_scaling_factor()` and `set_objective_offset()` methods.
  */
 template <typename i_t, typename f_t>
 class optimization_problem_t {

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

@@ -214,7 +214,7 @@ class pdlp_solver_settings_t {
    *
    * @note By default there is no time limit.
    * For performance reasons, cuOpt's does not constantly checks for time limit, thus, the solver
-   * might run a few milliseconds over the limit.
+   * might run slightly over the limit.
    * If set along iteration limit, the first limit reached will exit.
    *
    * @param time_limit Time limit to set in seconds

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

@@ -93,7 +93,7 @@ class optimization_problem_solution_t : public base_solution_t {
     /** Objective value for the extreme dual ray */
     f_t dual_ray_linear_objective;
 
-    /** Solve time in milliseconds */
+    /** Solve time in seconds */
     double solve_time;
   };
 
@@ -163,7 +163,7 @@ class optimization_problem_solution_t : public base_solution_t {
                                   bool deep_copy);
 
   /**
-   * @brief Set the solve time in milliseconds
+   * @brief Set the solve time in seconds
    *
    * @param ms Time in ms
    */
@@ -177,9 +177,9 @@ class optimization_problem_solution_t : public base_solution_t {
   void set_termination_status(pdlp_termination_status_t termination_status);
 
   /**
-   * @brief Get the solve time in milliseconds
+   * @brief Get the solve time in seconds
    *
-   * @return Time in ms
+   * @return Time in seconds
    */
   double get_solve_time() const;
 

cpp/include/cuopt/linear_programming/solve.hpp

@@ -38,7 +38,7 @@ namespace cuopt::linear_programming {
  * @tparam f_t Data type of the variables and their weights in the equations
  *
  * @param[in] op_problem  An optimization_problem_t<i_t, f_t> object with a
- * representation of a linear program on standard form.
+ * representation of a linear program
  * @param[in] settings  A pdlp_solver_settings_t<i_t, f_t> object with the settings for the PDLP
  * solver.
  * @param[in] problem_checking  If true, the problem is checked for consistency.
@@ -63,7 +63,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
  *
  * @param[in] handle_ptr  A raft::handle_t object with its corresponding CUDA stream.
  * @param[in] mps_data_model  An optimization_problem_t<i_t, f_t> object with a
- * representation of a linear program on standard form.
+ * representation of a linear program
  * @param[in] settings  A pdlp_solver_settings_t<i_t, f_t> object with the settings for the PDLP
  * solver.
  * @param[in] problem_checking  If true, the problem is checked for consistency.
@@ -86,7 +86,7 @@ optimization_problem_solution_t<i_t, f_t> solve_lp(
  * @tparam f_t Data type of the variables and their weights in the equations
  *
  * @param[in] op_problem  An optimization_problem_t<i_t, f_t> object with a
- * representation of a linear program on standard form.
+ * representation of a linear program
  * @return optimization_problem_solution_t<i_t, f_t> owning container for the solver solution
  */
 template <typename i_t, typename f_t>
@@ -101,7 +101,7 @@ mip_solution_t<i_t, f_t> solve_mip(
  * @tparam f_t Data type of the variables and their weights in the equations
  *
  * @param[in] mps_data_model  An optimization_problem_t<i_t, f_t> object with a
- * representation of a linear program on standard form.
+ * representation of a linear program
  * @return optimization_problem_solution_t<i_t, f_t> owning container for the solver solution
  */
 template <typename i_t, typename f_t>

cpp/include/cuopt/linear_programming/utilities/cython_solve.hpp

@@ -84,6 +84,8 @@ struct mip_ret_t {
   double max_constraint_violation_;
   double max_int_violation_;
   double max_variable_bound_violation_;
+  int nodes_;
+  int simplex_iterations_;
 };
 
 struct solver_ret_t {

cpp/libmps_parser/include/mps_parser/data_model_view.hpp

@@ -27,14 +27,11 @@
 namespace cuopt::mps_parser {
 
 /**
- * @brief Standard form representation of a Linear Program with a twist from the standard form for
- * Quadratic Programs.
+ * @brief A representation of a linear programming (LP) optimization problem
  *
  * @tparam f_t  Data type of the variables and their weights in the equations
  *
- * Standard form representation follows the description from the wiki article here:
- * https://en.wikipedia.org/wiki/Linear_programming#Standard_form. In other words,
- * this structure stores all information used to represent the following LP equation:
+ * A linear programming optimization problem is defined as follows:
  * <pre>
  * Minimize:
  *   dot(c, x)

cpp/src/dual_simplex/branch_and_bound.cpp

@@ -222,13 +222,21 @@ f_t sgn(f_t x)
   return x < 0 ? -1 : 1;
 }
 
+template <typename f_t>
+f_t relative_gap(f_t obj_value, f_t lower_bound)
+{
+  f_t user_mip_gap = obj_value == 0.0
+                       ? (lower_bound == 0.0 ? 0.0 : std::numeric_limits<f_t>::infinity())
+                       : std::abs(obj_value - lower_bound) / std::abs(obj_value);
+  if (user_mip_gap != user_mip_gap) { return std::numeric_limits<f_t>::infinity(); }
+  return user_mip_gap;
+}
+
 template <typename f_t>
 std::string user_mip_gap(f_t obj_value, f_t lower_bound)
 {
-  const f_t user_mip_gap = obj_value == 0.0
-                             ? (lower_bound == 0.0 ? 0.0 : std::numeric_limits<f_t>::infinity())
-                             : std::abs(obj_value - lower_bound) / std::abs(obj_value);
-  if (user_mip_gap != user_mip_gap) {
+  const f_t user_mip_gap = relative_gap(obj_value, lower_bound);
+  if (user_mip_gap == std::numeric_limits<f_t>::infinity()) {
     return "  -  ";
   } else {
     constexpr int BUFFER_LEN = 32;
@@ -545,7 +553,9 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
 
   f_t total_lp_iters = 0.0;
   f_t last_log       = 0;
-  while (gap > settings.mip_gap_tol && heap.size() > 0) {
+  while (gap > settings.absolute_mip_gap_tol &&
+         relative_gap(get_upper_bound<f_t>(), lower_bound) > settings.relative_mip_gap_tol &&
+         heap.size() > 0) {
     // Check if there are any solutions to repair
     std::vector<std::vector<f_t>> to_repair;
     global_variables::mutex_repair.lock();
@@ -607,7 +617,8 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     const i_t leaf_depth = node_ptr->depth;
     f_t now              = toc(start_time);
     f_t time_since_log   = last_log == 0 ? 1.0 : toc(last_log);
-    if ((nodes_explored % 1000 == 0 || gap < 10 * settings.mip_gap_tol || nodes_explored < 1000) &&
+    if ((nodes_explored % 1000 == 0 || gap < 10 * settings.absolute_mip_gap_tol ||
+         nodes_explored < 1000) &&
           (time_since_log >= 1) ||
         (time_since_log > 60) || now > settings.time_limit) {
       settings.log.printf(" %8d %8lu       %+13.6e  %+10.6e   %4d   %7.1e     %s %9.2f\n",
@@ -789,11 +800,16 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     compute_user_objective(original_lp, get_upper_bound<f_t>()),
     compute_user_objective(original_lp, lower_bound));
 
-  if (gap < settings.mip_gap_tol) {
+  if (gap <= settings.absolute_mip_gap_tol ||
+      relative_gap(get_upper_bound<f_t>(), lower_bound) <= settings.relative_mip_gap_tol) {
     status = mip_status_t::OPTIMAL;
-    if (gap > 0) {
-      settings.log.printf("Optimal solution found within MIP gap tolerance (%.1e)\n",
-                          settings.mip_gap_tol);
+    if (gap > 0 && gap <= settings.absolute_mip_gap_tol) {
+      settings.log.printf("Optimal solution found within absolute MIP gap tolerance (%.1e)\n",
+                          settings.absolute_mip_gap_tol);
+    } else if (gap > 0 &&
+               relative_gap(get_upper_bound<f_t>(), lower_bound) <= settings.relative_mip_gap_tol) {
+      settings.log.printf("Optimal solution found within relative MIP gap tolerance (%.1e)\n",
+                          settings.relative_mip_gap_tol);
     } else {
       settings.log.printf("Optimal solution found.\n");
     }
@@ -803,9 +819,10 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   uncrush_primal_solution(original_problem, original_lp, incumbent.x, solution.x);
-  solution.objective   = incumbent.objective;
-  solution.lower_bound = lower_bound;
-
+  solution.objective          = incumbent.objective;
+  solution.lower_bound        = lower_bound;
+  solution.nodes_explored     = nodes_explored;
+  solution.simplex_iterations = total_lp_iters;
   return status;
 }
 

cpp/src/dual_simplex/phase2.cpp

@@ -325,6 +325,7 @@ i_t phase2_ratio_test(const lp_problem_t<i_t, f_t>& lp,
 template <typename i_t, typename f_t>
 i_t bound_flipping_ratio_test(const lp_problem_t<i_t, f_t>& lp,
                               const simplex_solver_settings_t<i_t, f_t>& settings,
+                              f_t start_time,
                               const std::vector<variable_status_t>& vstatus,
                               const std::vector<i_t>& nonbasic_list,
                               const std::vector<f_t>& x,
@@ -421,6 +422,12 @@ i_t bound_flipping_ratio_test(const lp_problem_t<i_t, f_t>& lp,
       // we hit a variable that is not bounded. Exit
       break;
     }
+
+    if (toc(start_time) > settings.time_limit) { return -2; }
+    if (settings.concurrent_halt != nullptr &&
+        settings.concurrent_halt->load(std::memory_order_acquire) == 1) {
+      return -3;
+    }
   }
   // step_length, nonbasic_entering, and entering_index are defined after the
   // while loop
@@ -1327,6 +1334,7 @@ dual::status_t dual_phase2(i_t phase,
     } else if (bound_flip_ratio) {
       entering_index = phase2::bound_flipping_ratio_test(lp,
                                                          settings,
+                                                         start_time,
                                                          vstatus,
                                                          nonbasic_list,
                                                          x,
@@ -1340,6 +1348,8 @@ dual::status_t dual_phase2(i_t phase,
       entering_index = phase2::phase2_ratio_test(
         lp, settings, vstatus, nonbasic_list, z, delta_z, step_length, nonbasic_entering_index);
     }
+    if (entering_index == -2) { return dual::status_t::TIME_LIMIT; }
+    if (entering_index == -3) { return dual::status_t::CONCURRENT_LIMIT; }
     if (entering_index == -1) {
       if (primal_infeasibility > settings.primal_tol &&
           max_val < settings.steepest_edge_primal_tol) {

cpp/src/dual_simplex/simplex_solver_settings.hpp

@@ -34,7 +34,8 @@ struct simplex_solver_settings_t {
   simplex_solver_settings_t()
     : iteration_limit(std::numeric_limits<i_t>::max()),
       time_limit(std::numeric_limits<f_t>::infinity()),
-      mip_gap_tol(1e-3),
+      absolute_mip_gap_tol(0.0),
+      relative_mip_gap_tol(1e-3),
       integer_tol(1e-5),
       primal_tol(1e-6),
       dual_tol(1e-6),
@@ -73,15 +74,16 @@ struct simplex_solver_settings_t {
   void close_log_file() { log.close_log_file(); }
   i_t iteration_limit;
   f_t time_limit;
-  f_t mip_gap_tol;  // Tolerance on mip gap to declare optimal
-  f_t integer_tol;  // Tolerance on integralitiy violation
-  f_t primal_tol;   // Absolute primal infeasibility tolerance
-  f_t dual_tol;     // Absolute dual infeasibility tolerance
-  f_t pivot_tol;    // Simplex pivot tolerance
-  f_t tight_tol;    // A tight tolerance used to check for infeasibility
-  f_t fixed_tol;    // If l <= x <= u with u - l < fixed_tol a variable is consider fixed
-  f_t zero_tol;     // Values below this tolerance are considered numerically zero
-  f_t cut_off;      // If the dual objective is greater than the cutoff we stop
+  f_t absolute_mip_gap_tol;  // Tolerance on mip gap to declare optimal
+  f_t relative_mip_gap_tol;  // Tolerance on mip gap to declare optimal
+  f_t integer_tol;           // Tolerance on integralitiy violation
+  f_t primal_tol;            // Absolute primal infeasibility tolerance
+  f_t dual_tol;              // Absolute dual infeasibility tolerance
+  f_t pivot_tol;             // Simplex pivot tolerance
+  f_t tight_tol;             // A tight tolerance used to check for infeasibility
+  f_t fixed_tol;             // If l <= x <= u with u - l < fixed_tol a variable is consider fixed
+  f_t zero_tol;              // Values below this tolerance are considered numerically zero
+  f_t cut_off;               // If the dual objective is greater than the cutoff we stop
   f_t
     steepest_edge_ratio;  // the ratio of computed steepest edge mismatch from updated steepest edge
   f_t steepest_edge_primal_tol;    // Primal tolerance divided by steepest edge norm

cpp/src/dual_simplex/solution.hpp

@@ -68,6 +68,8 @@ class mip_solution_t {
   std::vector<f_t> x;
   f_t objective;
   f_t lower_bound;
+  i_t nodes_explored;
+  i_t simplex_iterations;
 };
 
 }  // namespace cuopt::linear_programming::dual_simplex

cpp/src/linear_programming/pdlp.cuh

@@ -61,7 +61,7 @@ class pdlp_solver_t {
    * For full description of algorithm, see https://arxiv.org/abs/2106.04756
    *
    * @param[in] op_problem An problem_t<i_t, f_t> object with a
-   * representation of a linear program on standard form.
+   * representation of a linear program
    */
   pdlp_solver_t(
     problem_t<i_t, f_t>& op_problem,