New heuristic improvements

benchmarks/linear_programming/cuopt/initial_problem_check.hpp

@@ -39,11 +39,12 @@ struct violation {
   }
 };
 
-void test_constraint_sanity(const cuopt::mps_parser::mps_data_model_t<int, double>& op_problem,
-                            const std::vector<double>& primal_vars,
-                            double abs_tol,
-                            double rel_tol,
-                            double int_tol = 1e-5)
+bool test_constraint_and_variable_sanity(
+  const cuopt::mps_parser::mps_data_model_t<int, double>& op_problem,
+  const std::vector<double>& primal_vars,
+  double abs_tol,
+  double rel_tol,
+  double int_tol = 1e-5)
 {
   const std::vector<double>& values                  = op_problem.get_constraint_matrix_values();
   const std::vector<int>& indices                    = op_problem.get_constraint_matrix_indices();
@@ -52,6 +53,7 @@ void test_constraint_sanity(const cuopt::mps_parser::mps_data_model_t<int, doubl
   const std::vector<double>& constraint_upper_bounds = op_problem.get_constraint_upper_bounds();
   const std::vector<double>& variable_lower_bounds   = op_problem.get_variable_lower_bounds();
   const std::vector<double>& variable_upper_bounds   = op_problem.get_variable_upper_bounds();
+  const std::vector<char>& variable_types            = op_problem.get_variable_types();
   std::vector<double> residual(constraint_lower_bounds.size(), 0.0);
   std::vector<double> viol(constraint_lower_bounds.size(), 0.0);
 
@@ -64,32 +66,41 @@ void test_constraint_sanity(const cuopt::mps_parser::mps_data_model_t<int, doubl
 
   auto functor = violation<double>{};
 
+  bool feasible = true;
   // Compute violation to lower/upper bound
   for (size_t i = 0; i < residual.size(); ++i) {
     double tolerance = abs_tol + combine_finite_abs_bounds<double>(constraint_lower_bounds[i],
                                                                    constraint_upper_bounds[i]) *
                                    rel_tol;
     double viol = functor(residual[i], constraint_lower_bounds[i], constraint_upper_bounds[i]);
     if (viol > tolerance) {
-      printf("error feasibility violation %f at cstr %d is more than total tolerance %f \n",
-             viol,
-             i,
-             tolerance);
-      exit(1);
+      feasible = false;
+      CUOPT_LOG_ERROR(
+        "feasibility violation %f at cstr %d is more than total tolerance %f lb %f ub %f \n",
+        viol,
+        i,
+        tolerance,
+        constraint_lower_bounds[i],
+        constraint_upper_bounds[i]);
     }
   }
-
+  bool feasible_variables = true;
   for (size_t i = 0; i < primal_vars.size(); ++i) {
+    if (variable_types[i] == 'I' && abs(primal_vars[i] - round(primal_vars[i])) > int_tol) {
+      feasible_variables = false;
+    }
     // Not always stricly true because we apply variable bound clamping on the scaled problem
     // After unscaling it, the variables might not respect exactly (this adding an epsilon)
     if (!(primal_vars[i] >= variable_lower_bounds[i] - int_tol &&
           primal_vars[i] <= variable_upper_bounds[i] + int_tol)) {
-      printf("error at bounds var %d lb %f ub %f val %f\n",
-             i,
-             variable_lower_bounds[i],
-             variable_upper_bounds[i],
-             primal_vars[i]);
-      exit(1);
+      CUOPT_LOG_ERROR("error at bounds var %d lb %f ub %f val %f\n",
+                      i,
+                      variable_lower_bounds[i],
+                      variable_upper_bounds[i],
+                      primal_vars[i]);
+      feasible_variables = false;
     }
   }
+  if (!feasible || !feasible_variables) { CUOPT_LOG_ERROR("Initial solution is infeasible"); }
+  return feasible_variables;
 }

benchmarks/linear_programming/cuopt/initial_solution_reader.hpp

@@ -14,6 +14,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
+
 #include <fstream>
 #include <iostream>
 #include <sstream>
@@ -24,7 +25,7 @@ class solution_reader_t {
  public:
   std::unordered_map<std::string, double> data_map;
 
-  bool readFromCsv(const std::string& filepath)
+  bool read_from_sol(const std::string& filepath)
   {
     std::ifstream file(filepath);
     if (!file.is_open()) {
@@ -37,11 +38,8 @@ class solution_reader_t {
       std::stringstream ss(line);
       std::string var_name;
       std::string value_str;
-
-      // Read var_name
-      if (!std::getline(ss, var_name, ',')) continue;
-      // Read value
-      if (!std::getline(ss, value_str)) continue;
+      ss >> var_name >> value_str;
+      if (var_name == "=obj=") continue;
 
       try {
         double value       = std::stod(value_str);

benchmarks/linear_programming/cuopt/run_mip.cpp

@@ -92,18 +92,17 @@ void write_to_output_file(const std::string& out_dir,
 
 inline auto make_async() { return std::make_shared<rmm::mr::cuda_async_memory_resource>(); }
 
-// reads a solution from an input file. The input file needs to be csv formatted
-// var_name,val
-std::vector<double> read_solution_from_file(const std::string file_path,
-                                            const std::vector<std::string>& var_names)
+void read_single_solution_from_path(const std::string& path,
+                                    const std::vector<std::string>& var_names,
+                                    std::vector<std::vector<double>>& solutions)
 {
   solution_reader_t reader;
-  bool success = reader.readFromCsv(file_path);
+  bool success = reader.read_from_sol(path);
   if (!success) {
-    CUOPT_LOG_INFO("Initial solution reading error!");
-    exit(-1);
+    CUOPT_LOG_ERROR("Initial solution reading error!");
   } else {
-    CUOPT_LOG_INFO("Success reading csv! Number of var vals %lu", reader.data_map.size());
+    CUOPT_LOG_INFO(
+      "Success reading file %s Number of var vals %lu", path.c_str(), reader.data_map.size());
   }
   std::vector<double> assignment;
   for (auto name : var_names) {
@@ -112,20 +111,42 @@ std::vector<double> read_solution_from_file(const std::string file_path,
     if ((it != reader.data_map.end())) {
       val = it->second;
     } else {
-      CUOPT_LOG_INFO("Variable %s has no input value ", name.c_str());
+      CUOPT_LOG_TRACE("Variable %s has no input value ", name.c_str());
       val = 0.;
     }
     assignment.push_back(val);
   }
-  CUOPT_LOG_INFO("Adding a solution with size %lu ", assignment.size());
-  return assignment;
+  if (assignment.size() > 0) {
+    CUOPT_LOG_INFO("Adding a solution with size %lu ", assignment.size());
+    solutions.push_back(assignment);
+  }
+}
+
+// reads a solution from an input file. The input file needs to be csv formatted
+// var_name,val
+std::vector<std::vector<double>> read_solution_from_dir(const std::string file_path,
+                                                        const std::string& mps_file_name,
+                                                        const std::vector<std::string>& var_names)
+{
+  std::vector<std::vector<double>> initial_solutions;
+  std::string mps_file_name_no_ext = mps_file_name.substr(0, mps_file_name.find_last_of("."));
+  // check if a directory with the given mps file exists
+  std::string initial_solution_dir = file_path + "/" + mps_file_name_no_ext;
+  if (std::filesystem::exists(initial_solution_dir)) {
+    for (const auto& entry : std::filesystem::directory_iterator(initial_solution_dir)) {
+      read_single_solution_from_path(entry.path(), var_names, initial_solutions);
+    }
+  } else {
+    read_single_solution_from_path(file_path, var_names, initial_solutions);
+  }
+  return initial_solutions;
 }
 
 int run_single_file(std::string file_path,
                     int device,
                     int batch_id,
                     std::string out_dir,
-                    std::optional<std::string> input_file_dir,
+                    std::optional<std::string> initial_solution_dir,
                     bool heuristics_only,
                     int num_cpu_threads,
                     bool write_log_file,
@@ -163,23 +184,36 @@ int run_single_file(std::string file_path,
     CUOPT_LOG_ERROR("Parsing MPS failed exiting!");
     return -1;
   }
-  if (input_file_dir.has_value()) {
-    auto initial_solution =
-      read_solution_from_file(input_file_dir.value(), mps_data_model.get_variable_names());
-    settings.set_initial_solution(initial_solution.data(), initial_solution.size());
-    test_constraint_sanity(mps_data_model,
-                           initial_solution,
-                           settings.tolerances.absolute_tolerance,
-                           settings.tolerances.relative_tolerance,
-                           settings.tolerances.integrality_tolerance);
+  if (initial_solution_dir.has_value()) {
+    auto initial_solutions = read_solution_from_dir(
+      initial_solution_dir.value(), base_filename, mps_data_model.get_variable_names());
+    for (auto& initial_solution : initial_solutions) {
+      bool feasible_variables =
+        test_constraint_and_variable_sanity(mps_data_model,
+                                            initial_solution,
+                                            settings.tolerances.absolute_tolerance,
+                                            settings.tolerances.relative_tolerance,
+                                            settings.tolerances.integrality_tolerance);
+      if (feasible_variables) {
+        settings.add_initial_solution(
+          initial_solution.data(), initial_solution.size(), handle_.get_stream());
+      }
+    }
   }
 
   settings.time_limit      = time_limit;
   settings.heuristics_only = heuristics_only;
   settings.num_cpu_threads = num_cpu_threads;
   settings.log_to_console  = log_to_console;
-  auto start_run_solver    = std::chrono::high_resolution_clock::now();
+  cuopt::linear_programming::benchmark_info_t benchmark_info;
+  settings.benchmark_info_ptr = &benchmark_info;
+  auto start_run_solver       = std::chrono::high_resolution_clock::now();
   auto solution = cuopt::linear_programming::solve_mip(&handle_, mps_data_model, settings);
+  CUOPT_LOG_INFO(
+    "first obj: %f last improvement of best feasible: %f last improvement after recombination: %f",
+    benchmark_info.objective_of_initial_population,
+    benchmark_info.last_improvement_of_best_feasible,
+    benchmark_info.last_improvement_after_recombination);
   // solution.write_to_sol_file(base_filename + ".sol", handle_.get_stream());
   std::chrono::milliseconds duration;
   auto end = std::chrono::high_resolution_clock::now();
@@ -199,7 +233,9 @@ int run_single_file(std::string file_path,
   std::stringstream ss;
   int decimal_places = 2;
   ss << std::fixed << std::setprecision(decimal_places) << base_filename << "," << sol_found << ","
-     << obj_val << "\n";
+     << obj_val << "," << benchmark_info.objective_of_initial_population << ","
+     << benchmark_info.last_improvement_of_best_feasible << ","
+     << benchmark_info.last_improvement_after_recombination << "\n";
   write_to_output_file(out_dir, base_filename, device, batch_id, ss.str());
   CUOPT_LOG_INFO("Results written to the file %s", base_filename.c_str());
   return sol_found;

cpp/cuopt_cli.cpp

@@ -135,7 +135,7 @@ int run_single_file(const std::string& file_path,
     if (is_mip && !solve_relaxation) {
       auto& mip_settings = settings.get_mip_settings();
       if (initial_solution.size() > 0) {
-        mip_settings.set_initial_solution(initial_solution.data(), initial_solution.size());
+        mip_settings.add_initial_solution(initial_solution.data(), initial_solution.size());
       }
       auto solution = cuopt::linear_programming::solve_mip(op_problem, mip_settings);
     } else {

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

@@ -29,6 +29,12 @@
 
 namespace cuopt::linear_programming {
 
+struct benchmark_info_t {
+  double last_improvement_of_best_feasible    = 0;
+  double last_improvement_after_recombination = 0;
+  double objective_of_initial_population      = std::numeric_limits<double>::max();
+};
+
 // Forward declare solver_settings_t for friend class
 template <typename i_t, typename f_t>
 class solver_settings_t;
@@ -44,35 +50,26 @@ class mip_solver_settings_t {
   void set_mip_callback(internals::base_solution_callback_t* callback = nullptr);
 
   /**
-   * @brief Set an primal solution.
+   * @brief Add an primal solution.
    *
-   * @note Default value is all 0 or the LP optimal point.
+   * @note This function can be called multiple times to add more solutions.
    *
-   * @param[in] initial_solution Device or host memory pointer to a floating
-   * point array of size size. cuOpt copies this data. Copy happens on the
-   * stream of the raft:handler passed to the problem.
+   * @param[in] initial_solution Device or host memory pointer to a floating point array of
+   * size size.
+   * cuOpt copies this data. Copy happens on the stream of the raft:handler passed to the problem.
    * @param size Size of the initial_solution array.
    */
-  void set_initial_solution(const f_t* initial_solution,
+  void add_initial_solution(const f_t* initial_solution,
                             i_t size,
                             rmm::cuda_stream_view stream = rmm::cuda_stream_default);
 
-  /**
-   * @brief Get the initial solution.
-   *
-   * @return Initial solution as a rmm::device_uvector<f_t>
-   */
-  rmm::device_uvector<f_t>& get_initial_solution() const;
-
   /**
    * @brief Get the callback for the user solution
    *
    * @return callback pointer
    */
   const std::vector<internals::base_solution_callback_t*> get_mip_callbacks() const;
 
-  bool has_initial_solution() const;
-
   struct tolerances_t {
     f_t presolve_absolute_tolerance = 1.0e-6;
     f_t absolute_tolerance          = 1.0e-4;
@@ -99,9 +96,11 @@ class mip_solver_settings_t {
   std::string sol_file;
   std::string user_problem_file;
 
-  /** Initial primal solution */
-  std::shared_ptr<rmm::device_uvector<f_t>> initial_solution_;
+  /** Initial primal solutions */
+  std::vector<std::shared_ptr<rmm::device_uvector<f_t>>> initial_solutions;
   bool mip_scaling = true;
+  // this is for extracting info from different places of the solver during benchmarks
+  benchmark_info_t* benchmark_info_ptr = nullptr;
 
  private:
   std::vector<internals::base_solution_callback_t*> mip_callbacks_;

cpp/include/cuopt/linear_programming/solver_settings.hpp

@@ -88,10 +88,11 @@ class solver_settings_t {
   const rmm::device_uvector<f_t>& get_initial_pdlp_dual_solution() const;
 
   // MIP Settings
-  void set_initial_mip_solution(const f_t* initial_solution, i_t size);
+  void add_initial_mip_solution(const f_t* initial_solution,
+                                i_t size,
+                                rmm::cuda_stream_view stream = rmm::cuda_stream_default);
   void set_mip_callback(internals::base_solution_callback_t* callback = nullptr);
 
-  const rmm::device_uvector<f_t>& get_initial_mip_solution() const;
   const pdlp_warm_start_data_view_t<i_t, f_t>& get_pdlp_warm_start_data_view() const noexcept;
   const std::vector<internals::base_solution_callback_t*> get_mip_callbacks() const;
 

cpp/src/linear_programming/utilities/problem_checking.cu

@@ -113,8 +113,8 @@ void problem_checking_t<i_t, f_t>::check_initial_solution_representation(
   const optimization_problem_t<i_t, f_t>& op_problem,
   const mip_solver_settings_t<i_t, f_t>& settings)
 {
-  if (settings.initial_solution_.get() != nullptr) {
-    check_initial_primal_representation(op_problem, settings.get_initial_solution());
+  for (const auto& initial_solution : settings.initial_solutions) {
+    check_initial_primal_representation(op_problem, *initial_solution);
   }
 }
 

cpp/src/math_optimization/solver_settings.cu

@@ -371,9 +371,11 @@ const rmm::device_uvector<f_t>& solver_settings_t<i_t, f_t>::get_initial_pdlp_du
 }
 
 template <typename i_t, typename f_t>
-void solver_settings_t<i_t, f_t>::set_initial_mip_solution(const f_t* solution, i_t size)
+void solver_settings_t<i_t, f_t>::add_initial_mip_solution(const f_t* solution,
+                                                           i_t size,
+                                                           rmm::cuda_stream_view stream)
 {
-  mip_settings.set_initial_solution(solution, size);
+  mip_settings.add_initial_solution(solution, size, stream);
 }
 
 template <typename i_t, typename f_t>
@@ -382,12 +384,6 @@ void solver_settings_t<i_t, f_t>::set_mip_callback(internals::base_solution_call
   mip_settings.set_mip_callback(callback);
 }
 
-template <typename i_t, typename f_t>
-const rmm::device_uvector<f_t>& solver_settings_t<i_t, f_t>::get_initial_mip_solution() const
-{
-  return mip_settings.get_initial_solution();
-}
-
 template <typename i_t, typename f_t>
 const std::vector<internals::base_solution_callback_t*>
 solver_settings_t<i_t, f_t>::get_mip_callbacks() const