Fix several bugs appeared in unit testing of JuMP interface (#149)

While testing the JuMP interface, there were several corner cases showed up. Especially for tiny problems and problems without constraints. This PR has several such fixes.  

In addition, following things are added: 

- Implement `cuOptGetVersion()` in C-API
- Change the max allowed absolute gap tolerance to INFINITY 

## Issue

closes #148 
closes #147 
closes #166

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

Approvers:
  - Nicolas Blin (https://github.com/Kh4ster)
  - Ramakrishnap (https://github.com/rgsl888prabhu)
  - Chris Maes (https://github.com/chris-maes)
  - Akif ÇÖRDÜK (https://github.com/akifcorduk)

URL: #149

Author: rg20

cpp/include/cuopt/linear_programming/cuopt_c.h

@@ -93,6 +93,22 @@ int8_t cuOptGetFloatSize();
  */
 int8_t cuOptGetIntSize();
 
+/**
+ * @brief Get the version of the library.
+ *
+ * @param[out] version_major - A pointer to a cuopt_int_t that will contain the major version
+ * number.
+ * @param[out] version_minor - A pointer to a cuopt_int_t that will contain the minor version
+ * number.
+ * @param[out] version_patch - A pointer to a cuopt_int_t that will contain the patch version
+ * number.
+ *
+ * @return A status code indicating success or failure.
+ */
+cuopt_int_t cuOptGetVersion(cuopt_int_t* version_major,
+                            cuopt_int_t* version_minor,
+                            cuopt_int_t* version_patch);
+
 /**
  * @brief Read an optimization problem from an MPS file.
  *
@@ -656,6 +672,18 @@ cuopt_int_t cuOptGetSolutionBound(cuOptSolution solution, cuopt_float_t* solutio
  */
 cuopt_int_t cuOptGetDualSolution(cuOptSolution solution, cuopt_float_t* dual_solution_ptr);
 
+/** @brief Get the dual objective value of an optimization problem.
+ *
+ * @param[in] solution - The solution object.
+ *
+ * @param[in, out] dual_objective_value_ptr - A pointer to a cuopt_float_t that will contain the
+ * dual objective value.
+ *
+ * @return A status code indicating success or failure.
+ */
+cuopt_int_t cuOptGetDualObjectiveValue(cuOptSolution solution,
+                                       cuopt_float_t* dual_objective_value_ptr);
+
 /** @brief Get the reduced costs of an optimization problem.
  *
  * @param[in] solution - The solution object.

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

@@ -214,6 +214,12 @@ class optimization_problem_solution_t : public base_solution_t {
    */
   f_t get_objective_value() const;
 
+  /**
+   * @brief Returns the dual objective value of the solution as a `f_t`.
+   * @return objective value of the dual problem
+   */
+  f_t get_dual_objective_value() const;
+
   /**
    * @brief Returns the solution for the values of the primal variables as a vector of `f_t`.
    *

cpp/src/dual_simplex/solve.cpp

@@ -96,14 +96,14 @@ template <typename i_t, typename f_t>
 f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, const std::vector<f_t>& x)
 {
   const f_t obj      = compute_objective(lp, x);
-  const f_t user_obj = obj * lp.obj_scale + lp.obj_constant;
+  const f_t user_obj = compute_user_objective(lp, obj);
   return user_obj;
 }
 
 template <typename i_t, typename f_t>
 f_t compute_user_objective(const lp_problem_t<i_t, f_t>& lp, f_t obj)
 {
-  const f_t user_obj = obj * lp.obj_scale + lp.obj_constant;
+  const f_t user_obj = lp.obj_scale * (obj + lp.obj_constant);
   return user_obj;
 }
 

cpp/src/linear_programming/cuopt_c.cpp

@@ -24,6 +24,8 @@
 
 #include <mps_parser/parser.hpp>
 
+#include <cuopt/version_config.hpp>
+
 #include <memory>
 #include <string>
 
@@ -59,6 +61,19 @@ int8_t cuOptGetFloatSize() { return sizeof(cuopt_float_t); }
 
 int8_t cuOptGetIntSize() { return sizeof(cuopt_int_t); }
 
+cuopt_int_t cuOptGetVersion(cuopt_int_t* version_major,
+                            cuopt_int_t* version_minor,
+                            cuopt_int_t* version_patch)
+{
+  if (version_major == nullptr || version_minor == nullptr || version_patch == nullptr) {
+    return CUOPT_INVALID_ARGUMENT;
+  }
+  *version_major = CUOPT_VERSION_MAJOR;
+  *version_minor = CUOPT_VERSION_MINOR;
+  *version_patch = CUOPT_VERSION_PATCH;
+  return CUOPT_SUCCESS;
+}
+
 cuopt_int_t cuOptReadProblem(const char* filename, cuOptOptimizationProblem* problem_ptr)
 {
   problem_and_stream_view_t* problem_and_stream = new problem_and_stream_view_t();
@@ -825,6 +840,24 @@ cuopt_int_t cuOptGetDualSolution(cuOptSolution solution, cuopt_float_t* dual_sol
   }
 }
 
+cuopt_int_t cuOptGetDualObjectiveValue(cuOptSolution solution,
+                                       cuopt_float_t* dual_objective_value_ptr)
+{
+  if (solution == nullptr) { return CUOPT_INVALID_ARGUMENT; }
+  if (dual_objective_value_ptr == nullptr) { return CUOPT_INVALID_ARGUMENT; }
+  solution_and_stream_view_t* solution_and_stream_view =
+    static_cast<solution_and_stream_view_t*>(solution);
+  if (solution_and_stream_view->is_mip) {
+    return CUOPT_INVALID_ARGUMENT;
+  } else {
+    optimization_problem_solution_t<cuopt_int_t, cuopt_float_t>* optimization_problem_solution =
+      static_cast<optimization_problem_solution_t<cuopt_int_t, cuopt_float_t>*>(
+        solution_and_stream_view->lp_solution_ptr);
+    *dual_objective_value_ptr = optimization_problem_solution->get_dual_objective_value();
+    return CUOPT_SUCCESS;
+  }
+}
+
 cuopt_int_t cuOptGetReducedCosts(cuOptSolution solution, cuopt_float_t* reduced_cost_ptr)
 {
   if (solution == nullptr) { return CUOPT_INVALID_ARGUMENT; }

cpp/src/linear_programming/solver_solution.cu

@@ -309,6 +309,12 @@ f_t optimization_problem_solution_t<i_t, f_t>::get_objective_value() const
   return termination_stats_.primal_objective;
 }
 
+template <typename i_t, typename f_t>
+f_t optimization_problem_solution_t<i_t, f_t>::get_dual_objective_value() const
+{
+  return termination_stats_.dual_objective;
+}
+
 template <typename i_t, typename f_t>
 rmm::device_uvector<f_t>& optimization_problem_solution_t<i_t, f_t>::get_primal_solution()
 {

cpp/src/linear_programming/termination_strategy/convergence_information.cu

@@ -254,7 +254,7 @@ __global__ void apply_objective_scaling_and_offset(f_t* objective,
 {
   if (threadIdx.x + blockIdx.x * blockDim.x > 0) { return; }
 
-  *objective = (objective_scaling_factor * *objective) + objective_offset;
+  *objective = objective_scaling_factor * (*objective + objective_offset);
 }
 
 template <typename i_t, typename f_t>

cpp/src/math_optimization/solver_settings.cu

@@ -79,7 +79,7 @@ solver_settings_t<i_t, f_t>::solver_settings_t() : pdlp_settings(), mip_settings
     {CUOPT_MIP_ABSOLUTE_TOLERANCE, &mip_settings.tolerances.absolute_tolerance, 0.0, 1e-1, 1e-4},
     {CUOPT_MIP_RELATIVE_TOLERANCE, &mip_settings.tolerances.relative_tolerance, 0.0, 1e-1, 1e-4},
     {CUOPT_MIP_INTEGRALITY_TOLERANCE, &mip_settings.tolerances.integrality_tolerance, 0.0, 1e-1, 1e-5},
-    {CUOPT_MIP_ABSOLUTE_GAP, &mip_settings.tolerances.absolute_mip_gap, 0.0, 1e-1, 1e-10},
+    {CUOPT_MIP_ABSOLUTE_GAP, &mip_settings.tolerances.absolute_mip_gap, 0.0, CUOPT_INFINITY, 1e-10},
     {CUOPT_MIP_RELATIVE_GAP, &mip_settings.tolerances.relative_mip_gap, 0.0, 1e-1, 1e-4},
     {CUOPT_PRIMAL_INFEASIBLE_TOLERANCE, &pdlp_settings.tolerances.primal_infeasible_tolerance, 0.0, 1e-1, 1e-8},
     {CUOPT_DUAL_INFEASIBLE_TOLERANCE, &pdlp_settings.tolerances.dual_infeasible_tolerance, 0.0, 1e-1, 1e-8}

cpp/src/mip/presolve/trivial_presolve.cuh

@@ -234,7 +234,6 @@ void update_from_csr(problem_t<i_t, f_t>& pb)
 
   //  update objective_offset
   pb.presolve_data.objective_offset +=
-    pb.presolve_data.objective_scaling_factor *
     thrust::transform_reduce(handle_ptr->get_thrust_policy(),
                              thrust::counting_iterator<i_t>(0),
                              thrust::counting_iterator<i_t>(pb.n_variables),

cpp/src/mip/problem/problem.cu

@@ -67,19 +67,26 @@ void problem_t<i_t, f_t>::op_problem_cstr_body(const optimization_problem_t<i_t,
 
   // Set variables bounds to default if not set and constraints bounds if user has set a row type
   set_bounds_if_not_set(*this);
+
+  const bool is_mip = original_problem_ptr->get_problem_category() != problem_category_t::LP;
+  if (is_mip) {
+    variable_types =
+      rmm::device_uvector<var_t>(problem_.get_variable_types(), handle_ptr->get_stream());
+    // round bounds to integer for integer variables, note: do this before checking sanity
+    round_bounds(*this);
+  }
+
   // check bounds sanity before, so that we can throw exceptions before going into asserts
   check_bounds_sanity(*this);
+
   // Check before any modifications
   check_problem_representation(false, false);
   // If maximization problem, convert the problem
   if (maximize) convert_to_maximization_problem(*this);
 
-  const bool is_mip = original_problem_ptr->get_problem_category() != problem_category_t::LP;
   if (is_mip) {
     // Resize what is needed for MIP
     raft::common::nvtx::range scope("trivial_presolve");
-    variable_types =
-      rmm::device_uvector<var_t>(problem_.get_variable_types(), handle_ptr->get_stream());
     integer_indices.resize(n_variables, handle_ptr->get_stream());
     is_binary_variable.resize(n_variables, handle_ptr->get_stream());
     compute_n_integer_vars();
@@ -93,7 +100,9 @@ void problem_t<i_t, f_t>::op_problem_cstr_body(const optimization_problem_t<i_t,
 }
 
 template <typename i_t, typename f_t>
-problem_t<i_t, f_t>::problem_t(const optimization_problem_t<i_t, f_t>& problem_)
+problem_t<i_t, f_t>::problem_t(
+  const optimization_problem_t<i_t, f_t>& problem_,
+  const typename mip_solver_settings_t<i_t, f_t>::tolerances_t tolerances_)
   : original_problem_ptr(&problem_),
     handle_ptr(problem_.get_handle_ptr()),
     n_variables(problem_.get_n_variables()),
@@ -129,7 +138,8 @@ problem_t<i_t, f_t>::problem_t(const optimization_problem_t<i_t, f_t>& problem_)
     related_variables_offsets(n_variables, problem_.get_handle_ptr()->get_stream()),
     var_names(problem_.get_variable_names()),
     row_names(problem_.get_row_names()),
-    objective_name(problem_.get_objective_name())
+    objective_name(problem_.get_objective_name()),
+    tolerances(tolerances_)
 {
   op_problem_cstr_body(problem_);
   branch_and_bound_callback = nullptr;
@@ -1430,7 +1440,7 @@ void problem_t<i_t, f_t>::get_host_user_problem(
 template <typename i_t, typename f_t>
 f_t problem_t<i_t, f_t>::get_user_obj_from_solver_obj(f_t solver_obj)
 {
-  return solver_obj * presolve_data.objective_scaling_factor + presolve_data.objective_offset;
+  return presolve_data.objective_scaling_factor * (solver_obj + presolve_data.objective_offset);
 }
 
 #if MIP_INSTANTIATE_FLOAT

cpp/src/mip/problem/problem.cuh

@@ -54,7 +54,8 @@ constexpr bool USE_REL_TOLERANCE   = true;
 template <typename i_t, typename f_t>
 class problem_t {
  public:
-  problem_t(const optimization_problem_t<i_t, f_t>& problem);
+  problem_t(const optimization_problem_t<i_t, f_t>& problem,
+            const typename mip_solver_settings_t<i_t, f_t>::tolerances_t tolerances_ = {});
   problem_t() = delete;
   // copy constructor
   problem_t(const problem_t<i_t, f_t>& problem);